search
LoginSignup
7

More than 3 years have passed since last update.

posted at

updated at

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:1つの枝を交換する(最頻用?)

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 をつかう

最初の2項目が"ここから""ここまで"線を引けという感じ(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

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
What you can do with signing up
7