R中常见的基因组区间答案

【问题标题】：Common genomic intervals in RR中常见的基因组区间
【发布时间】：2014-04-15 12:53:08
【问题描述】：

我想推断不同样本之间共享的基因组间隔。

我的意见：

sample    chr start end
NE001      1   100  200
NE001      2   100  200
NE002      1   50   150
NE002      2   50   150
NE003      2   250  300

我的预期输出：

chr start end  freq
1    100  150   2
2    100  150   2

其中“频率”是有多少样本有助于推断共享区域。在上面的例子中，freq = 2（NE001 和 NE002）。

干杯！

【问题讨论】：

标签： r overlap overlapping genome

【解决方案1】：

如果您的数据在 data.frame 中（见下文），我使用 Bioconductor GenomicRanges 包创建一个 Granges 实例，同时保留非范围列

library(GenomicRanges)
gr <- makeGRangesFromDataFrame(df, TRUE)

数据所代表的离散范围由disjoin函数给出，不相交范围（'query'）和你原来的（'subject'）之间的重叠是

d <- disjoin(gr)
olaps <- findOverlaps(d, gr)

将与每个重叠主题关联的样本信息与对应的查询分开，并将其与不相交的GRanges关联为

mcols(d) <- splitAsList(gr$sample[subjectHits(olaps)], queryHits(olaps))

导致例如

> d[elementLengths(d$value) > 1]
GRanges with 2 ranges and 1 metadata column:
      seqnames     ranges strand |           value
         <Rle>  <IRanges>  <Rle> | <CharacterList>
  [1]        1 [100, 150]      * |     NE001,NE002
  [2]        2 [100, 150]      * |     NE001,NE002
  ---
  seqlengths:
    1  2
   NA NA

我是这样输入您的数据的：

txt <- "sample    chr start end
NE001      1   100  200
NE001      2   100  200
NE002      1   50   150
NE002      2   50   150
NE003      2   250  300"
df <- read.table(textConnection(txt), header=TRUE, stringsAsFactors=FALSE)

【讨论】：

这是一种比我更好的方法（忘记了disjoin 函数）——并且可以正确处理高阶交叉点。

【解决方案2】：

鉴于这个问题背后的背景，我怀疑你值得学习 Bioconductor 的 GenomicRanges 软件包。

library(GenomicRanges)
gr <- GRanges(seqnames=df$chr, ranges=IRanges(start=df$start, end=df$end))
ov <- findOverlaps(gr,gr, type="any")
ov <- ov[queryHits(ov) != subjectHits(ov)]
between <- pintersect(gr[subjectHits(ov)], gr[queryHits(ov)])

方法是：找到所有自重叠，删除与自身进行比较的区间（第 4 行），然后找到每对剩余区间之间的交点。然后，您可以根据需要将结果制成表格。

【讨论】：

这将适用于引用的案例，但如果您希望以“嵌套”方式处理超过 2 个区域之间的交叉点，那么@martin 的方法会更好。

【解决方案3】：

这肯定很长（考虑到 expand.grid.df，在大型 data.frames 上可能效率很低，但是，我希望它能给你一个起点。需要注意的是，我没有基因组学背景（我'肯定会通过）所以不知道通用的包。当然这些是最好的方法。我只是觉得尝试解决方案会很有趣。

s<-"sample    chr start end
NE001      1   100  200
NE001      2   100  200
NE002      1   50   150
NE002      2   50   150
NE003      2   250  300"

dat<-read.table(text=s, header=T)

library(plyr)
between<-function(x,y,z) x<=y & y<=z
dat$id<-seq_along(dat[,1])
expand.grid.df <- function(...) Reduce(function(...) merge(..., by=NULL), list(...))
expdat<-ddply(dat, .(chr), function(x) expand.grid.df(x,x))
expdat<-subset(expdat, id.x!=id.y)
expdat$betweenL<-with(expdat, between(start.y, start.x, end.y))
expdat$betweenR<-with(expdat, between(start.x, start.y, end.x))
expdat<-subset(expdat, betweenL | betweenR)
expdat$commonstart<-with(expdat,ifelse(betweenL,start.x,start.y))
expdat$commonend<-with(expdat, ifelse(betweenL, end.y, end.x))
res<-ddply(expdat, .(chr, commonstart, commonend),summarize, freq=length(sample.x))
> res
  chr commonstart commonend freq
1   1         100       150    2
2   2         100       150    2

【讨论】：