r - 通过使用 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      +

200万行以上的小区间是这样的

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(from)和第 3(to)列。

情况类似

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

1. 我想知道每个较大间隔中有多少被较小间隔覆盖。
2. 此外，我会为每个大间隔关联相交的 $PhyloP 值，以便稍后访问以进行绘图。

Answer 1

Bioconductor IRanges 和 GenomicRanges 包具有 findOverlaps、countOverlaps、coverage 和其他旨在执行这些操作的基于区间的函数。您将使用 GRanges 函数来表示上面的每个 subject(“较大间隔”)和 query(“较小间隔”)对象.查看installation instructions然后是小插图，例如 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