3

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

我的输入：

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）。

干杯!

4

3 回答 3

6

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

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

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

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

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

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)
于 2014-04-15T13:35:57.693 回答
3

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

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 行），然后找到每对剩余间隔之间的交点。然后，您可以根据需要将结果制成表格。

于 2014-04-15T13:30:31.077 回答
1

这肯定很长（考虑到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
于 2014-04-15T13:39:59.177 回答