r - R中的常见基因组区间

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

我的输入:

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

干杯!

最佳答案

如果您的数据在 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)

关于r - R中的常见基因组区间，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/23084322/