上周的推文[计算基因组范围内邻近DSB的间距]介绍了邻近DSB
间距的计算方法。虽然结果没有问题,但总觉得哪里有点怪怪的。思来想去,是不是应该去除掉不关心的噪音,即区域间的gap
。例如,下图中DSB1
与DSB2
互为最邻近,DSB3
、DSB4
的最邻近分别为DSB2
、DSB5
,DSB5
与DSB6
互为最邻近。探索关心的肯定是最邻近DSB
之间的距离,而类似DSB3
与DSB4
之间的距离即为gap
,这样的gap
并不需要关心,留下来就属于噪音,为了更加精准的探索距离情况应该想办法去除掉。
这个问题整体诉求并不难,为了解决该问题第一时间还是想到了借助distanceToNearest
函数。先看一下这个函数的用法:
suppressPackageStartupMessages(library(GenomicRanges))
query <- GRanges(c("chr1", "chr2"), IRanges(c(1, 5), width=1))
subject <- GRanges("chr1", IRanges(c(6, 5, 13), c(10, 10, 15)))
d <- distanceToNearest(query, subject)
d
Hits object with 1 hit and 1 metadata column:
queryHits subjectHits | distance
<integer> <integer> | <integer>
[1] 1 2 | 3
-------
queryLength: 2 / subjectLength: 3
可以看到很好用,只要输入查询和搜索集即可,这样查询集里面的每一个区域都会在搜索集找到一个最近的区域。函数返回的结果是一个Hits
对象,里面记录了哪两个区域最靠近以及距离。有了这个函数,便可以很方便地计算DSB
与其最近DSB
之间的距离。接下来,问题就变成了如何去计算每一个DSB
最邻近的距离。是不是每个DSB
循环一遍就可以了呢?问题看似简单,但有些细节还是需要留心。毕竟,不是程序跑起来有结果就是对的。
肯定是要利用循环,每次取出一个作为查询区域,数据集剩余的区域作为搜索集,但与此同时会产生一个重复计算的问题。如第一次以DSB1
作为查询时,最邻近区域为DSB2
,第二次以DSB2
查询时,最邻近区域为DSB1
,如此一来这两者的距离就计算两次了。这样肯定是不合理的,计算的时候需要避免。
这里,为避免重复计算,本人采用的策略是构造一个新的id
,这样互为最邻近的区域就会形成一个id
组合,最后根据这个组合来去除重复。下面来看看具体的实现代码,为方便使用定义成一个函数:
calcdist <- function(gr, sname){
chr <- table(gr@seqnames) > 1
gr <- gr[gr@seqnames %in% names(chr[chr!=F])]
num <- length(gr)
gr$id <- 1:num
out <- data.frame()
for(idx in 1:num){
sgr <- gr[-idx]
hit <- distanceToNearest(gr[idx], sgr)
idv <- paste0(sort(c(gr[idx]$id, sgr[subjectHits(hit)]$id)), collapse='_')
dis <- mcols(hit)$distance
tmp <- data.frame(chr=as.vector(gr[idx]@seqnames), dist=dis, sample=sname, id=idv, stringsAsFactors=F)
out <- rbind(out, tmp)
}
out <- subset(out[!duplicated(out$id),], select=-id)
return(out)
}
函数的魅力在于便于重复使用,以后再有这样的需求调用一下函数就可以了。下面用样本测试一下:
gr <- ChIPpeakAnno::toGRanges('sample_peaks.bed')
gr
GRanges object with 8187 ranges and 1 metadata column:
seqnames ranges strand | score
<Rle> <IRanges> <Rle> | <numeric>
MACS_peak_1 GL456210.1 [ 21815, 21815] * | 424
MACS_peak_2 GL456211.1 [ 97142, 97142] * | 330
MACS_peak_3 GL456211.1 [174111, 174111] * | 146
MACS_peak_4 GL456212.1 [ 85136, 85136] * | 827
MACS_peak_7 GL456221.1 [128963, 128963] * | 491
... ... ... ... . ...
MACS_peak_8183 chrY [ 1404762, 1404762] * | 457
MACS_peak_8184 chrY [ 4036698, 4036698] * | 8868
MACS_peak_8185 chrY [11799470, 11799470] * | 623
MACS_peak_8186 chrY [16484365, 16484365] * | 611
MACS_peak_8187 chrY [90763211, 90763211] * | 891
system.time(out <- calcdist(gr, 'ctrl'))
user system elapsed
432.552 0.135 433.711
head(out)
chr dist sample
1 chr1 1484377 ctrl
2 chr1 201371 ctrl
3 chr1 50155 ctrl
4 chr1 66122 ctrl
5 chr1 59074 ctrl
虽然函数可以正常运行结果也没有问题,可是这个运行时间着实没有料到,8千多行循环一遍居然用了七八分钟,只能说还是尽量不要用R写for
循环的好,不然数据多了这耗时真的很严重。不过,话说回来,如果不着急的话这样也完全可以接受。
时间和效率有时候还是很重要的,如果想要节约时间,很快拿到结果,就得想想办法加速了,最简单快捷的方法自然是多线程并发。可以知道计算过程是染色体独立的,这样可以利用多线程同时计算多个染色体,如此便可以大大提高计算速度。下面用多线程来计算一遍:
suppressPackageStartupMessages(library(foreach))
suppressPackageStartupMessages(library(doParallel))
chr <- table(gr@seqnames) > 1
chr <- names(chr[chr!=F])
cl <- makeCluster(10)
registerDoParallel(cl)
system.time(
out2 <- foreach(i = chr, .combine = rbind) %dopar% {
suppressPackageStartupMessages(library(GenomicRanges))
calcdist(gr[gr@seqnames == i], 'ctrl')}
)
user system elapsed
0.073 0.004 44.498
stopCluster(cl)
head(out2)
chr dist sample
1 chr1 1484377 ctrl
2 chr1 201371 ctrl
3 chr1 50155 ctrl
4 chr1 66122 ctrl
5 chr1 59074 ctrl
可以看到相同的数据,用10个线程同时计算,全部时间消耗也只有半分钟左右,这速度可以说是所见即所得了。