Rで系統樹を扱う(ape, ggtree)

Rで系統樹を扱うときに必要なapeやggtreeの使い方がどうしても覚えられないので軽くまとめました。
内容については随時追加予定です。

参考にしたサイト

• ggtree: Elegant Graphics for Phylogenetic Tree Visualization and Annotation
• facet_plotの記事
• Data Integration, Manipulation and Visualization of Phylogenetic Trees 4, 7章あたり (そのまま4章に飛びます)
• Visualizing and Annotating Phylogenetic Trees with R+ggtree

系統樹データについて

treeオブジェクトの中身

treeの中には、大体以下の5つのデータが格納されています。
edge：どのノードとノードが繋がっているか（matrix）
edge.length：エッジの長さ (vector)
Nnode：ノード数(int)
node.label：ノードのラベル(vector)
tip.label：葉のラベル(vector)

読み込み・書き出し

apeの関数を使います。

sample.R

###ライブラリのダウンロード###
library(tidyverse)
library(ape)
library(ggtree)

#読み込み
tree <- read.tree('path/to/tree')
#書き出し
write.tree(phy=tree, file='path/to/out/tree')

系統樹を作る

適当な木を作りたい時は、rtree()を使います。

##乱数を指定
set.seed(12345678)
##引数は葉の数
tree <- rtree(10)

系統樹データを図示する

表示はggtreeを使います。
ggtreeを使うと系統樹だけではなく、ヒートマップなどのアノテーションを付けられるのもいいところです（後述）。
操作の基本はggplotに沿っていますが、独特な表記もありややこしいなと思います。。

sample.R

#とりあえずこれで系統樹は出てくる
ggtree(tree)

###ノード番号を表示させたい場合###
ggtree(tree) + geom_text(aes(label=node), hjust=-.2)

###葉のラベルを表示させたい場合###
ggtree(tree) + geom_tiplab() 

形状も様々な形に指定することができます。

sample.R

###いろいろな形状###
ggtree(tree, layout="slanted") 
ggtree(tree, layout="circular")
ggtree(tree, layout="fan", open.angle=120)
ggtree(tree, layout="equal_angle")
ggtree(tree, layout="daylight")
ggtree(tree, branch.length='none')
ggtree(tree, branch.length='none', layout='circular')
ggtree(tree, layout="daylight", branch.length = 'none')

その他のデザイン変更は以下のコマンドでできます。

sample.R

###系統樹全体のデザイン変更###
##color: 色
##size: 線の太さ
##linetype: 線の種類
ggtree(tree, color="firebrick", size=1, linetype="dotted") #1


###ノード全体にしるし###
ggtree(tree) + geom_nodepoint(color="pink", alpha=0.7, size=5) #2

###葉全体にしるし###
ggtree(tree) + geom_tippoint(color="darkgreen", shape=8, size=5) #3

###指定ノードにしるし###
##subsetを使う
ggtree(tree) + geom_point2(aes(subset=(node==15)), shape=21, size=5, fill='green') #4

##isTip(TRUE/FALSE)を使えばノード/葉の区別ができる
ggtree(tree) + geom_point(aes(shape=isTip, color=isTip), size=3) #5

###色分け###
##groupClade()を使ってクレードを分ける
tree2 <- groupClade(tree, c(15,17))
g6 <- ggtree(tree2, aes(color=group)) + 
  scale_color_manual(values=c("black", "firebrick", "steelblue")) #6

shapeの番号指定はggplot2と同じです。他の形を使いたい時はggplot2 point shapesを確認してください。
色の名前もここでチェックできます。

クレードの色分けはgroupClade()を使わない方法もあります。以下に２種類の方法を載せます。

sample.R

###横にラベリング
ggtree(tree) + 
  geom_tiplab(offset =.05) + 
  geom_cladelabel(node=15, label="Clade A", 
                  color="red2", offset=.3, align=TRUE) + 
  geom_cladelabel(node=13, label="Clade B", 
                  color="blue", offset=.3, align=TRUE) +
  theme_tree2() + 
  xlim(0,5) +
  theme_tree()

最後の3行では、右が切れてしまうことへの対処法が記されています(試しに３行を入れずに動かしてみると良いと思います)。
詳しい解説はVisualizing and Annotating Phylogenetic Trees with R+ggtreeにあります。

sample.R

###クレードを塗りつぶす###
ggtree(tree) + 
  geom_tiplab() + 
  geom_hilight(node=15, fill="gold", alpha = 0.5, extend = .1) + 
  geom_hilight(node=13, fill="purple", alpha = 0.5, extend = .2)

アノテーションをつける

treeのtip.labelを含むアノテーションデータを用いて、系統樹の横にtipデータのグラフを付け加えることができます。
ここが一番ややこしいです

sample.R

##アノテーションデータの生成
set.seed(1234)
d <- data.frame(label=tree$tip.label, var1=abs(rnorm(10)), var2=as.factor(rbinom(10,3,0.4)))
##アノテーションデータを組み込む
tree <- full_join(tree, d, by='label')
ggtree(tree) +  
  geom_tippoint(aes(colour=var1),size=5) + 
  scale_colour_gradient(low='blue', high='red') 

ggtree(tree) +  
  geom_tippoint(aes(shape=var2),size=5) 

これだと表現できる内容・数が限られてしまうので、横にどんどんアノテーションしていくためにfacet_plot()を使います。
facet_plot()ではtreeデータに情報を追加しなくても、アノテーションデータのtip.lable情報をもとにうまく当てはめてくれます。

sample.R

library(ggstance) #geom_barh用
tree3 <- rtree(15)
d <- tibble(label=tree3$tip.label, var1=runif(15, min=0, max=1) , var2=as.factor(rbinom(15,3,0.4)))
p <- ggtree(tree3) + geom_tiplab(size=3, align=T, linesize=0.3) + xlim_tree(0.01)
p <- facet_plot(p, panel="rate", data=d, geom=geom_barh, 
                mapping=aes(x=var1, fill=var2),stat = "identity")
p

ヒートマップを付け加える場合にはgheatmap()を、配列のアライメントを図示したい場合にはmsaplot()を使うことができます。
力尽きたので詳しい説明を見てください。上のfacet_plotについても説明があります。

系統樹の加工

　apeの関数を用いて、tipを指定し、そこだけを残したり、消したりできます。

sample.R

###tipを残す###
tree_filtered <- keep.tip(tree, tip.label)
###tipを消す###
tree_filtered <- drop.tip(tree, tip.label)

#tip(ラベル)はあらかじめ指定
#vector/numericどちらも可能

###指定したnode以下のtreeを取り出す###
tree_selected <- extract.clade(tree, node)

祖先や子孫を調べる

様々なライブラリを用いて系統関係を調べることができます。

sample.R

###node(番号)の祖先を調べる###
library(phangorn)
Ancestors(tree, node.num, type = c("all", "parent"))
Descendants(tree, node.num, type = c("tips", "children", "all"))

##node(番号)の祖先を調べる
#nodeはvectorでもできる
#all(すべての祖先)かparent(1つ上の親)を選択可能
#返り値はリスト

###nodes(tips)の共通祖先ノードを返す###
library(ape)
getMRCA(tree, nodes)