r - 通过使用 R 将较小的基因组间隔数据与较大的基因组间隔相交来覆盖

我想在 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(从)列和第 3(到)列中。

情况类似于

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

我想知道每个较大的间隔中有多少被较小的间隔覆盖。

此外，我会将每个大间隔的相交 $PhyloP 值关联起来，以便以后访问以进行绘图。

最佳答案

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

关于r - 通过使用 R 将较小的基因组间隔数据与较大的基因组间隔相交来覆盖，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/5200946/