R语言 重叠基因组范围

lhcgjxsq  于 2023-01-06  发布在  其他
关注(0)|答案(1)|浏览(171)

我有两个文件
编码

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

encode.total

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

我想比较两个文件,每个条目由chrstartstop组成。如果第一个文件中的start和stop值福尔斯同一染色体的第二个文件的start和stop之间,那么第一个文件中的start & stop值应该被第二个文件的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包.我的 Dataframe 很大,并且合并函数创建了一个非常大的 Dataframe (〉20亿行,这是不允许的),尽管我最终将 Dataframe 分成了一个较小的子集.但是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),]
djp7away

djp7away1#

使用BioconductorGenomicRanges软件包。

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("GenomicRanges")

假设有两个 Dataframe x0x1,就像原始示例中的encodeencode.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命令“手工”创建一个GRANGE示例,这里我们使用with(),这样我们就可以写GRanges(chr, IRanges(start=start, end=stop))来代替GRanges(x0$chr, IRanges(start=x0$start, end=x0$stop))
下一步是查找查询(gr 0)和主题(gr 1)之间的重叠

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

这将是快速上升到数百万的范围。

相关问题