слияние файлов генома на основе перекрытия

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 регионов, которые перекрываются ...

Так что мне интересно, может ли кто-нибудь помочь мне с Rs 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) 

Функция, которая принимает Тиббл с начала и конца, искать индекс в 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 и поставить запуск конец в одном Тиббл:

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]> 

пройти Тиббл в каждой строке функционировать 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 на основе смещения от перевода 0 на основе, полусегмент КОЙКИ файла к более дружественной/Bioconductor стандарт 1 на основе, отрезки.