5

添付のプロット (マンハッタン プロット) には、x 軸にゲノムからの染色体位置、Y 軸に -log(p) が含まれます。p は、その特定の位置からのポイント (バリアント) に関連付けられた p 値です。 ここに画像の説明を入力

次のRコードを使用して生成しました（ギャップパッケージから）：

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)) を超えるプロットのピークの座標を取得することに興味があります。

4

2 に答える 2

1

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

これは、任意のしきい値で使用できます。一次導関数を計算していないことに注意してください。いくつかのポインタについては、この質問を参照してください。

時系列の一次導関数の計算方法

一次導関数を計算した後、一次導関数が (ほぼ) ゼロである時系列のポイントを調べることで、ピークを見つけることができます。これらのピークを特定したら、どのピークがしきい値を超えているかを確認できます。

于 2013-02-25T13:52:52.317 に答える
0

グラフをプロットした後の私の経験に基づいて、次のRコードを使用してピーク座標を見つけることができます

plot(x[,1], x[,2])

Identify(x[,1], x[,2], labels=row.names(x))

ここで x[,1] は x 座標を指します (ゲノム座標と x[,2] は #your -log10P 値になります)

この時点で、マウスのポイントを使用してポイントを選択し、Enter キーを押します。これにより、ピークの位置が表示されます。次のコードを入力して、#座標を取得します。

座標 <- ロケーター(type="l")

座標

于 2016-09-30T22:38:50.943 に答える