重疊基因組範圍

Question

我有兩個文件

編碼

X.pattern.name  chr  start stop strand score p.value q.value matched.sequence 
1 V_CETS1P54_01 chr1 98769545 98769554 + 11.42280 8.89e-05 NA TCAGGATGTA 
2 V_CETS1P54_01 chr1 152013037 152013046 + 11.98020 4.74e-05 NA ACAGGAAGTT 
3 V_CETS1P54_01 chr1 168932563 168932572 + 11.60860 7.59e-05 NA ACCGGATGCT 

encode.total

chr  start   stop 
1 chr1 58708485 58708713 
2 chr1 58709084 58710538 
3 chr1 98766295 98766639 
4 chr1 98766902 98770338 
5 chr1 107885506 107889414 
6 chr1 138589531 138590856 
7 chr1 138591180 138593378 
8 chr1 152011746 152013185 
9 chr1 152014263 152014695 
10 chr1 168930561 168933076 
11 chr1 181808064 181808906 
12 chr1 184609002 184611519 
13 chr1 193288453 193289567 
14 chr1 193290105 193290490 
15 chr1 193290744 193291092 
16 chr1 196801920 196804954 

我想這兩個文件，通過CHR每個條目比較，開始和停止。如果第一個文件中的開始和結束值落在同一染色體的第二個文件的開始和結束之間，那麼在第一個文件中的開始&停止值應該被替換爲第二個文件的開始&停止值。我爲此寫了一個for循環，但它耗時太長。有什麼選擇？

代碼：

for(i in 1:nrow(encode)) 
{ 
    for(j in 1:nrow(encode.total)) 
    { 
      if(encode[i,2]==encode.total[j,1]) 
      { 
       if((encode[i,3]>=encode.total[j,2]) & (encode[i,4]<=encode.total[j,3])) 
       { 
        encode[i,3]=encode.total[j,2] 
        encode[i,4]=encode.total[j,3] 
       } 
      } 
    } 
} 

出於同樣的目的，我也試過GenomicRanges包像下面。我的數據框的大小是巨大的，它們上的合併函數創建了一個非常大的數據框（大於20億行，這是不允許的），儘管我最終將數據集劃分爲更小的數據框。但合併（）正在採取的記憶很多，終止R.

gr1<-GRanges(seqnames=encode$chr,IRanges(start=encode$start,end=encode$end)) 
gr2<-GRanges(seqnames=encode.total$chr, IRanges(start=encode.total$start,end=encode.total$end)) 

ranges <- merge(as.data.frame(gr1),as.data.frame(gr2),by="seqnames",suffixes=c("A","B")) 
ranges <- ranges[with(ranges, startB <= startA & endB >= endA),] 

Answer 1

使用BioconductorGenomicRanges包。

source("http://bioconductor.org/biocLite.R") 
biocLite("GenomicRanges") 

假設有兩個數據幀x0和x1，像在原始的示例encode和encode.total。我們希望將這些變成一個GRanges實例。我這樣做是與

library(GenomicRanges) 
gr0 = with(x0, GRanges(chr, IRanges(start=start, end=stop))) 
gr1 = with(x1, GRanges(chr, IRanges(start=start, end=stop))) 

它往往可以簡單地說makeGRangesFromDataFrame(x0)，或者使用標準的R命令「手工」打造一個農莊實例。這裏我們使用with()，這樣我們可以寫GRanges(chr, IRanges(start=start, end=stop))而不是GRanges(x0$chr, IRanges(start=x0$start, end=x0$stop))。

下一步是找到查詢（GR0）和主題（GR1）

hits = findOverlaps(gr0, gr1) 

這導致

> hits 
Hits of length 3 
queryLength: 3 
subjectLength: 16 
    queryHits subjectHits 
    <integer> <integer> 
1   1   4 
2   2   8 
3   3   10 

之間的重疊。然後更新相關的開始/結束座標

ranges(gr0)[queryHits(hits)] = ranges(gr1)[subjectHits(hits)] 

給予

> gr0 
GRanges with 3 ranges and 0 metadata columns: 
     seqnames     ranges strand 
     <Rle>    <IRanges> <Rle> 
    [1]  chr1 [ 98766902, 98770338]  * 
    [2]  chr1 [152011746, 152013185]  * 
    [3]  chr1 [168930561, 168933076]  * 
    --- 
    seqlengths: 
    chr1 
    NA 

這將快速上升到數百萬個範圍。