[連載] フリーソフトによるデータ解析・マイニング 第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、radial.phylog など) がある。パッケージ ade4 は CRAN ミラーサイトからダウンロードできる。
  関数 newick2phylog、hclust2phylog、taxo2phylog はそれぞれクラス newick、hclust、taxo を phylog クラスに変換する関数である (表3を参照)。関数 plot、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(tith.phy,clabel.n= 0.6, f = 0.75) #図7を作成
> radial.phylog(tith.phy) #図8を作成

関数plot.phylogの例

図7 関数 plot の例

関数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()

[1] "genbank" "embl" "emblwgs" "swissprot"
[5] "ensembl" "hogenom" "hogenomdna" "hovergendna"
[9] "hovergen" "hogenom5" "hogenom5dna" "hogenom4"
[13] "hogenom4dna" "homolens" "homolensdna" "hobacnucl"
[17] "hobacprot" "phever2" "phever2dna" "refseq"
[21] "greviews" "bacterial" "archaeal" "protozoan"
[25] "ensprotists" "ensfungi" "ensmetazoa" "ensplants"
[29] "ensemblbacteria" "mito" "polymorphix" "emglib"
[33] "refseqViruses" "ribodb" "taxodb"

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

> 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://www.ddbj.nig.ac.jp/wp-content/downloads/mag/No30.pdf ) に載せられているゲノムデータを GenBank や EMBL などでダウンロードできる。このサイトに関しては次のサイトに日本語による説明がある (http://wing.jst.go.jp/detail.php?wid=129&sw=#keitai ) 。 インターネットデータベースを利用する例として、図12に示されているデータ Archaea を GenBank からダウンロードするコマンド(データベースへのアクセス・操作、名前の確認、データのダウンロードなど)を次に示す。

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

図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 など)。