ggtreeで外部情報のついた系統樹を書く

系統樹を書く必要が出てきたので練習します！

1. まずは系統樹そのものを書いてみる

Netwik形式の系統樹でテストしたいので、まずはtree.dndというタイトルのテキストファイルに以下の文字列を入れて保存（練習なので枝が5本だけの超シンプルな系統樹にしました）

(((genome1:0.3,genome2:0.2):0.4,genome3:0.6):0.6,((genome4:0.3,genome5:0.2):0.3):0.2);

今回は円形 (circular) の系統樹を想定しているので、で系統樹を書けるかチェックします

> library(tidyverse)
> library(ggtree)
> tree <- read.tree("./tree.dnd")
> ggtree(tree, layout="circular")

geom_tiplab() を付けると各枝にdndファイルに入っている名前がつきます

> ggtree(tree, layout="circular") + geom_tiplab()

> ggtree(tree, layout="circular") + geom_tiplab() + geom_treescale()

必須ではありませんが、後半の処理のためにノード番号を可視化しておきます

> ggtree(tree, layout="circular") + geom_tiplab() + geom_treescale() + geom_text(aes(label=node), hjust = -1.5)

2. 系統樹の枝に色をつける

ノード指定で枝に色を付けてみようと思います

groupClade() で10番ノード以下をグルーピングして系統樹を書きます

> tree <- groupClade(tree, .node=c(10))
> ggtree(tree, layout="circular", aes(color=group)) + geom_tiplab() + geom_treescale() + geom_text(aes(label=node), hjust = -1.5)

色を指定することもできます。青と黒で指定してみます

> ggtree(tree, layout="circular", aes(color=group)) + geom_tiplab() + geom_treescale() + geom_text(aes(label=node), hjust = -1.5) + scale_color_manual(values=c("black", "blue"))

3. 系統樹をグラデーションにする

ゲノムの存在比によって色付けしたいので、新しく存在比データを読み込みます

まず練習用に2列のcsvファイルを作成 (rate.csv)

label,rate
genome1,0.28
genome2,0.50
genome3,0.20
genome4,0.01
genome5,0.01

ではこのデータを系統樹にのせていきます！

以下リンクで紹介されているツールfsatAnc() を用いる方法で描いてみます (リンク：http://www.phytools.org/eqg2015/asr.html)

> library(phytools)
> rate <- read.csv("./rate.csv", row.names = 1)
> rate <- as.matrix(rate)[,1]
> fit<-fastAnc(tree,rate,vars=TRUE,CI=TRUE)
> obj<-contMap(tree,rate, plot=TRUE,type="fan")

4. 系統樹にヒートマップをアノテーションする

遺伝子発現/有無情報を系統樹にアノテーションします

↓これもcsvファイルで作成

label,geneA,geneB,geneC
genome1,200,0,500
genome2,200,0,30
genome3,200,0,500
genome4,0,400,300
genome5,0,300,800

まずは系統樹、次に枝名がインデックスになるようにcsvを読み込みます

> p <- ggtree(tree, layout="circular", aes(color=group)) + geom_tiplab() + geom_treescale() + scale_color_manual(values=c("black", "blue"))
> h <- read.csv("./genes.csv", row.names = 1)
> gheatmap(p, h)