HWE2015年度の修士向け講義は、「基礎の基礎」に戻ります(こちら)日本人類遺伝学会の教育講演も、それに沿った話にしますそのためのメモこちらに図などがうまくいっていないかもしれないepub(文書の散逸防止のためにkindleにも置いています。ためのp値の話』をまとめたものです。2x2分割表で確認する、p値、オッズ比とその信頼区間こんな2x2分割表があったとする。検定する・推定するカイ二乗検定をして、縦軸と横軸が無関係かどうかの検定をする。```{r}chisq.test(tab)```p値が帰ってきて、棄却するべきかどうかの情報が得られる。p値が小さいほど、無関係という仮説(帰無仮説)を棄却することになる。棄却する、ということにすれば、棄却するかどうかの決断ができる。「その基準を有意水準として、帰無仮説を棄却する」と説明し、「関連があった」と主張する。その逆なら「関連がなかった」と主張する。その心は、「帰無仮説が正しいと仮定したときに、このテーブルと同じかそれよりもっと偏ったテーブルを観察する確率はpくらい小さいから(pくらい大きいから)」。オッズ比は、相対危険度の推定値として使えるから。オッズ比〜相対危険度が1であるときが、「無関係」に相当する。推定値なので、「この値らしい」という点推定値があり、また、推定には幅があるから、だいたいこのくらいの範囲という意味で信頼区間がある(細かいことはここでは省略する)。「だいたいこのくらいの範囲」の中に、「無関係の相対危険度〜1」が入っているかどうかで、およその帰無仮説を棄却するかどうかのめども立つ。oddsratio(tab,log=FALSE)orconfint(or)```このようにデータが得られたら、たいていの場合、「検定」と「推定」をする。「検定」するのは、「帰無仮説」を捨てるため。「推定」するのは、相対危険度など、「因子の強さ」の程度に見当をつけるため。場合は、「検定」をして、「関連あり」を選んだ後、「その関連の強さを数値として理解する」というプロセスをとる。p値とオッズ比の信頼区間との大まかな関係たくさんの2x2表を適当に作り、カイ二乗検定p値とオッズ比の95%信頼区間を出して、その関係をプロットしてみる。横軸にp値、縦軸に95%信頼区間をとると、オッズ比の推定値とその信頼区間は1より大きい側になることもあれば、小さい側になることもあり、また、p値の大小の順番がオッズ比の信頼区間の上下限の値の大小の順番と完璧に一致するわけではないので、多少のでこぼこはあるが、総じて、p値の大小とオッズ比の信頼区間の大小には滑らかな関係があり、p値が0.05のあたりで、オッズ比の95%信頼区間に1が入るかどうかが分離できることが見て取れる。p値と一様分布、QQプロット起きやすさを0から1までの数値で表したもの。色々な検定方法などでは、異なる名前の統計量が計算されるけれど、最終的には、p値に変換することで、解釈がしやすくなる。p値が小さければ小さいほど、観測データと同じかそれよりも起こりにくいデータは生じにくいことを意味する。したがって、p値が小さいとき、その観測データを引き起こす前提として設定した仮説(帰無仮説)は信用しにくいことになる。p値が小さければ小さいほど、「観測データと同じかそれよりも起こりにくいデータが得られる確率の総和」は小さくなる。「観測データそのものが起きる確率」自体が小さいわけではない。実際、p値は0から1のどの値も同程度に起きやすく起きにくい。それを『帰無仮説が成り立つとき、p値は一様分布に従う』と言う。p値が仮説を信じるかどうかにおいて、わかりやすいのは、この『一様分布に従う』という性質のお蔭。一様分布とは、次のヒストグラムに示すように、等確率の分布のこと。ヒストグラムを描いてみる。value")```分布が一様であることを示すには、小さい順に並べてプロットすることでも示せる。```{r}plot(sort(ps))```このように、「観察されたp値は一様分布である」ということを、「グラフにすると対角線上に並ぶ」ということを利用したのがQQプロット。p値・一様分布の場合には、ただ、小さい順に並べてプロットすればよいけれど、それ以外(たとえばカイ二乗値)の場合には、小さい順に並べてプロットすると直線ではなくてカーブを描く。「たくさんの観測されたカイ二乗値が、カイ二乗分布になるはず」と思っているときは、横軸に「カイ二乗分布からの値が同じ数だけランダムに得られたとして、それを小さい順に並べたときの理論的な値」を、縦軸に「観察された値」を取ってプロットしたものが、QQプロット。対角線上に並んでいれば、「理論通り」と言える。QQプロットは、「帰無仮説が成り立つときに、たくさんのp値が得られたとしたときに、得られるべきp値の理論値を横軸に、得られたp値の観測値を縦軸にしたプロット」。観測が理論通りなら、QQプロットは対角線上、に一直線になる。2x3分割表とはSNVのほとんどは、2アレル(M,m)。常染色体上のSNVは3ジェノタイプ(MM,Mm,mm)。非病気のとき、フェノタイプは0/1の2値。トレンド検定3つのジェノタイプはあるけれど、アレルがリスクを生み出しており、アレル本数が0,1,2本と増えるにしたがってリスクが増えると考えるのは、よい考え方だ。その考え方に沿った検定がトレント検定とロジスティック回帰検定です。テーブルを作ってみる。3ジェのタイプのコントロールの数は等しく、ケースの数が増えている。matrix(c(40,80,50,50,50,50),2,3)*3tab```トレンド検定して、オッズ比も計算する。oddsratio(tab[,c(1,3)],log=FALSE)or3vs1#confint(or3vs1)```オッズ比は、1.5、2と計算できましたが、『トレンド検定をするという気持ち』に照らすと、このオッズ比はちょっと違います。トレンド検定をするというのは、線形回帰をしているのと同じです。線形回帰では、ジェノタイプごとのケースの割合が直線に乗るように推定します。(1:3)*a[2]points(cbind(1:3,p.prop),col=2,pch=20,cex=3)abline(h=c(0,1))abline(v=c(1,2,3))```このように、ジェノタイプごとのフェノタイプ率が推定できるので、基準とするジェノタイプを決めれば、それに対するオッズ比も計算できる。トレンド検定・線形回帰では、このようにして計算したオッズ比を想定していたことになる。トレンド検定のp値を用いて「関連あり」と信じるのであれば、トレンド検定が支持しているオッズ比は、1.5,2ではなく、線形回帰から計算されるオッズ比になる。ロジスティック回帰とその検定ロジスティック回帰もトレンド検定と同様によく用いられる。glm(p~g,family=binomial)summary(glm.out)$coefficients[2,4]```2ジェノタイプのオッズ比が得られる。```{r,warning=FALSE}exp(summary(glm(p~g,family=binomial))$coefficients[2,1])^(1:2)#summary(glm(p~g,family=binomial))$coefficients```ロジスティック回帰によって推定されるオッズ比の特徴は、ホモ接合体のそれがヘテロ接合体のそれの二乗になっていることである。トレンド検定では直線を引いたが、ロジスティック回帰では、シグモイドカーブを引く。tmp/(1+tmp)points(x,y,col=3,type="l")#abline(lm.out$coefficients[1],lm.out$coefficients[2],col=2)abline(h=c(0,1))abline(v=c(1,2,3))```トレンド検定とロジスティック回帰とを併せて描いてみる。直線と曲線はそれなりに違うけれど、3ジェノタイプに対応するx軸上の位置はほとんど同じであることが分かる。c(-2,6),ylim=c(-1,2),xlab="genotype",ylab="phenotype",main="トレンド検定とロジスティック回帰検定")abline(lm.out$coefficients[1],lm.out$coefficients[2],col=2)points(x,y,col=3,type="l")abline(h=c(0,1))abline(v=c(1,2,3))```式に表して整理してみる。\text{ジェタノタイプ}}$$となるように線を引いている。遺伝形式と検定法トレンド検定では、ジェノタイプ別の発病率の上昇分がリスクアレルの本数に比例していると仮定していたので相加モデルと呼ばれ、ロジスティック回帰では、ジェノタイプ別のオッズ比がリスクアレルの本数のべき乗になっていると仮定したので相乗モデルと呼ばれます。これらはいずれも、リスクアレル本数が1のとき(ヘテロ接合体)のリスクを二つのホモ接合体のリスクの中間的な値とするモデルです。これ以外に、単一遺伝子病などに認められてきた遺伝形式に優性形式と劣性形式があります。合わせた検定というものがあります。優性形式では、リスクアレルを0本持つか、1本以上持つかの2つに分けて検定しますし、劣性形式では、その逆で、2本持つか、そうでないかに分けて検定します。2の両方の検定を実施して、結果として、優性形式と劣性形式とについて検定することになります。検定をするのかしないのか」です。検定をしなかったとして、結果をどう解釈するか」です。1つの2x3表に複数の検定をすると、複数のp値が得られる。本当は「優性、または、劣性」な影響があるとすると、偽陰性が増える。ジェノタイプ別のオッズ比を計算して、「優性形式っぽい」ときには、「優性なのかも」と思い、「劣性形式っぽい」ときには、「劣性なのかも」と思う、という考え方もある。トレンド検定・ロジスティック回帰をしたときには、「そのモデルに基づいたジェノタイプ別のオッズ比は『すでに出ている』」わけで、ジェノタイプ別にオッズ比を出しなおして、遺伝形式を問題にしないほうがよい、という考え方もある。相加・相乗モデルでの検定で、ひとまず「関連があるかないか」の選別をして、そのうえで、「相加・相乗モデル」と「優性(劣性)モデル」のどちらが、よりもっともらしいのか、を数字に変えることは可能だけれど、通常、そのようなことはやらない。戻って、「優性形式が真」のときに、「相加モデルだけ」で検定すると、どれくらい偽陰性になるのかを見てみることにする。優性形式のときに相加モデルで検定するとたとえば、頻度が0.3のアレルがリスクアレルであって、優性形式でフェノタイプに影響しているとする。今、リスクアレルを1本または2本持つと、リスクが1.2倍になるとし、ケース・コントロール、ともに1000人ずつで、関連を調べるとする。相加モデルで検定するか、優性モデルで検定するかで、p値はどのようになるか、その違いが問題となる。優性モデルで検定すると有意になるのは30%程度なのに対して、相加モデルで検定すると、23%程度となる。このパワーの差が、真のモデルに沿って検定するかどうかの違いである。どちらの検定でもp値が0.01未満になる場合、どちらの検定でもならない場合、どちらか片方の検定でのみp値が0.01未満になる場合の4通りのうち、優性モデルで検定すれば、有意になるのに、相加モデルで検定したがために有意にならなかった場合が問題である。図ではそれを赤点で表している。先ほどの例で、相加モデルと優性モデルとの両方で検定をして、どちらでもよいから、有意なp値が得られたら、「関連あり」とみなすことにすれば、パワーが上がってよいだろう、ということになる。図で示せば、色のついたところのすべてを関連ありとする、ということである。2plot(log10(p.add),log10(p.dom),pch=20,cex=0.1,col=col)abline(0,1,col=3)abline(h=-2,col=3)abline(v=-2,col=3)```しかし、残念なことに、今度は、パワーではなくて、偽陽性の方に注意しなくてはいけなくなる。リスクがないSNVで同じことをやってみる。「関連あり」となってしまうことがわかる。「関連あり」となることを許容しているのに、1.4%も「関連あり」としてしまうのは、偽陽性が1.4倍に増えるから問題である。偽陽性と偽陰性に折り合いをつける〜マルチプルテスティング対策〜優性モデルにもパワーを出しつつ、偽陽性を増やさないようにしたい。そのためには、優性モデル検定と相加モデル検定の両方を実施して、そのどちらかで小さいp値が出たら「関連あり」とすることにしつつ、その「関連あり」とみなすための判定基準のp値を厳しめにする、という作戦をとるとよい。図で示せば、基準となる線を引きなおして、偽陽性が1%になるようなところを見つければよい。p値基準を見つけてみる〜マルチプルテスティング対策〜力技で見つけてみる。min(abs(log10(f)+2)))[1]abline(v=log10(p.check[clst]),col=2)print(p.check[clst])```このように、「帰無仮説」のときに、複数のp値がどういう組み合わせで出るのかがわかっていれば、「関連あり」とみなすp値の基準を決めて、偽陽性が0.01になるようにできることがわかります。問題は、「複数のp値がどういう組み合わせで出るのか」がわかっていないことです。素朴なマルチプルテスティング対策2つの検定のp値を縦軸と横軸に取ります。2つの検定のp値がお互いに関係がないものとします。優性モデルと相加モデルのときとは違います。p=0.01はわかりにくいのでp=0.1で考えてみます。ボンフェロニ法、Sidak法2つの検定を実施したら、p値を2倍して解釈する。n個の検定を実施したら、p値をn倍して解釈する。黒が90%くらいになるとよいのですが、90%より多いです。「関連あり」は少なすぎるということです。ボンフェロニ法は「関連を少なくし過ぎる」ので「保守的すぎる」と言います。「2つの独立な検定」の場合で言えば、黒い点がある部分の面積が、$$0.9^2=0.81$$と、0.9より小さすぎたので、$$0.95^2=0.9025$$にしてみたら、今度は、0.9より大きくなってしまった、ということです。(0.9487)^2$$にすると、ちょうどよいはずなので、そうすることもできます。Sidakの補正方法です。すれば、2つの検定が独立であれば、この方法でうまくいきます。どちらか片方でp<0.1であるものを「関連あり」とすると、10%だとよいところが15%を「関連あり」とみなしてしまって、多すぎます。ボンフェロニ法にすると、7.7%まで減ってしまって、保守的にすぎます。2つの検定が独立だったときは、Sidakはよかったのですが、優性と相加のように相互に近い関係だろ、Sidak法でもボンフェロニ法とほとんど変わらず(7.7%)、まだまだ保守的です。分布がわからないので、どうするか分布がわかれば、補正の仕方はわかるということはわかりました。実際には分布がわからないから困っているわけで、分布がわからないときにどうするかが問題です。わからない分布を調べてから、それに基づいて補正するの2つです。分布がわからないままに、補正する〜ストイックであれば大丈夫〜補正をしないままに「関連あり」と主張すると、「それは言い過ぎ」と反論されます。逆に、ボンフェロニ法・Sidak法で補正したp値を使って「関連あり」と主張するのは、2つの検定が独立でも「言い過ぎ」ではないし、まして、2つの検定が独立でないときには、「もっと安心」なので、「関連あり」と言ってもよい、と言うことになります。逆に「関連なし〜関連ないと信じても悪くない」と主張する場合には、ボンフェロニ法・Sidak法で補正したp値を使った場合には、「それは言い過ぎ」ということになりますし、補正しないままのp値で判断するなら、「関連なし」と主張する場合には、「安心」ということになります。場合は「関連あり」と主張するので、「ボンフェロニ法・Sidak法で補正しているから安心」ということになります。『こうだったら、論文がアクセプトになって好都合なのに』と思ったら、『自分に不都合に』判断することです。『ストイックに・禁欲的に』ですね。分布がわからなければ、分布を調べればよい〜正確確率法、ランダマイゼーション・パーミュテーション法〜わからないなら、調べればよいです。こんな2x3表が得られたとします。今、知りたいのは、「関連がない」との仮定の下で、ケースとコントロールを、それぞれ55人、30人のサンプリングして、『このSNV』のタイピングをした結果の表をたくさん作ること、そして、その表に優性・相加の検定を実施してp値を求めることです。言っても、『このSNV』の集団でのジェノタイプ頻度の正確なところはわかっていませんから、仮想の表をたくさん作るのは難しいです。いうときは、分割表の周辺度数を固定して、周辺度数が同じ表をたくさん作るという方法もあります。いわゆる、分割表の正確確率検定で、『周辺度数を等しくする表のすべての生起確率を計算』して、検定p値を決めるのと同じ方法です。可能なすべての表を全部列挙してその生起確率を計算することは、できなくはないですが、結構、大変です。その大変さは、表の自由度に応じてどんどん大きくなります。2x3表の自由度は2なので、それほど大変ではないので、正確に確率を計算することはできますが、一般的なマルチプルテスティング対策としては、現実的ではありません。そこで乱数を発生させて、正確確率法と同じような分布を得る方法がランダマイゼーション法・(乱数を使った)パーミュテーション法です。matrix(c(10,5,10,20,10,30),2,3)tab```ランダマイゼーション法・パーミュテーション法は簡単です。今、2x3表では、サンプル全員のフェノタイプとジェノタイプが観察されていますが、ランダマイゼーション法では、誰にどのフェノタイプが観測されたかを、『なかったこと』にして、ランダムに割り付けなおします。すると、フェノタイプが観察と変わりますから、ジェノタイプと突き合わせてできる2x3表が変わります。変わりますが、ケース・コントロールの人数、ジェノタイプ別の人数は変わりません。周辺度数は変わっていないということです。この割り付けを$$n!$$通り、すべて実施すれば、正確確率法と同じになりますし、$$n!$$は多すぎるから、適当な値でやめにすれば、適当な数の表がランダムに作れたことになります。作った表で優性・相加の検定をしてp値を出してやれば、目的は達成されます。分割表の周辺度数が等しいという条件の下、ある分割表を観察する確率を正確に計算したり、それに準ずる方法としてランダマイゼーション・パーミュテーション法を用いることで、互いに独立でない複数の検定のp値のマルチプルテスティング補正ができることがわかりました。ここで、正確確率検定に話を戻します。分割表があったときに、トレンド検定を実施することもできますが、正確確率検定を実施することもできるからです。2x3表だと話が少し複雑になるので、2x2表に戻って考えることにします。>「分割表のセルの期待値が小さめのときには、表にカイ二乗検定をするのではなく、フィッシャーの正確確率検定をする方がよい」と言われています。正確確率検定の方が『保守的』だから例を見ておきます。matrix(c(18,2,4,10),2,2)tabchisq.test(tab,correct=FALSE)$p.valuefisher.test(tab)$p.value```これらは正しいので、表を観察したら、この原則に従うのが良いのです。良いのですが、GWASなどでは、山ほど検定を行います。実際にp<0.05となる割合は、普通のカイ二乗検定ならは0.05くらいですが、フィッシャーの正確確率検定をすると0.012くらいになります。偽陽性を出さないためには、よいことですが、パワーのことを考えると、随分、かい離した値と言えます。これは、表の期待値が小さめのときだけに正確確率検定を適用した場合でも、ほとんど変わりません。p値の分布も、『一様』からは程遠いです。このような条件でのp値はp<0.05の割合が少ないだけでなく、p値の平均も随分と違います。カイ二乗検定ならばp値の平均は0.5程度となり、一様分布らしさがありますが、正確検定を用いると、平均が0.7くらいになります。たくさんの検定の結果からp値の分布と取り出し、それが一様分布になっているとよい、という考えの下、たくさんのp値の処理をするようなときに、この正確検定の『保守性』が邪魔になることもあります。正確検定ならいつもよい、というわけではない、という話としてとらえておくことにします。一つのSNVに1つの検定をするか複数の検定をするか、という話をしてきました。一つの遺伝子に複数のSNVがあり、その遺伝子とフェノタイプとに関連があるのかを調べることと似ています。結局2x9表を観察することになります。今、2個のSNVについて、それぞれ検定をしてフェノタイプとの関係を調べるとします。matrix(c(1,2,3,1,2,3,1,2,3),3,3))```という重みを付けることだし、もう一つのSNVについて相加モデルで検定するというのは```{r}t(m33)```という重みを付けることですから、2x9表に2検定のマルチプルテスティングをしていることになって、先ほどの議論と同じです。2x9表の2つのSNVのそれぞれに相加モデルと優性モデルとをそれぞれやるのであれば、次のような4つの重みづけでマルチプルテスティングをしていることになります。matrix(c(0,1,1,0,1,1,0,1,1),3,3))t(m33.2)```検定間が独立であれば、Sidak法が有用で、そうでなければ、補正をしないなら、『保守的に』、補正をするなら『正確に確率を計算』するか『ランダマイゼーション法・パーミュテーション法』が有用なのでした。100人のケースコントロール関連検定を実施するとします。個々のSNVについて、相加モデル検定をすることにすれば、n.snv=20個の検定となります。今、有意水準0.01での棄却をするには、20個のp値の中の最小のp値の分布を調べて、ある確率0.01未満でしか観測されないくらい小さいp値とはどれほどかを調べることになります。LDのプロットではこれくらい強い領域であるとすると、それは約0.001と10分の1くらいに小さくする必要があることがランダマイゼーション法・パーミュテーション法から示されます。独立検定何個分に相当するか。これを実際に実施した検定数と比較してみるprint(correct.rate[1])```ここでは、ある特定のLD領域にある複数のSNVのことを扱いましたが、これをゲノムワイドに広げても考え方は同じです。?というのが、ゲノムワイドな補正p値の目安、と言うことです。百万個の場合にカバーしなかった場所をカバーすることにしたのなら、マルチプルテスティングで目指すべきp値は増やした分だけ下げる(十分の1にする)ことになります。「真の遺伝因子」があって、その近傍に「真の遺伝因子と連鎖不平衡関係にある代用マーカー」があるので、代用マーカーを調べることで、「真の遺伝子の位置」の見当をつけます。連鎖不平衡とは、A/aというバリアントとB/bというバリアントがあったときに、A-B,A-b,a-B,a-bの4つのペアの割合について、2x2分割表を作り、それについての独立性検定をするようなものです。サンプル総数が1の2x2表のカイ二乗値の最大値は1で、最小値は0となり、それが連鎖不平衡係数$$r^$$に相当します。念のためやっておきます。カイ二乗テストをすると、そのカイ二乗値はr^2r^2chisq.test(pH.d,correct=FALSE)$statistic```では、2SNVがあって、その片方が真のリスクを有しているとして、真のリスクSNVとLD関係にあるSNVとで相加モデル検定をしたときに、p値がどんな風になるかを見てみます。低頻度のアレルをリスクアレルとし、リスクアレルの保有コピー数が0,1,2のときに有病率を6%,8%,10%として見ます。このようなときに、ケース・コントロール1000人数ずつで解析をすると、2SNVのp値はいずれも低めに出て、両者には相関があります。真のリスクSNVのp値の方が少し低めです。hist(log10(c(p.out$p1,p.out$p2)),col=gray(0.99),border=FALSE)hist(log10(p.out$p1),density=17,add=TRUE,breaks=h$breaks)hist(log10(p.out$p2),col=2,density=20,add=TRUE,breaks=h$breaks)par(mfcol=c(1,1))cor(p.out$p1,p.out$p2)print(r.squared)mean(p.out$p1)mean(p.out$p2)```では同じ条件で、$$r^2$$だけを変化させるとどうなるかを見てみます。2SNVのp値の相関が$$r^2$$に応じて変化します。相関係数と$$r^2$$とが同じ値です。真のリスクSNVのp値がLD関係にあるSNVのp値以下になる頻度も$$r^2$$に応じて高くなることがわかります。$$r^2=0.9$$くらいでは、半分くらいがそうなっています。このように、$$r^2$$で計られたLDの強さは、関連検定のp値の相関係数と関連があり、また、強いLDにあれば、真のSNVよりp値が小さくなる確率も上がることがわかります。アレルで検定する・ハプロタイプで検定するあるアレル・あるハプロタイプに着目すると、個人は、そのアレルを0本持つか、1本持つか、2本持つかの3通りあります。1個のSNVについてこれを調べるには、相加モデル検定(かロジスティック回帰検定)をすればよいです。SNVを組み合わせてハプロタイプに分解するには、SNPごとのジェノタイプを実験で調べるだけでは無理なので、ハプロタイプ推定をすることもあります。この推定の大規模な場合がHLA領域のハプロタイプ推定になります。シークエンシングを用いて、うまくハプロタイプごと検定することもできるかもしれません。その場合は「推定」は不要です。matrix(c(m[1,1]*2+m[1,2],m[1,2]+2*m[1,3],m[2,1]*2+m[2,2],m[2,2]+2*m[2,3]),byrow=TRUE,nrow=2)a```というような2x2表にして検定する、ということです。カイ二乗検定してみます。```{r}prop.trend.test(m[1,],apply(m,2,sum))chisq.test(a,correct=FALSE)```同じです。トレンド検定なんていう面倒くさいことしなければよいということになりそうです。matrix(c(m[1,1]*2+m[1,2],m[1,2]+2*m[1,3],m[2,1]*2+m[2,2],m[2,2]+2*m[2,3]),byrow=TRUE,nrow=2)aprop.trend.test(m[1,],apply(m,2,sum))chisq.test(a,correct=FALSE)```また、同じです。matrix(c(m[1,1]*2+m[1,2],m[1,2]+2*m[1,3],m[2,1]*2+m[2,2],m[2,2]+2*m[2,3]),byrow=TRUE,nrow=2)aprop.trend.test(m[1,],apply(m,2,sum))chisq.test(a,correct=FALSE)```今度は違います。検定手法の名前が違うということは、トレンド検定と2x2表化したカイ二乗検定は『同じなわけではない』ということですね。同じになった場合と同じにならなかった場合との違いはなんなのでしょう。matrix(c(23,26,24,18,6,3),byrow=FALSE,nrow=2)apply(m,2,sum)```となっています。0.42$$$$0.3^2=0.09$$という比率になっています。これを遺伝学分野ではハーディ・ワインバーグ平衡(HWE)にある、と言いますが、HWEにあるかそうでないかが影響します。ケース・コントロールを合わせたときのジェノタイプカウントがHWEに近ければどちらの方法でやっても結果に大差はないが、HWEからの逸脱が大きいときには、2つの手法の結果は違う、違うからにはどちらかを採用するのがよいことになる、それはトレンド検定である、ということになり、いつもトレンド検定をできるのならばそうしておくのが良いでしょう。ハプロタイプの場合は、アレルの数がハプロタイプの種類数となりますが、考え方は同じです。ハプロタイプの本数をケース・コントロール別に数えて検定することは、ハプロタイプ全体でHWEになっていそうかどうかの影響を受けるということになります。実際のところ、ハプロタイプの本数推定しかできないこともあり、その場合には、ハプロタイプのペアがどうなっているのかまでわからないこともありますから、ハプロタイプ本数での検定をせざるを得ないのが現実ですし、それが許容されてもいます。ハプロタイプ本数をケース・コントロール間で比較するときも同様です。ハプロタイプ本数での検定をするときには、「ハプロタイプの種類数だけ検定を繰り返し」ますから、その分のマルチプルテスティングのことを考慮する必要があることは、もうわかると思います。rmultinom(1,n.co,prob=c(g12))matrix(c(g.ca,g.co),byrow=TRUE,nrow=2)```では、そのHWEについて考えてみます。HWEとはHWEとは、SNVの3ジェノタイプの割合の偏りがあるかどうかの検定です。偏りがないというのは、2アレルの組合わせのうちホモ接合体(MM,mm)が妙に多すぎたり、少なすぎたりしないということです。いうときに偏るかと言えば、Mアレルが多い集団とmアレルが多い集団があって、2つの集団が別々である場合です。極端な例では、Mアレルしかない地域とmアレルしかない地域では、それぞれ、全員がMM,mmですから、この2地域からのサンプルを合わせると、ヘテロ接合体Mmが一人もいないことになります。実際にはこれほど極端なことはないですが、地域差があって、サンプルがその混合であるときにはホモ接合体の割合が多くなります。偏りがないときというのは、地域差がなく、よく混ざり合っているときのことで、ランダム・メイティングしているとき、と言い換えることもできます。0.09$$となっているのが、偏りがない状態(平衡状態)です。HWE検定で実験の失敗を疑うHWE検定では、確かに、集団の偏り具合(集団の構造化)の有無がわかるのですが、たくさんのSNVを使ったGWASのようなスタディでは、集団の構造化については、個々のHWE検定では判断しません。あるSNVのデータを見ると、構造化していて、別のSNVのデータを見ると構造化していない、ということになったら、個人の集まりとしての構造化はどうなっているのかに対して、答えが山ほど出てしまって判断できなくなるからです。HWE検定では、もう一つ別のことを調べます。3ジェノタイプの観察人数をHWE検定したところ、HWEから外れていると判断できたとします。ように、ホモ接合体が多い方に外れていたり、v.hwD.2のようにホモ接合体が少ない方に外れていたりした場合です。c(44,52,4)```検定をすると、確かにp値が小さめです。pchisq(sum((v.hwD.2-v.hwE)^2/v.hwE),1,lower.tail=FALSE))```集団構造化によるHWEからの逸脱であれば、ホモ接合体が多めになるはずですが、ホモ接合体が少な目になることの説明はつきません。この場合は、「本当のジェノタイプはHWEなのだけれど、実験したら、何かの理由でホモ接合体のサンプルがヘテロ接合体とコールされてしまった」からかもしれません。実験エラーによるHWEからのずれは、ホモ接合体のコール数が多い方への偏りも起こし得るので、ホモ接合体が多く観察された場合には、集団構造化と実験エラーの両方を疑う必要があります。このようにHWE検定では集団構造化ではなく、実験エラーを見つけるために使うこともできます。2x3表でケース・コントロールの関連を検定したときには、p値のほかに、検定の強さを表すオッズ比のことを気にしました。HWE検定でもHWEからの偏りの強さのことを気にします。HWEに合致しているときが0、ヘテロ接合体が0人のときが1になるような指標があるとわかりやすいです。ホモ接合体が多い方に偏っているのか、少ない方に偏っているかもわかると便利です。\frac{\text{観測されたヘテロの割合}}{\text{アレル頻度から期待されるヘテロの割合}}$$のような式が好都合です。Fixation指標です。HWE検定のp値が小さいときHWE検定では『サンプルは集団構造化がない母集団を代表しており、実験もうまく行っている』というのが帰無仮説ですから、もしHWE検定のp値が小さかったときには、*『サンプルは集団構造化がある母集団を代表している』か*『実験がうまく行っていない』かのいずれかの可能性を否定しにくいことがわかります。後述しますが、ある程度の集団構造化があってもそれに対する方策はありますので、集団構造化だけを理由に特定のSNVを排除しなくてもよいです。HWE検定を使って排除したいのは、『実験がうまくいっていない」場合です。HWE検定もマルチプルテスティングの影響を受けますから、100万マーカーを検定すれば、HWE検定のp値が100万分の1であっても、それだけで問題にはなりません。実験の失敗による異常は、HWE検定p値が、そのほかの大多数の「成功していると思われるマーカーのp値の分布」に照らして、「外れている」ことから見当をつけるのがよいでしょう。p値が一様ではないp値は、一様分布に従っているから、その値を0.01と聞けば、「あー、0.01的に珍しいことなんだ」とわかるわけですから、p値の本領は一様分布になっていることです。実際にGWASを実施して、数十万個のp値を算出して、その分布を見てやると、一様分布になっていない。このときに大きく分けて2つのアプローチがあります。ほとんどすべてのp値は、帰無仮説が正しいはずだから、「本当は一様分布」なはず。それなのに、一様になっていないのだったら、「一様分布」になるように修正してしまおう、という作戦。一つは、たくさんのp値があるけれど、結構の数の検定が、実は対立仮説が真なのだから、全部のp値が作る分布が一様分布になるわけがない。真に関連がある結果はなるべく拾いつつ、本当は帰無仮説が正しいものが偽陽性にならないように工夫する。一様分布に修正する作戦〜ジェノミック・コントロールケース・コントロール関連をたくさんのマーカーについて検定を繰り返すとき(GWASのように)、ほぼすべてのマーカーが疾患と無関係なので、理想的にはp値は一様分布になります。それはカイ二乗検定をしている場合で言えば、たくさんの検定統計量がカイ二乗分布に従っている、ということです。ケース群とコントロール群とで「違いがない」のに、カイ二乗値のすべてが0にならず、大きめの値が出たりするか、と言えば、サンプリングのばらつきがあって、「たまたま、ケースとコントロールとで違いが出る」からです。この、「たまたま出る違い」は、「ケースとコントロールをサンプリングする母集団」が「きれいに均一」な場合には、カイ二乗分布そのものになるのですが、「母集団が、似ている人、似ていない人が混ざっている」と、たとえ、「母集団からケースとコントロールとをランダムにサンプリングし」ても、「たまたま」起きる「ケースサンプルとコントロールサンプル」の違いは少し大きめになります。hist(c(k,k.),col=gray(0.99),border=FALSE)hist(k,density=17,add=TRUE,breaks=h$breaks)hist(k.,col=2,density=20,add=TRUE,breaks=h$breaks)```p値で見ると、小さいp値が多めになって、このようになります。hist(c(ps,ps.),col=gray(0.99),border=FALSE)hist(ps,density=17,add=TRUE,breaks=h$breaks)hist(ps.,col=2,density=20,add=TRUE,breaks=h$breaks)```一様分布ではないp値はそのまま信じるわけにいかないので、補正が必要です。今の状況は、「ケースとコントロールとは、もともと違いがな」く、「サンプリングも母集団からうまくされている」けれど、「母集団が均質ではなくてばらつきがある」という単純な状態なので、補正も単純にできることが知られています。カイ二乗分布が大きめにシフトしているので、すべてのカイ二乗値をある値で割ってやって、うまく、p値が一様になるようにできることが知られています。このうまく補正する値が、ジェノミックコントロールのλ値です。1は補正の必要がないことを意味し、1より大きければ、p値が小さめに出ているので、大きめに補正してやる必要があることを意味します。補正前と補正後のp値をQQプロットした図を提示すれば次のようになります。補正前の赤いラインはカーブしていますが、補正後の黒い線は対角線になっていて、うまく補正で来ていることになります。pchisq(k.corr,df,lower.tail=FALSE)hist(ps)hist(ps.corr)qqplot(qunif(ppoints(n)),ps.corr,xlab="Theoretical",ylab="Observed",type="l")points(qunif(ppoints(n)),sort(ps),col=2,type="l")```ジェノミック・コントロールのようなタイプの補正では、小さいp値ほど補正によるp値の増大が大きくなります。元のp値(p.nominal),補正後p値(p.corr),両値の比(fold)としてそれを見てみます。たくさんの本物がある場合〜FDR〜たしかに、GWASでゲノム全体のどこに関連のある遺伝子があるかわからないけれど、あったとしても、ごく一部の遺伝子しか、フェノタイプと関係ないだろう、と考えるときには、得られたp値が一様分布をとると仮定するのが妥当です。そのうえで、ほとんどすべてのマーカーを関連なしとみなしたいので、「厳しく」棄却するのがよいです。その目的にかなっているのが、ボンフェロニ法やSidak法でした。2万数千個の遺伝子の発現量を組織Xと組織Yで比べたり、癌部と非癌部とで比べたりする場合には、そもそも、同じわけがない、と考えていますし、実際に、調べると、かなりの数の遺伝子の発現量が明らかに違うことがわかります。いうときに、厳しく切り捨てると、「かなりの『本物の関連』」も切り捨てられてしまうので、目的に適わないという問題が生じます。Rateと呼ばれる方法がそれです。ボンフェロニ法はSidak法では、「棄却基準」をすべてのp値に等しく適用しました。たくさんの検定をしたのだから、それ全体を使って判断しよう、という考え方です。次のQQプロットは、帰無仮説が真である検定200個と対立仮説が真である検定40個が混ざっている場合で、黒が帰無仮説を真とする結果、赤が対立仮説を真とする結果を表しています。order(k,decreasing=TRUE)qqplot(qunif(ppoints(length(k))),p,col=col[ord],pch=20,xlab="Theoretical",ylab="Observed")```赤のp値は小さい値が多いので、左下に多いです。黒のp値は一様に分布していますが、赤よりも相対的に大きいので、右に偏って分布しています。このようなたくさんのp値が得られたときに、赤・黒の情報はないわけですが、赤は、「拾っ」て、黒は「捨て」たいので、そのような方法を作ろう、という考え方です。ボンフェロニ法やSidak法では、以下のように水平線を引いて、その線より下に来た点を「拾う」方法です。```{r}qqplot(qunif(ppoints(length(k))),p,col=col[ord],pch=20,xlab="Theoretical",ylab="Observed")abline(h=0.01,col=4)```FDRは、水平線ではない基準を入れる方法で、たとえば、原点を通る次のような直線を引きます。```{r}qqplot(qunif(ppoints(length(k))),p,col=col[ord],pch=20,xlab="Theoretical",ylab="Observed")abline(0,0.1,col=4)```p値が小さい部分を拡大してみます。```{r}qqplot(qunif(ppoints(length(k))),p,col=col[ord],pch=20,xlab="Theoretical",ylab="Observed",xlim=c(0,0.2))abline(0,0.1,col=4)```なぜ、水平線にしないかを考えてみます。たくさんのp値があって、それが一様分布に従っているときに、一番小さいp値は、小さい値を中心に少しばらつきますし、一番大きいp値は、大きい値を中心にばらつきます。得られたp値を小さい順に並べたとき、一番小さいp値の場合にそれが、「帰無仮説を信じるには小さすぎる」という基準の値は、二番目に小さいp値のための基準値よりも小さくしてもよさそうです。三番目に小さいp値の場合はさらに緩めにすることができるでしょう。このように何番目に小さいかによって基準値を変えてやります。言っても2番目に小さいp値が2番目用の基準値をクリアしなかったけれど、3番目に小さいp値が3番目用の基準値をクリアしたからと言って、2番目に小さいp値の結果は、関連なしとし、3番目は関連ありとするのは、いかにもおかしいです。したがって、何番目に小さいかによって基準を定め、そういう基準をクリアしたp値を見つけたら、それ以下のp値はすべて「関連あり」側にする、という約束も付け加えます。この直線の引き方としては、傾きが0.05(たとえばの水準)を引くことです。たくさんの検定のうち、最も大きいp値を出したものについて、それが帰無仮説を棄却する場合を考えるとします。大きいp値を棄却することとしたら、それより小さいすべてのp値の検定の帰無仮説を棄却することとなります。すると、今、「帰無仮説が真」かもしれない仮説は、問題にしている検定の対象だけです。したがって、その検定を棄却するかどうかは、単独の検定について考えればよいから、通常の0.05や0.01を使えばよいです。それ以外のp値(二番目、三番目...に大きいp値)については、別途基準を決める必要がありますが、思い切って、原点からの直線を引いてしまえ、というやり方をここでは示しました。他にも複数のやり方が提唱されていますが、大まかな発想を理解するのには、この方法で十分でしょう。GWAS/GWASメタアナリシスとその周辺解析用のアプリケーションへのリンク複合遺伝性疾患とは----------News----------20万アクセス2008/07/2040万アクセス2010/09/12-------------------------はてなブログをはじめよう!ryamada22さんは、はてなブログを使っています。ブログを報告するスターをつけました引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。