配列解析4 補足資料

IT
配列解析4
補足資料
藤　博幸
BIO
最尤法(MaximumLikelihoodMethod)
現在観察される配列（アライメント）が実現する確からしさ
（尤度）をそれぞれの樹形(topology)ごとに計算し、最大尤度を
与える樹形を選択する方法。
→　どの樹形を仮定すれば、現在の状態（アライメント）が最も
良く説明できるかを計算する。
塩基・アミノ酸置換の確率モデル(Poisson,JTTなど)を使用する
完全探索と発見的探索
系統樹法では、ある基準値（step数、
枝の総和、尤度など）について、nOTU
における全ての樹形を評価し、最良の
ものを選ぶ。これを完全探索
(exhausNvesearch)とよぶ。
しかし、この方法ではOTU数が増える
と爆発的に樹形数が増加するため、
実用的には段階的に高い値をもつ
樹形を見出していく発見的探索
(heurisNcsearch)を使わざるを得ない。
ＯＴＵ数とラベル付き無根
系統樹の種類数との関係
3: 1
4: 3
5: 15
6: 105
7: 945
8: 10395
9: 135135
10: 2027025
11: 34459425
12: 654729075 n: (2n-5)!/{2n-3(n-3)!}
①　PhylogenyTest
Bootstrapを100回
②　SubsNtuNonModel
LGmodel
③　RatesandPaZerns
HasInvariantsites(I)
他はデフォルトで
発見的探索
発見的探索では、基準値（スコア）が計算されている樹形を微小変形し、
スコアがどう変化するかを観察する。スコアが上昇すればその樹形を次回の変形操作
の候補とする（下降しても候補として取り込む方法もある）。
1
3
1
2
7
3
6
4 5
Score=50
7
Score=55
6
（上昇）
2
4 5
3
2
7
1
6
4 5
Score=40
（下降）
よく使用される系統樹の変形操作
•  NNI(NearestNeighborInterchange)
•  SPR(SubtreePruningandRegra_ing)
•  TBR(TreeBisecNonandReconnecNon)
•  TS,TXS(TaXon/TaXaSwap)
NNI(NearestNeighborInterchange)
1
1
2
3
3
2
すぐ隣のOTUとの間の交換。
SPR(SubtreePruningandRegra_ing)
11
1
2
10
3
4
5
6
7
1
2
4
9
8
11
5 6
10
7
9
3
8
TBR(TreeBisecNonandReconnecNon)
11
1
2
5
1
10
9
9
6
7
8
5
6
7
11
3
3
5
4
2
10
4
6
4
8
7
TS,TXS(TaXon[TaXa]Swap)
1
8
2
7
8
2
3
7
4
5
6
3
1
4
5
6
Felsenstein(1981)による塩基配列の最尤法
得られた樹形（トポロジー）のそれぞれについて、
パラメータの最尤推定と最大対数尤度を求める
(1)塩基配列アラインメント
Seq1
Seq2
Seq3
Seq4
C
G
C
T
T
G
T
T
G
G
G
G
G
G
G
G
G
A
G
G
T
T
T
T
G
G
G
G
G
G
C
C
T
T
G
T
1番目のサイトについてとりあえず考える
Felsenstein(1981)による塩基配列の最尤法
(2)次の樹形（無根系統樹）を考える
Seq1
n1
Seq2 Seq3 Seq4
n2 n3
n4 n1-n6はそのノードにおける塩基
(A,C,G,T)
n5
5
6
n6
アライメント中の１つの座位を考える
前ページのアラインメントのサイト１を考える
Seq1 Seq2 Seq3 Seq4
C
G
C
T
d1 d2
d3 d4
n5
5
d5
6
n6
Felsenstein(1981)による塩基配列の最尤法
(3)樹形に基づき尤度 Lを考える
n1-n6はそのノードにおける塩基
Seq1 Seq2 Seq3 Seq4 (A,C,G,T)
C
G
C
T
#は各枝の長さ
d1 d2
d3 d4
n5
5
d5
6
n6
Pij(t)は塩基 iであったものが長さ
tの枝を経て塩基jとなる確率（こ
こで、置換モデルが必要となる）
L1=ΣΣπn5Pn5,C(d1)Pn5,G(d2)Pn5,n6(d5)Pn6,C(d3)Pn5,T(d4)
n5n6
n5,n6は,祖先の塩基で不明なので４種の塩基全ての可能性に
ついて和をとる。 πは塩基組成（アラインメントから推測）
Felsenstein(1981)による塩基配列の最尤法
(4)アライメントの各座位ごとに与えられた樹形の確率を計算。
ある樹形のアライメント全体の尤度Lは、全ての座位（1からn）
の確率の積
L=L1×L2×…×Ln
になるが、通常は尤度の対数をとり、総和を求める。
lnL=∑lnLi　(i=1,n)
今回の例の場合、
L=L1×L2×L3×L4×L5×L6×L7×L8×L9　積による表現は、全てのサイトは独立であることを仮定している。
Pij(t)の表現
Pij(t)をi行、j列の要素とする行列 P(t)を考える
dP(t)
= QP(t)
dt
P(t) = e
Qt
Qは、置換モデル
行列Aの指数関数は次のように実数の場合と同じように
定義される。
€
1 2 1 3
A
e = I + A + A + A + ...
2!
3!
QがJC69の場合
$−3λ λ
λ
λ '
&
)
λ
−3
λ
λ
λ
)
Q = [qij ] = &
& λ
λ −3λ λ )
&
)
λ
λ
λ
−3
λ
%
(
" p0 (t)
€
$
p1 (t)
Qt
$
P(t) = e =
$ p1 (t)
$
# p1 (t)
p1 (t)
p0 (t)
p1 (t)
p1 (t)
p1 (t)
p1 (t)
p0 (t)
p1 (t)
p1 (t)%
'
p1 (t)'
p1 (t)'
'
p0 (t)&
1 3 −4 λt
p0 (t) = + e
4 4
1 1 −4 λt
p1 (t) = − e
4 4
Qでも、P(t)でも、行、列はT,C,A,Gを表す
QがK80の場合
TCAG
T
C
A
G
€
%−(α + 2β )
(
α
β
β
'
*
α
−
α
+
2
β
β
β
(
)
'
*
'
*
β
β
−(α + 2β )
α
'
*
β
β
α
−
α
+
2
β
(
)
&
)
トランジションの置換速度 α
トランスバージョンの置換速度 β
€
K80モデル
" p0 (t)
$
p1 (t)
$
P(t) =
$ p2 (t)
$
# p2 (t)
p1 (t)
p0 (t)
p2 (t)
p2 (t)
p2 (t)
p2 (t)
p0 (t)
p1 (t)
p2 (t)%
'
p2 (t)' dとκの２つのパラメータ
p1 (t)' dは距離、κはトランジション
'
p0 (t)& とトランスバージョンの速度比
1 1 −2 βt−2αt 1 −4 βt 1 1
+ e
+ e
= + e
4 2
4
4 2
-2
d (κ +1)
κ +2
1 1 −2 βt−2αt 1 −4 βtt
− e
+ e
4 2
4
1 1 -2
= − e
4 2
d (κ +1)
κ +2
1 −4 κ d+2
+ e
= p1 (t)
4
+
1
e
4
−4
d
κ +2
= p0 (t)
1 1 −4 β t
− e
4 4
1 1 −4 κ d+2
= − e
= p2 (t)
4 4
Seq1 Seq2 Seq3 Seq4
C
G
C
T
d1 d2
d3 d4
n5
5
d5
6
n6
1 3 −4 λt
p0 (t) = + e
4 4
1 1 −4 λt
p1 (t) = − e
4 4
同じ塩基
異なる塩基
L1=n5n6
ΣΣπn5Pn5,C(d1)Pn5,G(d2)Pn5,n6(d5)Pn6,C(d3)Pn5,T(d4)
lnL=∑lnLi　(i=1~9)
€
λは置換速度、tは時間だが、分離して推定できないので、
λt=dとしてパラメータをLについて最尤法で推定する
また、その時の最大対数尤度 lnLを計算する
通常、最尤法というと
与えられたデータに対して、モデルのパラメータの最適化を、
尤度を最大化することで行う。
ex)100回のコイン投げ　表が47回出た。
コインを１回投げた時に、表が出る確率pを求めたい。
n回中のコイン投げで、表 がm回出る確率は、二項分布
から次のように表現できる
n!
L=
p m (1− p)n−m
m!(n − m)!
"
%
n!
ln L = ln $
' + m ln p + (n − m)ln(1− p)
# m!(n − m)! &
&
∂ln L ∂ #
n!
∂
∂
= ln %
( + m ln p + (n − m) ln(1− p)
∂p
∂p $ m!(n − m)! '
∂p
∂p
m n−m
= −
=0
p 1− p
m n−m
=
p 1− p
(1− p)m = (n − m)p
m = np
m
p=
n
よって例ではpは47/100=0.47
系統樹の場合　コイン投げの場合
観察データ　アライメント　コイン投げの全数　表が出た回数
尤度　与えられた　二項確率
樹形のもと
でアライメント
が観察される
確率
最尤法　尤度に含まれる　二項確率に
枝の長さや　含まれる
使用した進化
パラメータp
モデルなどの
パラメータ
モデル選択　（※　置換モデルではなく樹形のこと）
系統樹の最尤推定では、パラメータの最尤推定
に加えて、様々なトポロジー（系統樹の樹形）に
ついて、最大尤度（あるいは最大対数尤度）が
もとまる。
様々な樹形の中で、最大対数尤度が最大の樹形を
最尤系統樹として選択する
ここで、異なる系統樹を異なるモデルとして
対数尤度を指標として、モデル選択が行われている
最尤法といっているが、系統樹の場合、これが通常の
最尤法と異なっている
樹形1
樹形2
最大対数尤度1 最大対数尤度2
...
樹形n
最大対数尤度n
尤度最大の樹形を最尤推定系統樹として選択
樹形m
※　実際には、選択された
最大対数尤度max 最大対数尤度より低いが、
有意に差のない樹形が複数
選択される場合がある
モデル選択の指標
ln L
パラメータ数の効果を考慮
AIC = − 2 ln L + 2k
2k(k +1)
AICc = −2 ln L +
n − k +1
BIC = - 2 ln L + k ln n
k:パラメータ数
n:サンプル数
最尤法の長所と短所
長所
•  シミュレーションの結果では、多くの場合他の方法より
正しい樹形の推定がなされる
•  統計的基盤があるため、最尤系統樹がそれ以外の樹
形よりどれだけ確からしいかを知ることができる
短所
•  計算時間が非常にかかる
•  置換の確率モデルを変えると結果が変わることがある
サイト間の置換速度の変動
1.2の全てのモデル：配列中の異なるサイトは同じ速度である
ことを仮定
現実的ではない
(1)  各サイトの突然変異率が違う
(2)  各サイトの選択圧の違いから固定される速度が違う
置換のホットスポット、不変サイト
２本の配列のアラインメントはnサイトよりなり、不変サイトを
mサイト含むとする。
2本の配列で差異のあるサイト数をxとする。
p-distance=(x/n)しかし、有効なサイト数は(n-m)なので
x/(n-m)と計算されるべき。　(x/n)　<x/(n-m)
※　速度の違いを考慮しないと配列間距離が過小評価される
ガンマ距離(gammadistane)
任意のサイトの置換速度 r
＝ある確率分布から抽出された確率変数と仮定することで
サイト間の速度の変動を考慮できる。
確率分布としてよく使用されているものは、ガンマ分布
※　Felsenstein(2004)によれば、ガンマ分布を用いる
理由は数学的取り扱いやすさ以上の理由はない
ガンマ分布から得られるモデルは、添字”+Γ”を使って、
JC69+Γ、K80+Γと表現される。また、この時の距離を
Γ距離(gammadistance)とよぶ
※　MEGAでは＋G
Rのコンソール画面が立ち上がる
プロンプロメッセージ“>”の右に処理を書き込む
（１）　apeパッケージの読み込み
library(ape)
（２）　Newick形式のファイルの読み込み
tree <- read.tree(”Z:aaseq2.nwk")
ただし、aaseq2.nwkは、ドキュメントフォルダにおいておくこと
（３）以下を一行ずつ入力してどのような系統樹が描画されるか
を確認せよ
plot(tree)
plot(tree,type="c")
plot(tree,type="r")
plot(tree,type="u")
ファイルへの保存
次の処理を一行ずつ順番に入力して実行せよ
jpeg("aaseq2.jpg")
jpeg形式で保存するファイル名の指定
plot(tree,type="r")
系統樹の描画　typeは他のものを指定しても良い
dev.off()
ファイルへの描画の終了
ドキュメントフォルダにaaseq2.jpgというファイルが
作られている。
ダブルクリックするとファイルが開く
ブートストラップ解析
母集団の性質を、全数調査の代わりに無作為な標本抽出から推定
(randomsampling)
母集団
例えば
母集団：人間の成人男子
身長の平均（母平均）
無作為標本の平均値
によって推定
標本平均は、母集団の真の平均には
一致するとは限らないが、標本平均に
ついてはそれが従う確率分布（正規分布）
がわかっているので、推定値の信頼性
を評価できる
母集団が正規分布　に従っている場合
N( , 2 )
µσ
μの推定値である xは、次の正規分布に従う
2
x ~ N(µ, σ / n)
nが増えれば、推定値の分散は小さくなることがわかる
実際に試してみよう
平均50, 標準偏差10
cat(">> n = 100")
の正規分布に従う
x <- c()
乱数を100個生成
for (i in 1:1000) {
x <- c(x, mean(rnorm(100, mean=50, sd=10)))
}
cat("m = ", 50)
平均値をとる
理論値
cat("sd^2/n = ",100/100)
cat("sm = ",mean(x))
1000個の標本平均xの平均と分散
cat("sv = ", sd(x)^2)
cat("//")
cat(">> n = 100")
平均50, 標準偏差10
の正規分布に従う
乱数を1000個生成
y <- c()
for (i in 1:1000) {
y <- c(y, mean(rnorm(1000, mean=50, sd=10)))
}
cat("m = ", 50)
理論値
cat("sd^2/n = ",100/1000)
cat("sm = ",mean(y))
1000個の標本平均xの平均と分散
cat("sv = ", sd(y)^2)
> cat(">> n = 100")
>> n = 100
> x <- c()
> for (i in 1:1000) {
+
x <- c(x, mean(rnorm(100, mean=50, sd=10)))
+ }
> cat("m = ", 50)
m = 50
> cat("sd^2/n = ",100/100)
sd^2/n = 1
> cat("sm = ",mean(x))
sm = 49.99424
> cat("sv = ", sd(x)^2)
sv = 1.022322
>
> cat(">> n = 100")
>> n = 100
>
> y <- c()
> for (i in 1:1000) {
+
y <- c(y, mean(rnorm(1000, mean=50, sd=10)))
+ }
> cat("m = ", 50)
m = 50
> cat("sd^2/n = ",100/1000)
sd^2/n = 0.1
> cat("sm = ",mean(y))
sm = 49.99273
> cat("sv = ", sd(y)^2)
sv = 0.09781444
推定値は、その確率分布が仮定できるものばかりではない
確率分布を仮定できない推定値に対して、その信頼性を評価したい
ブートストラップ(bootstrap)法
1970年代にEfronらによって提案された方法
確率分布によらない
コンピュータの利用が必要(それまでの数理統計的なアプローチとは異なる方法）
1つの標本から復元抽出を繰り返して
大量の標本を生成し、それらの標本
から推定値 を計算し、母集団の性質
やモデルの推測の誤差などを分析する
無作為標本
母集団
母数θ
サイズn
X=(x1,x2,…,xn)
ブートストラップ標本
サイズn
XB1=(xB11,xB12,…,xB1n)
サイズn
XB2=(xB21,xB22,…,xB2n)
.
.
.
サイズn
^
推定値θ(X)
^
θ(XB1)
XBM=(xBM1,xBM2,…,xBMn)
^
θ(XB2)
.
.
.
^
θ(XBM)
無作為標本からブートストラップ標本の作成
無作為標本から同じサイズを復元抽出
(無作為標本のサイズがnであれば、n個を抽出)
サイズ6
無作為抽出
a
c
b
f d
c
e
a
復元
c
b
f d
e
a
b
f d
c
e
c
6回この処理を繰り返す
無作為標本からブートストラップ標本の作成
無作為標本から同じサイズを復元抽出
(無作為標本のサイズがnであれば、n個を抽出)
無作為標本
ブートストラップ標本
a
c
b
f d
サイズn
c
e
b
f b
c
a
サイズn
ブートストラップ標本
サイズn
^
θ(XB1)
XB1=(xB11,xB12,…,xB1n)
無作為標本
母集団
母数θ
サイズn
サイズn
X=(x1,x2,…,xn)
XB2=(xB21,xB22,…,xB2n)
.
.
.
サイズn
^
推定値θ(X)
^
θ(XB2)
.
.
.
^
θ(XBM)
XBM=(xBM1,xBM2,…,xBMn)
M個のブートストラップ標本それぞれから
推定値を求め、それによりオリジナルの推定値の
誤差や信頼区間を求める
M
θˆB = ∑θˆ ( X Bi )
i=1
1
sB =
θˆ ( X Bi ) − θˆB
M −1
(
)
2
M個の推定値を昇順に並べたM×α番目の値を100α%点とする
ブートストラップで行う意味はないが、先の正規分布(平均50,標準偏差10)から
のサンプルからブートストラップ解析を行ってみる
data <- rnorm(100, mean=50, sd=10)
z <- c()
for (i in 1:1000) {
z <- c(z, mean(sample(data, 100, replace=T)))
}
cat("sm = ", mean(z))
cat("sv = ",var(z))
quantile(data,p=c(0.025, 0.975))
data <- rnorm(1000, mean=50, sd=10)
w <- c()
for (i in 1:1000) {
w <- c(w, mean(sample(data, 1000, replace=T)))
}
cat("sm = ", mean(w))
cat("sv = ",var(w))
quantile(data,p=c(0.025, 0.975))
Rで無作為抽出をするには sample関数を使う
sample関数の引数
sample(x,size,replace=FALSE,prob=NULL)
X
無作為抽出されるデータのベクトル
size
無作為抽出する個数
replace
復元抽出を行うか否か
prob
ベクトルxの要素のそれぞれが抽出される確率
hZp://rplus.wb-nahce.info/rsemi_stat_basic/r_musakuityuusyutu.html　より
> data <- rnorm(100, mean=50, sd=10)
> z <- c()
> for (i in 1:1000) {
+
z <- c(z, mean(sample(data, 100, replace=T)))
+ }
> cat("sm = ", mean(z))
sm = 49.25009
> cat("sv = ",var(z))
sv = 1.115145
> quantile(data,p=c(0.025, 0.975))
2.5%
97.5%
29.29625 69.97043　95% bootstrap信頼区間
>
> data <- rnorm(1000, mean=50, sd=10)
> w <- c()
> for (i in 1:1000) {
+
w <- c(w, mean(sample(data, 1000, replace=T)))
+ }
> cat("sm = ", mean(w))
sm = 50.00343
> cat("sv = ",var(w))
sv = 0.1043676
> quantile(data,p=c(0.025, 0.975))
2.5%
97.5%
30.90386 70.41030
95% bootstrap信頼区間
Rにはbootstrapのためのパッケージbootがある。
bootでブートストラップ解析をする
library(boot)
data <- rnorm(100, mean=50, sd=10)
boot(data, function(x, i)mean(x[i]),R=1000)
data <- rnorm(1000, mean=50, sd=10)
boot(data, function(x, i)mean(x[i]),R=1000)
boot(data, function(x, i)mean(x[i]), R=1000)
# ブーツストラップの設定
１）「data」：
再抽出対象となるデータ名．ここでは"data"．
２）「function(x, i)mean(x[i])」：
誤差評価を行なう統計量の関数指定．ここでは標本平均．
３）「R」：
ブーツストラップ反復回数．ここでは1000回．
> library(boot)
> data <- rnorm(100, mean=50, sd=10)
> boot(data, function(x, i)mean(x[i]),R=1000)
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = data, statistic = function(x, i) mean(x[i]), R = 1000)
Bootstrap Statistics :
original
bias
std. error
t1* 49.06145 -0.02574031
0.9885968
> 0.9885968^2
[1] 0.9773236
> data <- rnorm(1000, mean=50, sd=10)
> boot(data, function(x, i)mean(x[i]),R=1000)
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = data, statistic = function(x, i) mean(x[i]), R = 1000)
Bootstrap Statistics :
original
bias
t1* 49.66522 0.01122756
> 0.3194149^2
[1] 0.1020259
std. error
0.3194149
「ブーツストラップ統計量」として出力されるのは，
元データの平均（original）
ブーツストラップ反復から得られた平均の originalからのバイアス（bias）
ブーツストラップから得られた標本平均の標準誤差（std.error）
である．
hZp://cse.naro.affrc.go.jp/minaka/R/R-resampling.pdfより
Bootstrap解析は、確率分布が仮定できない場合に有効
例えば、
パーセンタイル点
割合
オッズ比
相関係数
などに適用できる
x <- runif(100)
y <- c()
for (i in 1:100) {
y <- c(y, rnorm(1, mean=x[i], sd=0.5))
}
plot(x, y)
相関のある
データの生成
ブートストラップ標本の生成
data <- 1:length(x)
z <- c()
for (i in 1:1000) {
smp <- sample(data, length(x), replace=T)
z <- c(z, cor(x[smp], y[smp]))
}
ブートストラップ標本
cat("orignal =", cor(x, y))
についての相関係数
cat("sm = ", mean(z))
の計算
cat("sv = ", var(z))
quantile(z, p=c(0.025, 0.975))
>x <- runif(100)
> y <- c()
> for (i in 1:100) {
+
y <- c(y, rnorm(1, mean=x[i], sd=0.5))
+ }
> plot(x, y)
>
> data <- 1:length(x)
> z <- c()
> for (i in 1:1000) {
+ smp <- sample(data, length(x), replace=T)
+ z <- c(z, cor(x[smp], y[smp]))
+ }
> cat("orignal =", cor(x, y))
orignal = 0.506939
> cat("sm = ", mean(z))
sm = 0.5008048
> cat("sv = ", var(z))
sv = 0.00513788
> quantile(z, p=c(0.025, 0.975))
2.5%
97.5%
0.3486814 0.6270728