3

I am working with data from a meta-analysis of studies, where not every study had data for every genetic variant. I am trying to solve the following problem:

  • I have a data frame ("studies") with information about 10 studies, which includes the number of cases and controls in each study and the total study sample size (N = cases + controls).
  • I also have a data frame ("rawresults") with approximately 10 million rows, each of which was generated by summing some unknown subset of the sample sizes (N) from the table of studies.
  • I need to figure out which studies were added up to create each total

Here are examples of what the data looks like:

> head(studies,2)
  study cases controls     N
1     A  3747     8024 11771
2     B  5367     5780 11147

> head(rawresults,2)
          ID      N
1 rs58241367  65280
2 rs85436064 107624

Here is an example of what I want to make from it (the contributing_studies column in optional, I can get rid of it if doing so allows for a better solution):

> head(final,2)
          ID contributing_studies cases controls     N
1 rs10685984            B,C,D,F,G 26221    19987 46208
2 rs12123751            A,C,D,G,J 25631    23509 49140

So far my best idea for how to approach this problem is to brute-force it. Each of the ten studies has two possible states (contributing vs. not contributing) for a given total, so that's 2^10 = 1024 possible sums. Some sums may not be unique (there may be more than one combination of study Ns that could create that sum) and I would plan to exclude those as ambiguous. I have included code below with an example of that solution as an answer.

What I want to ask: Does R have a better way to do this? Perhaps some library or function that exists for handling this kind of problem? Or is there something else I could do to make it faster and more efficient?

Here is code to simulate data for the scenario:

set.seed(1)

# Make "studies"
studies <- data.frame(toupper(letters[1:10]),round(rnorm(10,5000,2000)),round(rnorm(10,5000,2000)),stringsAsFactors=F)
colnames(studies) <- c('study','cases','controls')
studies$N <- studies$cases + studies$controls

# Make "rawresults"
rawresults <- data.frame(character(length=50),numeric(length=50),stringsAsFactors=F)
colnames(rawresults) <- c('ID','N')
for(i in seq(1,50)) {
  numstudies <- sample(seq(5,10),1)
  rawresults[i,'N'] <- sum(sample(studies$N,numstudies))
  rawresults[i,'ID'] <- paste0('rs',sample(seq(1,99999999),1))
}

Edit: Faster code to simulate data for the scenario, so it will be possible to simulate millions of rawresults rows. Inspired by Allan Cameron's solution below, which also uses combn.

set.seed(1)

# Make "studies"
studies <- data.frame(toupper(letters[1:10]),round(rnorm(10,5000,2000)),round(rnorm(10,5000,2000)),stringsAsFactors=F)
colnames(studies) <- c('study','cases','controls')
studies$N <- studies$cases + studies$controls

# Make "rawresults"
num_results <- 50 # Number of results to simulate
possible_ns <- unlist(sapply(1:10,combn,x=studies$N,sum))
rawresults <- data.frame(paste0('rs',sample(1:99999999,num_results)),sample(possible_ns,num_results,rep=T),stringsAsFactors=F)
colnames(rawresults) <- c('ID','N')
7
  • If you the contribution for each row is different.. it is not possible. If it is the same, then you take some rows, and do a non-negative least square.. r-bloggers.com/non-negative-least-squares
    – StupidWolf
    Commented Jun 29, 2020 at 22:24
  • @MrFlick I have edited it to move my brute-force solution to an answer as you suggested - thanks! Commented Jun 29, 2020 at 22:31
  • @StupidWolf Thanks for the link! When you say if the contribution for each row is different, what do you mean? About half of my rows have all studies. The rest have different totals made up of some studies, but the totals repeat, none are unique to a single row. Commented Jun 29, 2020 at 22:33
  • let's say every row is has a total as response, and the total has x1,x2,x3 contribution from study1, study2,study3 that makes up the total.. what i meant is, is x1,x2,x3 the same for every row. maybe I am complicating things.
    – StupidWolf
    Commented Jun 29, 2020 at 22:38
  • 1

3 Answers 3

5

I think there is certainly a way of making the code a bit shorter by using the combn function:

getSum <- function(x) colSums(matrix(studies$N[combn(1:10, x)], nrow = x))
getInd <- function(x) apply(combn(1:10, x), 2, function(y) paste(y, collapse = ", "))
all_sums <- do.call(c, lapply(1:10, getSum))
all_inds <- do.call(c, lapply(1:10, getInd))
rawresults$studies <- all_inds[match(rawresults$N, all_sums)]

This yields:

rawresults
#>            ID      N                       studies
#> 1  rs58241367  65280              2, 4, 6, 7, 8, 9
#> 2  rs85436064 107624 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
#> 3  rs64407295 107624 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
#> 4  rs78593369  83683       1, 3, 4, 5, 6, 7, 8, 10
#> 5  rs18774630 107624 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
#> 6  rs99681114  49670                3, 6, 7, 9, 10
#> 7   rs8426283  69694          2, 3, 4, 5, 6, 7, 10
#> 8  rs81116968  75972          2, 4, 5, 6, 7, 8, 10
#> 9  rs54871836  55138                1, 3, 5, 9, 10
#> 10 rs21386862  87919       1, 2, 3, 5, 6, 8, 9, 10
#> 11 rs16179951  73195           1, 2, 3, 4, 6, 8, 9
#> 12  rs8492848  74843          1, 3, 4, 5, 7, 9, 10
#> 13 rs81342050  56555                 2, 4, 5, 7, 9
#> 14 rs44945811  59794             1, 2, 3, 6, 7, 10
#> 15 rs43997413  94715    1, 2, 3, 4, 6, 7, 8, 9, 10
#> 16 rs97055078  94715    1, 2, 3, 4, 6, 7, 8, 9, 10
#> 17 rs19941161  51106                1, 3, 4, 5, 10
#> 18 rs30748262  94259    1, 2, 3, 4, 5, 6, 7, 9, 10
#> 19 rs10959349  64914             2, 4, 6, 8, 9, 10
#> 20   rs982457  99355    1, 2, 3, 4, 5, 7, 8, 9, 10
#> 21  rs2007022 107624 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
#> 22 rs21833759 100537    1, 2, 4, 5, 6, 7, 8, 9, 10
#> 23 rs74715222  87172       1, 2, 4, 5, 6, 7, 9, 10
#> 24 rs49182929  54248                 4, 5, 6, 7, 8
#> 25 rs95501056  64548              1, 2, 3, 5, 6, 8
#> 26 rs57556390  53270                 2, 3, 4, 5, 8
#> 27 rs98150573 107624 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
#> 28 rs12123751  49140                1, 3, 4, 7, 10
#> 29 rs97971902 107624 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
#> 30 rs89722202  50848                 2, 3, 4, 5, 7
#> 31 rs98543720  65904              1, 4, 6, 7, 8, 9
#> 32 rs31961269  70318          1, 3, 4, 5, 6, 7, 10
#> 33  rs9764985 107624 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
#> 34 rs10685984  46208                 2, 3, 4, 6, 7
#> 35 rs34912847  77631           1, 3, 4, 5, 7, 8, 9
#> 36 rs94227949  76514           2, 3, 5, 6, 7, 8, 9
#> 37 rs74751789  67467             1, 2, 5, 6, 9, 10
#> 38 rs47441925  66565             1, 2, 4, 7, 8, 10
#> 39  rs4502074  69920              2, 4, 5, 7, 8, 9
#> 40  rs2741222  57384                1, 4, 5, 8, 10
#> 41 rs11561555  79017           1, 2, 4, 5, 6, 8, 9
#> 42 rs96740802  85900        1, 3, 4, 5, 6, 7, 8, 9
#> 43  rs7081648  96681    1, 2, 3, 4, 5, 6, 8, 9, 10
#> 44 rs68842412  74957           1, 3, 4, 5, 6, 8, 9
#> 45 rs78557028 107624 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
#> 46 rs75435758  84706       3, 4, 5, 6, 7, 8, 9, 10
#> 47 rs58940608  58751             2, 3, 4, 5, 6, 10
#> 48 rs25060056  66639             2, 5, 6, 7, 9, 10
#> 49 rs98833089  99355    1, 2, 3, 4, 5, 7, 8, 9, 10
#> 50 rs70239332  94259    1, 2, 3, 4, 5, 6, 7, 9, 10

This takes under 10 milliseconds to find all the combinations of N and generate the index strings. It only does this once, then the efficiency is just down to the efficiency of the match function, which is about as good as R offers for this kind of algorithm. The whole thing runs in about 11ms on my modest machine.

Created on 2020-06-29 by the reprex package (v0.3.0)

1
  • Fantastic! combn() makes it more compact, and match() makes it much, much faster. I scaled up the simulation to 2 million rows, and even when I have extra data (cases, controls, etc.) that I add to rawresults the same way, using match() to add the data takes 1/3 of the time it would take using merge()! Commented Jun 30, 2020 at 0:55
1

Here is the solution I came up with using a brute-force approach.

I am the person who asked a question, and I am hoping someone will have a better solution than mine, eg. something in R meant for deconvolution (I think that's the right word?) of sums that could be applied to this problem.

#### Bruce-force generate all possible combinations of studies ####

sumstudies <- function(whichstudies,whichcolumn) {
  # Convert integer "whichstudies" to binary, then use the binary digits to decide which studies are included or excluded for this combination
  in_or_out <- as.logical(intToBits(whichstudies)[1:10])
  # Return appropriate combination of data from included studies (sum if numeric, paste otherwise)
  if(is.numeric(studies[,whichcolumn])) {
    return(sum(studies[in_or_out,whichcolumn]))
  } else {
    return(paste(studies[in_or_out,whichcolumn],collapse=','))
  }
}

# Create a data frame with all 1024 possible combinations of studies
allcombos <- data.frame(matrix(nrow=1024,ncol=4))
colnames(allcombos) <- c('contributing_studies','cases','controls','N')
allcombos$contributing_studies <- sapply(seq(1,1024),sumstudies,'study')
allcombos$N <- sapply(seq(1,1024),sumstudies,'N')
allcombos$cases <- sapply(seq(1,1024),sumstudies,'cases')
allcombos$controls <- sapply(seq(1,1024),sumstudies,'controls')

# Get rid of Ns that can be made by summing more than one different combination of studies, since we wouldn't know which solution was correct
duplicates <- duplicated(allcombos$N) | duplicated(allcombos$N,fromLast=T)
allcombos[duplicates,] <- NA # Set all affected rows to NA

#### Match the data about all possible combinations to the Ns in rawresults ####

final <- merge(rawresults,allcombos,by='N',all.x=T,all.y=F)
final <- final[,c('ID','contributing_studies','cases','controls','N')]
1

For anyone who Googles across this later, here's the code I wrote (based on Allan Cameron's solution) that adds all the fields you might want when working with genetic meta-analysis data. It took only 16 seconds to do 2 million rows.

allcomb <- function(studies,excl_ambig=T) {
  n_studies <- nrow(studies)
  N <- unlist(sapply(1:n_studies,combn,x=studies$N,sum))
  cases <- unlist(sapply(1:n_studies,combn,x=studies$cases,sum))
  controls <- unlist(sapply(1:n_studies,combn,x=studies$controls,sum))
  contrib_studies <- unlist(sapply(1:n_studies,combn,x=studies$study,paste,collapse=','))
  combined <- data.frame(contrib_studies,cases,controls,N,stringsAsFactors=F)
  # Flag ambiguous lines where more than one possible combination of studies exists to produce that sum
  combined$ambig <- duplicated(combined$N) | duplicated(combined$N,fromLast=T)
  if(excl_ambig) {
    combined <- combined[!combined$ambiguous,]
  }
  return(combined)
}
allcombos <- allcomb(studies)
rawresults[,c('contrib_studies','cases','controls','N')] <- allcombos[match(rawresults$N, allcombos$N),]

Accepting Allan's solution as the answer since this was his idea, I just built on it!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.