使用ggtree實現進化樹的可視化和注釋

本文作者：餘光創，目前就讀於香港大學公共衛生系，開發過多個R/Bioconductor包，包括ChIPseeker, clusterProfiler, DOSE,ggtree，GOSemSim和ReactomePA。

進化樹看起來和層次聚類很像。有必要解釋一下兩者的一些區別。

層次聚類的側重點在於分類，把距離近的聚在一起。而進化樹的構建雖然也可以說是一個聚類過程，但側重點在於推測進化關係和進化距離(evolutionary distance)。

層次聚類的輸入是距離，比如euclidean或manhattan距離。把距離近的聚在一起。而進化樹推斷是從生物序列（DNA或胺基酸）的比對開始。最簡單的方法是計算一下序列中不匹配的數目，稱之為hamming distance（通常用序列長度做歸一化），使用距離當然也可以應用層次聚類的方法。進化樹的構建最簡單的方法是非加權配對平均法（Unweighted Pair Group Method with Arithmetic Mean, UPGMA），這其實是使用average linkage的層次聚類。這種方法在進化樹推斷上現在基本沒人用。更為常用的是鄰接法（neighbor joining），兩個節點距離其它節點都比較遠，而這兩個節點又比較近，它們就是neighbor，可以看出neighbor不一定是距離最近的兩個節點。真正做進化的人，基本不用這些基於距離的方法。現在主流的方法是最大似然法(Maximum likelihood, ML)，通過進化模型（evolutionary model)估計拓樸結構和分支長度，估計的結果具有最高的概率能夠產生觀測數據（多序列比對）。另外還有最大簡約法和貝葉斯推斷等方法用於構建進化樹。

Newick是最常用的存儲進化樹的文件格式，如上面這個樹，拓樸結構用newick格式可以表示為：

(B,(A,C,E),D);

括號最外層是根節點，它有三個子節點，B, (A,C,E)和D，而節點(A,C,E)也有三個子節點A，C和E。

加上分支長度，使用 : 來分隔：

(B:6.0,(A:5.0,C:3.0,E:4.0):5.0,D:11.0);

比如A:5.0代表的是A與其父節點的距離是5.0。

內部節點也可以有label，寫在相應的括號外面，如下所示：

(B:6.0,(A:5.0,C:3.0,E:4.0)Ancestor1:5.0,D:11.0);

這是最為廣泛支持的文件格式，很多進化樹可視軟體只支持newick格式。

ggtree的開發源自於我需要在樹上做注釋，發現並沒有軟體可以很容易地實現，通常情況下我們把統計信息加到節點的label上來展示，比如CodeML的dN/dS分析，輸出文件裡就給用戶準備了newick樹文本，把dN/dS (ω) 加於節點label之上:

codeml_file <- system.file("extdata/PAML_Codeml/mlc", package="ggtree")tree_text <- readLines(codeml_file)[375:376]tree_text

## [1] "w ratios as labels for TreeView:" ## [2] "(K #0.0224 , N #0.0095 , (D #0.0385 , (L #0.0001 , (J #0.0457 , (G #0.1621 , ((C #0.0461 , (E #0.0641 , O #0.0538 ) #0.0001 ) #0.0395 , (H #0.1028 , (I #0.0001 , (B #0.0001 , (A #0.0646 , (F #0.2980 , M #0.0738 ) #0.0453 ) #0.0863 ) #1.5591 ) #0.0001 ) #0.0001 ) #0.0549 ) #0.0419 ) #0.0001 ) #0.0964 ) #0.0129 );"

這種做法只能展示一元信息，而且修改節點label真心是個髒活，滿滿的都是不爽，我心中理想的方式是樹與注釋信息分開，注釋信息可以方便地通過圖層加上去，而且可以自由組合。於是著手開發ggtree。ggtree是個簡單易用的R包，一行代碼ggtree(read.tree(file))即可實現樹的可視化。而注釋通過圖層來實現，多個圖層可以完成複雜的注釋，這得力於ggtree的設計。其中最重要的一點是如何來解析進化樹。

除了ggtree之外，我所了解到的其它畫樹軟體在畫樹的時候都把樹當成是線條的集合。很明顯畫出來的進化樹就是一堆線條，但是線條表示的是父節點和子節點的關係，除此之外沒有任何意義，而節點在進化樹上代表物種，葉子節點是我們構建進化樹的物種，內部節點是根據葉子節點推斷的共同祖先。我們所有的進化分析、推斷、實驗都是針對節點，節點才是進化樹上有意義的實體。這是ggtree設計的基礎，ggtree只映射節點到坐標系統中，而線條在geom_tree圖層中計算並畫出來。這是與其它軟體最根本的不同，也是ggtree能夠簡單地用圖層加注釋信息的基礎。

有很多可視化包基於ggplot2實現，包括各種gg打頭的，號稱擴展了ggplot2，支持圖形語法(grammar of graphics)，我並不認同。雖然基於ggplot2產生的圖，我們可以用theme來進一步調整細節，用scale_系列函數來調整顏色和標尺的映射，但這些不足以稱之為』支持圖形語法』，圖形語法最關鍵核心的部分我認為是圖層和映射。

像ggphylo, OutbreakTools和phyloseq這幾個包都有基於ggplot2的畫樹函數，但其實都不支持圖形語法，它們所實現的是複雜的函數，畫完就完事了，用戶並不能使用圖層來添加相關的信息。

如果在OutbreakTools這個包中：

if (show.tip.label) { p <- p + geom_text(data = df.tip, aes(x = x, y = y, label = label), hjust = 0, size = tip.label.size)}

如果show.tip.label=FALSE，當函數返回p時df.tip就被扔掉，用戶想要再加tip.label就不可能了。ggphylo和phyloseq都是類似的實現，這些包把樹解析為線條，所以節點相關的信息需要額外的data.frame來存儲，並且只有極少數的預設參數，比如上面例子中的tip.label。在上面的例子中，用戶連更改tip.label的顏色都不可能，更別說使用額外的注釋信息了。

這幾個包所實現的畫圖函數，都可以很容易地用ggtree實現，並且經過測試，ggtree運行速度比這幾個包都要快。更多信息請參考ggtree的wiki頁面。

ggtree是真正擴展ggplot2，支持圖形語法的包。我們首先擴展ggplot支持tree object做為輸入，並實現geom_tree圖層來畫線條。

library(ggplot2)
library(ggtree)
set.seed(2015-11-26)
tree <- rtree(30)
ggplot(tree, aes(x, y)) + geom_tree()

ggtree函數是ggplot() + geom_tree() + xlab(NA) + ylab(NA) + theme_tree()的簡單組合。

ggtree(tree)

想要加tip.label，用geom_tiplab圖層，並且ggplot2的圖層都可以直接應用於ggtree。

ggtree(tree) + geom_tiplab() + geom_point(color='firebrick')

ggtree提供了多個函數可以把clade放大縮小(scaleClade)，摺疊(collapse)和展開(expand)，位置調換(flip)和旋轉(rotate)，以及分類(groupOTU, groupClade)。

nwk <- system.file("extdata", "sample.nwk", package="ggtree")
tree <- read.tree(nwk)
p <- ggtree(tree)
cp <- ggtree(tree) %>% collapse(node=21) + ggtitle('collapse')
ep <- cp %>% expand(node=21) + ggtitle('expand')
hp <- p %>% hilight(node=21) + ggtitle('hilight')
rp <- hp %>% rotate(node=21) + ggtitle('rotate')
library(gridExtra)
grid.arrange(cp, ep, hp, rp, ncol=2)

ggtree支持的文件格式包括Newick, Nexus, NHX和jplace。

上面已經展示了Newick格式，下面的例子是NHX格式：

nhxfile = system.file("extdata", "ADH.nhx", package="ggtree")
nhx <- read.nhx(nhxfile)
ggtree(nhx, ladderize=F) + geom_tiplab() + geom_point(aes(color=S), size=8, alpha=.3) + theme(legend.position="right") + geom_text(aes(label=branch.length, x=branch), vjust=-.5) + xlim(NA, 0.3)

我們知道FigTree是針對BEAST的輸出設計的，可以把BEAST的統計推斷拿來給樹做注釋，但很多的進化分析軟體並沒有相應的畫樹軟體支持，用戶很難把信息展示出來。
ggtree支持ape, phangorn, r8s, RAxML, PAML, HYPHY, EPA, pplacer和BEAST的輸出。相應的統計分析結果可以應用於樹的注釋。可以說ggtree把這些軟體分析的結果帶給了R用戶，通過ggtree的解析， 這些進化分析結果可以進一步在R裡進行處理和統計分析，並不單單是在ggtree中展示而已。

raxml_file <- system.file("extdata/RAxML", "RAxML_bipartitionsBranchLabels.H3", package="ggtree")
raxml <- read.raxml(raxml_file)ggtree(raxml) + geom_text(aes(label=bootstrap, color=bootstrap)) + scale_color_gradient(high='red', low='darkgreen') + theme(legend.position='right')

multiPhylo也是支持的，所以100顆bootstrap樹可以同時用一行代碼展示出來。

btree_file <- system.file("extdata/RAxML", "RAxML_bootstrap.H3", package="ggtree")
btree = read.tree(btree_file)
ggtree(btree) + facet_wrap(~.id, ncol=10)

如果不分面，這100顆樹會重疊畫在一起，這也能很好地展示bootstrap分析的結果，bootstrap值低的clade，線條會比較亂，而bootstrap值高的地方，線條一致性比較好。

使用BaseML預測的祖先序列，ggtree解析結果的同時會把父節點到子節點的subsitution給統計出來，可以直接在樹上注釋：

rstfile <- system.file("extdata/PAML_Baseml", "rst", package="ggtree")
rst <- read.paml_rst(rstfile)
p <- ggtree(rst) + geom_text(aes(label=marginal_AA_subs, x=branch), vjust=-.5)
print(p)

不同於BaseML以鹼基為單位，CodeML預測祖先序列，以密碼子為單位。ggtree定義了一個操作符%<%，如果有相同的注釋信息要展示，可以用tree object來更新tree view。

rstfile <- system.file("extdata/PAML_Codeml", "rst", package="ggtree")
crst <- read.paml_rst(rstfile)
p %<% crst

像上面的例子，用crst來更新p，就是用crst畫出來的樹+注釋。對比兩圖，可以發現BaseML和CodeML推測的祖先序列是稍有不同的。

CodeML的dN/dS分析，我們可以直接把數據拿來給樹上色。同樣道理分類數據也可以拿來上色。

mlc_file <- system.file("examples/mlc", package="ggtree")
mlc <- read.codeml_mlc(mlc_file)
ggtree(mlc, aes(color=dN_vs_dS)) + scale_color_continuous(limits=c(0, 1.5), high='red', low='green', oob=scales::squish, name='dN/dS') + theme(legend.position='right')

進化樹已經被廣泛應用於各種跨學科的研究中，隨著實驗技術的發展，各種數據也更易於獲得，使用用戶數據注釋進化樹，也是ggtree所支持的。

nwk <- system.file("extdata", "sample.nwk", package="ggtree")
tree <- read.tree(nwk)
p <- ggtree(tree)
dd <- data.frame(taxa = LETTERS[1:13], place = c(rep("GZ", 5), rep("HK", 3), rep("CZ", 4), NA), value = round(abs(rnorm(13, mean=70, sd=10)), digits=1))
dd <- dd[sample(1:13, 13), ]
row.names(dd) <- NULL
print(dd)

## taxa place value## 1 A GZ 68.8## 2 J CZ 56.8## 3 L CZ 74.7## 4 C GZ 53.3## 5 F HK 62.8## 6 B GZ 60.8## 7 E GZ 87.1## 8 M <NA> 70.9## 9 H HK 67.0## 10 G HK 59.8## 11 I CZ 77.7## 12 K CZ 69.8## 13 D GZ 66.3

在上面的例子中，使用一個分類數據和一個連續型數據，輸入的唯一要求是第一列是taxon label。ggtree中定義了操作符%<+%，來添加數據。添加之後，用戶的數據對ggplot是可見的。可以用於樹的注釋。

p <- p %<+% dd + geom_text(aes(color=place, label=label), hjust=-0.5) + geom_tippoint(aes(size=value, shape=place, color=place), alpha=0.25)p+theme(legend.position="right")

ggtree還支持用戶把自己的數據和樹保存為jplace格式。

更多的實例請參考vignette。

ggtree允許把不同軟體的分析結果整合在一起，同時在樹上展示或者比較結果。在我們提交的論文中，使用了整合BEAST和CodeML的例子，在樹上展示dN/dS、時間軸、胺基酸替換、clade support values、物種和基因型 (genotype）等多維信息，6種不同的信息同時展示在一顆進化樹上，這是個複雜的例子，我們在附件1中展示了可重複的代碼。如果有興趣，可以留意一下我們的文章。 :)

我們把樹當成節點的集合，而不是線條的集合，這一點回歸到了進化樹的本質意義上，使這一實現成為可能。而支持圖形語法，與ggplot2的無縫銜接又讓注釋變得更加容易。ggtree為我們打開了各種注釋和操作的可能性。甚至於可以創造出好玩的圖，比如使用showtext來加載圖形化的字體、用emoji來畫樹、使用圖片來注釋樹等等。

一個比較正經又好玩的是使用PhyloPic資料庫上的圖形。

pp <- ggtree(tree) %>% phylopic("79ad5f09-cf21-4c89-8e7d-0c82a00ce728", color="steelblue", alpha = .3)
pp + geom_tiplab(align=T, linetype='dashed', linesize=.5) + geom_tippoint(color='firebrick', size=2)

另一個好玩又為我們展現各種可能性的是 subview 函數，它使得圖上加小圖變得特別容易。並且已經被應用於地圖上加餅圖。解決這個問題的初衷在於，想要給節點加餅圖注釋。有了subview函數之後，這會變得很容易，當然我還沒有寫出給節點加餅圖的函數，因為我還沒有這個需求，得有一些實際的數據做參考，這樣才能夠設計出更易用的函數呈現給用戶。

很多的功能還在開發之中，有問題/建議請及時在Github上報告(中英文都可以)。

本文微信編輯：馮璟爍

統計之都：專業、人本、正直的中國統計學門戶網站。

關注方式：掃描下圖二維碼。或查找公眾帳號，搜索 統計之都 或 CapStat 即可。

往期推送：進入統計之都會話窗口，點擊右上角小人圖標，查看歷史消息即可。