Охват путем пересечения меньших данных геномного интервала по большим геномным интервалам с использованием R

Question

Я хочу пересечь два геномных интервала в R. И я хочу получить статистику покрытия меньшего интервала за больший интервал.

Данные с большим интервалом представляют собой фрейм данных, подобный этому ...

Chr  Start     End       Name         Val    Strand
chr7 145444998 146102295 CCDS5889.1   0      +
chr7 146102406 146167735 CCDS5889.1   0      +
chr7 146167929 146371931 CCDS5889.1   0      +

меньший интервал с более чем 2 миллионами строк выглядит следующим образом.

Chr  Start     End       Name         Val    Strand PhyloP   
chr7 145444386 145444387 CCDS5889.1   0      +      0.684764
chr7 145444387 145444388 CCDS5889.1   0      +      0.684764
chr7 145444388 145444389 CCDS5889.1   0      +      0.684764
chr7 145444389 145444390 CCDS5889.1   0      +      0.684764

Данные интервала находятся во 2-м (с) и 3-м (до) столбцах в обоих фреймах данных.

Ситуация похожа на

Large Interval:    [-----]   [-----]     [--------------]   [-------------------]
Small Interval: |||  ||||  |||||||||||  ||||||||   ||||  || |||||||||   ||    ||||||||

1. Я хочу знать, сколько из каждого из больших интервалов охватывается меньшими интервалами.
2. ТакжеЯ хотел бы связать пересекающиеся значения $ PhyloP для каждого из больших интервалов для последующего доступа для построения графиков.

Answer 1

Пакеты Bioconductor IRanges и GenomicRanges имеют findOverlaps, countOverlaps, coverage и другие интервальные функции, предназначенные для выполнения этих операций.Вы бы использовали функцию GRanges для представления каждого из объектов subject ("большие интервалы") и query ("меньшие интервалы") выше.См. инструкции по установке и затем виньетки, например, browseVignettes("GenomicRanges")

Более подробно, вот ваши данные

sdf <- read.table(textConnection(
"Chr  Start     End       Name         Val    Strand
chr7 145444998 146102295 CCDS5889.1   0      +
chr7 146102406 146167735 CCDS5889.1   0      +
chr7 146167929 146371931 CCDS5889.1   0      +"), header=TRUE)

qdf <- read.table(textConnection(
"Chr  Start     End       Name         Val    Strand PhyloP   
chr7 145444386 145444387 CCDS5889.1   0      +      0.684764
chr7 145444387 145444388 CCDS5889.1   0      +      0.684764
chr7 145444388 145444389 CCDS5889.1   0      +      0.684764
chr7 145444389 145444390 CCDS5889.1   0      +      0.684764"), header=TRUE)

и здесь мы конвертируем их в GRanges и найдите пересечение между запросом и темой

library(GenomicRanges)
subj <-
    with(sdf, GRanges(Chr, IRanges(Start, End), Strand, Val=Val))
query <-
    with(qdf, GRanges(Chr, IRanges(Start, End), Strand, Val=Val,
                      PhyloP=PhyloP, names=Name))
intersect(query, subj)

Ответ здесь не очень интересен

> intersect(query, subj)
GRanges with 0 ranges and 0 elementMetadata values
     seqnames ranges strand |

seqlengths
 chr7
   NA

Чтобы быть немного более полезным, вот запрос, который разбивает по вашимобщая область

tiles <- successiveIRanges(rep(100, 950), 900, 145444998)
query <- GRanges("chr7", tiles, "+")

Мы находим пересекающиеся диапазоны, находим, какой предмет перекрывается каждым пересекающимся диапазоном, и суммируем ширину перекрывающихся диапазонов

int <- intersect(query, subj)
tapply(int, subjectHits(findOverlaps(int, subj)),
       function(x) sum(width(x)))

, ведущую к

    1     2     3 
65800  6500 20400