原文鏈接： 2 Manipulating Tree with Data
2.1使用Tidy接口操作樹文件
所有被treeio解析/合并的樹數(shù)據(jù) treeio (Wang et al., 2020)都可以使用 tidytree包轉(zhuǎn)換成tidy格式數(shù)據(jù)。tidytree包提供了tidy接口來操作帶有關(guān)聯(lián)數(shù)據(jù)的樹。例如恃锉，可以將外部數(shù)據(jù)與系統(tǒng)發(fā)育聯(lián)系起來盯荤，也可以使用tidyverse verbs合并從不同來源獲得的進(jìn)化數(shù)據(jù)橄杨。在處理完樹數(shù)據(jù)后秦踪，可以將其轉(zhuǎn)換回一個樹數(shù)據(jù)對象并導(dǎo)出到單個樹文件中懂鸵，在R中進(jìn)一步分析或使用ggtree (Yu et al., 2017)和 ggtreeExtra (Xu, Dai, et al., 2021)進(jìn)行可視化。
2.1.1 phylo對象
phylo是R語言中系統(tǒng)發(fā)育分析包包 ape 中的基礎(chǔ)命令 (Paradis et al., 2004)痪宰。該領(lǐng)域的大部分R包廣泛依賴于phylo對象叼架。tidytree 包提供了 as_tibble 方法,將 phylo 對象轉(zhuǎn)換為 tidy 數(shù)據(jù)框架,一個tbl_tree對象畔裕。
library(ape)
set.seed(2017)
tree <- rtree(4)
tree
##
## Phylogenetic tree with 4 tips and 3 internal nodes.
##
## Tip labels:
## t4, t1, t3, t2
##
## Rooted; includes branch lengths.
x <- as_tibble(tree)
x
## # A tibble: 7 × 4
## parent node branch.length label
## <int> <int> <dbl> <chr>
## 1 5 1 0.435 t4
## 2 7 2 0.674 t1
## 3 7 3 0.00202 t3
## 4 6 4 0.0251 t2
## 5 5 5 NA <NA>
## 6 5 6 0.472 <NA>
## 7 6 7 0.274 <NA>
可以使用 as.phylo() 方法將 tbl_tree 對象轉(zhuǎn)換回 phylo 對象。
as.phylo(x)
##
## Phylogenetic tree with 4 tips and 3 internal nodes.
##
## Tip labels:
## t4, t1, t3, t2
##
## Rooted; includes branch lengths.
x <- as_tibble(tree)
x
## # A tibble: 7 × 4
## parent node branch.length label
## <int> <int> <dbl> <chr>
## 1 5 1 0.435 t4
## 2 7 2 0.674 t1
## 3 7 3 0.00202 t3
## 4 6 4 0.0251 t2
## 5 5 5 NA <NA>
## 6 5 6 0.472 <NA>
## 7 6 7 0.274 <NA>
tbl_tree對象可以使用as.phylo()方法轉(zhuǎn)換回一個phylo對象碉碉。
as.phylo(x)
##
## Phylogenetic tree with 4 tips and 3 internal nodes.
##
## Tip labels:
## t4, t1, t3, t2
##
## Rooted; includes branch lengths.
使用tbl_tree對象使樹和數(shù)據(jù)操作更有效柴钻、更容易(參見FAQ中的示例)。例如,我們可以使用dplyr命令full_join()將進(jìn)化信息與phylogeny聯(lián)系起來:
d <- tibble(label = paste0('t', 1:4),
trait = rnorm(4))
y <- full_join(x, d, by = 'label')
y
## # A tibble: 7 × 5
## parent node branch.length label trait
## <int> <int> <dbl> <chr> <dbl>
## 1 5 1 0.435 t4 0.943
## 2 7 2 0.674 t1 -0.171
## 3 7 3 0.00202 t3 0.570
## 4 6 4 0.0251 t2 -0.283
## 5 5 5 NA <NA> NA
## 6 5 6 0.472 <NA> NA
## 7 6 7 0.274 <NA> NA
2.1.2 treedata對象
tidytree包定義了treedata類來存儲帶有相關(guān)數(shù)據(jù)的系統(tǒng)發(fā)育樹垢粮。在將外部數(shù)據(jù)映射到樹結(jié)構(gòu)之后贴届，tbl_tree對象可以轉(zhuǎn)換為一個樹數(shù)據(jù)對象。
as.treedata(y)
## 'treedata' S4 object'.
##
## ...@ phylo:
##
## Phylogenetic tree with 4 tips and 3 internal nodes.
##
## Tip labels:
## t4, t1, t3, t2
##
## Rooted; includes branch lengths.
##
## with the following features available:
## 'trait'.
treeio包(Wang et al.蜡吧， 2020)中使用treedata類來存儲毫蚓，由常用軟件(BEAST、EPA昔善、HyPhy元潘、MrBayes、pml君仆、PHYLDOG翩概、PPLACER、r8s返咱、RAxML和RevBayes等)推斷的進(jìn)化證據(jù)(詳見第一章)钥庇。
tidytree包還提供了as_tibble()方法來將treedata對象轉(zhuǎn)換為tidy數(shù)據(jù)框架。系統(tǒng)發(fā)生樹結(jié)構(gòu)和進(jìn)化信息存儲在tbl_tree對象中咖摹，使得操作不同軟件推斷的進(jìn)化統(tǒng)計數(shù)據(jù)和將外部數(shù)據(jù)連接到相同的樹結(jié)構(gòu)變得一致和容易评姨。
y %>% as.treedata %>% as_tibble
## # A tibble: 7 × 5
## parent node branch.length label trait
## <int> <int> <dbl> <chr> <dbl>
## 1 5 1 0.435 t4 0.943
## 2 7 2 0.674 t1 -0.171
## 3 7 3 0.00202 t3 0.570
## 4 6 4 0.0251 t2 -0.283
## 5 5 5 NA <NA> NA
## 6 5 6 0.472 <NA> NA
## 7 6 7 0.274 <NA> NA
2.1.3 訪問相關(guān)的節(jié)點
dplyr語句可以直接應(yīng)用于tbl_tree來操作樹數(shù)據(jù)。此外萤晴，tidytree還提供了幾個語句來過濾相關(guān)節(jié)點吐句，包括child()、parent()店读、offspring()嗦枢、ancestor()、sibling()和MRCA()屯断。
這些語句接受一個tbl_tree對象和一個可以是節(jié)點號或標(biāo)簽的選定節(jié)點文虏。
child(y, 5)
## # A tibble: 2 × 5
## parent node branch.length label trait
## <int> <int> <dbl> <chr> <dbl>
## 1 5 1 0.435 t4 0.943
## 2 5 6 0.472 <NA> NA
parent(y, 2)
## # A tibble: 1 × 5
## parent node branch.length label trait
## <int> <int> <dbl> <chr> <dbl>
## 1 6 7 0.274 <NA> NA
offspring(y, 5)
## # A tibble: 6 × 5
## parent node branch.length label trait
## <int> <int> <dbl> <chr> <dbl>
## 1 5 1 0.435 t4 0.943
## 2 7 2 0.674 t1 -0.171
## 3 7 3 0.00202 t3 0.570
## 4 6 4 0.0251 t2 -0.283
## 5 5 6 0.472 <NA> NA
## 6 6 7 0.274 <NA> NA
ancestor(y, 2)
## # A tibble: 3 × 5
## parent node branch.length label trait
## <int> <int> <dbl> <chr> <dbl>
## 1 5 5 NA <NA> NA
## 2 5 6 0.472 <NA> NA
## 3 6 7 0.274 <NA> NA
sibling(y, 2)
## # A tibble: 1 × 5
## parent node branch.length label trait
## <int> <int> <dbl> <chr> <dbl>
## 1 7 3 0.00202 t3 0.570
MRCA(y, 2, 3)
## # A tibble: 1 × 5
## parent node branch.length label trait
## <int> <int> <dbl> <chr> <dbl>
## 1 6 7 0.274 <NA> NA
所有這些方法也是在treeio中實現(xiàn)的，用于處理phylo和treedata對象裹纳。您可以嘗試使用樹對象訪問相關(guān)節(jié)點。例如紧武，下面的命令將輸出所選內(nèi)部節(jié)點5的子節(jié)點：
child(tree, 5)
## [1] 1 6
注意剃氧，樹對象的方法輸出相關(guān)的節(jié)點號，而tbl_tree對象的方法輸出一個包含相關(guān)信息的tibble對象阻星。
2.2 數(shù)據(jù)集成
2.2.1 結(jié)合樹的數(shù)據(jù)
treeio包(Wang et al.朋鞍， 2020)作為一種基礎(chǔ)工具已添，可將從通用分析程序推斷出的各種類型的系統(tǒng)發(fā)育數(shù)據(jù)導(dǎo)入r中并使用。例如滥酥，由 CODEML估計的dN/dS或祖先序列更舞，以及由BEAST/MrBayes推斷出的進(jìn)化支持值(后向)。此外坎吻，treeio支持將外部數(shù)據(jù)與系統(tǒng)發(fā)育聯(lián)系起來缆蝉。它將這些外部系統(tǒng)發(fā)育數(shù)據(jù)(來自軟件輸出或外部來源)帶到R社區(qū)，并使其可用于R中的進(jìn)一步分析瘦真。此外刊头，treeio還可以將多個系統(tǒng)發(fā)育樹與它們的節(jié)點/分支特定的屬性數(shù)據(jù)組合成一個。因此诸尽，從本質(zhì)上講原杂，一個這樣的屬性(如替換率)可以映射到同一節(jié)點/分支的另一個屬性(如dN/dS)，以便進(jìn)行比較和進(jìn)一步計算(Yu et al., 2017; Yu et al., 2018))您机。
之前發(fā)表的一組數(shù)據(jù)穿肄，76個H3血凝素基因序列包含豬和人甲型流感病毒(Liang et al.， 2014)际看，在這里被用來證明比較不同軟件推斷的進(jìn)化統(tǒng)計數(shù)據(jù)的實用性咸产。利用BEAST對數(shù)據(jù)集進(jìn)行時間尺度估計，CODEML對數(shù)據(jù)集進(jìn)行同義和非同義替換估計仿村。在這個例子中锐朴，我們首先使用read.beast()函數(shù)將BEAST的輸出解析為兩個樹數(shù)據(jù)對象得哆，然后使用read.codeml()函數(shù)將CODEML的輸出解析為兩個樹數(shù)據(jù)對象馋评。然后扇谣，通過merge_tree()函數(shù)將這兩個包含獨立的特定于節(jié)點/分支的數(shù)據(jù)集的對象合并盛霎。
beast_file <- system.file("examples/MCC_FluA_H3.tree", package="ggtree")
rst_file <- system.file("examples/rst", package="ggtree")
mlc_file <- system.file("examples/mlc", package="ggtree")
beast_tree <- read.beast(beast_file)
codeml_tree <- read.codeml(rst_file, mlc_file)
merged_tree <- merge_tree(beast_tree, codeml_tree)
merged_tree
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/ggtree/examples/MCC_FluA_H3.tree',
## '/home/ygc/R/library/ggtree/examples/rst',
## '/home/ygc/R/library/ggtree/examples/mlc'.
##
## ...@ phylo:
##
## Phylogenetic tree with 76 tips and 75 internal nodes.
##
## Tip labels:
## A/Hokkaido/30-1-a/2013, A/New_York/334/2004,
## A/New_York/463/2005, A/New_York/452/1999,
## A/New_York/238/2005, A/New_York/523/1998, ...
##
## Rooted; includes branch lengths.
##
## with the following features available:
## 'height', 'height_0.95_HPD', 'height_median',
## 'height_range', 'length', 'length_0.95_HPD',
## 'length_median', 'length_range', 'posterior', 'rate',
## 'rate_0.95_HPD', 'rate_median', 'rate_range', 'subs',
## 'AA_subs', 't', 'N', 'S', 'dN_vs_dS', 'dN', 'dS',
## 'N_x_dN', 'S_x_dS'.
合并beast_tree和codeml_tree對象后告喊，從BEAST和CODEML輸出文件導(dǎo)入的所有特定于節(jié)點/分支的數(shù)據(jù)都可以在merged_tree對象中使用嘿悬。使用整潔樹包將樹對象轉(zhuǎn)換為整潔的數(shù)據(jù)幀齿兔，并將其可視化為CODEML推斷的dN/dS仆救、dN和dS與BEAST推斷的相同分支上的速率(以替換單位/站點/年計算的替換率)的hexbin散點圖云矫。
library(dplyr)
df <- merged_tree %>%
as_tibble() %>%
select(dN_vs_dS, dN, dS, rate) %>%
subset(dN_vs_dS >=0 & dN_vs_dS <= 1.5) %>%
tidyr::gather(type, value, dN_vs_dS:dS)
df$type[df$type == 'dN_vs_dS'] <- 'dN/dS'
df$type <- factor(df$type, levels=c("dN/dS", "dN", "dS"))
ggplot(df, aes(rate, value)) + geom_hex() +
facet_wrap(~type, scale='free_y')
合并beast_tree和codeml_tree對象后膳沽，從BEAST和CODEML輸出文件導(dǎo)入的所有特定于節(jié)點/分支的數(shù)據(jù)都可以在merged_tree對象中使用。使用整潔樹包將樹對象轉(zhuǎn)換為整潔的數(shù)據(jù)幀让禀，并將其可視化為CODEML推斷的dN/dS挑社、dN和dS與BEAST推斷的相同分支上的速率(以替換單位/站點/年計算的替換率)的hexbin散點圖。
輸出如圖2.1所示巡揍。然后痛阻，我們可以使用Pearson相關(guān)性測試這些節(jié)點/分支特定數(shù)據(jù)的相關(guān)性，在本例中腮敌，這表明dN和dS阱当，而不是dN/dS與速率顯著(p值)相關(guān)俏扩。
使用merge_tree()函數(shù)，我們能夠比較使用來自不同軟件包的相同模型的分析結(jié)果弊添，或者使用不同或相同軟件的不同模型的分析結(jié)果录淡。它還允許用戶集成來自不同軟件包的不同分析結(jié)果。合并樹數(shù)據(jù)并不局限于軟件發(fā)現(xiàn)油坝，也允許將外部數(shù)據(jù)與分析發(fā)現(xiàn)聯(lián)系起來嫉戚。merge_tree()函數(shù)是可鏈接的，允許幾個樹對象合并成一個免钻。
phylo <- as.phylo(beast_tree)
N <- Nnode2(phylo)
d <- tibble(node = 1:N, fake_trait = rnorm(N), another_trait = runif(N))
fake_tree <- treedata(phylo = phylo, data = d)
triple_tree <- merge_tree(merged_tree, fake_tree)
triple_tree
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/ggtree/examples/MCC_FluA_H3.tree',
## '/home/ygc/R/library/ggtree/examples/rst',
## '/home/ygc/R/library/ggtree/examples/mlc'.
##
## ...@ phylo:
##
## Phylogenetic tree with 76 tips and 75 internal nodes.
##
## Tip labels:
## A/Hokkaido/30-1-a/2013, A/New_York/334/2004,
## A/New_York/463/2005, A/New_York/452/1999,
## A/New_York/238/2005, A/New_York/523/1998, ...
##
## Rooted; includes branch lengths.
##
## with the following features available:
## 'height', 'height_0.95_HPD', 'height_median',
## 'height_range', 'length', 'length_0.95_HPD',
## 'length_median', 'length_range', 'posterior', 'rate',
## 'rate_0.95_HPD', 'rate_median', 'rate_range', 'subs',
## 'AA_subs', 't', 'N', 'S', 'dN_vs_dS', 'dN', 'dS',
## 'N_x_dN', 'S_x_dS', 'fake_trait', 'another_trait'.
##
## # The associated data tibble abstraction: 151 × 28
## # The 'node', 'label' and 'isTip' are from the phylo tree.
## node label isTip height height_0.95_HPD
## <int> <chr> <lgl> <dbl> <list>
## 1 1 A/Hokkaido/30-1-… TRUE 0 <lgl [1]>
## 2 2 A/New_York/334/2… TRUE 9 <dbl [2]>
## 3 3 A/New_York/463/2… TRUE 8 <dbl [2]>
## 4 4 A/New_York/452/1… TRUE 14 <dbl [2]>
## 5 5 A/New_York/238/2… TRUE 8 <dbl [2]>
## 6 6 A/New_York/523/1… TRUE 15 <dbl [2]>
## 7 7 A/New_York/435/2… TRUE 13 <dbl [2]>
## 8 8 A/New_York/582/1… TRUE 17 <dbl [2]>
## 9 9 A/New_York/603/1… TRUE 17 <dbl [2]>
## 10 10 A/New_York/657/1… TRUE 19 <dbl [2]>
## # … with 141 more rows, and 23 more variables:
## # height_median <dbl>, height_range <list>,
## # length <dbl>, length_0.95_HPD <list>,
## # length_median <dbl>, length_range <list>,
## # posterior <dbl>, rate <dbl>, rate_0.95_HPD <list>,
## # rate_median <dbl>, rate_range <list>, subs <chr>,
## # AA_subs <chr>, t <dbl>, N <dbl>, S <dbl>, …
上面顯示的triple_tree對象包含了從BEAST和CODEML獲得的分析結(jié)果彼水，以及從外部來源獲得的進(jìn)化特征。所有這些信息都可以使用ggtree (Yu et al.极舔， 2017)和ggtreeExtra (Xu, Dai, et al.凤覆， 2021)來注釋樹。
2.2.2 連接外部數(shù)據(jù)與系統(tǒng)發(fā)育
除上述與樹相關(guān)的分析結(jié)果外拆魏，還有廣泛的異質(zhì)性數(shù)據(jù)盯桦，包括表型數(shù)據(jù)、實驗數(shù)據(jù)和臨床數(shù)據(jù)等渤刃，需要整合并與系統(tǒng)發(fā)育聯(lián)系起來拥峦。例如，在病毒進(jìn)化的研究中卖子，樹節(jié)點可能與流行病學(xué)信息相關(guān)略号，如位置、年齡和亞型洋闽。功能注釋可能需要映射到基因樹玄柠，以進(jìn)行比較基因組學(xué)研究诫舅。為了便于數(shù)據(jù)集成虚汛，treeio提供了full_join()方法來將外部數(shù)據(jù)鏈接到系統(tǒng)發(fā)育，并將其存儲在phylo或treedata對象中瓢娜。請注意，將外部數(shù)據(jù)鏈接到一個門系對象將產(chǎn)生一個樹數(shù)據(jù)對象來存儲具有關(guān)聯(lián)數(shù)據(jù)的輸入門系谷丸。full_join方法也可以在整潔的數(shù)據(jù)幀級(即前面描述的tbl_tree對象)和ggtree級(在第7章中描述)使用 (Yu et al., 2018)揩慕。
下面的示例計算引導(dǎo)值蜗搔，并通過匹配它們的節(jié)點號將這些值與樹(一個門系對象)合并厉膀。
library(ape)
data(woodmouse)
d <- dist.dna(woodmouse)
tr <- nj(d)
bp <- boot.phylo(tr, woodmouse, function(x) nj(dist.dna(x)))
## 'treedata' S4 object'.
##
## ...@ phylo:
##
## Phylogenetic tree with 15 tips and 13 internal nodes.
##
## Tip labels:
## No305, No304, No306, No0906S, No0908S, No0909S, ...
##
## Unrooted; includes branch lengths.
##
## with the following features available:
## 'bootstrap'.
##
## # The associated data tibble abstraction: 28 × 4
## # The 'node', 'label' and 'isTip' are from the phylo tree.
## node label isTip bootstrap
## <int> <chr> <lgl> <int>
## 1 1 No305 TRUE NA
## 2 2 No304 TRUE NA
## 3 3 No306 TRUE NA
## 4 4 No0906S TRUE NA
## 5 5 No0908S TRUE NA
## 6 6 No0909S TRUE NA
## 7 7 No0910S TRUE NA
## 8 8 No0912S TRUE NA
## 9 9 No0913S TRUE NA
## 10 10 No1103S TRUE NA
## # … with 18 more rows
另一個例子演示了通過匹配它們的提示標(biāo)簽來將進(jìn)化特征與樹(一個treedata對象)合并。
file <- system.file("extdata/BEAST", "beast_mcc.tree", package="treeio")
beast <- read.beast(file)
x <- tibble(label = as.phylo(beast)$tip.label, trait = rnorm(Ntip(beast)))
full_join(beast, x, by="label")
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/BEAST/beast_mcc.tree'.
##
## ...@ phylo:
##
## Phylogenetic tree with 15 tips and 14 internal nodes.
##
## Tip labels:
## A_1995, B_1996, C_1995, D_1987, E_1996, F_1997, ...
##
## Rooted; includes branch lengths.
##
## with the following features available:
## 'height', 'height_0.95_HPD', 'height_median',
## 'height_range', 'length', 'length_0.95_HPD',
## 'length_median', 'length_range', 'posterior', 'rate',
## 'rate_0.95_HPD', 'rate_median', 'rate_range', 'trait'.
##
## # The associated data tibble abstraction: 29 × 17
## # The 'node', 'label' and 'isTip' are from the phylo tree.
## node label isTip height height_0.95_HPD
## <int> <chr> <lgl> <dbl> <list>
## 1 1 A_1995 TRUE 18 <dbl [2]>
## 2 2 B_1996 TRUE 17 <dbl [2]>
## 3 3 C_1995 TRUE 18 <dbl [2]>
## 4 4 D_1987 TRUE 26 <dbl [2]>
## 5 5 E_1996 TRUE 17 <dbl [2]>
## 6 6 F_1997 TRUE 16 <dbl [2]>
## 7 7 G_1992 TRUE 21 <dbl [2]>
## 8 8 H_1992 TRUE 21 <dbl [2]>
## 9 9 I_1994 TRUE 19 <dbl [2]>
## 10 10 J_1983 TRUE 30 <dbl [2]>
## # … with 19 more rows, and 12 more variables:
## # height_median <dbl>, height_range <list>,
## # length <dbl>, length_0.95_HPD <list>,
## # length_median <dbl>, length_range <list>,
## # posterior <dbl>, rate <dbl>, rate_0.95_HPD <list>,
## # rate_median <dbl>, rate_range <list>, trait <dbl>
操縱樹對象會因為用于處理對象格式的碎片化函數(shù)而難以處理，更不用說將外部數(shù)據(jù)連接到系統(tǒng)發(fā)育結(jié)構(gòu)了旨椒。使用treeio包 (Wang et al., 2020)，可以很容易地組合來自不同來源的樹數(shù)據(jù)阔拳。此外嗤形，使用tidytree包，使用整潔的數(shù)據(jù)原理操作樹更容易吭产，并且與已經(jīng)廣泛使用的工具(包括dplyr、tidyr臣淤、ggplot2和ggtree)一致 (Yu et al., 2017)达吞。
2.2.3分組分類單元
groupOTU()和groupClade()方法被設(shè)計用來向輸入樹對象添加分類單元分組信息。這些方法分別在tidytree荒典、treeio和ggtree中實現(xiàn)，以支持為tbl_tree吞鸭、phylo和tredata以及ggtree對象添加分組信息寺董。這個分組信息可以在樹的可視化中直接使用ggtree(圖6.5)(e.g., coloring a tree based on grouping information)。
2.2.3.1 groupClade
groupClade()方法接受一個內(nèi)部節(jié)點或內(nèi)部節(jié)點的向量刻剥，以添加所選clade/clades的分組信息遮咖。
nwk <- '(((((((A:4,B:4):6,C:5):8,D:6):3,E:21):10,((F:4,G:12):14,H:8):13):
13,((I:5,J:2):30,(K:11,L:11):2):17):4,M:56);'
tree <- read.tree(text=nwk)
groupClade(as_tibble(tree), c(17, 21))
## # A tibble: 25 × 5
## parent node branch.length label group
## <int> <int> <dbl> <chr> <fct>
## 1 20 1 4 A 1
## 2 20 2 4 B 1
## 3 19 3 5 C 1
## 4 18 4 6 D 1
## 5 17 5 21 E 1
## 6 22 6 4 F 2
## 7 22 7 12 G 2
## 8 21 8 8 H 2
## 9 24 9 5 I 0
## 10 24 10 2 J 0
## # … with 15 more rows
2.2.3.2 groupOTU
set.seed(2017)
tr <- rtree(4)
x <- as_tibble(tr)
## the input nodes can be node ID or label
groupOTU(x, c('t1', 't4'), group_name = "fake_group")
## # A tibble: 7 × 5
## parent node branch.length label fake_group
## <int> <int> <dbl> <chr> <fct>
## 1 5 1 0.435 t4 1
## 2 7 2 0.674 t1 1
## 3 7 3 0.00202 t3 0
## 4 6 4 0.0251 t2 0
## 5 5 5 NA <NA> 1
## 6 5 6 0.472 <NA> 1
## 7 6 7 0.274 <NA> 1
groupClade()和groupOTU()都可以使用tbl_tree、phylo和tredata以及ggtree對象造虏。下面是一個將groupOTU()與一個 phylo tree 對象一起使用的例子御吞。
groupOTU(tr, c('t2', 't4'), group_name = "fake_group") %>%
as_tibble
## # A tibble: 7 × 5
## parent node branch.length label fake_group
## <int> <int> <dbl> <chr> <fct>
## 1 5 1 0.435 t4 1
## 2 7 2 0.674 t1 0
## 3 7 3 0.00202 t3 0
## 4 6 4 0.0251 t2 1
## 5 5 5 NA <NA> 1
## 6 5 6 0.472 <NA> 1
## 7 6 7 0.274 <NA> 0
另一個使用ggtree對象的例子可以在6.4章節(jié)中找到。
groupOTU將從輸入節(jié)點回溯到最近的公共祖先漓藕。在這個例子中陶珠，節(jié)點1、4享钞、5和6被分組在一起( (4 (t2) -> 6 -> 5 and 1 (t4) -> 5).
相關(guān)的操作分類單位(OTUs)是分組的揍诽，它們不一定在一個進(jìn)化支內(nèi)。它們可以是單系的(枝)栗竖，多系的或副系的暑脆。
cls <- list(c1=c("A", "B", "C", "D", "E"),
c2=c("F", "G", "H"),
c3=c("L", "K", "I", "J"),
c4="M")
as_tibble(tree) %>% groupOTU(cls)
## # A tibble: 25 × 5
## parent node branch.length label group
## <int> <int> <dbl> <chr> <fct>
## 1 20 1 4 A c1
## 2 20 2 4 B c1
## 3 19 3 5 C c1
## 4 18 4 6 D c1
## 5 17 5 21 E c1
## 6 22 6 4 F c2
## 7 22 7 12 G c2
## 8 21 8 8 H c2
## 9 24 9 5 I c3
## 10 24 10 2 J c3
## # … with 15 more rows
如果在回溯到最近的共同祖先時出現(xiàn)沖突，用戶可以將 overlap 參數(shù)設(shè)置為origin(第一個計數(shù))狐肢、overwrite(默認(rèn)值添吗，最后一個計數(shù))或abandon(取消選擇以便分組)。
2.3 樹的置根
系統(tǒng)發(fā)生樹可以用指定的外群重新根化份名。ape包實現(xiàn)了一個root()方法來重新根化存儲在phylo對象中的樹碟联，而treeio包為treedata對象提供了root()方法。此方法設(shè)計用于使用與指定外組相關(guān)的數(shù)據(jù)或基于在 ape 包中實現(xiàn)的root()在指定節(jié)點上重新建立系統(tǒng)發(fā)育樹的根僵腺。
我們首先使用left_join()將外部數(shù)據(jù)鏈接到樹玄帕，并將所有信息存儲在一個樹數(shù)據(jù)對象trda中。
library(ggtree)
library(treeio)
library(tidytree)
library(TDbook)
# load `tree_boots`, `df_tip_data`, and `df_inode_data` from 'TDbook'
trda <- tree_boots %>%
left_join(df_tip_data, by=c("label" = "Newick_label")) %>%
left_join(df_inode_data, by=c("label" = "newick_label"))
trda
## 'treedata' S4 object'.
##
## ...@ phylo:
##
## Phylogenetic tree with 7 tips and 6 internal nodes.
##
## Tip labels:
## Rangifer_tarandus, Cervus_elaphus, Bos_taurus,
## Ovis_orientalis, Suricata_suricatta,
## Cystophora_cristata, ...
## Node labels:
## Mammalia, Artiodactyla, Cervidae, Bovidae,
## Carnivora, Caniformia
##
## Rooted; includes branch lengths.
##
## with the following features available:
## '', 'vernacularName', 'imageURL', 'imageLicense',
## 'imageAuthor', 'infoURL', 'mass_in_kg',
## 'trophic_habit', 'ncbi_taxid', 'rank',
## 'vernacularName.y', 'infoURL.y', 'rank.y',
## 'bootstrap', 'posterior'.
##
## # The associated data tibble abstraction: 13 × 17
## # The 'node', 'label' and 'isTip' are from the phylo tree.
## node label isTip vernacularName imageURL
## <int> <chr> <lgl> <chr> <chr>
## 1 1 Rangifer_tarand… TRUE Reindeer http://…
## 2 2 Cervus_elaphus TRUE Red deer http://…
## 3 3 Bos_taurus TRUE Cattle https:/…
## 4 4 Ovis_orientalis TRUE Asiatic moufl… http://…
## 5 5 Suricata_surica… TRUE Meerkat http://…
## 6 6 Cystophora_cris… TRUE Hooded seal http://…
## 7 7 Mephitis_mephit… TRUE Striped skunk http://…
## 8 8 Mammalia FALSE <NA> <NA>
## 9 9 Artiodactyla FALSE <NA> <NA>
## 10 10 Cervidae FALSE <NA> <NA>
## # … with 3 more rows, and 12 more variables:
## # imageLicense <chr>, imageAuthor <chr>,
## # infoURL <chr>, mass_in_kg <dbl>,
## # trophic_habit <chr>, ncbi_taxid <int>, rank <chr>,
## # vernacularName.y <chr>, infoURL.y <chr>,
## # rank.y <chr>, bootstrap <int>, posterior <dbl>
然后想邦，我們可以用關(guān)聯(lián)到分支和節(jié)點的數(shù)據(jù)映射來重新根化樹裤纹，如圖2.2所示。這個圖是使用ggtree實現(xiàn)的(參見第4章和第5章)。
# reroot
trda2 <- root(trda, outgroup = "Suricata_suricatta", edgelabel = TRUE)
# The original tree
p1 <- trda %>%
ggtree() +
geom_nodelab(
mapping = aes(
x = branch,
label = bootstrap
),
nudge_y = 0.36
) +
xlim(-.1, 4) +
geom_tippoint(
mapping = aes(
shape = trophic_habit,
color = trophic_habit,
size = mass_in_kg
)
) +
scale_size_continuous(range = c(3, 10)) +
geom_tiplab(
offset = .14,
) +
geom_nodelab(
mapping = aes(
label = vernacularName.y,
fill = posterior
),
geom = "label"
) +
scale_fill_gradientn(colors = RColorBrewer::brewer.pal(3, "YlGnBu")) +
theme(legend.position = "right")
# after reroot
p2 <- trda2 %>%
ggtree() +
geom_nodelab(
mapping = aes(
x = branch,
label = bootstrap
),
nudge_y = 0.36
) +
xlim(-.1, 5) +
geom_tippoint(
mapping = aes(
shape = trophic_habit,
color = trophic_habit,
size = mass_in_kg
)
) +
scale_size_continuous(range = c(3, 10)) +
geom_tiplab(
offset = .14,
) +
geom_nodelab(
mapping = aes(
label = vernacularName.y,
fill = posterior
),
geom = "label"
) +
scale_fill_gradientn(colors = RColorBrewer::brewer.pal(3, "YlGnBu")) +
theme(legend.position = "right")
plot_list(p1, p2, tag_levels='A', ncol=2)
圖2.2:對一個樹進(jìn)行置根鹰椒。原始樹(A)和重根樹(B)與相關(guān)的數(shù)據(jù)正確地映射到樹的分支或節(jié)點锡移。(A)和(B)分別出現(xiàn)在通向根尖節(jié)Suricata_suricatta的分枝生根前和生根后。
outgroup參數(shù)表示特定的新outgroup漆际，它可以是節(jié)點標(biāo)簽(字符)或節(jié)點編號淆珊。如果它是一個值，這意味著使用本技巧下面的節(jié)點作為新根;如果它有多個值奸汇，這意味著最近的常見值將被用作新根施符。注意，如果節(jié)點標(biāo)簽應(yīng)該被視為邊緣標(biāo)簽擂找，則應(yīng)該將edgelabel設(shè)置為TRUE戳吝，以返回節(jié)點和相關(guān)數(shù)據(jù)之間的正確關(guān)系。有關(guān)re-root的更多細(xì)節(jié)贯涎，包括注意事項和陷阱听哭，請參閱綜述文章 (Czech et al., 2017)
。
2.4 調(diào)節(jié)分支
系統(tǒng)發(fā)育數(shù)據(jù)可以合并進(jìn)行聯(lián)合分析(圖2.1)塘雳。它們可以作為更復(fù)雜的注釋顯示在相同的樹結(jié)構(gòu)上陆盘，以幫助可視化地檢查它們的進(jìn)化模式。存儲在樹數(shù)據(jù)對象中的所有數(shù)值數(shù)據(jù)都可以用于重新縮放樹分支败明。例如隘马，CodeML推斷dN/dS、dN和dS妻顶，所有這些統(tǒng)計數(shù)據(jù)都可以用作分支長度(圖2.3)祟霍。所有這些值也可以用來為樹著色(會話4.3.4)，并且可以投影到垂直維度以創(chuàng)建一個二維樹或表現(xiàn)型圖(會話4.2.2和圖4.5和4.11)盈包。
p1 <- ggtree(merged_tree) + theme_tree2()
p2 <- ggtree(rescale_tree(merged_tree, 'dN')) + theme_tree2()
p3 <- ggtree(rescale_tree(merged_tree, 'rate')) + theme_tree2()
plot_list(p1, p2, p3, ncol=3, tag_levels='A')
圖2.3:重新縮放樹枝樹的分支按時間(從根開始的年份)縮放(A)沸呐。使用dN作為分支長度(B)重新縮放樹。使用替換率(C)重新縮放樹呢燥。
2.5用數(shù)據(jù)劃分樹的子集
2.5.1在系統(tǒng)發(fā)育樹中去除tips
有時崭添，我們希望從系統(tǒng)發(fā)育樹中刪除所選擇的tps。這是由于序列質(zhì)量低叛氨、序列組裝錯誤呼渣、部分序列比對錯誤、系統(tǒng)發(fā)育推斷錯誤等原因造成的寞埠。
假設(shè)我們想從樹中刪除三個提示(紅色)(圖2.4A)屁置， drop.tip()方法刪除指定的提示并更新樹(圖2.4B)。所有相關(guān)的數(shù)據(jù)將在更新后的樹中進(jìn)行維護(hù)仁连。
f <- system.file("extdata/NHX", "phyldog.nhx", package="treeio")
nhx <- read.nhx(f)
to_drop <- c("Physonect_sp_@2066767",
"Lychnagalma_utricularia@2253871",
"Kephyes_ovata@2606431")
p1 <- ggtree(nhx) + geom_tiplab(aes(color = label %in% to_drop)) +
scale_color_manual(values=c("black", "red")) + xlim(0, 0.8)
nhx_reduced <- drop.tip(nhx, to_drop)
p2 <- ggtree(nhx_reduced) + geom_tiplab() + xlim(0, 0.8)
plot_list(p1, p2, ncol=2, tag_levels = "A")
圖2.4:移除樹中的tips梢蓝角。原始樹有三個tips(紅色)要刪除(A)阱穗。更新后的樹刪除選中的tips(B)。
2.5.2根據(jù)tip laber 劃分樹的子集
有時一棵樹可能很大使鹅，很難只看到感興趣的部分揪阶。tree_子集()函數(shù)是在treeio包中創(chuàng)建的(Wang et al.， 2020)患朱，以提取樹部分的子集鲁僚，同時仍然保持樹部分的結(jié)構(gòu)。圖2.5A中的beast_tree有點擁擠裁厅。顯然冰沙，我們可以讓圖形變高以給標(biāo)簽留出更多的空間(類似于在FigTree中使用擴展滑塊)，或者我們可以讓文本變小执虹。然而拓挥，當(dāng)你有很多 tips 時(例如，成百上千的tips)声畏，這些解決方案并不總是適用的。特別是姻成，當(dāng)您只對樹的特定tips周圍的部分感興趣時插龄。
假設(shè)你對樹中的tip /Swine/HK/168/2012 感興趣(圖2.5A)，你想看看這個tips的直系親屬科展。
默認(rèn)情況下均牢，tree_subset()函數(shù)將在內(nèi)部調(diào)用groupOTU()來從其他提示中分配指定的組提示(圖2.5B)。此外才睹，分支長度和相關(guān)的數(shù)據(jù)在子集設(shè)置后進(jìn)行維護(hù)(圖2.5C)徘跪。默認(rèn)情況下，樹的根總是錨定在子集樹的零的位置琅攘，所有的距離都是相對于這個根的垮庐。如果您希望所有的距離都相對于原始根，您可以指定根位置(by root.position parameter)坞琴。位置參數(shù))到子集樹根邊哨查，即從原始根到子集樹根的分支長度之和(圖2.5D和E)。
beast_file <- system.file("examples/MCC_FluA_H3.tree", package="ggtree")
beast_tree <- read.beast(beast_file)
p1 = ggtree(beast_tree) +
geom_tiplab() + xlim(0, 40) + theme_tree2()
tree2 = tree_subset(beast_tree, "A/Swine/HK/168/2012", levels_back=4)
p2 <- ggtree(tree2, aes(color=group)) +
scale_color_manual(values = c("black", "red"), guide = 'none') +
geom_tiplab() + xlim(0, 4) + theme_tree2()
p3 <- p2 +
geom_point(aes(fill = rate), shape = 21, size = 4) +
scale_fill_continuous(low = 'blue', high = 'red') +
xlim(0,5) + theme(legend.position = 'right')
p4 <- ggtree(tree2, aes(color=group),
root.position = as.phylo(tree2)$root.edge) +
geom_tiplab() + xlim(18, 24) +
scale_color_manual(values = c("black", "red"), guide = 'none') +
theme_tree2()
p5 <- p4 +
geom_rootedge() + xlim(0, 40)
plot_list(p1, p2, p3, p4, p5,
design="AABBCC\nAADDEE", tag_levels='A')
圖2.5:特定提示的子集樹剧辐。原始樹(A).子集樹(B).帶數(shù)據(jù)的子集樹(C).將子集樹相對于原始位置可視化寒亥，(D)不帶根的樹，(E)有根樹荧关。
2.5.3根據(jù)內(nèi)部節(jié)點編號劃分子樹
tree_subset()函數(shù)將整個進(jìn)化枝也追溯到特定的水平看進(jìn)化枝的相關(guān)類群(圖2.6 a和B)。我們可以使用tree_subset()函數(shù)來放大所選部分和整個樹的一部分,它類似于ape::zoom()函數(shù)來查看一個非常大的樹(圖2.6C和D)。用戶也可以使用viewClade()函數(shù)來限制樹在特定演化枝上的可視化柒室，如session 6.1所示胎食。
clade <- tree_subset(beast_tree, node=121, levels_back=0)
clade2 <- tree_subset(beast_tree, node=121, levels_back=2)
p1 <- ggtree(clade) + geom_tiplab() + xlim(0, 5)
p2 <- ggtree(clade2, aes(color=group)) + geom_tiplab() +
xlim(0, 8) + scale_color_manual(values=c("black", "red"))
library(ape)
library(tidytree)
library(treeio)
data(chiroptera)
nodes <- grep("Plecotus", chiroptera$tip.label)
chiroptera <- groupOTU(chiroptera, nodes)
clade <- MRCA(chiroptera, nodes)
x <- tree_subset(chiroptera, clade, levels_back = 0)
p3 <- ggtree(chiroptera, aes(colour = group)) +
scale_color_manual(values=c("black", "red")) +
theme(legend.position = "none")
p4 <- ggtree(x) + geom_tiplab() + xlim(0, 5)
plot_list(p1, p2, p3, p4,
ncol=2, tag_levels = 'A')
提取一個枝支(A)。提取一個枝支并追溯它以查看它的近親(B)。查看一棵非常大的樹(C)和它的選定部分(D)胸竞。
2.6 操縱樹進(jìn)行可視化
盡管ggtree實現(xiàn)了幾種可視化的數(shù)據(jù)樹探索方法欺嗤，但您可能希望執(zhí)行一些不直接支持的操作。在本例中卫枝，您需要使用用于可視化的節(jié)點協(xié)調(diào)位置來操作樹數(shù)據(jù)煎饼。這對于ggtree來說非常簡單。用戶可以使用內(nèi)部調(diào)用 tidytree::as_tibble() 的 fortify() 方法將樹轉(zhuǎn)換為整潔的數(shù)據(jù)幀校赤，并添加用于繪制樹的坐標(biāo)位置列(例如吆玖，x、y马篮、分支和角度)沾乘。您也可以通過ggtree(tree)$data訪問數(shù)據(jù)。
下面是一個繪制兩棵面對面樹的示例浑测，類似于由ape::cophyloplot()函數(shù)生成的圖(圖2.7)翅阵。
library(dplyr)
library(ggtree)
set.seed(1024)
x <- rtree(30)
y <- rtree(30)
p1 <- ggtree(x, layout='roundrect') +
geom_hilight(
mapping=aes(subset = node %in% c(38, 48, 58, 36),
node = node,
fill = as.factor(node)
)
) +
labs(fill = "clades for tree in left" )
p2 <- ggtree(y)
d1 <- p1$data
d2 <- p2$data
## 以x軸翻轉(zhuǎn)，并設(shè)置偏移量迁央，使第二顆樹位于第一棵樹的右側(cè)
d2$x <- max(d2$x) - d2$x + max(d1$x) + 1
pp <- p1 + geom_tree(data=d2, layout='ellipse') +
ggnewscale::new_scale_fill() +
geom_hilight(
data = d2,
mapping = aes(
subset = node %in% c(38, 48, 58),
node=node,
fill=as.factor(node))
) +
labs(fill = "clades for tree in right" )
dd <- bind_rows(d1, d2) %>%
filter(!is.na(label))
pp + geom_line(aes(x, y, group=label), data=dd, color='grey') +
geom_tiplab(geom = 'shadowtext', bg.colour = alpha('firebrick', .5)) +
geom_tiplab(data = d2, hjust=1, geom = 'shadowtext',
bg.colour = alpha('firebrick', .5))
圖2.7:面對面地繪制兩棵系統(tǒng)發(fā)育樹掷匠。使用ggtree()繪制樹(左邊)，然后通過geom_tree()添加樹的另一層(右邊)岖圈。繪制的樹的相對位置可以手動調(diào)整讹语，并且向每棵樹添加圖層(例如，提示標(biāo)簽和突出顯示分支)是獨立的蜂科。
在一個圖中連接多棵樹中對應(yīng)的分類單元也比較容易; 例如顽决，繪制由流感病毒所有內(nèi)部基因片段構(gòu)建的樹形圖，并在樹形上連接等效菌株 (Venkatesh et al., 2018))导匣。圖2.8演示了通過組合多個geom_tree()層來繪制多棵樹的用法才菠。
z <- rtree(30)
d2 <- fortify(y)
d3 <- fortify(z)
d2$x <- d2$x + max(d1$x) + 1
d3$x <- d3$x + max(d2$x) + 1
dd = bind_rows(d1, d2, d3) %>%
filter(!is.na(label))
p1 + geom_tree(data = d2) + geom_tree(data = d3) + geom_tiplab(data=d3) +
geom_line(aes(x, y, group=label, color=node < 15), data=dd, alpha=.3)
圖2.8:并排繪制多個系統(tǒng)發(fā)育樹。使用ggtree()繪制樹贡定，然后使用geom_tree()添加多層樹鸠儿。
2.7 總結(jié)
treeio包允許我們將不同的系統(tǒng)發(fā)育相關(guān)數(shù)據(jù)導(dǎo)入R，然而厕氨，系統(tǒng)發(fā)育樹的存儲方式是為了便于計算處理进每，需要專門知識來操作和探索樹數(shù)據(jù)。tidytree包為探索樹數(shù)據(jù)提供了一個tidy界面命斧，而ggtree提供了一組實用工具田晚，可以使用圖形語法可視化和探索樹數(shù)據(jù)。這種完整的軟件包使得普通用戶可以很容易地與樹狀數(shù)據(jù)進(jìn)行交互国葬，并允許我們整合來自不同來源的系統(tǒng)發(fā)育相關(guān)數(shù)據(jù)(例如贤徒，實驗結(jié)果或分析結(jié)果)芹壕，這為整合和比較研究創(chuàng)造了可能性。此外接奈，這個軟件包套件將系統(tǒng)發(fā)育分析帶入了tidyverse界面踢涌，并將我們帶到處理系統(tǒng)發(fā)育數(shù)據(jù)的新的層次。
