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

Rと系統樹(2)

1. 系統樹のデザインと操作

(1) 先端のデザイン

 パッケージapeには系統樹の先端(tip)のデザインのための関数tiplabelsがある。関数tiplabelsの書き式とデフォルトの設定を次に示し、用いる主な引数の説明を表1に示す。

tiplabels(text, tip, adj = c(0.5, 0.5), frame = "rect",pch = NULL, pie = NULL, thermo = NULL, piecol = NULL, col = "black", bg = "yellow", ...)

表1  関数tiplabelsに用いる主な引数
text 表示する文字列、あるいはシンボルのベクトル。省略することが可能。
tip 表示するノードの番号のベクトル。省略することが可能。
adj 表示する文字列、あるいはシンボルの座標(x,y)。xは左右、yは上下。
frame 表示する文字列を四角("rect"デフォルト値)、あるいは円("circle")で囲む
pch 表示するシンボルのベクトル
pie 円グラフを表示する数値ベクトル。
thermo 温度計型のグラフを表示する数値ベクトル
piecol 引数pie、thermoに用いる色
col 用いるテキストあるいはシンボルの色

 関数tiplabelsの使用例として、データwoodmouesを用いて、系統樹の先端に異なる種類のシンボルを示すコマンドとその例を次に示す。

> library(ape);data(woodmouse);
> wood.dist<-dist.dna(woodmouse)
> wood.tr<-nj(wood.dist)
> lab<-c(rep(10,3),rep(11,2),rep(12,3), rep(13,7)) #印の番号を作成する
> plot(wood.tr, "c", FALSE, font = 1, label.offset = 2,x.lim = 20, no.margin = TRUE)
> tiplabels(pch =lab,col =lab, adj = 1.5, cex = 2)

関数tiplabelsの例

図1 関数tiplabelsの例

 引数pie、thermoに各個体の特徴となる数値ベクトルを用いると円グラフ、温度計グラフを付け加えることができる。

(2) ノードのデザイン

 パッケージapeには系統樹のノードのデザインのための関数nodelabelsがある。関数nodelabelsの書き式とデフォルトの設定を次に示す。用いる引数の項目は関数tiplabelsとほぼ同じである。

nodelabels(text, node, adj = c(0.5, 0.5), frame = "rect",pch = NULL, thermo = NULL, pie = NULL, piecol = NULL,col = "black", bg = "lightblue", ...)

 まず系統樹のノードに番号を付ける例を次に示す。

> plot(wood.tr,no.margin = TRUE)
> nodelabels(wood.tr$node.label,adj=0, cex=1.3)

関数nodelabelsの例1

図2 関数nodelabelsの例1

 ノードには複数の値を表示することもできる。その例として3セットの乱数(rs,rs2,rs3)をノードに表示するコマンドを次に示す。文字を表示する位置は引数adjで指定する。

> plot(wood.tr, use.edge.length = FALSE)
> rs1 <- round(runif(13, 60, 100), 0)
> rs2 <- round(runif(13, 60, 100), 0)
> rs3 <- round(runif(13, 60, 100), 0)
> nodelabels(rs1, adj = -0.2, cex = 0.8, font = 2)
> nodelabels(rs2, adj = c(1.2, 1), frame = "n", cex = 0.8, font = 3)
> nodelabels(rs3, adj = c(1.2, -0.2), frame = "n", cex = 0.8)

関数nodelabelsの例2

図3 関数nodelabelsの例2

 系統樹のノードに何らかの値を円グラフや温度計グラフなどで示すこともできる。その例を次に示す。

> plot(wood.tr, "c",use.edge.length = FALSE)
> nodelabels(thermo = runif(13), cex = .8)

関数nodelabelsの例3

図4 関数nodelabelsの例3

(3) 系統樹の分割とズーム

 系統樹の、ある部分を異なる色で示すには関数plotの引数edge.colorを用いる。その例を次に示す。

> col<-rep("black",nrow(wood.tr$edge))
> col[1:5]<-"lightblue"
> plot(wood.tr,edge.color=col, cex=1.3, edge.width=2)

系統樹の分割表示の例

図5 系統樹の分割表示の例

 系統樹が大きいと、詳細の部分が読み取れない。パッケージapeには大きい系統樹の一部分を切り取って表示する関数zoomがある。ここではパッケージapeに用意されている系統樹のデータchiropteraを用いることにする。データchiropteraは幾つかのソースを複合したコウモリの系統樹である。各個体のラベルはchiroptera$tip.labelで呼び出すことができる。
 関数zoomを用いて巨大な系統樹の一部分を切り取って表示するためには個体の番号、あるいは個体のラベルを指定することが必要である。データchiropteraの中のラベルに文字列"Plecotus"と"Dobsonia"が含まれている部分を切り取って拡大し、表示する例を示す。

> data(chiroptera)
> zoom.p<- list(grep("Plecotus", chiroptera$tip.label), grep("Dobsonia", chiroptera$tip.label)
> zoom(chiroptera, zoom.p, cex=1.5)

関数zoomの使用例

図6 関数zoomの使用例

2. パッケージade4

 パッケージade4には系統樹に関連する幾つかの関数(newick2phylog, hclust2phylog, taxo2phylog, plot.phylog, radial.phylogなど)がある。パッケージade4はCRANミラーサイトからダウンロードできる。
 関数newick2phylog, hclust2phylog, taxo2phylogはそれぞれクラスnewick、hclust、taxoをphylogクラスに変換する関数である(表3を参照)。関数plot.phylog, radial.phylog, dotchart.phylogは異なる系統樹のグラフを作成する関数である。
 パッケージade4の中に用意されているデータセットtithoniaを用いた例を示す。まず次のようなコマンドでデータ属性および構造を確認する。

> library(ade4);data(tithonia)
> class(tithonia);length(tithonia)

[1] "list"
[1] 2

> tithonia[1];tithonia[2]
<結果は省略>

 返されたデータから分かるように、tithoniaは2つの形式で記録されたデータリストである。データtithonia[2]を用いた系統樹のグラフを作成する例を次に示す。

> tith2<- as.data.frame(tithonia[2])
> tith.dis<-dist(tith2)
> tith.hc<-hclust(tith.dis,"ave")
> tith.phy<-hclust2phylog(tith.hc))
> plot.phylog(tith.phy,clabel.n= 0.6, f = 0.75) #図7を作成
> radial.phylog(tith.phy) #図8を作成

関数plot.phylogの例

図7 関数plot.phylogの例

関数radial.phylogの例

図8 関数radial.phylogの例

3.データの形式と利用

(1) 結果のオブジェクトのクラス

 Rには遺伝子データ処理に関連するパッケージが幾つかある。パッケージによって返す結果のオブジェクトのクラスが異なる場合がある。表2に幾つかのパッケージが返す結果のオブジェクトのクラスと主な項目の名称を示す。

表2 オブジェクトのクラス
パッケージ クラス 主な項目
ape phylo $edg, $edg.length, $tip.label, $node.label$rood.edg
ape matching $matching,$edg.length$edge.tip.label,$node.label
apTreeshape treeshape $merge,$names
ade4 phylog $tre, $leaves, $nodes, $parts, $paths, $droot, $call, $Wmat, $Wdist, $Wvalues, $Wscores, $Amat, $Avalues, $Adim, $Ascores, $Aparam, $Bindica, $Bscores, $Bvalues, $Blabels

 また、表2に示すクラスを互いに変換する幾つかの関数を表3に示す。

表3 クラスの変換関数
変換前 変換後 関数
phylo phylog newick2phylog()
matching as.matching()
treeshape as.treeshape()
phylog phylo as.phylo()
matching phylo as.phylo()
treeshape phylo as.phylo()
hclust phylo as.phylo()
phylog hclust2phylog()

(2) 系統樹データ

 系統樹をテキスト形式で記録したデータとしてNewick フォーマットがある。Newickは 系統樹を次のように、階層的に丸括弧で区切ったテキスト形式のデータである。枝の長さはコロンの右に数値で示す。

((Human:0.1,Gorilla:0.1):0.4,(Mouse:0.2,Rat:0.2):0.3);

 パッケージapeでは関数read.tree, パッケージade4では関数newick2phylogで読み込み・変換を行う。その使用例を次に示し、その結果を用いて作成したグラフを図9に示す。

> a<-"((Human:0.1,Gorilla:0.1):0.4,(Mouse: 0.2,Rat:0.2):0.3);"
> par(mfrow=c(2,1))
> plot(read.tree(text=a))
> plot(newick2phylog(a))

 第2節で用いたデータtithoniaのリスト1 (tithonia[1])はNewickフォーマットである。

関数plot.phylogの例

図9 関数plot.phylogの例

(3)  遺伝子配列データ

 遺伝子配列データの記録にも幾つかのフォーマットがある。パッケージapeの中の関数read.dnaは3種類のフォーマット("interleaved", "sequential", "fasta")を読み込むことができる。
 パッケージseqinrの中には関数read.fasta, read.alignmentがある。関数read.alignmentは"mase", "clustal", "msf", "phylip", "fasta"フォーマットのアレンジメントのデータを読み込むことができる。

(4) インターネットデータベースの利用

 パッケージape、apTreeshape, seqinrにはインターネットデータベースを利用する関数が用意されている。パッケージseqinrはCRANミラーサイトからダウンロードすることができる。パッケージapeにはデータベースGenBank(http://www.ncbi.nlm.nih.gov/Genbank/)を利用する関数read.GenBankがある。
 GenBankにはフウキンチョウ(Tanager)の遺伝子データがある。その一部分をダウンロードして用いる例を示す。
 まず次のように遺伝子データのIDのリストを作成する。

> tan.id<-c("AF006213", "AF006219", "AF006233","U15717", "U15718", "U15719", "U15720","U15721")

 これらのデータの詳細はGenBankのホームページから検索できる。関数read.GenBankとIDを用いてデータをダウンロードし、その系統樹を作成する例を次に示す。

> tan8<-read.GenBank(tan.id)
> ta.tree<-nj(dist.dna(tan8))
> plot(ta.tree,use.edge.length =F)

フウキンチョウの系統樹の例

図10 フウキンチョウの系統樹の例

 パッケージapTreeshapeの中の関数panditはインターネットデータベースPANDIT(http://www.ebi.ac.uk/goldman-srv/pandit/)を利用することができる。
 パッケージseqinrはより多くのデータベースを利用する環境を提供する。使用可能なデータベースは関数choosebank()で確認し、アクセスすることができる。

> library(seqinr)
> choosebank()

choosebank実行結果

 これらのデータベースで利用可能なデータの量に関しては、次のコマンドでグラフに表示することができる。

> banks <- c(choosebank())
> ntaxa <- numeric(0)
> for (i in banks) {
  ntaxa[i]<-as.numeric(choosebank(i)$totspecs);closebank()
 }
> dotchart(log10(ntaxa[order(ntaxa)]), pch = 19, main = "Number of taxa in available databases",xlab = "Log10(number of taxa)")

データベースにおける利用可能なデータの量

図11 データベースにおける利用可能なデータの量

 図12に示す国立遺伝学研究所のサイト(http://gib.genes.nig.ac.jp/)に載せられているゲノムデータをGenBankやEMBLなどでダウンロードできる。このサイトに関しては次のサイトに日本語による説明がある(http://wing.jst.go.jp/detail.php?wid=129&sw= #keitai)。 インターネットデータベースを利用する例として、図12に示されているデータArchaeaをGeaBankからダウンロードするコマンド(データベースへのアクセス・操作、名前の確認、データのダウンロードなど)を次に示す。

国立遺伝学研究所のデータ公開サイト

図12 国立遺伝学研究所のデータ公開サイト

> choosebank("genbank")
> query("AA", "sp=Archaea")
> sapply(AA$req[1:4], getName)

[1] "A75589" "A80402" "A80403" "A80404"

> sq<-lapply(bb$req[1:4], getSequence)
> sq[1]

[[1]]
[1] "a" "t" "c" "t" "g" "g" "t" "g" "t" "t"
<中略>
[3101] "a" "a" "a" "t" "c" "g" "t" "a" "t" "a"
[3111] "t" "t"

 本稿で用いたパッケージ以外にもゲノムデータに関する幾つかのパッケージがある(GeneR, genArise, GeneSpringなど)。