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

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

## 参考にしたサイト
* [ggtree: Elegant Graphics for Phylogenetic Tree Visualization and Annotation](https://guangchuangyu.github.io/ggtree-book/chapter-ggtree.html#introduction-1)
* [facet_plotの記事](https://guangchuangyu.github.io/2016/10/facet_plot-a-general-solution-to-associate-data-with-phylogenetic-tree/)
* [Data Integration, Manipulation and Visualization of Phylogenetic Trees](http://yulab-smu.top/treedata-book/chapter4.html) 4, 7章あたり (そのまま4章に飛びます)
* [Visualizing and Annotating Phylogenetic Trees with R+ggtree](https://bioconnector.github.io/workshops/r-ggtree.html)

## 系統樹データについて
#### 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() 

```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/605546/4035ad1a-e113-8ab9-9ed0-acde4cc75395.png)

形状も様々な形に指定することができます。
```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')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/605546/47140250-9e34-5156-e034-47fa0b315582.png)

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

```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
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/605546/ec62ccb7-b2bb-0640-19ed-bc2207f67147.png)

shapeの番号指定はggplot2と同じです。他の形を使いたい時は[ggplot2 point shapes](http://www.sthda.com/english/wiki/ggplot2-point-shapes)を確認してください。
色の名前も[ここ](http://sape.inf.usi.ch/sites/default/files/ggplot2-colour-names.png)でチェックできます。

クレードの色分けは`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()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/605546/a25c6258-c4f7-6767-67c2-b7aef97f5c22.png)

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



```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)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/605546/fef2d50c-e5e8-5acb-a848-4bf3ae8e0a93.png)



#### アノテーションをつける
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) 
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/605546/58fa9540-ddaa-699e-0a4e-11dcebffbf8e.png)

これだと表現できる内容・数が限られてしまうので、横にどんどんアノテーションしていくために`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
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/605546/94a9e2ab-3ced-e8be-f323-4ce407d1f150.png)

ヒートマップを付け加える場合には`gheatmap()`を、配列のアライメントを図示したい場合には`msaplot()`を使うことができます。
力尽きたので[詳しい説明](https://yulab-smu.top/treedata-book/chapter7.html)を見てください。上の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)

```