日本人類遺伝学会 2014/11/20 京都大学 医学研究科 統計遺伝学分野 山田 亮

Presentation on theme: "日本人類遺伝学会 2014/11/20 京都大学 医学研究科 統計遺伝学分野 山田 亮"— Presentation transcript:

1 日本人類遺伝学会 2014/11/20 京都大学 医学研究科 統計遺伝学分野 山田 亮
P値から考えるゲノム疫学解析 日本人類遺伝学会 2014/11/20 京都大学　医学研究科　統計遺伝学分野 山田　亮

2 SNPを用いたGWASなどのスタディに関わっている・かかわろうとしている大学院修士の学生さん、くらいを主なターゲットとします
ご参加のみなさんには、 SNPって何？　検定とp値？　そもそもたくさん検定をするって何？　という方や あそこの分割表の検定は、厳密には、●●しないといけないのでは という方が混ざっている・・・と思いますが、せっかくですので、何かしら得るものがありますように 今日のスライドの大部分のデータ処理はフリーソフトRで実施しています。そのコードを含めた資料は、分野facebookから入手可能です。聴講しながら、叩いていただいても結構です。

3 2x2分割表で確認する、p値、オッズ比とその信頼区間

4 2x2分割表で確認する、p値、オッズ比とその信頼区間
あり なし コントロール 20 27 47 ケース 28 25 53 48 52 100

5 検定する・推定する

6 検定する・推定する カイ二乗統計量 自由度 1 p-値

7 検定 p値 棄却するべきかどうかの情報 p値が小さいほど、無関係という仮説（帰無仮説）を棄却する

8 あり なし コ 20 27 47 ケ 28 25 53 48 52 100 推定 オッズ比　0.6614 95%信頼区間　0.3 ～1.458

9 推定 オッズ比を計算する。 オッズ比は、相対危険度の推定値 オッズ比～相対危険度が１であるときが、「無関係」に相当する 推定値は 点推定値
信頼区間 「だいたいこのくらいの範囲」の中に、「無関係の相対危険度～１」が入っているかどうかで、およその帰無仮説を棄却するかどうかのめども立つ

10 p値と CIが1をまたぐか の関係 緩く正しく 微妙に違う OR 95% CI 1.0 p値

11 p値はなんのため？ 一様分布 小さいp値も大きいp値も同じくらい出やすい 0.5が出やすいわけではない 1が出やすいわけではない

12 p値はなんのため？ 一様分布 小さいp値も大きいp値も同じくらい出やすい 0.5が出やすいわけではない 1が出やすいわけではない

13 一様分布のヒストグラム 1 p値

14 順番に並べてプロットすると 直線になる １番 １万番 小さい順

15 小さい順に並べて 対角線を描くのは QQプロット
対数を取ってもOK p値だけではなく カイ二乗値でもＯＫ Lanktree M B et al. Stroke. 2010;41: Copyright © American Heart Association, Inc. All rights reserved.

16 QQプロット 理論に合っているかを確認する
観測値 1 期待値・理論値

17 p値の基礎、終了 ぴーち

18 SNV解析の基本、2x3分割表 MM Mm mm コントロール 10 30 ケース 5 20 55 15 40 85

19 コントロールに比べて、ケースが： 0.5、2、3倍 MM Mm mm コントロール 10 30 ケース 5 20 55 15 40 85

20 さっそく「関連検定」しよう

21 さっそく「関連検定」しよう 面倒くささの元は

22 さっそく「関連検定」しよう 面倒くささの元は いくつも検定法があること

23 2x3分割表の検定法 ２つ トレンド検定とロジスティック回帰検定

24 トレンド検定する p 値 MM Mm mm コントロール 10 30 ケース 5 20 55 15 40 85

25 MM Mm mm コントロール 10 30 ケース 5 20 55 15 40 85 オッズ比は？
p 値 , (10x20)/(10x5) = 4 MM Mm mm コントロール 10 30 ケース 5 20 55 15 40 85

26 MM Mm mm コントロール 10 30 ケース 5 20 55 15 40 85 オッズ比その２
p 値 , (10x30)/(10x5) = 6 MM Mm mm コントロール 10 30 ケース 5 20 55 15 40 85

27 トレンド検定する p 値 , オッズ比 4 (Mm vs. MM) 6 (mm vs. MM)

28 トレンド検定のオッズ比… トレンド検定をするというのは、線形回帰をしているのと同じです。 線形回帰では、ジェノタイプごとのケースの割合が直線に乗るように推定します この『トレンド検定をするという気持ち』に照らすと、このオッズ比はちょっと違います どう違う？

29 トレンド検定という線形回帰 フェノタイプの現れる確率 1.0 0.0 ジェノタイプ

30 線形回帰は「傾き」が大事 フェノタイプの現れる確率 (p2/1-p2) 1-p2 (p1/1-p1) 1.0 1-p1 p2 p1 0.0
ジェノタイプ

31 線形回帰は「傾き」が大事 オッズ比 2.1 4 (Mm vs. MM) 5.1 6 (mm vs. MM) フェノタイプの現れる確率
(p2/1-p2) 1-p2 (p1/1-p1) 1.0 1-p1 p2 p1 0.0 オッズ比 2.1 4 (Mm vs. MM) 5.1 6 (mm vs. MM) ジェノタイプ

32 線形回帰は「傾き」が大事 オッズ比 2.1 4 (Mm vs. MM) 5.1 6 (mm vs. MM) フェノタイプの現れる確率
(p2/1-p2) 1-p2 (p1/1-p1) 1.0 1-p1 p2 p1 0.0 オッズ比 2.1 4 (Mm vs. MM) 5.1 6 (mm vs. MM) ジェノタイプ

33 ロジスティック回帰する p 値 (トレンド検定) p 値 (ロジスティック回帰検定) 大差ない

34 ロジスティック回帰する フェノタイプの現れる確率 1.0 0.0 ジェノタイプ

35 ロジスティック回帰する オッズ比 2.1 → 2.3 (Mm vs. MM) 5.1 → 5.2 (mm vs. MM)
フェノタイプの現れる確率 1-p2 1.0 1-p1 p2 p1 0.0 オッズ比 2.1 → 2.3 (Mm vs. MM) 5.1 → 5.2 (mm vs. MM) ジェノタイプ

36 ロジスティック回帰する オッズ比 2.1 → 2.3 5.1 → 5.2 = 2.3 x 2.3 フェノタイプの現れる確率 1-p2 1.0
0.0 オッズ比 2.1 → 2.3 5.1 → 5.2 = 2.3 x 2.3 ジェノタイプ

37 トレンド検定とロジスティック回帰

38 トレンド検定とロジスティック回帰 直線と曲線の 違いはあるけれど 大差ない

39 大差がない ロジスティック　log10(p) トレンド　log10(p)

40 違いは トレンド検定は ロジスティック回帰検定は 相加モデルに相当する 計算が簡単
正確検定もできる(低アレル頻度・少サンプル数の場合の対処法がある) 個人別のデータが不要 「線形回帰」に基づくオッズ比が結果に付いていないことが多い 共変量を組み込めない(線形回帰にすればよいが、それならロジスティック回帰をすればよい) ロジスティック回帰検定は 相乗モデルに相当する 計算が面倒ではある(けれど計算機がやってくれるので問題はないが、検定個数が大量になるとそれなりに影響してくる) オッズ比が結果についてくることが多い(ただし、係数は対数で返ってくることが通例) 年齢・性別などの共変量を組み込みやすい

41 遺伝形式 相加　相乗　優性　劣性 面倒くささの元は いくつも検定法があること

42 優性

43 優性

44 劣性

45 優性・劣性形式の検定 するのか、しないか することのメリットとデメリット しないことのメリットとデメリット
優性・劣性形式に照らした結果がわかる 1つの2x3表に複数の検定をすると、複数のp値が得られる。複数のp値が出たら、そのp値はには「補正」をしないといけない しないことのメリットとデメリット 1つのp値しか出ていなければ、複数のp値の補正について悩む必要はない 優性・劣性形式に照らしての判断になっていない

46 優性・劣性形式の検定をしなかったら 本当は「優性、または、劣性」な影響があるとすると、偽陰性が増える。
では、「優性形式が真」のときに、「相加モデルだけ」で検定すると、どれくらい偽陰性になるのかを見てみることにする。

47 アレル頻度 0.3 RR(優性) RR=1.2 1000人 vs. 1000人 優性モデルlog10(p) 相加モデルlog10(p)

48 アレル頻度 0.3 RR(優性) RR=1.2 1000人 vs. 1000人 優性モデルlog10(p) 30% 相加モデルlog10(p)

49 アレル頻度 0.3 RR(優性) RR=1.2 1000人 vs. 1000人 優性モデルlog10(p) 30% 24% 相加モデルlog10(p)

50 アレル頻度 0.3 RR(優性) RR=1.2 1000人 vs. 1000人 優性モデルlog10(p) 偽陰性 30% 24% 相加モデルlog10(p)

51 優性モデルでは拾えずに相加モデルで「たまたま拾う」こともある
アレル頻度 0.3 RR(優性) RR=1.2 1000人 vs. 1000人 優性モデルlog10(p) 偽陰性 相加モデルlog10(p)

52 優性モデルlog10(p) 相加　優性　両方を併せると パワーが上がる 相加モデルlog10(p)

53 パワーが上がるのはよいことだ

54 相加・優性どちらも『あり』にしたら 偽陽性が1.8倍になった
真の優性座位 無関係の座位

55 パワーが上がると 偽陽性が増える いいことがあると 悪いこともある
パワーが上がると 偽陽性が増える いいことがあると 悪いこともある

56 パワーを上げつつ、偽陽性を増やさない 0.01より小さい『新たなp値基準』を作る
1%

57 稼ぐ 失う

58 『新たなp値基準』を探す この図があれば、できる けれど、この図はない(すぐには手に入らない) 1%

59 素朴なマルチプルテスティング対策

60 素朴なマルチプルテスティング対策 黒、赤、緑、青、の比率がわかれば、 『新たなp値基準』はわかる

61 素朴なマルチプルテスティング対策 黒、赤、緑、青、の比率がわかれば、 『新たなp値基準』はわかる

62 素朴なマルチプルテスティング対策 黒、赤、緑、青、の比率がわかれば、 『新たなp値基準』はわかる

63 素朴なマルチプルテスティング対策 黒、赤、緑、青、の比率は 　わかる　　　　一様分布なら

64 １辺の長さが１の正方形 ?の長さを求めなさい 0.99 ? ?

65 １辺の長さが１の正方形 ?の長さを求めなさい 0.99 Sidak法 ? ?

66 細長い白枠長方形2個の面積が0.01になるとき ? の長さはいくつか？
0.99より大きくする ボンフェリニ法 ? ?

67 問題は、偏りがあること 優性・相加のp値には相関がある

68 分布がわからないので、どうするか 分布がわからないままに、補正する わからない分布を調べてから、それに基づいて補正する

69 分布がわからないままに、補正する

70 分布が違うのに、それでよいのか？ ボンフェロニ・Sidakを使うと 偽陽性が少なくなる パワーが弱くなる ～ストイックであれば大丈夫～ ～保守的であれば大丈夫～

71 分布がわからなければ 分布を調べればよい 正確確率法 ランダマイゼーション・ パーミュテーション法

72 本当に知りたいこと 本当は「関連がない」ときに 相互に相関のある複数の検定を実施したときに
最も小さいp値は、どれくらい小さければ0.01並みに珍しいか

73 本当に知りたいこと 本当は「関連がない」ときに 相互に相関のある複数の検定を実施したときに
最も小さいp値は、どれくらい小さければ0.01並みに珍しいか これは、ちょっと面倒なので、少し変えます

74 本当に知りたいこと 本当は「関連がない」ときに 相互に相関のある複数の検定を実施したときに
『今、観測された分割表』の周辺度数を満足する場合のすべてを考慮して 最も小さいp値は、どれくらい小さければ0.01並みに珍しいか

75 『今、観測された分割表』の周辺度数を満足する場合のすべてを考慮する
２つの方法 本当に「すべての場合」を考慮する 正確確率法 乱数を使って「一部の場合」を考慮して代用する モンテカルロ・ランダマイゼーション法、パーミュテーション法

76 正確確率法とランダマイゼーション法 の違い
『正確』 すべての場合を扱えるのは、自由度2くらいまで。それは、2x3表が１個ある場合。 ランダマイゼーション法 『推定値』 試行ごとに少し違う 1000の場合をやれば、最小p値は0.001、10000回やれば、最小p値は0.0001

77 GWAS基準の有意p値はとても小さいけれど、 それはどうするの？
正確確率法 『正確』 すべての場合を扱えるのは、自由度2くらいまで。それは、2x3表が１個ある場合。 ランダマイゼーション法 『推定値』 試行ごとに少し違う 1000の場合をやれば、最小p値は0.001、10000回やれば、最小p値は0.0001

78 GWAS基準は「デフォルト推奨値」

79 正確確率法って？

80 正確確率法って？ サンプル数が少ないときにカイ二乗検定の代わりに使う方法 サンプル数が少ない、というより、分割表のセルの値が小さいとき…
分割表のセルの値が小さいとき、というより、セルの期待値が小さいとき…

81 どうしてか？ セルの期待値が小さめのときには、カイ二乗検定のp値は『不正確』だから 正確確率検定の方が『保守的』だから
『保守的』なことは、『よいこと』だから

82 たくさんのp値を正確法で得ると… 正確確率検定 カイ二乗検定

83 たくさんのp値を正確法で得ると… 正確検定は、「保守的」なので、その結果をたくさん集めると、一様分布からは随分ずれる

84 連鎖不平衡とp値

85 連鎖不平衡とマルチプルテスティング 1つのSNP 複数の遺伝的モデル 複数の検定 1つの遺伝子 複数のSNP 個々のSNPに１つの検定

86 相互に独立ではない 複数の検定 １つのSNP 複数のSNP 相互に連鎖不平衡にある 相加・優性・劣性の３検定は、相互に独立ではない

87 相互に独立ではない 複数の検定 １つのSNP 複数のSNP 相互に連鎖不平衡にある 相加・優性・劣性の３検定は、相互に独立ではない
マルチプルテスティング補正が必要

88 連鎖不平衡領域 SNPごとに相加検定 20個のp値

89 ケース・コントロールスタディを実施 有意水準0.01で関連ありとするには、どれくらい小さいp値が適当？

90 クイズ、1-6のどれ？ 0.01/20 = 0.0005 0.000502 = 1-(1-0.01)^(1/2) 0.00045
= 1-(1-0.01)^(1/2)

91 クイズ、1-6のどれ？ 0.01/20 = 0.0005 ボンフェロニ 0.000502 = 1-(1-0.01)^(1/2) Sidak

92 クイズ、1-6のどれ？ 0.01/20 = 0.0005 ボンフェロニ 0.000502 = 1-(1-0.01)^(1/2) Sidak

93 ボンフェロニやSidakより小さいわけがない
クイズ、1-6のどれ？ 0.01/20 = ボンフェロニ Sidak ボンフェロニやSidakより小さいわけがない

94 ３つの数字、３つのLD図 0.01/20 = 0.0005 ボンフェロニ 0.0005023906 Sidak 0.00045

95 強 弱 中 中 弱 強 ３つの数字、３つのLD図 0.01/20 = 0.0005 ボンフェロニ 0.0005023906 Sidak
強 弱 中 中 弱 強

96 GWAS基準は「デフォルト推奨値」

97 SNPの数を十倍の 一千万個に増やしたら？

98 連鎖不平衡にあるマーカーで 代用する LDマッピングの原理そのもの SNP 1 : A / a の２アレル
SNP 2 : B / b の２アレル ハプロタイプは４種類 AB Ab aB

99 2SNPの4ハプロタイプは 2x2 分割表 B b A 0.78 0.02 0.8 a 0.18 0.2 1

100 2x2 分割表なら カイ二乗検定しよう カイ二乗値 = 0.81 = r2
B b A 0.78 0.02 0.8 a 0.18 0.2 1

101 LD関係にある 2 SNPのカイ二乗値の相関の良さとLDインデックス r2
カイ二乗値の相関係数 LDのr2

102 LD関係にあるSNPで 代用したときのパワー (r2 = 0.81)
本体 代用

103 LD関係にあるSNPによるパワー 代用SNPの場合 真のリスクSNPの場合

104 p値の高低、どちらが小さい？ 代用SNPの場合 代用SNPのp値の方が小さい 真のリスクSNPの場合

105 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 r2 =1 0.9

106 p値の高低、どちらが「本物」？ 代用SNPの場合 代用SNPのp値の方が小さい 真のリスクSNPの場合

107 p値の高低、どちらが「本物」？ 1 代用のp値が本物以下の割合 LDのr2 0.1 1

108 アレルで、ハプロタイプで 検定する MM Mm mm コントロール 10 30 ケース 5 20 55 15 40 85

109 アレル本数で考える 2x3表の相加モデル(トレンド)検定 2x2表を作って普通にカイ二乗検定(独立性検定)

110 MM Mm mm 10 30 5 20 55 15 40 85 M m 30 60 80 110 170

111 MM Mm mm 10 30 5 20 55 15 40 85 M m 30 60 80 110 170

112 20 24 6 50 29 18 3 49 42 9 100 クイズ ３つの分割表 相加と2x2とが同じ表が2個ある。どれ？ 23 24 6 50 26 18 3 49 42 9 100 25 22 6 50 26 18 3 51 9 100

113 20 24 6 50 29 18 3 49 42 9 100 クイズ ３つの分割表 相加と2x2とが同じ表が2個ある。どれ？ 23 24 6 50 26 18 3 49 42 9 100 25 22 6 50 26 18 3 51 9 100

114 20 24 6 50 29 18 3 49 42 9 100 クイズ ３つの分割表 相加と2x2とが同じ表が2個ある。どれ？ 23 24 6 50 26 18 3 49 42 9 100 25 22 6 50 26 18 3 51 9 100

115 20 24 6 50 29 18 3 49 42 9 100 ハーディ・ワインバーグ平衡かどうか 23 24 6 50 26 18 3 49 42 9 100 25 22 6 50 26 18 3 51 9 100

116 SNPのアレル単位でかんがえるのも ハプロタイプで考えるのも 基本は同じ
SNPのアレルの場合は2x3表の相加モデル(トレンド)検定がある ハプロタイプの方は、ディプロタイプがわからないことが多く、やりようがないかもしれない

117 ハーディ・ワインバーグ平衡のp値

118 HWE検定p値が小さいとき 『サンプルは集団構造化がある母集団を代表している』 『母集団を代表していない』 『実験がうまく行っていない』

119 HWE検定p値が小さいとき 『サンプルは集団構造化がある母集団を代表している』 『母集団を代表していない』 『実験がうまく行っていない』
GWASならば補正方法がある 『母集団を代表していない』 GWASならば個々のHWE検定を問題にする必要はない 『実験がうまく行っていない』 GWASにおいて、個々のHWE検定p値を利用するべきは、これ

120 HWE検定で実験の失敗を疑う HWE検定p値が小さいとき 『サンプルは集団構造化がある母集団を代表している』 『母集団を代表していない』
『実験がうまく行っていない』

121 「ずれ」を見るなら、QQプロット

122 p値が一様じゃない p値は、一様分布に従っているから、その値を0.01と聞けば、「あー、0.01的に珍しいことなんだ」とわかるわけですから、p値の本領は一様分布になっていることです。 しかしながら、実際にGWASを実施して、数十万個のp値を算出して、その分布を見てやると、一様分布になっていない。

123 ２つのアプローチ 「本当は一様分布」なはず。「一様分布」になるように修正してしまおう、という作戦。
一様分布になるわけがない。個々の検定結果のp値を見て、対立仮説が真なのか、帰無仮説が真なのかを選別する情報が得られればよい、という作戦

124 一様分布に修正する作戦 ジェノミック・コントロール

125 単純な１要因 集団が完全には均一でないときに、帰無仮説の検定結果(カイ二乗値、p値)が理論的分布から外れる
その外れ方は、「うまく混ざっていない」という単純な要因で説明できると仮定すると 観測されたカイ二乗値の中央値が理論的な中央値になるように、割り算補正すると解決することが知られている じゃあ、そうしてしまおう、というのがジェノミックコントロール法

126 中央値が揃うように補正

127 たくさんの本物がある場合～FDR～

128 「合否」の基準を 一律にせず 何番目に小さいかで 手加減する

129 色々方法はあるが、基礎的なFDRは 直線に照らして「合否判定」
0.05

130 まとめ 京大 統計遺伝学 p値は、判断するための値 パワーと偽陽性とは、お互い様 たくさんのp値があったら、その特性に応じて補正する
0 から 1 、一様分布 一様分布であることを使って判断したい パワーと偽陽性とは、お互い様 たくさんのp値があったら、その特性に応じて補正する マルチプルテスティング補正 相互に関連しあう検定があったら補正は少し甘くする FDRを使うことも 京大 統計遺伝学