用R语言找曲线的峰的位置
在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），那么我们也可用这个函数来看在哪些位置会成峰。
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2019-10-29，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读