Take the 2-minute tour ×
Code Review Stack Exchange is a question and answer site for peer programmer code reviews. It's 100% free, no registration required.

I wrote an R code that basically performs 2 operations:

  1. For each segment in file A, find all segments in file B that lie in that segment.

  2. Find the percentage of overlap for each case in previous point.

The code that I have written works fine but given I have 74 files of type A and the rows in file B are about 200K, the approximate run time for this on my cluster is about ~18 days! I am writing the code and input files below and would really appreciate your comments on how I can improve the run time, with code if possible.

File A

Chromosome    Start       End
1         0           2420037
1         2420037     2522634
1         2522634     10794763
1         10794763    10925915
1         10925915    11280057
...

File B

chr       start   end     variantaccession    variantsubtype
1     10001   127330  nsv7879             Gain+Loss
1     10001   846808  dgv2n71             Gain
1     10377   177417  dgv1e1              Complex
1     10377   177417  nsv428112           Gain
1     10377   1018704 dgv3e1              Complex
1     10499   177368  esv27265            Gain+Loss
...

A lot of rows from file B will map to each row of file A. Once my code finds indices of mapped rows from file B it also extracts certain other values part of additional columns of file B.

setwd("/abc/xyz")

filenames <- basename(list.files(path=getwd())) # reading list of all files of type A

# Reading File B
ref <- read.table("xyz.txt", header=TRUE,sep="\t")

for (i in 1:length(filenames))
{
  # Reading file
  file <- read.table(filenames[i],header=TRUE,sep="\t")
  cat("Working with ",filenames[i], " i.e. ",i,"/",length(filenames),"\n")
  # Start the clock!
  ptm <- proc.time()
  for(j in 1:nrow(file))
  {
    cat(j,"\n")
    ind <- which(file$Chromosome[j] == ref$chr)
    ref2 <- ref[ind,]
    B <- ""
    C <- ""
    D <- ""
    count <- 0
    for(k in 1:length(ind))
    {
      if((ref2$start[k]>=file$Start[j] & ref2$start[k]<=file$End[j]) | (ref2$end[k]>=file$Start[j] & ref2$end[k]<=file$End[j]))
      {
        B[k] <- as.character(ref2$variantaccession[k])
        C <- as.integer(append(C,ref2$start[k]:ref2$end[k]))
        D[k] <- as.character(ref2$variantsubtype[k])
        count <- count+1
      }
    }
    file$CNVcount[j] <- count
    C <- C[2:length(C)]
    E <- file$Start[j]:file$End[j]
    file$PerCNVoverlap[j] <- (length(E)/length(which(C %in% E)))*100
    B2 <- paste(unique(B[complete.cases(B)]),collapse=",")
    B3 <- paste(B2,",",sep="")
    file$CNVs[j] <- B3
    D2 <- paste(unique(D[complete.cases(D)]),collapse=",")
    D3 <- paste(D2,",",sep="")
    file$CNVtypes <- D3
  }
  # Stop the clock
  proc.time() - ptm
  file <- file[,c(1,2,3,4,5,6,8,7)]
  write.table(file,filenames[i],col.names=TRUE,row.names=FALSE,sep="\t")
}

How I am calculating overlap:

  File A, row 01: 1 2 3 4 5 6 7 8 9 10
  File B, seg 01:   2 3 4
  File B, seg 32:         5 6 7 
  File B, seg 12:     3 4 5 6 7 8 9
Consensus File B:   2 3 4 5 6 7 8 9
Percentage overlap = 8/10 -> 80%

Required Output (i.e add additional columns to File A)

Chromosome    Start       End    Overlap_Count     Overlap_Percentage     Overlaps_Names     Overlap_subtype
1         0           2420037    800             22.12%            nsv7879,dgv1e1...    Gain+Loss,Complex     
1         2420037     2522634    626             35.12%            nsv7879,dgv1e1...    Gain+Loss,Complex     
1         2522634     10794763    200             17.12%            nsv7879,dgv1e1...    Gain+Loss,Complex     
1         10794763    10925915    75             42.12%            nsv7879,dgv1e1...    Gain+Loss,Complex     
1         10925915    11280057    800             22.12%            nsv7879,dgv1e1...    Gain+Loss,Complex
share|improve this question

migrated from stackoverflow.com Jul 23 at 11:45

This question came from our site for professional and enthusiast programmers.

    
Have you tried using bedtools instead of R? –  zx8754 Jun 10 at 11:14
    
@zx8754 No, but I have heard about it and I can surely give it a go if I am sure it can calculate the % overlap for me. –  Jason Jun 10 at 11:18
    
If you use intersect of 2 bed files using bedtools, then you can count just count of rows on output file wc -l output.txt. –  zx8754 Jun 10 at 11:19
    
If nothing else, initialize outside the loop, e.g. D<-vector(length=length(ind)) rather than rebuilding each time thru the loop. I suspect you could do a histogram on the start and end values in "File B" with the bins defined by the values in "File A" to see which entries fall into a "File A-bin" . It's not clear what you mean by "Percentage overlap"- I would call that "Coverage", which is simply a count of how many elements in each row of A exist somewhere in B; a much simpler problem to solve. –  Carl Witthoft Jun 10 at 11:31

1 Answer 1

I read the data in

a <- read.delim(textConnection("Chromosome\tStart\tEnd
1\t0\t2420037
1\t2420037\t2522634
1\t2522634\t10794763
1\t10794763\t10925915
1\t10925915\t11280057"))

b <- read.delim(textConnection("chr\tstart\tend\tvariantaccession\tvariantsubtype
1\t10001\t127330\tnsv7879\tGain+Loss
1\t10001\t846808\tdgv2n71\tGain
1\t10377\t177417\tdgv1e1\tComplex
1\t10377\t177417\tnsv428112\tGain
1\t10377\t1018704\tdgv3e1\tComplex
1\t10499\t177368\tesv27265\tGain+Loss"))

and then used the GenomicRanges package to efficiently find overlaps between the query "A" and subject "B".

library(GenomicRanges)
A <- makeGRangesFromDataFrame(a)
B <- makeGRangesFromDataFrame(b, keep.extra.columns=TRUE)
olaps <- findOverlaps(A, B)

olaps is basically a two-column matrix, indicating which element(s) of A overlap the corresponding elements of B. In the invocation above A and B enter symmetrically, but could be re-ordered if there were restrictions on the type of overlap (argument type, e.g., all of B strictly within A). To find the extent of overlap, I looked at the parallel (element wise) intersection of each overlap

> isect <- pintersect(A[queryHits(olaps)], B[subjectHits(olaps)])
> isect
GRanges with 6 ranges and 0 metadata columns:
      seqnames           ranges strand
         <Rle>        <IRanges>  <Rle>
  [1]        1 [10001,  127330]      *
  [2]        1 [10001,  846808]      *
  [3]        1 [10377,  177417]      *
  [4]        1 [10377,  177417]      *
  [5]        1 [10377, 1018704]      *
  [6]        1 [10499,  177368]      *
  ---
  seqlengths:
    1
   NA

I could then create a data.frame with whatever information I want, e.g.,

data.frame(query=queryHits(olaps), subject=subjectHits(olaps),
           olap_width=width(isect), 
           query_width=width(A)[queryHits(olaps)],
           variantaccession=B$variantaccession[subjectHits(olaps)])

This will be fast for millions of records in A and B.

makeGRangesFromDataFrame is available in the current release version of the package, installed by default when using R-3.1. It's easy to make a GRanges 'by hand', e.g.,

A <- with(a, GRanges(Chromosome, IRanges(Start, End)))
B <- with(b, GRanges(chr, IRanges(start, end), variantaccession=variantaccession))
share|improve this answer
    
Thanks. That looks impressive, somehow I can't find the function makeGRangesFromDataFrame() in the package GenomicRanges. It's also not listed in package functions. –  Jason Jun 10 at 12:04
    
@Jason this was introduced in the most recent version of GenomicRanges, which is installed by default if you are using R-3.1. It's easy to make a GRanges object and I'll illustrate in an updated answer. –  Martin Morgan Jun 10 at 12:07
    
Excellent -- always better to use easily-found tools than to reinvent code. –  Carl Witthoft Jun 10 at 13:34
    
@MartinMorgan Thanks Martin, that makes a lot of sense and GenomicRanges does help. However, since it will take me some time to get used to its abilities, could you please help tailor the output from your code to the output I require? I've updated my output in my question (which I forgot to put in initially) Your code sure helps find the overlapping segments but the variable I am most interested in is percentage overlap calculation of which I have show above. Thanks –  Jason Jun 10 at 13:51
    
@Jason -- can you clarify what you mean, for instance using the numerical values in your example above, by 'percentage overlap'? –  Martin Morgan Jun 10 at 22:33

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.