系統樹を最尤推定する

目次

・系統樹を推定する

・突然変異と系統分岐

・塩基置換モデル

　　・JC68における塩基置換関数の導出

・系統樹の尤度関数

・系統樹の最尤推定

　　・系統樹の最尤推定は時間がかかる

・Rで最尤推定してみる

・補足1：より一般的な塩基置換モデルをどう解くか

・補足2：推移行列の固有値と正規直交固有ベクトル

系統樹を推定する

　私たち生物学者たちにとって、生物の進化というのは興味深いトピックである。進化生物学の中にも、いろいろな研究の分野があるが、その一つのテーマが、過去から現在までにかけて、どのように種が分かたれてきたか、を推定するというものである。この種が分かたれてきた歴史を図として表したものを系統樹phylogenetic treeと呼ぶ。例えば、下記なら、種1～3の最も近い共通祖先Most recent common ancestorから、まず種1と2のグループと種3に分かれ、その後に種1と種2が分かれている。

では、このような系統樹はどのような計算で得られるのだろうか。系統樹はいうなれば、「近い」種間は短い枝で、「遠い」種間は長い枝でつないだ図である。ゆえに、種の「近さ・遠さ」を定義する必要がある。

　その定義の仕方は複数の種類がある。第1に「何を近さを測るためのデータ」とするか、によって複数の種類がある。使われるデータには、DNA配列、アミノ酸配列、表現型、などがあり、それぞれで何の家系かの意味合いの異なる。第2にデータをもとに「どんな距離」を定義するか、がある。これを距離行列法distance matrix methodと呼び、距離の定義の仕方によって、近隣結合法neighbor joining method（NJ法）や非加重結合法Unweighted Pair Group Method with Arithmetic mean（UPGMA）などがある。一方、系統の分岐に関する目的関数を定義し、これを最適化する手法がある。これは形質状態法character state methodと呼び、その中には最節約法maximum parsimony methodや最尤法maximum likelihood methodがある。

　本稿では、「DNA配列」を対象とした「最尤法」による系統樹の作成を紹介する。例えば、上図であれば、右に与えられた塩基配列からどうやって、系統樹を計算するか、を考える。いろいろな系統樹を描けること以上に、その背後にある方法論に注目し、その原理の理解を目指そう。様々な専門知識を要求するため、適宜、不明な用語は各自で調査してほしい。

突然変異と系統分岐

　生物のDNAの塩基配列はアデニン（A）、チミン（T）、グアニン（G）、シトシン（C）からなる。過去の生命体が有していたDNA配列は不変ではなく、時間経過とともにその配列に変化が生じ、現在の生命体に受け継がれている。例えば、配列中のある個所のGがAに変わることもある。このような塩基が置き換わることを塩基置換substitutionと呼ぶ。他には、TCCTと並んでいた配列中にTCACTみたいにAが挿入されることもある。あるいはTCTとなってCが抜け落ちることもある。もっと長い配列が挿入したり、欠失したりすることもある。これら塩基配列の挿入と欠失はまとめて挿入/欠失indel（インデル）と呼ぶ。

これらの変異はまとめて、突然変異mutationと呼ぶ。以降、突然変異のうち、塩基置換についてだけ議論する。インデルは挙動が複雑で、モデル化が難しいからである。

　ここで、系統が分岐する過程を想像しよう。いろんな過程があることが指摘されているが、今回はもっとも想像しやすいものを単純化して考える。例えば、下図のように左では、過去から現在にかけて任意交配可能なグループ＝集団が分かたれることなく現在に至っている。このとき、途中で突然変異が生じても、ある変異を持った個体は、別の変異を持った個体と交配可能で、様々な変異の組み合わせの遺伝子型が集団中に存在することになる。ゆえに遺伝子型は個体を識別することは可能でも、集団を分割できる指標にはならない。

　一方、右は地理的障壁によって右と左の集団の間で交配は不可能な状況である。すると、左集団で「TAGTGCT」という遺伝子型が突然変異で生じても、これは右集団では観察できない。逆に「TAATCCT」は左集団では観察できない。突然変異が長期にわたって各集団で独立に蓄積すると、特定の集団固有の変異＝遺伝子型、が生じるようになる。集団固有の変異を指標にすれば、集団を識別することができる。

　このような状態を、遺伝的分化genetic differentiationしているという。突然変異は時間と相関して蓄積し、さらにその多くは蓄積速度が一定であると考えられている（木村資生の分子進化の中立説neutral theory of molecular evolution）。ゆえに多くの突然変異に関しては、比べる集団間でどれだけ塩基配列が異なっているかを調べることで集団間の遺伝的分化を定量でき、さらに突然変異率を時計のように使って分岐年代を推定することも可能である。あたかも、塩基配列情報が突然変異という形で経過時間を記録している様から、これを分子進化時計molecular clockと呼んでいる。

塩基置換モデル

　過去から現在までどんな過程を経て塩基置換が蓄積してきたかをモデル化できれば、そのモデルをもとに生物の分岐の歴史を推定することができる。塩基置換に関するモデルは、そのまま塩基置換モデルと呼ばれ、いくつか提唱されている。代表的、かつシンプルな過程のモデルは下記に示すようなJC68モデルやK80モデルである。

　JC68モデルは、JunksとCantorが提唱したモデルで、例えばAはTGCのいずれにも同じ速度で変化する。TはT以外に、GはG以外に、CはC以外にすべて同じ速度で変化する。極めてシンプルなモデルであり、同時に現実の塩基置換ではこれにそぐわないことが多い。例えば、AとGはプリン塩基、TとCはピリミジン塩基と呼ばれ、AとG、TとCは共通の骨格を持っている。この時、共通の骨格を持っている塩基間は置換しやすいが、異なる骨格の塩基間では10倍程度置換しにくいことが知られている。共通骨格の塩基置換（A⇔G、T⇔C）をトランジションtransition、異なる骨格間の塩基置換（例えば、A⇔Tなど）をトランスヴァージョンtransversionと呼ぶ。

　K80モデルは、木村が提唱したモデルで、トランジションとトランスヴァージョンで異なる置換速度を考慮している。

　ほかにもK80モデルを発展させて、Hasegwa, Kishino, Yanoが提唱した、塩基組成も考慮に入れたHKY89も知られる。

・JC68における塩基置換関数の導出

　では一番シンプルなJC68に関して、時間tにおける塩基状態を表す関数を導出しよう（より一般の場合は追記）。下記のように、まずはAに注目し、時間tにおいてAである確率PA(t)を導出することを目指そう。それを使って、塩基iからjに変化する確率Pij(t)や変化しない確率Pii(t)を求める。αを塩基置換する確率とする。

するとある単位時間Δt後にAである確率は、漸化式を使って右のようにあらわせる。これを下記のように式変形していく。初期値t = 0の時に、注目する塩基がAであるか否かで場合分けが必要なのでそれに注意。

最終的に、上記の議論はA以外でも成り立つので同様の式がつかえて、以上から、全塩基置換パターンの関数を得られたことになる。

系統樹の尤度関数

　塩基置換モデルを導出できたことで、とうとう系統樹の尤度関数を定義できる。系統樹の尤度関数を定義するときは、まず、どのように種が分かたれたか＝樹形をこちらで設定してやる必要がある。この樹形のことをトポロジーtopologyと呼ぶ。

　例えば、下記のように真の樹形は不明だが、種3が時間t1の時点で分岐し、そのあと種1と種2が時間t2の時点でわかれた、というような樹形を仮定し、それの尤度関数を定義しよう。パラメータも、トポロジーもいろいろ設定できるが、ここでは最もシンプルなものを想定しておく。サイトaは塩基配列の左からa番目の塩基とする。サイトaにおける種1～3の塩基をia、ja、kaとする。塩基xとyは過去における塩基状態である。

観測値はia＝G、ja＝G、ka＝Aであり、このデータをもとにt1、t2および塩基置換速度λを最尤推定する。天下り的だが、尤度関数は下記のように定義される。Pij(t)は上記で計算した塩基置換モデルである。

いかつい式だが、恐れずに眺めてみるとサイトaにおける尤度関数は、仮定した樹形のときに観測値になる確率を愚直に計算しているに過ぎない。

ただし問題になるのが、過去の塩基の状態xやyを知らない、ということだ。そこで、ここでも愚直に、すべてのパターンを足すことにする。つまり、x = Aのとき、Tのとき、……、およびy = Aのとき、……、このすべてのパターンの確率を足し上げて尤度にするのである。これがサイトaにおける尤度関数の意味である。

　また、過去の塩基xの割合πxも同様にわからないのだが、もし、JCモデルのように本当にランダムに突然変異するなら、あるサイトにおける塩基の割合は過去の塩基の割合と（解析種数が大きければ）大体等しいだろう。ゆえに、現在の割合から過去の割合を推定したものをπxに代入して使う。いうなれば、πxが初期状態の確率を表していると考えればよいだろう。例えば、今回ならサイトaにおいて、Gは2個、Aは1個で、残りは0だからπA= 1/3、πG＝2/3、πT＝πC＝0である。

　サイトaにおける尤度関数に関して、この後の計算のために、展開して、行列のかたちに帰着させておく（見た目はいかついが、行列に帰着するほうが計算速度が速い）。

最後、全サイトにわたって尤度を掛け算すれば、塩基配列全体での尤度関数が完成である。塩基配列の長さがnなら、サイト1～nにわたって掛け算する。確率を表すから、通常は対数をとったもののほうがPC計算上の扱いは楽である。

系統樹の最尤推定

　さて、上記のように尤度関数を設定することができた。あとは、その関数を最大値を与えるt1、t2、λを探索をする。通常の回帰分析ならここで終わりだ。しかし、系統樹の作成においてはこれで終わらない。なぜなら、上記の尤度関数は種3が一番初めに分かれたことを仮定するモデルだからだ。このトポロジーが最適とは限らない。種1～3のうち、種3が一番初めに分かれたモデルを、

　　(種1, 種2), 種3

とあらわすとすると、他にも

　　(種1, 種3), 種2

　　(種2, 種1), 種1

を検証する必要がある。今回の場合、()内は入れ替えても同じトポロジーなので、何が1番始めに分かれた仮定するか、にトポロジーは依存している。ゆえに今回は3つのトポロジーについて、すべて最尤推定し、その中で最も尤度の大きいトポロジーを採用したうえで、最大尤度を与えるt1、t2、λを推定値とする必要がある。

・系統樹の最尤推定は時間がかかる

　3種の系統樹に関する尤度関数はそこそこ複雑だった。これがもし4種とか、5種になったらどうなるだろう。当然のごとく、尤度関数は複雑化する。そうなると、関数一つでの最尤推定には時間がかかる。さらに、トポロジーも考慮するべきものが爆発的に増加する。3種なら、トポロジーは3つだったが、4種になれば下図のように15種類まで膨れ上がる。ちなみに5種なら105種類だ。一般にn種いるとき、検討するべき系統樹の数は二重階乗を使って、(2n-3)!! = (2n-3)×(2n-5)×(2n-7)×……×5×3個である。

ゆえに、関数の複雑化と検討するべきトポロジーの個数の増加により、系統樹の推定は莫大な時間がかかる。

　当然ながら、nが多くなった時、そのすべてのトポロジーを検証すると、私たちが生きている間に計算が終わることはないだろう。そこで通常は、「だいたいこのトポロジーが正解に近いだろう」というものを用意して、そこから出発してちょっとずつトポロジーを変えて、尤度も最大になるように系統樹の探索を行う（局所探索と呼ばれる）。だいたい正解のトポロジーは、より計算の速いNJ法系統樹や最節約系統樹などが使われる。

Rで最尤推定してみる

　では、塩基配列が与えられた時、自分で最尤推定を試みてみよう。今回は3種、i、j、kの塩基配列を使って考えてみる。下記のように塩基配列をもとに系統樹を作成する。一つ一つのコマンドはここでは説明しない。ここでの目標は、系統樹を描くことよりも、系統樹の最尤推定を実感することである。

------------------------------------------------------

library(plotn)

library(viridis)

library(Biostrings)

library(msa)

library(phangorn)

library(ape)

i_DNA <- "AGGATAATCCATTTGATCTCCTGCTCCTGGTGTATTGTAC"

j_DNA <- "AGGATAATCCCTATGATATGATGCTCCTTGTGTATTGTAC"

k_DNA <- "ATGATAGTGCCTGTGATCTTATGATCCTGGTGCAGTGTAC"

sequences <- c(i_DNA, j_DNA, k_DNA)

biostrings_sequences <- DNAStringSet(sequences)

msa_result_sample <- msa(biostrings_sequences, method = "ClustalW")

print(biostrings_sequences)

## DNAStringSet object of length 3:

##     width seq

## [1]    40 AGGATAATCCATTTGATCTCCTGCTCCTGGTGTATTGTAC

## [2]    40 AGGATAATCCCTATGATATGATGCTCCTTGTGTATTGTAC

## [3]    40 ATGATAGTGCCTGTGATCTTATGATCCTGGTGCAGTGTAC

phyDat_msa_sample <- as.phyDat(msa_result_sample)

names(phyDat_msa_sample) <- c("Sp.i", "Sp.j", "Sp.k")

pdist <- dist.dna(as.DNAbin(phyDat_msa_sample), model = "JC69")

pdist

##           Sp.i      Sp.j

## Sp.j 0.1673577          

## Sp.k 0.3040988 0.3040988

nj_tree <- nj(pdist)

nj_tree$edge.length[which(nj_tree$edge.length < 0)] <- 0

nj_tree <- midpoint(multi2di(nj_tree))

plot(nj_tree)

tree_UPGMA <- upgma(pdist)

tree_UPGMA <- midpoint(tree_UPGMA)

plot(tree_UPGMA)#図1の描画

------------------------------------------------------

すると、種kが先に分かれて、そのあとに種iとjが分かれたことが推定された。では、下記のようにπxや塩基置換モデルを定義して、尤度関数を定義しよう。logL関数内の計算は、途中で紹介した尤度関数を行列に分解した時の方法で計算している。logL関数は、sequence引数には、塩基配列のベクトル、outgroup引数には、最も初めに分岐したと考える種の塩基配列が格納されている番号をとる。例えば、sequence引数に、種i、j、kの順で塩基配列が格納されていたとき、outgroup = 1と指定すれば、種iが最も初めに分かれたと仮定していることになる。param引数には、listオブジェクトをとり、list[[1]] = t1、list[[2]] =t2、list[[3]] =λをまとめたものをとる。

------------------------------------------------------

g <- function(sequences, site, x){#πxの計算

  tmp <- strsplit(sequences,"")

  b_site <- NULL

  for(i in 1:length(tmp)){

    b_site <- c(b_site, tmp[[i]][site])

  }

  sum(b_site == x)/length(b_site)

}

p <- function(x, y, t, lambda){#塩基置換モデル

  mapply(function(xt, yt){

      if(xt == yt){

    1/4 + 3/4 * exp(-4/3 * lambda * t)

  } else{

    1/4 - 1/4 * exp(-4/3 * lambda * t)

  }

  }, xt = x, yt = y)

}

logL <- function(sequences, outgroup, params){

  outseq <- strsplit(sequences[outgroup],"")[[1]]

  inseq <- strsplit(sequences[-outgroup],"")

  mapply(function(t_1, t_2, lambda_1){

    res <- sapply(1:length(outseq), function(kt){

      tmp1 <- c(g(sequences, kt, "A") * p("A", outseq[kt], t_1+t_2, lambda_1), 

                g(sequences, kt, "T") * p("T", outseq[kt], t_1+t_2, lambda_1), 

                g(sequences, kt, "G") * p("G", outseq[kt], t_1+t_2, lambda_1), 

                g(sequences, kt, "C") * p("C", outseq[kt], t_1+t_2, lambda_1))

      PM <- outer(c("A","T","G","C"), c("A","T","G","C"), FUN = p, t = t_1, lambda = lambda_1)

      Pv <- c(p("A", inseq[[1]][kt], t_2, lambda_1) * p("A", inseq[[2]][kt], t_2, lambda_1),

              p("T", inseq[[1]][kt], t_2, lambda_1) * p("T", inseq[[2]][kt], t_2, lambda_1),

              p("G", inseq[[1]][kt], t_2, lambda_1) * p("G", inseq[[2]][kt], t_2, lambda_1),

              p("C", inseq[[1]][kt], t_2, lambda_1) * p("C", inseq[[2]][kt], t_2, lambda_1))

      return(log(tmp1 %*% as.vector(PM %*% Pv)))

    })

    return(sum(res))

  }, t_1 = params[[1]], t_2 = params[[2]], lambda_1 = params[[3]])

}

------------------------------------------------------

このようにlogL関数を定義し、下記のようにoptim関数をcontrol = list(fnscale = -1)として実行することで、logL関数を最大化するt1、t2、λを得ることができる。outgroup引数を1、2、3としたとき、それぞれで最尤推定を行う。

------------------------------------------------------

optim(list(0.1, 0.1, 0.1), logL, sequences = sequences, outgroup = 1, 

      method = "L-BFGS-B", lower = 0, control = list(fnscale = -1))

## $par

## [1] 0.0000000 0.3731721 0.3456129

## 

## $value

## [1] -62.44325

## 

## $counts

## function gradient 

##        9        9 

## 

## $convergence

## [1] 0

## 

## $message

## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

optim(c(0.1, 0.1, 0.1), logL, sequences = sequences, outgroup = 2, 

      method = "L-BFGS-B", lower = 0, control = list(fnscale = -1))

## $par

## [1] 0.0000000 0.3731721 0.3456129

## 

## $value

## [1] -62.44325

## 

## $counts

## function gradient 

##        9        9 

## 

## $convergence

## [1] 0

## 

## $message

## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

optim(c(0.1, 0.1, 0.1), logL, sequences = sequences, outgroup = 3, 

      method = "L-BFGS-B", lower = 0, control = list(fnscale = -1))

## $par

## [1] 0.1123253 0.2970536 0.3437593

## 

## $value

## [1] -62.15676

## 

## $counts

## function gradient 

##       10       10 

## 

## $convergence

## [1] 0

## 

## $message

## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

------------------------------------------------------

得られた結果のうち、valueがlogLの最大値で、parがその時のt1、t2、λである。すると、outgroup = 3のとき、1や2の時と比べてlogLは最も大きな値を示し、その時、t1 = 0.112、t2 = 0.30、λ = 0.343であることがわかる。つまり、系統樹としては種kが0.412時点で種ijグループと分岐し、種iとjが0.112時点で分岐するモデルが最適であることがわかる。

　では最適なλの時、t1とt2が様々な値をとった時の尤度関数を図示してみよう。

------------------------------------------------------

res <- optim(c(0.1, 0.1, 0.1), logL, sequences = sequences, outgroup = 3, 

      method = "L-BFGS-B", lower = 0, control = list(fnscale = -1))

lambda <- res$par[3]

n <- 40

t1 <- seq(0.001,1,length=n)

t2 <- seq(0.001,1,length=n)

t <- data.frame(t1 = rep(t1, each = n), 

                t2 = rep(t2, times = n))

t$logL <- logL(sequences, 3, params = list(t$t1, t$t2, lambda))

lL <- matrix(t$logL, byrow = T, nrow = n)

lev <- seq(round(min(lL), digits = 1), round(max(lL), digits = 1), length = 10)

image(t1, t2, lL,

      col = magma(400), xlab = "t1", ylab = "t2", las = 1)#図2の図示

contour(t1, t2, lL, levels = lev, add = T)

abline(v = res$par[1])

abline(h = res$par[2])

------------------------------------------------------

このようにして系統樹の最尤推定は行われる。私は昔、なぜATGCの文字情報が、数値として扱えるのかが理解できていなかったが、裏ではこのように関数の設定が行われていたのだ。もちろん、ここでの例は極めて単純化した話であるが、この単純な例の発展であることには変わりない。系統樹を最尤推定する、というトピックに関して理解の助けになれたら幸いである。

補足1：より一般的な塩基置換モデルをどう解くか

　塩基置換モデルに関して、上記ではJC68モデルを例にとってきたが、それ以外の塩基置換モデルはどう解くべきだろうか。実は、行列の計算さえできれば、そう難しい話ではない。もう一度、JC68モデルを使って説明しよう。今度は、Aだけに注目するのではなく、4つの塩基を同時に考える。

一部の矢印しか表示していないが、各塩基の漸化式は上記のようになることがわかるだろう。これを以下のように計算を進めていく。

最終的にベクトルの微分のかたちに落ち着く。この時、右辺に登場する行列が推移行列transition matrixと呼ばれる（個体群生態学でも登場する）。1～4行目はそれぞれATGCに対応し、1～4列目も同様にATGCに対応する。そして、列の塩基から行の塩基にどれくらいの速度で変化するかを記述したもの、である。例えば、1行2列目の項はTからAにαの割合で加わることを表すし、3行3列目の項はGからGに-3αの割合で加わる＝3αの割合で減っていくことを示す。推移行列の各塩基置換モデルによって推移行列は異なっており、上記のかたちがJC68モデルの推移行列のかたちとなる。

　ここで線形代数の知識が必要だが、推移行列をRと表せば、Rの固有値、固有ベクトルを用いて、下記のようにあらわすことが可能である。

V^(-1)は行列で、P(t)は縦ベクトルだからV^(-1)P(t)は縦ベクトルで、これをひとまとめとしてみれば、上記の解は、積分定数を使って下記のとおり。

ゆえにさらに計算を続けると

推移行列Rは、各塩基から各塩基に、どのような割合で変化するかがわかれば記述でき、各モデルに対応するRがわかれば固有値、固有ベクトルも計算が可能である。そして、固有値と固有ベクトルをまとめた行列Vから、比較的簡単に各塩基から各塩基に変化する確率を計算が可能である。

　例えば、JC68モデルなら、Rは上記の通りで、またRの固有値とそれに対応する正規直交固有ベクトルをまとめた行列Vは以下の通り（固有ベクトルや正規直交ベクトルの求め方は、さらに追記する）。

K80モデルなら、以下の通りになる。トランジションの変化率はα、トランスヴァージョンの変化率はβで表されている。同じ塩基のまま変わらないときは、2塩基分のトランスヴァージョンと1円基分のトランジションの差し引きがあることも反映されている。例えば、Aから変化するとき、TとCはトランスヴァージョンで、Gはトランジションの変化であり、これは残りの塩基のいずれを開始としても関係は変わらない。

補足2：推移行列の固有値と正規直交固有ベクトル

　では、JC68やK80モデルについて、固有値や正規直交固有ベクトルを求めてみよう。以下のように計算する。JC68について、

K80について、