所附图(曼哈顿图)在基因组的x轴染色体上包含位置，在Y轴上包含-log(p)，其中p是与该特定位置上的点(变量)相关的p值。
我已经使用以下R代码(从gap包中)生成它:
require(gap)
affy <-c(40220, 41400, 33801, 32334, 32056, 31470, 25835, 27457, 22864, 28501, 26273,
24954, 19188, 15721, 14356, 15309, 11281, 14881, 6399, 12400, 7125, 6207)
CM <- cumsum(affy)
n.markers <- sum(affy)
n.chr <- length(affy)
test <- data.frame(chr=rep(1:n.chr,affy),pos=1:n.markers,p=runif(n.markers))
oldpar <- par()
par(cex=0.6)
colors <- c("red","blue","green","cyan","yellow","gray","magenta","red","blue","green", "cyan","yellow","gray","magenta","red","blue","green","cyan","yellow","gray","magenta","red")
mhtplot(test,control=mht.control(colors=colors),pch=19,bg=colors)
> head(test)
chr pos p
1 1 1 0.79296584
2 1 2 0.96675136
3 1 3 0.43870076
4 1 4 0.79825513
5 1 5 0.87554143
6 1 6 0.01207523
我对使图的峰坐标高于某个阈值(-log(p))感兴趣。
最佳答案
如果您希望值的索引高于第99个百分位，请执行以下操作:
# Add new column with log values
test = transform(test, log_p = -log10(test[["p"]]))
# Get the 99th percentile
pct99 = quantile(test[["log_p"]], 0.99)
...并从原始数据
test中获取值:peaks = test[test[["log_p"]] > pct99,]
> head(peaks)
chr pos p log_p
5 1 5 0.002798126 2.553133
135 1 135 0.003077302 2.511830
211 1 211 0.003174833 2.498279
586 1 586 0.005766859 2.239061
598 1 598 0.008864987 2.052322
790 1 790 0.001284629 2.891222
您可以将其与任何阈值一起使用。请注意，我尚未计算出一阶导数，请参见以下问题以获取一些指针:
How to calculate first derivative of time series
计算完一阶导数后，可以通过查看时间序列中一阶导数(几乎)为零的点来找到峰值。识别出这些峰值之后，您可以检查哪些峰值在阈值之上。
关于r - 曼哈顿图中的峰值检测，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/15068445/