Q-value(R)

FDR制御法の一つです。FDR制御法でよく出てくる、p値のアナロジーとして使われている小文字のq-valueとは違って、これはQ−valueという手法名(もしくはソフトウェア名)なのでご注意を。

BH法は、p値が一様分布している事を前提にFDRの期待値を計算しますが、実際の生物のデータは一様分布しているとは限りません。例えば、GEOに登録されているGSE22531でのSEW-2871を投与されたマウス(処置群)と、何も処置をしていないマウス(対照群)とでt検定を行った時の、p値の分布を見てみます。
このように実際のp値の分布は何故か0側に寄る事が多いです。BH法では(http://www.med.osaka-u.ac.jp/pub/kid/clinicaljournalclub1.html)の図にあるように、閾値をずらすと一定の割合で有意とするp値が検出するようになっていると私は理解していますが、上の図のような0側か1側に偏ったデータでBH法をかけると、閾値を大きい値にしないとまったく発現変動遺伝子が検出されなかったり、反対に小さな閾値で大量にとれすぎてしまうという問題があります。これはBH法が前提としている一様分布と実際の分布のズレが原因だと考えられます。例えばこの場合q値の分布は以下のようになり、閾値を0.2くらいまでひきあげないと何もとれてきません。


q<0.2で8遺伝子
[1] 0.1647064 0.1647064 0.1647064 0.1647064
[5] 0.1647064 0.1647064 0.1647064 0.1647064

論文を読んでいると、q値はp値での閾値よりも大きくても許されている印象がありますが(0.1〜0.3とかでも)、あまりに大きいq値だとやはりおかしいです。例えばうちのデータはもっとひどくて閾値を0.5以上に設定しないと何も検出されません。でもこれでは判定した発現変動遺伝子のうち半分以上は間違いですと言っているわけなので、とても使い物になりません。

この実際の分布が一様分布しない現象を、以下のように帰無仮説の分布と対立仮説の分布とが別々に分布しているからだと仮定したものがQ-valueです。というか最近の手法はみんなそのような仮定をしているっぽいです。

帰無仮説:発現変動していない遺伝子。平均値近辺で分布。
対立仮説:発現変動遺伝子。平均値から離れて分布。

帰無仮説の分布は、平均値近辺に分布しており、この分布から求められたp値は一様分布します。
一方、対立仮説の分布はこの帰無仮説の分布から離れたところで分布しているため、p値は小さくなり、0側に密集していると考えられます。
<イメージ>
この場合、実データの分布は帰無仮説の分布と対立仮説の分布の混合分布とされるので、0側に寄り気味な一様分布っぽい分布が観測されるのだと考えられます。

この帰無仮説と対立仮説の密度比が各々
π0:1 - π0
の比だったとします(足したら1として)。

この混合分布が一様分布する帰無仮説の分布に、0側によっている対立仮説が足されたのだとすると、1に近い方は殆ど対立仮説が含まれていないため、そこでのヒストグラムの高さは真の
π0に近い値をとっていると考えられます。ただし、1に近い程今度はπ0の予測値がばらつきやすいという問題もあります。

そのため、ヒストグラムのどの場所で高さをとるかという問題を考えたのがQ-valueです。
Q-valueでは、高さをとる場所をγ、γから推測される
π0の値をπ0(γ)としてプロットし、このデータを自然スプライン関数で回帰します。このスプライン関数をγ=0まで外推した値をπ0の予測値とします。

Q−valueのやっている事

実際にπ0の予測値を計算している様子

Q−valueはこの
π0をBH法のq値にかけます。それ以外の計算手順はBH法と同じなので、先にこのπ0の予測値をかけておいたp値をBH法に適用すれば、Q値を求める事ができます。

この手法により、Q<0.15でも発現変動遺伝子が検出されるようになりました。

Q<0.15で10遺伝子
 [1] 0.1119015 0.1414379 0.1119015 0.1119015
 [5] 0.1119015 0.1119015 0.1119015 0.1119015
 [9] 0.1119015 0.1414379

(補足)
Tcl/Tkとの依存関係でQVALUEのインストールがMacだとうまくいかなかったので、Q-valueの計算だけWindowsでやっちゃいました。みんなここではまってるっぽいです。

↓この方も
http://blog.hackingisbelieving.org/2010/09/r-qvalue-package-tcltk_10.html
Rをインストールする時にTcl/Tkの場所を指定すればいいらしいです。誰かやった人いたら教えてください。


ソースコード(ここから)

#ダウンロードしたGSE22531のzipを解凍
#GSE22531_RAWフォルダ内に作業ディレクトリを指定
#パッケージ読み込み
library("affy")
#そのディレクトリにあるCELファイルを認識
Data <- ReadAffy()
#MAS5で正規化
est <- mas5(Data)
#対数変換
express <- exprs(est)

#対照群
control <- express[,1:3]
#処置群
treatment <- express[,7:9]

#p値とq値
p <- c()
q <- c()

#全遺伝子のt検定のp値
for(i in 1:31099){
p[i] <- t.test(control[i,],treatment[i,])$p.value
}

#BH法のq値
q <- p.adjust(p,"BH")

#p値の分布
jpeg(file="pvalue.jpeg")
hist(p)
dev.off()

#q値の分布
jpeg(file="qvalue.jpeg")
hist(q)
dev.off()

#閾値を0.2に設定して8遺伝子
q[q<0.2]

#QVALUE計算用にp値を一度保存
write.table(p,"p-value.txt",quote=F,row.names=F,col.names=F)

#QVALUEのGUI操作でpi_0を概算(0.6794)
#新たにQ-Valueでのq値を計算するために
#このpi_0をあらかじめかけておいたp値を用意
p_Q <- p * 0.6794

#再びBH法
q_Q <- p.adjust(p_Q,"BH")

#閾値0.15に設定して10遺伝子(若干増加)
q_Q[q_Q<0.15]

ソースコード(ここまで)
(forループの部分はもう少し最適化したほうがいいかも)

参考文献
*データを正確に解釈するための6つのポイント、No.5 大規模データの解析における問題点、DNAマイクロアレイによる遺伝子発現量の測定を例として
*Q-valueの文献
Statistical significance for genomewide studies (2003)
*Q-valueを使った事例
Large-scale meta-analysis of cancer microarray data identifies common transcriptional profiles of neoplastic transfprmation and progression (2004)
*Q−valueのサイト
http://genomine.org/qvalue/
*hoxo_mさんのサイト
http://d.hatena.ne.jp/hoxo_m/
*ベイズウィキのサイト
http://hawaii.sys.i.kyoto-u.ac.jp/~oba/bayeswiki/index.php?QValue
Comments