重叠的基因组范围

Posted 2023-04-19

【中文标题】重叠的基因组范围

【英文标题】：Overlapping Genomic Ranges

【发布时间】：2013-09-30 19:06:07

【问题描述】：

我有两个文件

编码

X.pattern.name       chr       start    stop    strand  score   p.value q.value matched.sequence
1   V_CETS1P54_01   chr1    98769545    98769554    +   11.42280    8.89e-05    NA  TCAGGATGTA
2   V_CETS1P54_01   chr1    152013037   152013046   +   11.98020    4.74e-05    NA  ACAGGAAGTT
3   V_CETS1P54_01   chr1    168932563   168932572   +   11.60860    7.59e-05    NA  ACCGGATGCT

编码.总数

    chr     start           stop
1   chr1    58708485    58708713
2   chr1    58709084    58710538
3   chr1    98766295    98766639
4   chr1    98766902    98770338
5   chr1    107885506   107889414
6   chr1    138589531   138590856
7   chr1    138591180   138593378
8   chr1    152011746   152013185
9   chr1    152014263   152014695
10  chr1    168930561   168933076
11  chr1    181808064   181808906
12  chr1    184609002   184611519
13  chr1    193288453   193289567
14  chr1    193290105   193290490
15  chr1    193290744   193291092
16  chr1    196801920   196804954

我想比较两个文件，每个条目按 chr、start 和 stop。如果第一个文件中的起始值和终止值介于同一染色体的第二个文件的起始值和终止值之间，则第一个文件中的起始值和终止值应替换为第二个文件的起始值和终止值。我为此目的编写了一个 for 循环，但它花费的时间太长。有哪些替代方案？

代码：

for(i in 1:nrow(encode))

     for(j in 1:nrow(encode.total))
     
           if(encode[i,2]==encode.total[j,1])
           
               if((encode[i,3]>=encode.total[j,2]) & (encode[i,4]<=encode.total[j,3]))
               
                   encode[i,3]=encode.total[j,2]
                   encode[i,4]=encode.total[j,3]
               
           
     

出于同样的目的，我还尝试了如下所示的 GenomicRanges 包。我的数据帧的大小很大，它们上的合并功能会创建一个非常大的数据帧（不允许超过 20 亿行），尽管我最终将数据帧子集化为一个较小的数据帧。但是 merge() 会占用大量内存并终止 R。

gr1<-GRanges(seqnames=encode$chr,IRanges(start=encode$start,end=encode$end))
gr2<-GRanges(seqnames=encode.total$chr, IRanges(start=encode.total$start,end=encode.total$end))

ranges <- merge(as.data.frame(gr1),as.data.frame(gr2),by="seqnames",suffixes=c("A","B"))
ranges <- ranges[with(ranges, startB <= startA & endB >= endA),]

【问题讨论】：

在您的示例中，两个文件中的 start 和 stop 值是相同的，所以什么都不会改变...

我已经编辑了我的文件，请立即查看。

重复？ ***.com/questions/3916195/…

@Henrik 您指向的问题适用于 Ranges 对象，我有一个数据框对象。我不熟悉 R 中的 Ranges 对象。

问题中的 OP 以“两个 data.frames 每个包含三列：chrom、start & stop”开头。在我看来，它与您的数据非常相似。

【参考方案1】：

使用Bioconductor GenomicRanges 包。

source("http://bioconductor.org/biocLite.R")
biocLite("GenomicRanges")

假设有两个数据框x0 和x1，如原始示例中的encode 和encode.total。我们想把这些变成一个 Granges 实例。我这样做了

library(GenomicRanges)
gr0 = with(x0, GRanges(chr, IRanges(start=start, end=stop)))
gr1 = with(x1, GRanges(chr, IRanges(start=start, end=stop)))

通常可以简单地说makeGRangesFromDataFrame(x0)，或使用标准R 命令“手动”创建一个GRanges 实例。这里我们使用with()，这样我们就可以写GRanges(chr, IRanges(start=start, end=stop))而不是GRanges(x0$chr, IRanges(start=x0$start, end=x0$stop))。

下一步是查找查询（gr0）和主题（gr1）之间的重叠

hits = findOverlaps(gr0, gr1)

导致

> hits
Hits of length 3
queryLength: 3
subjectLength: 16
  queryHits subjectHits 
   <integer>   <integer> 
 1         1           4 
 2         2           8 
 3         3          10 

然后更新了相关的开始/结束坐标

ranges(gr0)[queryHits(hits)] = ranges(gr1)[subjectHits(hits)]

给予

> gr0
GRanges with 3 ranges and 0 metadata columns:
      seqnames                 ranges strand
         <Rle>              <IRanges>  <Rle>
  [1]     chr1 [ 98766902,  98770338]      *
  [2]     chr1 [152011746, 152013185]      *
  [3]     chr1 [168930561, 168933076]      *
  ---
  seqlengths:
   chr1
     NA

这将很快达到数百万个范围。

【讨论】：

这就像一个魅力！非常有用且非常高效

with(x0... 和 with(x1... 的意义何在。 x0 和 x1 到底是什么。

@by0 我在答案中添加了尝试解释 x0 和 x1（包含基因组范围信息的两个数据框）和 with()（一个基本 R 函数，可以更轻松地引用数据中的列.frame 是它的第一个参数）。

@MartinMorgan 你能回答以下问题吗？谢谢。 ***.com/questions/48187739/…