Go to Qiita Advent Calendar 2024 Top
LoginSignup
7
7

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@MaedaTaro_Umiushi

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

Last updated at Posted at 2018-06-29

#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

7
7
0

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
  3. You can use dark theme
What you can do with signing up
Sign upLogin
7
7

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?