Transcript スライド - 京都大学

生命情報科学人材養成コンソーシアム
バイオインフォマティクス実習コース
遺伝医学
2010年10月29日（金曜日）
京都大学大学院医学研究科
附属ゲノム医学センター
統計遺伝学分野
教授
山田 亮
© 2010 山田 亮
本日の構成
「遺伝統計学」
•
•
•
•
•
•
遺伝統計学概観
遺伝学基礎
遺伝統計学でのデータの扱い
分割表検定
連鎖解析
たくさんの検定
目標
• データ解析をする立場で、研究対象を捉える
ことにイメージを持つ
• データを眺めるスタイルを理解する
• 既存の手法や用語に対して、自分の理解の
役に立つのかどうかという視点で客観的(批
判的)に対する
内容は、受講者の希望に応じて変更しながら進みます
せっかくなので知りたいことを知って帰りましょう
受講の動機は？
個別研究？/研究全般のアシスト？/
”遺伝統計学”研究?/シリーズ受講なのでついでに?
出身・背景は？
生物・遺伝系？/統計・計算機系？
今日の教材の元ネタ
遺伝統計学の基礎
Rによる遺伝因子解析・遺伝子
機能解析
ISBN 978-4-274-06822-5
Download Sources of R
http://www.genome.med.kyoto-u.ac.jp/wiki_tokyo/index.php/StatGenetOhm
※ 使用するRコードの入手とその使用法
-"http://www.genome.med.kyotou.ac.jp/wiki_tokyo/index.php/StatGenetOhm"から、
-Rsrc.zipファイルをダウンロードして解凍
-Rの作業ディレクトリを解凍してできたRsrcディレクトリに指定
する
-"Rコードのソースを読み込み"で、"StatGenetDemo.R"を読み
込むと、デモ開始
-必要なパッケージのインストールが行われ、Rsrcディレクトリ
内のRソースを順に実行する
-ソースファイルを順送りするときと、表示図を確認するための
プロンプトが現れるので、"Enter"キー等でデモを進める
遺伝統計学概観
似ている 似ていない
Alike, not alike
離散的
４カテゴリ：Ａ、Ｔ、Ｇ、Ｃ
長さ 30億塩基対
tattgcccac ggccctccag ccaagaagaa atccacaggt tcctccacat ggcccctgga
ccctggggta gaggtgaccc tgacgatgaa agtggccagt ggtagcacag gcgaccagaa
ggttcagatt tcatactacg gacccaagac tccaccagtc aaagctctac tctacctcac
cggggtggaa atctccttgt gcgcagacat cacccgcacc ggcaaagtga agccaaccag
agctgtgaaa gatcagagga cctggacctg gggcccttgt ggacagggtg ccatcctgct
ggtgaactgt gacagagaca atctcgaatc ttctgccatg gactgcgagg atgatgaagt
gcttgacagc gaagacctgc aggacatgtc gctgatgacc ctgagcacga agacccccaa
ggacttcttc acaaaccata cactggtgct ccacgtggcc aggtctgaga tggacaaagt
gagggtgttt caggccacac ggggcaaact gtcctccaag tgcagcgtag tcttgggtcc
caagtggccc tctcactacc tgatggtccc cggtggaaag cacaacatgg acttctacgt
ggaggccctc gctttcccgg acaccgactt cccggggctc attaccctca ccatctccct
gctggacacg tccaacctgg agctccccga ggctgtggtg ttccaagaca gcgtggtctt
ccgcgtggcg ccctggatca tgacccccaa cacccagccc ccgcaggagg tgtacgcgtg
cagtattttt gaaaatgagg acttcctgaa gtcagtgact actctggcca tgaaagccaa
gtgcaagctg accatctgcc ctgaggagga gaacatggat gaccagtgga tgcaggatga
aatggagatc ggctacatcc aagccccaca caaaacgctg cccgtggtct tcgactctcc
aaggaacaga ggcctgaagg agtttcccat caaacgcgtg atgggtccag attttggcta
tgtaactcga gggccccaaa cagggggtat cagtggactg gactcctttg ggaacctgga
agtgagcccc ccagtcacag tcaggggcaa ggaatacccg ctgggcagga ttctcttcgg
ggacagctgt tatcccagca atgacagccg gcagatgcac caggccctgc aggacttcct
cagtgcccag caggtgcagg cccctgtgaa gctctattct gactggctgt ccgtgggcca
cgtggacgag ttcctgagct ttgtgccagc acccgacagg aagggcttcc ggctgctcct
ggccagcccc aggtcctgct acaaactgtt ccaggagcag cagaatgagg gccacgggga
ggccctgctg ttcgaaggga tcaagaaaaa aaaacagcag aaaataaaga acattctgtc
aaacaagaca ttgagagaac ataattcatt tgtggagaga tgcatcgact ggaaccgcga
染色体伝達図
系統樹
遺伝統計学の位置づけ
The Lady Tasting Tea: How Statistics
Revolutionized Science in the Twentieth
Century [ペーパーバック]
David Salsburg
ISBN-13: 978-0805071344
統計学を拓いた異才たち―経験則から科学
へ進展した一世紀 [単行本]
デイヴィッド サルツブルグ (著), David S.
Salsburg (原著), 竹内 惠行 (翻訳), 熊谷 悦
生 (翻訳)
ISBN-13: 978-4532351946
遺伝学基礎
•
•
•
•
血縁関係 → 似る
似る程度にばらつき
似る特徴もあれば似ない特徴もある
似るか似ないかには「理由」がある
遺伝的に伝わることと多様であること
Heredity and Variation
• 何が伝わり、何が多様か？
• “もの”と”こと”
• 遺伝子型と表現型
ジェノタイプとフェノタイプ
フェノタイプ分類のこつ
ジェノタイプ
染色体
遺伝的多様性の３要素
• 変異
• 交叉・組換え
• 遺伝的浮動
遺伝子多型・変異
遺伝子多型
構造分類
遺伝子多型
サイズ分類
遺伝子座 アレル
ハプロタイプ ディプロタイプ
交叉
組み換え
• 起源の一致(ＩＢＳ)は、２本の染色体ができる
までの交叉の回数の偶奇で決まる
交叉回数
ポアッソン分布
R2-6.R
#R2-6.R
# 可能箇所すべてで交叉がおきるかどうかを試す方法
RecombSim<-function(L=10000,r=0.001,Niter=1000){
# Lは配列長,rは箇所あたりの交叉確率,Niterはシミュレー
ション試行回数
# 行数Niter、列数L-1箇所の行列にする
m<-matrix(rbinom((L-1)*Niter,1,r),nrow=Niter)
apply(m,1,sum)
}
# ポアッソン分布を使う方法
RecombPois<-function(L=10000,r=0.001,Niter=1000){
rpois(Niter,(L-1)*r)
# rpois() 関数については付録A.5 確率分布関数・疑似乱数
列の発生を参照
}
#R2-6.R(続き)
L<-10000;r<-0.0001;Niter<-1000
NumSim<-RecombSim(L=L,r=r,Niter=Niter)
NumPois<-RecombPois(L=L,r=r,Niter=Niter)
ylim<-c(0,max(NumSim,NumPois))
plot(ppoints(Niter,a=0),sort(NumSim),ylim=ylim,col=gray(6/8))
par(new=T)
plot(ppoints(Niter,a=0),sort(NumPois),type="l",ylim=ylim)
交叉間距離
指数分布
R2-7.R
#R2-7.R
Niter<-1000 # シミュレーション回数
L<-1000000 #染色体の長さ
r<-0.0001 #塩基間あたりの交叉確率
# 交叉箇所数をポアッソン分布からの乱数で指定し、交叉箇所をsample関数で指定する
crosses<-sort(sample(1:(L-1),rpois(1,(L-1)*r),replace=FALSE))
# 交叉間距離のベクトルを作る
A<-c(0,crosses) # 染色体の始点と交叉箇所のベクトル
B<-c(crosses,L) # 交叉箇所と染色体の終点のベクトル
C<-B-A #交叉間距離のベクトル
# 平均がmean(C)の指数分布からの乱数をlength(C)の数だけ発生させてプロット
rexps<-rexp(length(C),1/mean(C))
# rexp() 関数については付録A.5 確率分布関数、疑似乱数列の発生を参照
# 交叉間距離をソートしてプロット
ylim<-c(0,max(C,rexps))
plot(sort(C),ylim=ylim,cex=0.5,pch=15) #交叉間距離の昇順プロット
par(new=T)
plot(sort(rexps),col="red",ylim=ylim,type="l") # 指数分布からの乱数の昇順プロット
遺伝的浮動
#”a”と”b”の集団を作る
Na<-1;Nb<-11;k<-4
Ns<-Na+Nb
A<-c(rep("a",Na),rep("b",Nb))
A
#k倍する
B<-rep(A,k)
#Ns個抜き取る
sample(B,Ns)
Niter<-1000
#Niter回、繰り返して
Numa<-rep(0,Niter)
for(i in 1:Niter){
S<-sample(B,Ns)
Numa[i]<-length(which(S=="a"))
多型性を失う
}
ドリフトアウト
#"a"の抜き取られ数の分布をみる
hist(Numa)
#"a"が２個抜き取られた確率は？
length(which(Numa==2))/Niter
#n!を使って計算してみる
gamma(12+1)/(gamma(2+1)*gamma(10+1))*gamma(36+1)/(gamma(2+1)*gamma(34
+1))*gamma(4+1)*gamma(44+1)/gamma(48+1)
#"a"の抜き取り数は、0,1,2,3,4のいずれか
#その確率を計算して足し合わせてみる
P0<gamma(12+1)/(gamma(0+1)*gamma(12+1))*gamma(36+1)/(gamma(4+1)*gamma(32+1)
)*gamma(4+1)*gamma(44+1)/gamma(48+1)
P1<gamma(12+1)/(gamma(1+1)*gamma(11+1))*gamma(36+1)/(gamma(3+1)*gamma(33+1)
)*gamma(4+1)*gamma(44+1)/gamma(48+1)
P2<gamma(12+1)/(gamma(2+1)*gamma(10+1))*gamma(36+1)/(gamma(2+1)*gamma(34+1)
)*gamma(4+1)*gamma(44+1)/gamma(48+1)
P3<gamma(12+1)/(gamma(3+1)*gamma(9+1))*gamma(36+1)/(gamma(1+1)*gamma(35+1))
*gamma(4+1)*gamma(44+1)/gamma(48+1)
P4<gamma(12+1)/(gamma(4+1)*gamma(8+1))*gamma(36+1)/(gamma(0+1)*gamma(36+1))
*gamma(4+1)*gamma(44+1)/gamma(48+1)
P0+P1+P2+P3+P4
酔歩 Random walk
R16-sup1.R
#R16-sup1.R
nstep<-100
rwalk<-matrix(0,nstep,2)
rtheta<-rnorm(nstep-1)
stepx<-cos(rtheta)
stepy<-sin(rtheta)
for(i in 1:(nstep-1)){
rwalk[i+1,1]<-rwalk[i,1]+stepx[i]
rwalk[i+1,2]<-rwalk[i,2]+stepy[i]
}
plot(rwalk,type="l")
拡散
差が大きいときには急速に
差が小さいときにはゆっくりと
差と「時間微分」が比例
拡散方程式
R6-4.R
差が大きいときには急速に
差が小さいときにはゆっくりと
差と「時間微分」が比例
拡散方程式
# R6-4.R
# pa,pb:２集団の人口,d:単位時間あたりの移住人数,t:世代
pa<-9000
pb<-1000
d<-100
t<-0:100
fa<-(pa+pb*exp(-2*d*(pa+pb)/(pa*pb)*t))/(pa+pb)
fb<-(pa*(1-exp(-2*d*(pa+pb)/(pa*pb)*t)))/(pa+pb)
plot(t,fa,ylim=c(0,1),type="l",xlab="t",ylab="frequency")
par(new=T)
plot(t,fb,ylim=c(0,1),type="l",xlab="t",ylab="frequency")
アレルの行く末
どちらかのアレルのみになる
アレル頻度
R2-8.R
片方の
アレル
に固定
時間
低頻度なが
ら、こちら側
にも固定
遺伝統計学でのデータの扱い
カテゴリカルデータ
• 全体を覆う
• 重なりがない
カテゴリと順序
量的形質とカテゴリ
カテゴリと次元
カテゴリカルデータの表現
分割表
自由度
自由度 記述するための変数の数
２ｘ２ table
いくつの数値で表を説明するか？
テーブルを説明するのに
必要な変数の数
(自由度)
自由度 １
有意性を判断する
• 説明をするのに、変数を増やすか増やさない
か、それが問題
珍しさで順序をつける
• δが大きさで順序をつける
• δが大きいほど「珍しい」
独立を仮定したとき
順列(パーミュテーション)
項目１
項目２
○１
■１
○２
■２
○３
□３
●４
項目１
項目２
○１
□５
○２
□４
○３
■１
□４
●４
□３
●５
□５
●５
■２
●６
■６
●６
■６
項目１と項目２
のラベル付け
は独立だと仮
定すれば
並べ替えない 並べ替える
順列(パーミュテーション)
N個の並べ替え：N!
独立を仮定したとき
順列(パーミュテーション)
項目１
項目２
○１
■１
○２
■２
○３
□３
●４
項目１
項目２
○１
□５
○２
□４
○３
■１
□４
●４
□３
●５
□５
●５
■２
●６
■６
●６
■６
項目１と項目２
のラベル付け
は独立だと仮
定すれば
順列(パーミュテーション)
N個の並べ替え：N!
並べ替えごとに、分割表ができる
δの分布がとれる
珍しいδとありきたりなδがわかる
観測データのδのありきたりな程度を
0-1の値で表わす
完全パーミュテーション
と
モンテカルロパーミュテーション
• 全順列(N!)は膨大
– しらみつぶしにするか
– 1部を使うか(乱数を使う：モンテカルロ)
パーミュテーションの網羅と
サンプリング
library(gtools)
n<-3
permutations(n,n)
permutations(n,n,repeats=TRUE)
sample(1:n,n)
パーミュテーション(順列)の計算
パーミュテーション(順列)の足し合わせ
• N!は膨大
• モンテカルロをするにしても膨大
– モンテカルロの回数とパーミュテーションテストの
ｐ値
• 便法はないか？
便法
• ないこともある
便法
• ないこともある
• 分割表にはある
分割表の正確生起確率
珍しさの計算
G1,G2をn..人に割り付ける場合は
n..!/(n1.! n2.!)通り
A,aをn..人に割り付ける場合は
n..!/(n.1! n.2!)
“G1-A”,”G1-a”,”G2-A”,”G2-a”をn..
人に割り付ける場合はn..!/(n11!
n12!n21!n22!)
正確確率と完全パーミュテーションの
結果は同じ
• ラベルの枚数を考慮
するかしないか
– しない
• パーミュテーション
– する
• 正確確率計算
δと生起確率
生起確率
δ
R13-2.R
正確確率の計算
正確確率の足し合わせ
• 計算自体が面倒臭い
• 足し合わせるためには、「膨大な」観測可能テ
ーブルのすべてについて計算が必要
– 面倒臭い
– 自由度が上がると非現実的
• 便法はないか？
正確確率のプロットを近似できる関数があれば・・・
生起確率
δ
R13-2.R
２つの便法
• 確率
– ピアソンの独立性検定のカイ自乗値
• 尤度
– 尤度比検定
正確生起確率とはなんだったか
• 「ある仮説(帰無仮説)」の下で
• 「あるデータ」が観察される
• 確率
正確生起確率はなんだったか
• 「ある仮説(帰無仮説)」の下で
•
•
•
•
「あるデータ」
「別のデータ」
…
「すべてのデータ」が観察される
• 確率
→ 確率分布に照らす(カイ自乗検定)
正確生起確率はなんだったか
•
•
•
•
「ある仮説(帰無仮説)」
「別のある仮説(対立仮説)」
…
「すべての仮説(対立仮説と帰無仮説)」の下
で
• 「あるデータ」が観察される
• 確率
→ 仮説に基づく確率(尤度)の比較
確率と尤度
Probability and likelihood
仮説を固定、観察を動かす
確率：G1,G2に差がないときに
n11=x (x=0,1,2,…)という観察を
する確率
仮説を動かす、観察を固定
尤度：G1ではAの割合がp1でG2
ではAの割合がp2であるという
仮定のもとでn11=n11という観察
をする確率(p1=0～1,p2=0～1)
observation
hypothesis
ピアソンの独立性検定のカイ自乗値
独立仮説でのデータ観測確率を知る
尤度による方法
• 確率pの事象がn回続けて起きる確率
– p^n
• 確率pの事象がn回続けて起きて、そのあとm
回続けて起きない確率
– p^n x (1-p)^m
• n+m回のうちn回起きてm回起きない確率
– C x p^n (1-p)^m
– C: (n+m)からnを取り出す場合の数
２ｘ２表の尤度
A
a
All
A
a
All
G1
N1A
N1a
N1
G1
P(1A)
P(1a)
1
G2
N2A
N2a
N2
G2
P(2A)
P(2a)
1
Gall
NA
Na
T
Gall
P(A)
P(a)
1
ピアソンのカイ自乗検定の自由度
尤度比検定の自由度
• ピアソン
– 表の中身を変更するときに、自由に動かせるセ
ルの数
• 尤度比検定
– 比較する２つの仮説で使う変数の差
• ２ｘ２表の場合はなぜ、自由度が１？
検定３種 Three types of tests
– 正確確率検定 Exact tests
• パーミュテーションテスト Permutation-based
• テーブルの正確生起確率による Exact Probability
based on table
– ピアソンの独立性検定 Pearson's independence
test
– 尤度比検定 Likelihood ratio test
だいたい同じ 少し違う
Similar each other but a bit different
漸近近似
Nが大きいとほとんど一緒
カテゴリカルデータの検定
• 検定の手法
– Ｐ値に持ち込む形式の違いでわける
• 正確確率に基づく検定
• 期待値表からの距離に基づく検定
• 尤度に基づく検定
カテゴリと順序
量的形質とカテゴリ
カテゴリと次元
ディプロタイプというカテゴリ
A
a
A
AA
Aa
a
aA
aa
• AaとaA
– Aを父からaを母から
– Aを母からaを父から
– アレルの並び順を区別するかしな
いか、と、エピジェネティック効果
• AA, Aa(=aA),aaの３つ
ディプロタイプというカテゴリ
• ＳＮＰの場合
– ３タイプが対等
– ３タイプに順序(アレルについて評価)
• アレルの本数の力
– 相乗的
– 相加的
– それ以外
» 優性 劣性 ヘテロが特別
多アレルの場合
• 順序あり
– リピート数など（タンデムリピート、ＣＮＶ)
• 順序なし
– ハプロタイプ
順序あり・多アレルの
ディプロタイプの場合
0
1
2
0
(0,0)
(0,1)
(0,2)
1
(1,0)
(1,1)
(1,2)
2
(2,0)
(2,1)
(2,2)
順序をつける
足し算・・・
(0,0)=0, (0,1)=(1,0)=1,…,(2,2)=2
順序なし(かもしれない)多アレル
ハプロタイプというカテゴリ
• ＳＮＰのアレルの組み合わせ
– 組み合わせは多次元
– 多次元のまま扱う
– 順序を入れられる？
• ディプロタイプにすると？
– Na+Na(Na-1)/2 カテゴリ
普通の表
サンプル
Ｇ
Ｐ
S1
G1
P1
S2
G1
P2
S3
G2
P2
S4
G1
P1
S5
G3
P1
…
…
…
Sn
G2
P2
P1
P2
G1
N11
N12
G2
N21
N22
G3
N31
N32
普通の表
• 2x2, 2x3, …, NxM 表
– どれも同じ
– 自由度は(N-1)(M-1)
• 自由度が大きくなると正確確率検定は大変すぎる
正確確率検定の重さ
#行列を作る
m<-matrix(c(10,20,30,40,50,60),nrow=2,nrow=3)
m
#作る行列の値の納め方を変える
m<-matrix(c(10,20,30,40,50,60),nrow=2,nrow=3,byrow=TRUE)
m
fisher.test(m)
#大きな表、大きな自由度
nrow=4;ncol=4
m<-matrix(round(runif(nrow*ncol)*10,0),nrow=nrow,ncol=ncol)
m
fisher.test(m)
普通の表の検定
• 検定の手法
– Ｐ値に持ち込む形式の違いでわける
• 正確検定・ピアソン・尤度比検定
普通の表の検定
• 検定の手法
– Ｐ値に持ち込む形式の違いでわける
• 正確検定・ピアソン・尤度比検定
– 表の形とカテゴリの順序関係の違いでわける
普通の表の検定
• 検定の手法
• 確率・尤度
– 生起確率を元に、それをまねる確率分布(カイ自乗のような)を使う方
法
• 傾向性検定、ＫＷ、ＪＴ
– 比較するべき仮説(帰無と対立)を変数でモデル化する方法(尤度比
検定のような)
• 線形回帰・ロジスティック・それ以外のもっといろいろな回帰(一般化線形回帰とか)
普通の表
サンプル
Ｇ
Ｐ
S1
G1
P1
S2
G1
P2
S3
G2
P2
S4
G1
P1
S5
G3
P1
…
…
…
Sn
G2
P2
P1
P2
G1
N11
N12
G2
N21
N22
G3
N31
N32
普通ではない表
• ハーディ・ワインバーグ平衡検定を例に
ハーディワインバーグ平衡(HWE)とは
AA
Aa
aa
ALL
N(AA)
N(Aa)
N(aa)
N(ALL)
ランダムメイティング
HWEの期待値表
AA
Aa
aa
ALL
N(AA)
N(Aa)
N(aa)
N(ALL)
E(AA)=P(A)^2 N(ALL)
E(Aa)=2P(A)P(a)
N(ALL)
E(aa)=P(a)^2 N(ALL)
N(ALL)
• P(AA)=P(A) ^2
• P(Aa)=2P(A)p(a)
• P(aa)=P(a)^2
ＨＷＥと近交係数
• N(AA)/N(ALL) = E(AA) + f p(A)p(a)
• N(Aa) /N(ALL)= E(Aa) – 2f p(A)p(a)
• N(aa)/N(ALL) = E(aa) + f p(A)p(a)
• f:近交係数 (HWEの指数)
– f=0 : HWE
– f=1 : 全員ホモ接合体
HWEの観察表からカイ二乗検定
AA
Aa
aa
ALL
N(AA)
N(Aa)
N(aa)
N(ALL)
• アレル本数
• M(A)=(2N(AA)+N(Aa))/2
• M(a)=(N(Aa)+2(aa))/2
• アレル頻度 (の推定値)
• P(A)=M(A)/(2 N(ALL))
• P(a)=M(a)/(2 N(ALL))
HWEの期待値からのずれ
AA
Aa
aa
ALL
N(AA)
N(Aa)
N(aa)
N(ALL)
E(AA)=P(A)^2 N(ALL)
E(Aa)=2P(A)P(a)
N(ALL)
E(aa)=P(a)^2 N(ALL)
N(ALL)
• カイ自乗値
• 自由度:1 (自由変数:f)
•N(AA)/N(ALL) = E(AA)+ f p(A)p(a)
•N(Aa) /N(ALL)= E(Aa) – 2f p(A)p(a)
•N(aa)/N(ALL) = E(aa) + f p(A)p(A)
ＨＷＥ検定カイ自乗値と近交係数
• N(AA)/N(ALL) = E(AA) + f p(A)p(a)
• N(Aa) /N(ALL)= E(Aa) – 2f p(A)p(a)
• N(aa)/N(ALL) = E(aa) + f p(A)p(a)
• f:近交係数 (HWEの指数)
– f=0 : HWE
– f=1 : 全員ホモ接合体
– f^2 x N(ALL)=カイ自乗値
Chi^2 = N f^2
普通ではない表 HWEの表
AA
Aa
aa
ALL
N(AA)
N(Aa)
N(aa)
N(ALL)
ＨＷＥ正確確率検定
• 正確率検定はできる
– データの生起確率を計算できて
– すべての場合を網羅できればよ
い
ＨＷＥ正確確率検定
• 正確率検定はできる
– データの生起確率を計算できて
R4-1.R
ＨＷＥ正確確率検定
• 正確率検定はできる
– すべての場合を網羅できればよ
い
• ヘテロの人数は奇数か偶数かの
どちらか
R4-1.R
ハーディ・ワインバーグ平衡
↓
連鎖平衡
アレル関連 連鎖不平衡 連鎖平衡
• ハプロタイプ
アレル関連
数を確認
• アレル数が２の多型
– k個
• ハプロタイプ種類数
– 2^k
• ディプロタイプ種類数
– 3^k
アレルの関係が独立とは
• 独立は直交
– 内積がゼロ～相関がゼロ
• 独立は掛け算
ハプロイド：２ローカスの独立
• ハプロイド
– 同一染色体上
– 同一配偶子(精子・卵子)内
LE状態
• 2SNP haplotype ((1)AB,(2)Ab,(3)aB,(4)ab)
– H1 = pA pB
– H2 = pA pb
– H3 = pa pB
– H4 = pa pb
２ローカスが独立でなくなるとき
• 同一分子上にあるから、自由でない
• 交差してよければ・・・
• 遠ければ遠いほど
– 完全な自由：無限遠：異なる染色体
２ローカスが独立でなくなるとき
• 「無限」に遠くても不自由
– メイティングが不自由(HWD,集団の構造化)
LE状態からのずれを表す
• 2SNP haplotype ((1)AB,(2)Ab,(3)aB,(4)ab)
– H1 = pA pB +d
– H2 = pA pb -d
– H3 = pa pB -d
– H4 = pa pb + d
– d=r SQRT(pA pa pB pb)
– r^2 : LD 指標
HWD と LD
MM : p(M)^2 +f p(M)p(m)
Mm : 2p(M)p(m) – fp(M)p(m)
mm : p(m)^2 +f p(M)p(m)
H1 = pA pB +d
H2 = pA pb -d
H3 = pa pB -d
H4 = pa pb + d
MM : p(M)p(M) + d
Mm : p(M)p(m) – d
mM : p(m)p(M) – d
mm : p(m)p(m) + d
d= f sqrt(p(M)p(m)p(M)p(m))
H1 = pA pB + d
H2 = pA pb – d
H3 = pa pB – d
H4 = pa pb + d
d=r sqrt(pA pa pB pb)
HWE LE
指数と統計量
MM : p(M)p(M) +d
Mm : p(M)p(m) – d
mM : p(m)p(M) – d
mm : p(m)p(m) + d
d= f sqrt(p(M)p(m)p(M)p(m))
H1 = pA pB +d
H2 = pA pb -d
H3 = pa pB –d
H4 = pa pb + d
d=r sqrt(pA pa pB pb)
r^2 : LD index
N : No. samples
Chi^2 = N f^2
Chi^2 = N r^2
r : 相関係数
LDは染色体上に広がる
染色体上に広がるLE/LDの具合
•
•
•
•
K個のマーカー
ペアワイズな比較
K^2 な分布
K(K-1)/2 ペア
• Mx とMy (x , yは異なる) との関係(K(K-1)/2 ペア)
• Mx とMx (自身) との関係
• 分散・共分散行列
(K^2の要素)
多数のＳＮＰ
ペアについてのr^2から絵を描く
R5-2.R
R5-3.R
m<-matrix(rbinom(120,1,0.5),20,6)
heatmap(m)
cormatrix<-cor(m);rsqmatrix<cormatrix^2
image(1:nrow(rsqmatrix),1:ncol(rsqmatri
x),rsqmatrix,col=gray((100:0)/100))
ペアワイズLDプロット
• ペアワイズ情報
を用いる
• マーカーの位置
を変えない
• パターンから染
色体に沿った交
叉・組み換えの
様子を読み取る
LDプロットができるまで
M
1
M
2
M
3
M
4
…
N
1
0
0
1
1
…
N
2
O
1
1
1
…
N
3
1
0
0
0
…
N
4
1
1
1
0
…
N
5
0
1
0
1
…
N
6
1
0
1
1
…
…
NxM
MxM
NxMのデータの処理
ペアワイズな関係
多数の要素のペアワイズ関係の図示
要素の順番を並びかえれば、それは
階層的クラスタリング
R5-1.R
階層的クラスタリングは
樹形図の１方法
クラスタリング
• 階層的クラスタリング
• 非階層的クラスタリング
• クラスタ化する方法の違い
• クラスタ化して残す情報の違い
『階層的』クラスタリング
•
•
•
•
樹形図化
ペアワイズ比較情報
５つの要素(葉)が
３つの結合点(節)を介し
て
• ７本の線(辺)で結びつけ
られている
• 葉への要素の割り当てパ
ターンと辺の長短がクラ
スタを定めている
R6-2.R
『非階層的』クラスタリング R6-3.R
• 要素の数だけ点があ
る
• 点が空間に配置され
ている
• 位置の情報が、帰属
クラスタを決めている
類別してプロットする
主成分分析 PCA
分散の分解
Decomposition of variance into pieces
R7-5.R
適切な軸
Appropriate axes
固有値分解・主成分分析
(PCA)
R7-5.R
• 正規直交基底
• どうして「直交」
• 分散が基底成分の分散に
分解できるから
エッセンス 連鎖解析
木の解析
尤度の解析
家系図 Pedigree
家系の観察
• マーカーＸに想定される木
• マーカーＹに想定される木
• 表現型に想定される木
木の解析
尤度の解析
• マーカーのジェノタイプの伝達木
• 個体のフェノタイプの伝達木
• 異なる伝達木を一致させるためのねじれ
– 交叉
• 交叉の起きやすさ(尤度)は距離の関数
• 面倒くさいのは、木の場合分けが多いから
伝達の木
たくさんあるけれど、数えられる
教科書 第10章
• マーカーＸの木(Tx)、表現
型の木(Tp)、マーカーＹの
木(Ty)が異なるときには、
間で『組換え』が起きたと
考える
• 距離が遠ければ、組換え
が起きやすい
• 木のパターンによって、
Tx,Tp,Tyの遠近関係として「
ありそうな位置」が推定で
きる
マルチプルテスティング
たくさんの検定 教科書 第17章
均一分布
•
•
•
•
•
N<-1000
X<-runif(N)
plot(X)
plot(sort(X))
plot(ppoints(N,a=0),sort(X))
均一分布からのk個の乱数の最小
値はどれくらい小さいか
•
•
•
•
•
•
N<-10000
k<-5
Xs<-matrix(runif(N*k),N,k)
Mins<-apply(Xs,1,min)
mean(Mins)
1/(k+1)
• kの値を振ってやってみる
１番小さいｐ値の期待値
• k個の独立なテスト
•
1/(k+1)
i番目に小さいｐ値の期待値?
•
•
•
•
N<-1000
k<-10
Xs<-matrix(runif(N*k),N,k)
Sorted<matrix(apply(Xs,1,sort),N,k,byrow=TRUE)
• plot(apply(Sorted,2,mean))
• kの値を振ってやってみる
期待値は1/(k+1)
分布は？
•
•
•
•
•
•
•
N<-1000
k<-100
Xs<-matrix(runif(N*k),N,k)
Mins<-apply(Xs,1,min)
mean(Mins)
plot(density(Mins))
abline(v=mean(Mins))
このグラフをどう読む
• plot(sort(Mins))
• abline(h=mean(Mins))
i番目に大きいｐ値の期待値はi/(k+1)
分布は？
• doubleSorted<-apply(Sorted,2,sort)
matplot(doubleSorted,type="l")
独立なp値
• Bonferroni’s method
– すべての検定を「排他的」に考える
– もっとも保守的
– p_corrected= kp
• Sidak’s family-wise error rate
– すべての検定を「独立」と考える
– p_corrected=1-(1-p)^k
k<-100000
plot(1:k,1/k*(1:k),ylim=c(0,1),type="l")
par(new=TRUE)
plot(1:k,1-(11/k)^(1:k),ylim=c(0,1),col="red",type="l")
実際には、すべての検定は排他
的でもなく、独立でもない
• パーミュテーションで、「全通り」を調べる
• モンテカルロ・パーミュテーションで「全通り」
に近い分布を調べる
対立仮説が正しいとき
• 非心カイ自乗分布
• カイ自乗分布に１変数追加
強弱いろいろな対立仮説が成り立つとき
• たとえば、集団構造化があるときのＧＷＡＳ
• 「集団構造化」を１追加変数でモデル化すれ
ば、観察データから、その変数の値を推定す
ることができて、その推定値に基づいて補正
ができる
たくさんの対立仮説が正しいとき
• 1変数でモデル化するほど単純ではない
大雑把に「正しい対立仮説の数」が分かって
いる(と仮定する)
• False Discovery Rate（FDR)