1

我尝试使用 findOverlap 来解决这个问题，但我只找到没有条件的重叠区域，所以如果我有一些条件来选择数据。我应该怎么做？

假设我有两个数据框，如下所示

数据框

Sample, start, stop, event, probe, length, length/probe, region
CNV1234,  2000,  3000,  CN gain,  23,  235, 9, intron
CNV1534,  1200,  1800,  CN loss,  60,  600  10, exon

数据框 b

Sample, start, stop, event, probe, length, length/probe, region
CNV234,  2500,  3500,  CN gain,  23,  235, 9, exon 
CNV34,  1200,  1800,  CN loss,  60,  600  10, intron

我有两个问题

首先，我想找到这两个数据帧之间的重叠，其中 CNV 的长度重叠超过 50%，并且这种重叠位于内含子区域

其次，我想知道重叠区域的长度。

我希望我的结果有一个像这样的数据框

Sample, start, stop, event, probe, length, length/probe, region, overlap, length of overlap
CNV1234,  2000,  3000,  CN gain,  23,  235, 9, intron, CNV234, 500
4

1 回答 1

1

这是你的数据

a <- read.csv(textConnection(
    "Sample, start, stop, event, probe, length, length/probe, region
     CNV1234,  2000,  3000,  CN gain,  23,  235, 9, intron
     CNV1534,  1200,  1800,  CN loss,  60,  600  10, exon"))

b <- read.csv(textConnection(
    "Sample, start, stop, event, probe, length, length/probe, region
     CNV234,  2500,  3500,  CN gain,  23,  235, 9, exon 
     CNV34,  1200,  1800,  CN loss,  60,  600  10, intron"))

加载GenomicRanges包（我假设您的数据实际上来自多个染色体，并且您希望对所有染色体都这样做；“A”是染色体名称）

library(GenomicRanges)
gr1 <- with(a, GRanges("A", IRanges(start, stop - 1L),
                       Sample=Sample, event=event))
gr2 <- with(b, GRanges("A", IRanges(start, stop - 1L, names=Sample),
                       Region=Sample))

请注意 GRanges 如何表示范围（从 1 开始，包括开始和结束坐标）。查找这些对象之间的所有重叠（您可以使用 min.overlaps 排除一些重叠，例如，那些小于最小宽度 1/2 的重叠）

h <- findOverlaps(gr1, gr2)

目前尚不清楚宽度的“50%”是什么——a 的宽度？乙？-- 所以我计算所有重叠的宽度，然后保留那些宽度大于“a”宽度的1/2

wd <- width(pintersect(gr1[queryHits(h)], gr2[subjectHits(h)]))
ok <- wd > width(gr1[queryHits(h)]) / 2
h <- h[ok]

最后，我通过选择满足重叠标准的查询并添加元数据列 ( mcols) 和它们重叠区域的重叠宽度来组装结果

result <- gr1[queryHits(h)]
mcols(result) <- cbind(mcols(result), mcols(gr2[subjectHits(h)]))
result$`width of overlap` <- wd[ok]

结果可能会被强制返回到带有 的数据框as.data.frame(result)，或者您的下游分析可能是使用 Granges 基础架构自然完成的？

最好在 Bioconductor邮件列表中询问有关Bioconductor软件包的问题（无需订阅）。可能有更有效的方法可以做到这一点，该邮件列表中的人员将提供这些解决方案。

于 2013-08-13T06:12:29.057 回答