基於重疊合並2個基因組文件

Question

我有兩個「BED」文件。每一個指定基因組中的一組區域（開始和結束列），並且這些文件中的每一個指定給定基因組區域的特徵（例如NRL和另一個返回這些區域的'可映射性'）

它們是安排如下：

head(file1) 
    chr start  end mappability 
    chr1 3000066 3000100  1.0000 
    chr1 3000100 3000130  0.5000 
    chr1 3000130 3000199  0.0625 
    chr1 3000199 3000277  0.0500 


head(file2) 
    chr start  end NRL 
    chr1 3000000 3000067 250 
    chr1 3000067 3000079 300 
    chr1 3000079 3000084 200 
    chr1 3000084 3000099 130 

的問題是，這些文件是不同的實驗結果，而不是所有的都是兩個文件之間記錄的區域重疊......所以我需要找出哪些地區重疊。 ..

我對此的嘗試迄今爲止：

file1-read.table("file1.txt", sep='\t', header = F) 
file2=read.table("file2.txt", sep='\t', header = F) 


overlapping_regions<-function(file1, file2){ 
    for(i in file1[,2]){ 
    x<-seq(file1[i,2], file1[i,3]) 
    for(j in file1[,2]){ 
     y<-seq(file2[j,2], file2[j,3]) 
     if(setequal(union(x, y), c(setdiff(x, y), intersect(x, y), setdiff(y, x)))){ 
     ####GET OVERLAP 
     } 
    } 
    } 
} 

與上述戰略的第一個問題是，我得到上述錯誤：

Error in seq.default(file1[i, 2], file1[i, 3]) : 

「從」

不能NA，NaN或無窮

其次，我不知道這策略是正確的，因爲我希望每個文件的每一行都能與另一個進行比較，找到ANY區域重疊...

所以我想知道如果有人可以幫我一個盧比cript解析這些文件，以便我可以創建一個新文件，其中包含每個開始和結束之間指定列的重疊區域，並保留與每個原始文件相關的功能...

所以我希望我的輸出是這樣的：

head(files_merged) 

chr overlap mappability  NRL GC_content more_features...... 
chr1 start-end  1.0000  250 
chr1 start-end  0.5000  300 
chr1 start-end  0.0625  200 

我問這與嘗試應用機器學習算法來嘗試預測基因組特徵。

我可以看到（很明顯）我的計劃如何存在缺陷，因爲一個文件中指定的區域可能比另一個文件中的區域小得多。因此，我也開始建議更好的方法來做到這一點？

Answer 1

這可能有點長，但你可以試試看。

我創建類似dataframes，但並不確切：使用地圖和嵌套函數需要

df1 <- data.frame(chr=rep("chr1",4), 
        start=c(100,200,300,400), 
        end=c(200,300,400,500), 
        mappability=c(1,0.5,0.0625,0.05)) 

df2 <- data.frame(chr=rep("chr1",4), 
        start=c(90,190,290,380), 
        end=c(120,220,320,390), 
        NRL=c(250,300,200,130)) 

加載庫：

library(purrr) 
library(tidyr) 

，需要一個tibble有開始和結束的函數，尋找df1中存在重疊並返回行號的索引。 您可以根據自己的邊界，約束或重疊的定義編輯條件：

xx <- function(x){ 
     y <- (x$start<df1$start & x$end<df1$end & x$end>df1$start) | (x$start>df1$start & x$start<df1$start & x$end>df1$end) 

     z <- which(y==TRUE) 

     ifelse((length(z)>0),z,0) %>% 
       as.integer() 
} 

巢DF2，把開始結束在一個tibble：

df2 <- df2 %>% 
     nest(start,end,.key=data.df2) 

# A tibble: 4 x 3 
    chr NRL   data.df2 
    <fctr> <dbl>   <list> 
1 chr1 250 <tibble [1 x 2]> 
2 chr1 300 <tibble [1 x 2]> 
3 chr1 200 <tibble [1 x 2]> 
4 chr1 130 <tibble [1 x 2]> 

通過tibble每一行中發揮作用xx這將返回重疊行（如果有情況下，將有多個條目，該功能可能需要改變，我們將使用地圖，而不是map_int）

df2 <- df2 %>% 
     mutate(idx=map_int(data.df2,xx)) %>% 
     unnest %>% 
     filter(idx!=0) 

解開並刪除沒有交集的行之後，我們將在df2中包含df1中具有重疊條目的條目。

# A tibble: 3 x 5 
    chr NRL idx start end 
    <fctr> <dbl> <int> <dbl> <dbl> 
1 chr1 250  1 90 120 
2 chr1 300  2 190 220 
3 chr1 200  3 290 320 

我們將增加一個IDX列DF1能夠合併：

DF1 < - DF1％>％ 變異（IDX = seq_along（DF1））

chr start end mappability idx 
1 chr1 100 200  1.0000 1 
2 chr1 200 300  0.5000 2 
3 chr1 300 400  0.0625 3 
4 chr1 400 500  0.0500 4 

現在根據索引合併df1和df2：

df_all <- merge(df1,df2,by=c("idx"), 
     all.x = FALSE, 
     all.y = TRUE 
    ) 

TOu會有類似的東西，在這裏你可以清潔並且在每行計算重疊：

idx chr.x start.x end.x mappability chr.y NRL start.y end.y 
1 1 chr1  100 200  1.0000 chr1 250  90 120 
2 2 chr1  200 300  0.5000 chr1 300  190 220 
3 3 chr1  300 400  0.0625 chr1 200  290 320 

Answer 2

的問題也被要求在Bioconductor support site，在這裏我提供一個同樣長的答案。用於通過@OmaymaS提供的數據的結果是

> olaps 
GRanges object with 6 ranges and 2 metadata columns: 
     seqnames  ranges strand | mappability  NRL 
     <Rle> <IRanges> <Rle> | <numeric> <numeric> 
    [1]  chr1 [101, 120]  * |   1  250 
    [2]  chr1 [191, 200]  * |   1  300 
    [3]  chr1 [201, 220]  * |   0.5  300 
    [4]  chr1 [291, 300]  * |   0.5  200 
    [5]  chr1 [301, 320]  * |  0.0625  200 
    [6]  chr1 [381, 390]  * |  0.0625  130 
    ------- 
    seqinfo: 1 sequence from an unspecified genome; no seqlengths 

與基於1從BED文件的基於0的，半開區間的平移偏移到更友好/基於1的Bioconductor的標準，閉區間。