ggtreeを使ってRで系統樹を扱う

#ggtreeとは
The University of Hong KongのYu & Lamによって作られたR上で系統樹を扱うパッケージ

G Yu, DK Smith, H Zhu, Y Guan, TTY Lam*. ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods in Ecology and Evolution. 2017, 8(1):28-36. doi:10.1111/2041-210X.12628.

以下の内容は
https://bioconductor.org/packages/release/bioc/vignettes/ggtree/inst/doc/treeAnnotation.html
をざっくり翻訳し適宜コメントを追加したもの

#ライブラリ
以下パッケージが入っていることを前提とする
library(ggtree)
library(ape)
library(Biostrings)
library(ggplot2)
library(magrittr)
library(colorspace)

#Newickファイルの読み込み

iq-treeが吐くnewickファイルは普通にapeライブラリのread.treeで読み込める

```
tree <- read.tree("iq_tree_run.nex.contree")
```

情報リッチなデータのときは以下など参考
https://bioconductor.org/packages/devel/bioc/vignettes/ggtree/inst/doc/treeVisualization.html#time-scaled-tree


#とりあえずプロット

```
ggtree(tree) 
```


ggtree()だけだとOTUラベルが表示されない。
ラベルを描画するには geom_tiplab()

```
ggtree(tree) +  geom_tiplab()
```

ラベルを文字列リストとして得るには
tree$tip.label


#樹形の変更

##ノード番号の調べ方
Nodeには番号が付いている、以下のコマンドで表示できる。この番号が樹をいじる際の重要情報

```
ggtree(tree) +geom_text2(aes(subset=!isTip, label=node), hjust=-.3)+  geom_tiplab()
```


MRCAを使えば、OTU間のノードIDを直接調べられる

```
MRCA(tree, tip=c('Halimeda_borneensis', 'Halimeda_discoidea'))
```


##クレードのグルーピング
viewCladeを使えば特定のクレードだけ描画することができる

```
p0 <- ggtree(tree) 
viewClade(p0 +geom_tiplab(), node=26)
```

groupCladeを使えば特定のクレードだけグルーピングすることができる

```
tree <- groupClade(tree, node=26)
```

colorやlinetypeに"group"と指定することで色や線形を変えられる
ggtree(tree, aes(color=group, linetype=group))

まとめて使うとこういう感じ

```
tree <- groupClade(tree, node=28)
ggtree(tree, aes(color=group, linetype=group))
```

複数のノードをグルーピングする場合はc()でID番号を指定する。

```
tree <- groupClade(tree, node=c(26, 28))
ggtree(tree, aes(color=group, linetype=group)) 
```

geom_tiplabと併用すれば、特定グループだけにOTUラベルを表示することもできる

```
ggtree(tree, aes(color=group, linetype=group)) + geom_tiplab(aes(subset=(group==1)))
```

groupOTUを使えばOTU名からグルーピングをすることもできる。

```
tree <- groupOTU(tree, focus=c("Bryopsis_hypnoides", "kpDNA_bry", "Lambia_antarctica"))
ggtree(tree, aes(color=group, linetype=group)) + geom_tiplab(aes(subset=(group==1)))
```

複数のグループに名前をつけて、色分けしたりもできる

＜グループを定義＞

```
cls <- list(Derbesiaceae=c("Derbesia_sp"),
            Bryopsidaceae=c("Bryopsis_plumosa", "Bryopsis_hypnoides", "kpDNA_bry", "Lambia_antarctica"),
            Codiaceae=c("Codium_decorticatum"),
            Ostreobiaceae=c("Ostreobium_quekettii"),
            Udoteaceae=c("Avrainvillea_mazei","Tydemania_expeditionis", "kpDNA_Poro", "kpDNA_Rhip"),
            Caulerpaceae=c("Caulerpa_lentillifera", "Caulerpa_racemosa", "Caulerpa_cliftonii"),
            Halimedaceae=c("Halimeda_discoidea", "Halimeda_borneensis"))
            
```

##色分け

```
tree <- groupOTU(tree, cls)
library("colorspace")
ggtree(tree, aes(color=group)) + geom_tiplab() +
  scale_color_manual(values=c("black", rainbow_hcl(7))) + theme(legend.position="right")
```

#################################################################################
#クレードの交換
##outgroupの指定

```
tre1 <- root(tree, outgroup = "Chlamydomonas_reinhardtii")
tre1 <- groupOTU(tre1, cls)

ggtree(tre1, aes(color=group)) + geom_tiplab() +
  scale_color_manual(values=c("black", rainbow_hcl(7))) + theme(legend.position="right")
```


##rotate：ノードを中心に枝を回転させる(2OTUをつなげているところしか働かない：バグ？）

```
p0 <- ggtree(tre1) +geom_text2(aes(subset=!isTip, label=node), hjust=-.3)+  geom_tiplab()　#node番号を確認
tree <- groupClade(tre1, c(24, 25))
p <- ggtree(tree, aes(color=group)) + scale_color_manual(values=c("black", "firebrick", "steelblue"))
tre2 <- rotate(tree, 24) %>% rotate(25)
p1 <- ggtree(tre2, aes(color=group)) + scale_color_manual(values=c("black", "firebrick", "steelblue"))+  geom_tiplab()
p2 <-ggtree(tree, aes(color=group)) + scale_color_manual(values=c("black", "firebrick", "steelblue"))+  geom_tiplab()
multiplot(p1, p2, ncol=2)
```

##rotate：１つの枝を交換する（最頻用？）

```
tree <- groupClade(tre1, c(22, 25))
p <- ggtree(tree, aes(color=group)) + scale_color_manual(values=c("black", "firebrick", "steelblue"))
multiplot(p, flip(p, 22, 25), ncol=2)
```

##collapseでノードをまとめる

```
ggtree(tre1) +geom_text2(aes(subset=!isTip, label=node), hjust=-.3)+  geom_tiplab()
tre2 <- ggtree(tre1) %>% collapse(node=32) #なぜかうまく動かない
```

##scaleCladeで枝の幅を変える

```
ggtree(tre1) %>% scaleClade(27, scale=0.5) + geom_hilight(27, "steelblue")+  geom_tiplab()
```
#################################################################################
#アノテーションをつける
##geom_cladelabel:OTUラベルの右側にバーでクレード名をつけられる

```
tree <- groupOTU(tree, cls)
library("colorspace")
p <- ggtree(tree, aes(color=group)) + geom_tiplab() +
  scale_color_manual(values=c("black", rainbow_hcl(7))) + theme(legend.position="right") + geom_text2(aes(subset=!isTip, label=node))

p+geom_cladelabel(node=21, label="Bryopsidineae") +  geom_cladelabel(node=19, label="Halimedineae")

```


###バーの位置がずれるときはalignで揃えるとよい
p+geom_cladelabel(node=21, label="Bryopsidineae", align=TRUE, offset=.5) +  geom_cladelabel(node=19, label="Halimedineae", align=TRUE, offset=.5)

###バーの色はcolorで変える
p+geom_cladelabel(node=21, label="Bryopsidineae", align=TRUE, offset=.5, color='darkgreen')+
  geom_cladelabel(node=19, label="Halimedineae", align=TRUE, offset=.5, color='green')

###文字の角度もangleで変えられる
p+geom_cladelabel(node=21, label="Bryopsidineae", align=TRUE, angle=270, hjust='center', offset.text=.02, offset=.5, color='darkgreen')+
  geom_cladelabel(node=19, label="Halimedineae", align=TRUE, angle=270, hjust='center', offset.text=.02, offset=.5, color='green')


###geomで文字を飾ることもできる
p+geom_cladelabel(node=21, label="Bryopsidineae", align=TRUE, angle=270, hjust='center', offset.text=.02, offset=.5, color='darkgreen', geom='label')+
  geom_cladelabel(node=19, label="Halimedineae", align=TRUE, angle=270, hjust='center', offset.text=.02, offset=.5, color='green')

##まとめたいOTUが単系統で無いときは、geom_strip をつかう
## 最初の２項目が"ここから""ここまで"線を引けという感じ(OTU名でもnode numberでも指定できるがOTUのときlabelを使うとおかしくなるバグがある）

p + geom_strip("kpDNA_bry", "Ostreobium_quekettii",
barsize=2, color='red', barextend=.5)



##geom_hilight でクレードを囲んで色付できる alphaで透明度指定
p+ geom_hilight(node=21, fill="steelblue", alpha=.3) +
  geom_hilight(node=19, fill="darkgreen", alpha=.3)

##geom_balanceでもクレードを囲んで色付できる。入れ子構造の描画に便利
p+ geom_balance(node=21, fill="steelblue", alpha=.3) +
  geom_balance(node=19, fill="darkgreen", alpha=.3)

##########################################################################################3
#外部データを取り込んで、描画する

##数量データ（ゲノムサイズなど）をドットの大きさで表示することが出来る

cpDNA_size <- read.table("./gene/cpDNA_size.txt",header=T)
row.names(cpDNA_size) <- NULL
p <- p %<+% cpDNA_size


p + geom_tiplab() +
  geom_tippoint(aes(size=Size), alpha=0.5,align=TRUE, offset=.5)+ 
  geom_text(aes(label=Size), hjust=1, vjust=1.4, size=3)+
  xlim(NA, 1)

p+ xlim(NA, 1)
p+theme(legend.position="right")



##gheatmap：質的データ（遺伝子の有無等）をboxカラーで示す(テーブルデータに注意）
genotype <- read.table("./gene/sum2_PS1.txt", sep="\t",header=T)

colnames(genotype) <- sub("\\.$", "", colnames(genotype))

p_f <- gheatmap(p, genotype, offset = .5, width=2, font.size=3, colnames_angle=-45, hjust=0)+
  scale_fill_brewer(genotype,  palette = "Set1", direction = 1)

  
  
   scale_fill_manual(breaks=c("Cp_My",
                             "Cp_Pre",
                             "NO_may_Nuc_Pre",
                             "Nuc_My",
                             "Nuc_Pre",
                             "Nuc_Pre_transcript",
                             "Slug",
                             "unkonwn"
                             ), values=c("palegreen", "darkgreen", "white", "orange", "darkorange", "magenta", "blue", "gray"))

p_f