最近客户想验证一下基因组内邻近DSB之间的距离分布在实验组和对照组是否有差别。不熟悉DSB的朋友可以先看一看这篇推文：END-seq数据分析。既然提了需求，那就用数据验证一下呗。
其实，这个任务并不难，无非就是将基因组内邻近的DSB位点都相互算一下距离，注意这里说的距离是线性距离。现在明确了问题，下面就是具体该怎么做的事了。当然，首先想到的是利用R里面现有的包来做这件事，比如利用IRanges包的distanceToNearest函数。
其实，这个问题不用现有包，只用原始R代码来做也不难，下面咱们利用原始代码来解决。首先，确定一下距离的计算方式，如下图所示：1、前一个DSB结束位置到后一个DSB的起点；2、两个DSB位点中心之间的距离，或者submit位点之间的距离。
两种方式计算过程大体相似，第二种更为简单一些。这里就展示一下稍微麻烦一点的计算过程，第一种方式：
peak <- read.table('sample_peaks.bed', header=F, stringsAsFactors=F, sep='\t')
peak[11:16,]
V1 V2 V3 V4 V5
11 chr1 3118979 3132796 MACS_peak_14 1405.41
12 chr1 4603853 4617340 MACS_peak_15 3100.00
13 chr1 4801376 4815127 MACS_peak_16 272.22
14 chr1 4851617 4864724 MACS_peak_17 483.81
15 chr1 4917189 4930952 MACS_peak_18 3100.00
16 chr1 5019922 5033338 MACS_peak_19 3100.00
peakstart <- peak[,c('V1', 'V2')]
colnames(peakstart) <- c('chr', 'pos')
peakend <- peak[,c('V1', 'V3')]
colnames(peakend) <- c('chr', 'pos')
peak <- rbind(peakstart, peakend)
peak <- peak[order(peak$chr, peak$pos),]
chrname <- unique(peak$chr)
out <- data.frame()
for(chr in chrname){
chrpeak <- peak[peak$chr == chr, 'pos']
if(length(chrpeak)<4){
next
}
tmp <- data.frame(chr=chr, dist=diff(chrpeak), sample='sample1', stringsAsFactors=F)
out <- rbind(out, tmp)
}
out <- out[seq(2, nrow(out), by=2), ]
head(out)
chr dist sample
2 chr1 1471057 sample1
4 chr1 184036 sample1
6 chr1 36490 sample1
8 chr1 52465 sample1
10 chr1 88970 sample1
12 chr1 45880 sample1
上面的计算过程有几个条件：一，DSB位点之间没有重叠，这个条件基本没有问题，使用macs2这样的peakcalling软件得到的peak之间不会有重叠；二，按染色体位置排序；三，各个染色体单独计算。
计算过程的原理也很简单，将染色体上DSB位点一字排开，依次用后面的位置坐标减去前面的，diff可以直接计算向量相邻两个元素的差值，然后去除DSB本身的宽度数据即可得到邻近DSB的间距。此时，咱们可将所有样本都循环计算一遍，得到数据后便可以比较样本间是否有差别了。
网友评论