[連載]フリーソフトによるデータ解析・マイニング第42回

Rと系統樹(1)

1. 系統樹

 同一祖先を持っていることを前提としている系統内の個体の子孫関係の推定結果を樹木の枝分岐の形式で示すグラフを系統樹と言う。その例としては、生物種や遺伝子などの進化関係を樹木状に表現したものがあげられる。系統樹の枝分かれは系統の分岐を示し、枝の長さは進化の程度や時間経過を表す。子孫は枝分かれする各ノードを最も近い共通祖先とする。系統樹の嚆矢はドイツの生物学者エルンスト・ヘッケル (Ernst Haeckel)であると言われている。
 系統樹は生物学で主に研究・応用されていると思われがちであるが、人文科学の分野でも研究・応用が進められている[1][2][3]。それは、多くの文化現象(言語、写本、民俗など)、あるいは文化遺物が、同一の祖先から時間の流れとともに進化・変化したと考えられるからである[4][5]。近年、進化の視点で系統樹の手法を用いた工業製品の分析に関する研究が見られるようになっている[6]。このようなことから、系統樹はもはや生物学研究のみの専用ツールではなくなった。
 今日の系統樹を作成する基本技法は、1960年代前後から統計的分類学と集団遺伝学分野からそれぞれ独立に始まった研究成果に基づいている。系統樹は図1に示す有根系統樹と無根系統樹に分けることができる。図1(a)のような樹根(root)がある系統樹を有根系統樹、図1(b)のような樹根がない系統樹を無根系統樹と呼ぶ。

系統樹の例

図1 系統樹の例

2. 系統樹の推定(距離法)

 研究対象の本当の系統関係が未知であり、その系統関係を明らかにするために系統樹というツールが用いられている。あるデータを基に何らかの方法で作成した系統樹は、本当にそうであるかどうかは誰も断言できない。あくまでも、ある理論に基づいたもっともらしい推定結果に過ぎない。
 系統樹の推定には、幾つかの方法が提案されているが、本稿では、R環境上での距離法を中心として簡潔に説明する。
 距離法は、個体の間の何らかの距離を求め、その距離行列を基に系統樹を作成する方法である。基本的な考え方およびプロセスは、階層的クラスタリンク法と同じである。

(1) 距離

 量的なデータにおける距離に関しては、多次元尺度法や階層的クラスター分析の回(本欄の第27回〜29回(2005年10~12号))で説明を行った。本稿の系統樹の説明には、遺伝子データを用いるので、まず遺伝子配列データ間の距離について触れておく。遺伝子データは、アルファベットの文字列で表記されている。
 パッケージapeの中には、woodmouseというデータがある。データwoodmouse はMichaux らによるアカネズミのミトコンドリアのデータ配列の一部分であり、15個体により構成されたリストデータである。それぞれは1個体のデータの配列であり、その長さは965文字列になっている。そのデータ概略は次のコマンドで読み取ることができる。

> library(ape)
> data(woodmouse)
> str(woodmouse)
strコマンドの結果

 データの中に現れた文字の相対頻度は関数base.freqを用いて求めることができる。

> base.freq(woodmouse)

a c g t
0.3065414 0.2613083 0.1260264 0.3061239

 返された結果から分かるように、データwoodmouseは文字"a","c","g","t"で表記されたDNA配列である。
 データ形式をデータフレーム形式に変換すると、次のように個体ごとの文字の相対頻度を求め、図示することができる。

> wm.f<-t(as.data.frame(woodmouse))
> barplot(apply(wm.f,1,base.freq),las=2)

個体ごとのゲノムの相対頻度

図2 個体ごとのゲノムの相対頻度

 長さが同じである文字列間の距離の測度としては、表記上の一致度を測るハミング距離がよく知られている。しかし、DNAの変化は置換、挿入、欠失、逆位があり、また塩基置換の比率なども考慮しなくてはならないので、ハミング距離では不十分である。DNAの距離計測に関しては、幾つかの方法が提案されている。
 パッケージapeには、関数dist.dnaがある。その書き式を次に示す。

dist.dna(x, model = "K80", variance = FALSE,gamma = FALSE, pairwise.deletion = FALSE, base.freq = NULL, as.matrix = FALSE)

 引数xは、DNA配列のマトリックス、データフレーム、リスト形式のデータオブジェクトである。引数modelは、距離を求める方法であり、表1に示す方法がオプションとして用意されている。DNAの距離を求める際には、塩基置換の比率に関して何らかの仮定を置かなければならない。オプションに用意された方法の塩基サイトの変化率に関する仮説とマトリックスを表1に示す。表の中ではDNAの表記に用いる文字は大文字A,G,C,Tを用いる。

表1 関数dist.dnaの引数model
JC69' JukesとCantorが1969年に提案したモデル。どの塩基サイト(A,G,C,T)でも同じの変化率αで互いに変化すると仮定している。
関数dist.dnaの引数modelJC69の場合
K80 Kimura が1980年に提案したモデル。サイト当たりの変化率を次のように仮定している。
関数dist.dnaの引数modelK80の場合
F81 Felsensteinが1981年にJC69を一般化して、次の変化率マトリックスを仮定したモデル。ただし、πA、πG、πC、πTは塩基頻度である。 関数dist.dnaの引数modelF81の場合
K81 Kimuraが1981年にK80を一般化したモデル。変化率A⇔C,G⇔TとA⇔T,C⇔Gが異なると仮定している。
関数dist.dnaの引数modelK81の場合
F84 Felsensteinが1984年に、K80を一般化したモデル。塩基頻度が不平衡状態にあると仮定している。
関数dist.dnaの引数modelF84の場合
T92 Tamuraが1992年にK80モデルを一般化したモデル。塩基頻度(GとCの頻度、あるいは含量)θ=πG+πCを用いる。
関数dist.dnaの引数modelT92の場合
TN93 Tamura と Neiが1993年に提案したモデル。A⇔C,C⇔Tが異なる比率αRγで置換すると仮定している。
関数dist.dnaの引数modelTN93の場合
αRγであるとF84、αRγRγであるとHasegawa, Kisino and Yanoが1985年に提案したモデルに等しい。
GG95 Galtier and Gouyが1995年に提案したモデル。G+Cの含量が時間とともに変化することを許す不平衡モデル。その変化率マトリックスはT92に似ている。

 表1の諸方法の距離の計算などに関しては、参考文献[7]が詳しい。
 関数dist.dnaの引数modelのデフォルトは "K80"になっている。データwoodmouseの一部分を用いて、デフォルト値と"T92"で求めた距離を次に示す。

>dist.dna(woodmouse[1:5])

No305 No304 No306 No0906S
No304 0.0160
No306 0.0138 0.0042
No0906S 0.0192 0.0138 0.0095
No0908S 0.0171 0.0117 0.0074 0.0127

>dist.dna(woodmouse[1:5],model="GG95")

No305 No304 No306 No0906S
No304 0.0000
No306 0.0000 0.0000
No0906S 0.0132 0.0131 0.0131
No0908S 0.0000 0.0000 0.0000 0.0131

 また、パッケージapeには遺伝子データの距離を求める関数dist.geneもある。

(2) UPGMA系統樹

 距離法の中で、最も簡単な方法は非加重平均結合(UPGMA: unweighted pair-group method using arithmetic average)法である。この方法は、階層的クラスタリンクの群平均法に等しいので、次のように階層的クラスタリンク法で系統樹を作成することができる。

> wm.d<-dist.dna(woodmouse)
> wm.hc<-hclust(wm.d,"average")
> wm.phy<-as.phylo(wm.hc)
> plot(wm.phy)

 用いた関数as.phyloは、オブジェクトをphyloクラスに変換する。関数plotは関数 plot.phyloの略である。

K80距離のUPGMA系統樹

図3 K80距離のUPGMA系統樹

(3) 関数plot.phylo

 関数plot.phyloは、異なるタイプの系統樹の作成やデザインを行う関数である。その書き式次に示し、用いる引数の説明を表2に示す。表2の中の"De"はデフォルトの設定を指す。

plot(x, type = "phylogram", …)

表2 関数plot.phyloの主な引数
引数 機能 オプション
x phyloクラスのオブジェクト
type 樹の種類 "p","u"(De),"c", "r"
use.edge.length 樹の枝の長さの表現 TRUE, FALSE(De)
node.pos ノードの接続位置 1,2,De=NULL
edge.color 枝の色 De=" black"
edge.width 枝の太さ De=1
font ラベルの文字フォント 1(plain text),2(bold),3(italic)(De), 4(bold italic)
cex 文字サイズ De=1
adj ラベルの文字列の位置 0は左寄せ(De),0.5はセンターリング,1は右寄せ
srt ラベルの回転角度 負数が時計回りの方向(De=0)
label.offset ラベルと樹の線との間隔 De=0
direction 樹の向き "r"(右、De), "l" (左), "u" (上), "d" (下)
lab4ut 無根樹のラベルの書き方 "horizontal"(De),"radial"

 また、no.margin、root.edge、underscoreなどの引数がある。次にデータwoodmouseの"K80"距離を用いた系統樹の作成例を示す。

> plot(wm.phy,type="c")  #図(a)
> plot(wm.phy,type="c",use.edge.length = FALSE)  #図(b)
> plot(wm.phy,type="r")  #図(c)
> plot(wm.phy,type="u",use.edge.length = FALSE,lab4ut="radial")  #図(d)


異なるタイプの系統樹

図4 異なるタイプの系統樹

(4) 近隣結合法(NJ: Neighbor-Joining)

 近隣結合法は、Saito and Nei が1987年に考案した系統樹の作成方法である。近隣結合法の特徴の1つは、計算量が少ないため、大量の個体の系統樹の作成に向いている。個体の数が少ない(4~5)ときは、最小進化法と等価である。
 パッケージapeには、近隣結合法の関数njがある。関数njに用いる引数は距離マトリックスのみである。関数njの有根系統樹と無根系統樹の作成の例として、データwoodmouseのK80距離を用いた結果を次に示す。

> wm.d<-dist.dna(woodmouse)
> wm.nj<-nj(wm.d)
> plot(wm.nj)

NJ法の有根系統樹

図5 NJ法の有根系統樹

> plot(wm.nj,type="u",use.edge.length = FALSE,lab4ut="radial")

NJ法の無根系統樹

図6 NJ法の無根系統樹

 また、パッケージapeには、距離マトリックスを基に最小スパニング樹 (spanning tree)を作成する関数mstがある。

参考文献
[1] Michael J. O'Brien and R. Lee Lyman (2003). Cladistics and Archaeology, The University of Utah Press.
[2] Ruth Mace, Clare J. Holden, and Stephen Shennan(2005). The Evolution of Cultural Diversity:
  A Phylogenetic Approach,UCL Press
[3] Carl P. Lipo, Michael J. O'Brien, Mark Collard, and Stephen J. Shennan (2005) Mapping Our Ancestors:
  Phylogenetic Approaches in Anthropology and Prehistory, Transaction Publishers
[4] 村上征勝編(2006). 文化情報学入門、勉誠出版(p.36-48、矢野環著、文化系統学−歴史を復元する−).
[5] 三中信宏(2006). 系統樹思考の世界:すべてはツリーとともに,講談社
[6] 石山智明,伊藤則人,柴田裕介,土松隆志*,池上高志(2005). 系統樹から迫る非生命進化:鳥居・雑煮・デジタルカメラ,
   第7 回日本進化学会国際シンポジウム
[7] 根井正利/S.クマー共著、根井正利監訳(2006). 分子の進化と系統学、培風館