在NGS科研领域,做ChIP-seq/CLIP-seq等研究蛋白与DNA/RNA结合规律的时候,经常会用到peak calling的算法。这个方法会在全基因组/转录组范围内找DNA结合位点,一般先通过确定测序数据的depth peak,然后用case vs control样本,看depth peak的改变的倍数来确定正真的peak的分布。
假如我们有一组数据,我们画它的分布曲线如下
aa=100:1
bb=sin(aa/3)
cc=aa*bb
plot(cc, type="l")
我们想找到那些峰的位置,那么我们可用R语言这样来实现
aa=100:1
bb=sin(aa/3)
cc=aa*bb
plot(cc, type="l")
find_peaks <- function (x, m = 3){
shape <- diff(sign(diff(x, na.pad = FALSE)))
pks <- sapply(which(shape < 0), FUN = function(i){
z <- i - m + 1
z <- ifelse(z > 0, z, 1)
w <- i + m + 1
w <- ifelse(w < length(x), w, length(x))
if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])) return(i + 1) else return(numeric(0))
})
pks <- unlist(pks)
pks
}
abline(v=find_peaks(cc))
生信的临床应用领域我们经常需要看血浆游离DNA的长度分布(胎儿游离DNA或循环肿瘤游离DNA),那么我们也可用这个函数来看在哪些位置会成峰。