上周的推文[计算基因组范围内邻近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个线程同时计算，全部时间消耗也只有半分钟左右，这速度可以说是所见即所得了。
