是否有可能使用R data.table函数foverlaps来查找两个表中的重叠范围的交集？
我想使用foverlaps来查找两个bed文件的相交范围，并将包含重叠范围的所有行合并成一行。 在下面的例子中，我有两个基因组范围的表。 这些表被称为“床”文件，其具有基于零的起始坐标和染色体中的特征的一个基于终点的位置。 例如，START = 9，STOP = 20被解释为跨越基数10到20（含）。 这些床文件可以包含数百万行。 解决scheme需要给出相同的结果，无论提供两个文件相交的顺序如何。
第一个表
> table1 CHROMOSOME START STOP 1: 1 1 10 2: 1 20 50 3: 1 70 130 4: X 1 20 5: Y 5 200 第二个表
> table2 CHROMOSOME START STOP 1: 1 5 12 2: 1 15 55 3: 1 60 65 4: 1 100 110 5: 1 130 131 6: X 60 80 7: Y 1 15 8: Y 10 50
我在想，新的foverlaps函数可能是一个非常快速的方法来find这两个表中的相交范围，以产生一个如下的表格：
结果表：
> resultTable CHROMOSOME START STOP 1: 1 5 10 2: 1 20 50 3: 1 100 110 4: Y 5 50
这是可能的，还是有更好的方式来做data.table吗？
我还想首先确认在一张表中，对于任何给定的染色体，STOP坐标不与下一行的起始坐标重叠。 例如，染色体Y：1-15和染色体Y：10-50将需要被坍缩成染色体Y：1-50（参见第二表7和8）。 这不应该是这样，但function应该可能检查。 下面是一个关于潜在重叠如何崩溃的真实例子：
CHROM START STOP 1：1 721281 721619 2：1 721430 721906 3：1 721751 722042
期望的输出：
CHROM START STOP 1：1 721281 722042
创build示例表的函数如下所示：
table1 <- data.table( CHROMOSOME = as.character(c("1","1","1","X","Y")) , START = c(1,20,70,1,5) , STOP = c(10,50,130,20,200) ) table2 <- data.table( CHROMOSOME = as.character(c("1","1","1","1","1","X","Y","Y")) , START = c(5,15,60,100,130,60,1,10) , STOP = c(12,55,65,110,131,80,15,50) )
@Seth提供了使用data.table foverlaps函数解决交叉点重叠问题的最快方法。 但是，这个解决scheme没有考虑到input床文件可能有重叠范围需要被缩减到单个区域的事实。 @Martin Morgan用他的解决scheme使用了GenomicRanges软件包，解决了这个问题。 但是，Martin的解决scheme并没有使用foverlapsfunction。 @Arun指出，表中不同行的重叠范围当前不可能使用foverlaps。 感谢提供的答案，以及一些关于计算器的额外研究，我想出了这个混合解决scheme。
在每个文件中创build没有重叠区域的示例BED文件。
chr <- c(1:22,"X","Y","MT") #bedA contains 5 million rows bedA <- data.table( CHROM = as.vector(sapply(chr, function(x) rep(x,200000))), START = rep(as.integer(seq(1,200000000,1000)),25), STOP = rep(as.integer(seq(500,200000000,1000)),25), key = c("CHROM","START","STOP") ) #bedB contains 500 thousand rows bedB <- data.table( CHROM = as.vector(sapply(chr, function(x) rep(x,20000))), START = rep(as.integer(seq(200,200000000,10000)),25), STOP = rep(as.integer(seq(600,200000000,10000)),25), key = c("CHROM","START","STOP") )
现在创build一个新的床文件，其中包含bedA和bedB中的相交区域。
#This solution uses foverlaps system.time(tmpA <- intersectBedFiles.foverlaps(bedA,bedB)) user system elapsed 1.25 0.02 1.37 #This solution uses GenomicRanges system.time(tmpB <- intersectBedFiles.GR(bedA,bedB)) user system elapsed 12.95 0.06 13.04 identical(tmpA,tmpB) [1] TRUE
现在，修改bedA和bedB，使它们包含重叠的区域：
#Create overlapping ranges makeOverlaps <- as.integer(c(0,0,600,0,0,0,600,0,0,0)) bedC <- bedA[, STOP := STOP + makeOverlaps, by=CHROM] bedD <- bedB[, STOP := STOP + makeOverlaps, by=CHROM]
使用Foverlaps或GenomicRanges函数testing重叠范围的床文件的时间。
#This solution uses foverlaps to find the intersection and then run GenomicRanges on the result system.time(tmpC <- intersectBedFiles.foverlaps(bedC,bedD)) user system elapsed 1.83 0.05 1.89 #This solution uses GenomicRanges system.time(tmpD <- intersectBedFiles.GR(bedC,bedD)) user system elapsed 12.95 0.04 12.99 identical(tmpC,tmpD) [1] TRUE
获胜者： foverlaps ！
function使用
这是基于foverlaps的函数，如果存在重叠范围（使用rowShift函数进行检查），则只会调用GenomicRanges函数（reduceBed.GenomicRanges）。
intersectBedFiles.foverlaps <- function(bed1,bed2) { require(data.table) bedKey <- c("CHROM","START","STOP") if(nrow(bed1)>nrow(bed2)) { bed <- foverlaps(bed1, bed2, nomatch = 0) } else { bed <- foverlaps(bed2, bed1, nomatch = 0) } bed[, START := pmax(START, i.START)] bed[, STOP := pmin(STOP, i.STOP)] bed[, `:=`(i.START = NULL, i.STOP = NULL)] if(!identical(key(bed),bedKey)) setkeyv(bed,bedKey) if(any(bed[, STOP+1 >= rowShift(START), by=CHROM][,V1], na.rm = T)) { bed <- reduceBed.GenomicRanges(bed) } return(bed) } rowShift <- function(x, shiftLen = 1L) { #Note this function was described in this thread: #http://stackoverflow.com/questions/14689424/use-a-value-from-the-previous-row-in-an-r-data-table-calculation r <- (1L + shiftLen):(length(x) + shiftLen) r[r<1] <- NA return(x[r]) } reduceBed.GenomicRanges <- function(bed) { setnames(bed,colnames(bed),bedKey) if(!identical(key(bed),bedKey)) setkeyv(bed,bedKey) grBed <- makeGRangesFromDataFrame(bed, seqnames.field = "CHROM",start.field="START",end.field="STOP") grBed <- reduce(grBed) grBed <- data.table( CHROM=as.character(seqnames(grBed)), START=start(grBed), STOP=end(grBed), key = c("CHROM","START","STOP")) return(grBed) }
该function严格使用GenomicRanges软件包，产生相同的结果，但速度比foverlapsfunction慢大约10倍。
intersectBedFiles.GR <- function(bed1,bed2) { require(data.table) require(GenomicRanges) bed1 <- makeGRangesFromDataFrame(bed1, seqnames.field = "CHROM",start.field="START",end.field="STOP") bed2 <- makeGRangesFromDataFrame(bed2, seqnames.field = "CHROM",start.field="START",end.field="STOP") grMerge <- suppressWarnings(intersect(bed1,bed2)) resultTable <- data.table( CHROM=as.character(seqnames(grMerge)), START=start(grMerge), STOP=end(grMerge), key = c("CHROM","START","STOP")) return(resultTable) }
使用IRanges进行额外比较
我发现一个解决scheme使用IRanges折叠重叠区域，但比GenomicRanges慢10倍以上。
reduceBed.IRanges <- function(bed) { bed.tmp <- bed bed.tmp[,group := { ir <- IRanges(START, STOP); subjectHits(findOverlaps(ir, reduce(ir))) }, by=CHROM] bed.tmp <- bed.tmp[, list(CHROM=unique(CHROM), START=min(START), STOP=max(STOP)), by=list(group,CHROM)] setkeyv(bed.tmp,bedKey) bed[,group := NULL] return(bed.tmp[,-c(1:2),with=F]) } system.time(bedC.reduced <- reduceBed.GenomicRanges(bedC)) user system elapsed 10.86 0.01 10.89 system.time(bedD.reduced <- reduceBed.IRanges(bedC)) user system elapsed 137.12 0.14 137.58 identical(bedC.reduced,bedD.reduced) [1] TRUE
foverlaps()将会很好地完成。
首先设置两个表的键：
setkey(table1, CHROMOSOME, START, STOP) setkey(table2, CHROMOSOME, START, STOP)
现在，使用不匹配nomatch = 0 foverlaps()将不匹配的行放入table2 。
resultTable <- foverlaps(table1, table2, nomatch = 0)
接下来为START和STOPselect适当的值，并删除多余的列。
resultTable[, START := pmax(START, i.START)] resultTable[, STOP := pmin(STOP, i.STOP)] resultTable[, `:=`(i.START = NULL, i.STOP = NULL)]
重叠的STOP到未来的START应该是一个不同的问题。 这实际上是我的一个，所以当我有一个很好的答案的时候，也许我会问这个问题，然后再回来。
如果你没有卡在data.table解决scheme， GenomicRanges
source("http://bioconductor.org/biocLite.R") biocLite("GenomicRanges")
给
> library(GenomicRanges) > intersect(makeGRangesFromDataFrame(table1), makeGRangesFromDataFrame(table2)) GRanges object with 5 ranges and 0 metadata columns: seqnames ranges strand <Rle> <IRanges> <Rle> [1] 1 [ 5, 10] * [2] 1 [ 20, 50] * [3] 1 [100, 110] * [4] 1 [130, 130] * [5] Y [ 5, 50] * ------- seqinfo: 3 sequences from an unspecified genome; no seqlengths
在基因组学中最重叠的范围问题中，我们有一个大数据集x （通常是测序读数）和另一个较小的数据集y （通常是基因模型，外显子，内含子等）。 我们的任务是找出x哪些区间与y哪个区间重叠，或者x多less个区间在每个y区间重叠。
在foverlaps() ，我们不必在更大的数据集x上设置setkey() ，这是一个相当昂贵的操作。 但是，你需要有一个关键的设置。 对于你的情况，从这个例子看来， table2是大= x ， table1 = y 。
require(data.table) setkey(table1) # key columns = chr, start, end ans = foverlaps(table2, table1, type="any", nomatch=0L) ans[, `:=`(i.START = pmax(START, i.START), i.STOP = pmin(STOP, i.STOP))] ans = ans[, .(i.START[1L], i.STOP[.N]), by=.(CHROMOSOME, START, STOP)] # CHROMOSOME START STOP V1 V2 # 1: 1 1 10 5 10 # 2: 1 20 50 20 50 # 3: 1 70 130 100 130 # 4: Y 5 200 5 50
但是我同意能够一步完成这件事情会很棒。 不知道现在怎么样，但也许使用额外的值reduce和mult=参数intersect 。
