找到点到范围的重叠

Question

我有一个数据帧df1:

df1 <- read.table(text=" Chr06  79641   
Chr06   82862   
Chr06   387314  
Chr06   656098  
Chr06   678491  
Chr06   1018696", header=FALSE, stringsAsFactors=FALSE)

我想检查df1中的每一行是否包含在df2的范围内.df2中的column2是范围的开头,column3是范围的结尾.范围之间没有重叠(行之间).df2中的数据按Column1和column2排序.我为此写了一个循环,但我对此并不满意,因为如果我在df1中有几千行,它会运行很长时间.我想找到一种更有效的方法来完成这项工作(更好的没有循环).谢谢.df2数据框:

df2 <- read.table(text=" Chr05  799 870
Chr06   77914   77942
Chr06   78233   78269
Chr06   78719   78836
Chr06   79720   87043
Chr06   87223   87305
Chr06   380020  380060
Chr06   387314  387371
Chr06   654907  654988
Chr06   657929  658057
Chr06   677198  677229
Chr06   679555  680170
Chr06   1015425 1015475
Chr06   1018676 1018736
Chr06   1020564 1020592", header=FALSE, stringsAsFactors=FALSE)

我的剧本:

df1$V3 <- FALSE
for (i in 1:dim(df1)[1]) {
  for (j in 1:dim(df2)[1]) {
    if ((df1[i,1] == df2[j,1]) && (df1[i,2]>=df2[j,2]) 
        && (df1[i,2]<=df2[j,3])) {
      df1[i,3]<-TRUE
      break;
    }
  }
}
df1

预期结果显示为df1.

Answer 1

使用GenomicRanges:

#Convert to Granges objects
gr1 <- GRanges(seqnames = df1$V1,
               ranges = IRanges(df1$V2, df1$V2))

gr2 <- GRanges(seqnames = df2$V1,
               ranges = IRanges(df2$V2, df2$V3))
#Subset gr1
subsetByOverlaps(gr1, gr2)

# GRanges object with 3 ranges and 0 metadata columns:
#       seqnames             ranges strand
#          <Rle>          <IRanges>  <Rle>
#  [1]    Chr06 [  82862,   82862]      *
#  [2]    Chr06 [ 387314,  387314]      *
#  [3]    Chr06 [1018696, 1018696]      *
#   -------
#   seqinfo: 1 sequence from an unspecified genome; no seqlengths

#Or we can use merge
mergeByOverlaps(gr1, gr2)

# DataFrame with 3 rows and 2 columns
#                          gr1                        gr2
#                    <GRanges>                  <GRanges>
# 1 Chr06:*:[  82862,   82862] Chr06:*:[  79720,   87043]
# 2 Chr06:*:[ 387314,  387314] Chr06:*:[ 387314,  387371]
# 3 Chr06:*:[1018696, 1018696] Chr06:*:[1018676, 1018736]

另外,看看bedtools:

总的来说,bedtools实用程序的工具,用于基因组学分析任务的宽范围内的瑞士军刀.最广泛使用的工具使基因组算法成为可能:即基因组上的集合论.例如,bedtools允许人们在广泛使用的基因组文件格式(例如BAM,BED,GFF/GTF,VCF)中交叉,合并,计数,补充和混洗来自多个文件的基因组间隔.虽然每个单独的工具被设计为执行相对简单的任务(例如,交叉两个间隔文件),但是可以通过在UNIX命令行上组合多个bedtools操作来进行非常复杂的分析.

Answer 2

以下是一种data.table解决方案GenomicRanges:

library(data.table)
dt1 <- data.table(df1)[, V3 := V2]
dt2 <- data.table(df2, key = c("V2", "V3"))
foverlaps(dt1, dt2)[V1 == i.V1][, -c(4, 6), with = F]
#       V1      V2      V3    i.V3
# 1: Chr06   79720   87043   82862
# 2: Chr06  387314  387371  387314
# 3: Chr06 1018676 1018736 1018696