研究背景

ape (Analyses of Phylogenetics and Evolution) 包是由Emmanuel Paradis等人开发的R包，于2004年发表在Bioinformatics上 (Paradis et al., 2004)。到目前为止，这篇论文在谷歌学术上引用超过10,000次。现在ape的最新版本为5.5 (Paradis and Schliep, 2019)，相应论文在谷歌学术上引用超过2,000次。此外，Paradis于2012年编写了《Analysis of Phylogenetics and Evolution with R》一书，介绍了ape 包的使用方法。不过对于R语言比较陌生的读者，使用ape包多少会遇到一些困难。本教程对ape包的常见函数进行了系统地整理，可供读者参考。

软件和相关软件包

1. R software version 4.0 (https://www.r-project.org/)，官网下载安装。
2. Rstudio version 1.4 (https://rstudio.com/products/rstudio/)，官网下载安装。
3. ape package version 5.5 (Paradis and Schliep, 2019)，ape 包的安装直接使用 install.packages(“ape”) 命令即可。
   

分析流程

ape 包涉及的函数很多，本教程重点介绍ape包在以下几个方面的应用：

1. 序列的获取，
2. 序列的读取和写入，
3. 序列特征分析，
4. 序列比对，
5. 序列的可视化，
6. 计算遗传距离，
7.  构建进化树，
8. 树的读取和写入，
9. 树的常见操作，
10. 树的可视化。

1. 随机生成序列
   ape 包的rDNAbin函数可以随机生成DNA序列。下面的代码可以随机生成10条序列长度为100 bp 的DNA序列。随机生成的序列一般用于数据模拟及测试。
   ### 加载 ape 包
   library(ape)
   ### 随机生成10条长度为100 bp 的序列
   DNAbin <- rDNAbin(rep(100, 10))
   ### 展示序列结构
   str(DNAbin)
   ## List of 10
   ## $ Ind_1 : raw [1:100] 48 28 48 88 ...
   ## $ Ind_2 : raw [1:100] 88 48 18 48 ...
   ## $ Ind_3 : raw [1:100] 88 28 88 48 ...
   ## $ Ind_4 : raw [1:100] 48 48 48 18 ...
   ## $ Ind_5 : raw [1:100] 88 18 48 88 ...
   ## $ Ind_6 : raw [1:100] 48 48 28 28 ...
   ## $ Ind_7 : raw [1:100] 28 88 48 88 ...
   ## $ Ind_8 : raw [1:100] 48 88 28 48 ...
   ## $ Ind_9 : raw [1:100] 88 48 18 88 ...
   ## $ Ind_10: raw [1:100] 88 28 48 88 ...
   ## - attr(*, "class")= chr "DNAbin"
   
   
   
2. 使用自带数据
   ape 包有一些自带数据，woodmouse 就是其中之一，它包含15条长度为 965bp 的 Cytb 基因序列。使用时可以直接用 data(woodmouse) 来导入。用 str 函数可以展示 woodmouse 的数据结构，用 alview 函数可以直接显示序列的碱基。
   ### 载入 woodmouse 数据
   data(woodmouse)
   ### 展示数据结构
   str(woodmouse)
   ## 'DNAbin' raw [1:15, 1:965] n a a a ...
   ## - attr(*, "dimnames")=List of 2
   ## ..$ : chr [1:15] "No305" "No304" "No306" "No0906S" ...
   ## ..$ : NULL
   ### 显示前 60bp 序列
   alview(woodmouse[, 1:60])
   ## 000000000111111111122222222223333333333444444444455555555556
   ## 123456789012345678901234567890123456789012345678901234567890
   ## No305 NTTCGAAAAACACACCCACTACTAAAANTTATCAGTCACTCCTTCATCGACTTACCAGCT
   ## No304 A..........................A......AC........................
   ## No306 A..........................A......A.........................
   ## No0906S A..........................A.C....A...............T.........
   ## No0908S A..........................A......A........................C
   ## No0909S A..........................A......A.........................
   ## No0910S A..........................A......A......T........T.........
   ## No0912S A..........................A......A.........................
   ## No0913S A..........................A......AC........................
   ## No1103S A..........................A....T.A.........................
   ## No1007S A..........................A......A.........................
   ## No1114S .NNNNNNNNNNNNNNNNNNNNNNNNNN.NNNNNNNNNNNNNNNNN...............
   ## No1202S A..........................A......A...............T.........
   ## No1206S A..........................A......A...............T..G......
   ## No1208S .NN........................A......A.........................
   通过观察我们可以发现woodmouse前60 bp 的序列中， No1114S 含有多个连续 “N” 碱基。
   
   
   
3. 下载GenBank的分子序列
   ape 包的自带数据 woodmouse，只是一部分数据。woodmouse 的完整数据集可以通过 read.GenBank 函数下载。read.GenBank 函数可以使用 GenBank 登录号（accession numbers）来下载 GenBank 中的数据， 此函数能否运行成功取决于网络情况。ape 包自带数据选自 Michaux et al. (2003) 发表的论文。文章给出了 woodmouse 完整数据的登录号：AJ511877-AJ511987。接下来可以通过以下步骤进行下载：
   ### 生成GenBank 登录号
   acc_no <- paste0("AJ", 511877:511987)
   ### 下载序列
   woodmouse_all <- read.GenBank(acc_no)
   ## Warning in read.GenBank(acc_no): cannot get the following sequence(s):
   ## AJ511888, AJ511894, AJ511895, AJ511927, AJ511933, AJ511934
   ### 获取物种名称 (前10个)
   attr(woodmouse_all, "species")[1:10]
   ## [1] "Apodemus_sylvaticus" "Apodemus_sylvaticus" "Apodemus_sylvaticus"
   ## [4] "Apodemus_sylvaticus" "Apodemus_sylvaticus" "Apodemus_sylvaticus"
   ## [7] "Apodemus_sylvaticus" "Apodemus_sylvaticus" "Apodemus_sylvaticus"
   ## [10] "Apodemus_sylvaticus"
   ### 得到序列的描述信息 (前10个)
   attr(woodmouse_all, "description")[1:10]
   ## [1] "AJ511877.1 Apodemus sylvaticus mitochondrial partial cytb gene for cytochrome b, tissue library JRM-104"
   ## [2] "AJ511878.1 Apodemus sylvaticus mitochondrial partial cytb gene for cytochrome b, tissue library JRM-105"
   ## [3] "AJ511879.1 Apodemus sylvaticus mitochondrial partial cytb gene for cytochrome b, tissue library JRM-107"
   ## [4] "AJ511880.1 Apodemus sylvaticus mitochondrial partial cytb gene for cytochrome b, tissue library JRM-131"
   ## [5] "AJ511881.1 Apodemus sylvaticus mitochondrial partial cytb gene for cytochrome b, tissue library JRM-132"
   ## [6] "AJ511882.1 Apodemus sylvaticus mitochondrial partial cytb gene for cytochrome b, tissue library JRM-133"
   ## [7] "AJ511883.1 Apodemus sylvaticus mitochondrial partial cytb gene for cytochrome b, tissue library JRM-144"
   ## [8] "AJ511884.1 Apodemus sylvaticus mitochondrial partial cytb gene for cytochrome b, tissue library JRM-145"
   ## [9] "AJ511885.1 Apodemus sylvaticus mitochondrial partial cytb gene for cytochrome b, tissue library JRM-270"
   ## [10] "AJ511886.1 Apodemus sylvaticus mitochondrial partial cytb gene for cytochrome b, tissue library JRM-271, 157BP"
   通过下载我们可以得到105条序列，在这里我们要筛选序列长度为 965bp 的序列，可以通过下列代码实现：
   ### 计算woodmouse_all序列的长度
   len <- sapply(woodmouse_all, length)
   ### 筛选序列长度为 965 bp 的序列
   woodmouse_len <- woodmouse_all[grep("965", len)]
   ### 随机选择15条序列
   set.seed(12345)
   woodmouse_sel <- sample(woodmouse_len, 15)
   woodmouse_sel
   ## 15 DNA sequences in binary format stored in a list.
   ##
   ## All sequences of same length: 965
   ##
   ## Labels:
   ## AJ511910
   ## AJ511912
   ## AJ511930
   ## AJ511935
   ## AJ511928
   ## AJ511940
   ## ...
   ##
   ## Base composition:
   ## a c g t
   ## 0.309 0.259 0.124 0.307
   ## (Total: 14.47 kb)
   data(woodmouse)
   identical(woodmouse_sel, woodmouse)
   ## [1] FALSE
   需要注意的是，由于是随机选取的15条序列，woodmouse_sel 和 woodmouse 数据并不完全一致。
   

二、序列的读取和写入
使用 write.FASTA，write.dna 或 write.nexus.data 函数可以写入DNA 序列。使用函数 read.FASTA，read.dna和 read.nexus.data 可以读取DNA序列。接下来我们可以把上一步获得的 woodmouse_sel数据保存到本地。
### FASTA 格式
write.FASTA(woodmouse_sel, "woodmouse_sel.fas")
### NEXUS 格式
write.nexus.data(woodmouse_sel, file = "woodmouse_sel.nex", format = "dna")
### 读取序列，FASTA 格式
woodmouse_sel <- read.FASTA("woodmouse_sel.fas")
woodmouse_sel
## 15 DNA sequences in binary format stored in a list.
##
## All sequences of same length: 965
##
## Labels:
## AJ511910
## AJ511912
## AJ511930
## AJ511935
## AJ511928
## AJ511940
## ...
##
## Base composition:
## a c g t
## 0.309 0.259 0.124 0.307
## (Total: 14.47 kb)
### 读取序列，NEXUS 格式
woodmouse_sel_nex <- read.nexus.data("woodmouse_sel.nex")
str(woodmouse_sel_nex)
## List of 15
## $ AJ511910: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511912: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511930: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511935: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511928: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511940: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511906: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511937: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511879: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511921: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511938: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511891: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511897: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511904: chr [1:965] "a" "t" "t" "c" ...
## $ AJ511914: chr [1:965] "a" "t" "t" "c" ...

1. 计算DNA序列的碱基频率
   接下来可以使用 base.freq 函数计算 woodmouse_sel 序列的碱基频率，使用 GC.content 函数计算 woodmouse_sel 序列的GC含量。
   ### 计算四种碱基的频率
   base.freq(woodmouse_sel, freq = FALSE)
   ## a c g t
   ## 0.3091526 0.2591977 0.1242292 0.3074205
   ### 显示四种碱基的频数
   base.freq(woodmouse_sel, freq = TRUE)
   ## a c g t
   ## 4462 3741 1793 4437
   ### 计算GC 含量
   GC.content(woodmouse_sel)
   ## [1] 0.3834269
   
   
2. 计算DNA序列的多态性位点
   同时使用 seg.sites 函数可以得到 woodmouse_sel 序列的多态性位点。
   seg.sites(woodmouse_sel)
   ## [1] 30 33 48 52 60 72 99 106 114 117 123 136 147 164 189 198 207 210 213
   ## [20] 225 261 282 300 309 312 314 316 336 343 345 349 360 364 384 386 405 435 436
   ## [39] 439 444 453 477 489 516 529 531 534 540 546 555 558 576 585 587 588 624 639
   ## [58] 657 666 669 697 738 739 744 747 750 756 759 761 762 768 780 789 858 873 874
   ## [77] 882 886 891 898 909 912 918 919 933 939 945 951 952 954 958
3. 计算DNA序列的dN/dS比值
   使用 dnds 函数可以计算 woodmouse_sel 序列的dN/dS比值，在这里要注意选择正确的遗传密码。
   res <- dnds(woodmouse_sel, code = 2, quiet = TRUE) # 使用正确的遗传密码
   ## Warning in dnds(woodmouse_sel, code = 2, quiet = TRUE): sequence length not a
   ## multiple of 3: 2 nucleotides dropped
   res
   ## AJ511910 AJ511912 AJ511930 AJ511935 AJ511928 AJ511940
   ## AJ511912 0.58275270
   ## AJ511930 0.02134700 0.01041933
   ## AJ511935 0.10969093 0.06522544 0.02795323
   ## AJ511928 0.04328522 0.03097372 0.04956510 0.04088191
   ## AJ511940 0.03083661 0.01917935 0.03455807 0.03723846 0.28891084
   ## AJ511906 0.18705497 0.08037003 0.02179929 0.09491645 0.04212085 0.03003039
   ## AJ511937 0.04809339 0.03434489 0.07878045 0.05446881 0.32376358 0.15412459
   ## AJ511879 0.11536818 0.00000000 0.01092493 0.07552502 0.03255668 0.02014204
   ## AJ511921 0.07444135 0.03482987 0.01935405 0.07971876 0.03909170 0.02788470
   ## AJ511938 0.04414976 0.03156118 0.06985973 0.05012969 0.17625775 0.17476286
   ## AJ511891 0.18714304 0.08055267 0.02284600 0.09534667 0.04656098 0.03308829
   ## AJ511897 0.22980609 0.09568942 0.02113103 0.10194056 0.04406411 0.03053050
   ## AJ511904 0.23111838 0.09623807 0.02087321 0.10239574 0.04118446 0.02935392
   ## AJ511914 0.10148947 0.04657559 0.02181212 0.06201525 0.04425439 0.03152534
   ## AJ511906 AJ511937 AJ511879 AJ511921 AJ511938 AJ511891
   ## AJ511912
   ## AJ511930
   ## AJ511935
   ## AJ511928
   ## AJ511940
   ## AJ511906
   ## AJ511937 0.04671703
   ## AJ511879 0.08023448 0.03627198
   ## AJ511921 0.06501695 0.04306317 0.03715269
   ## AJ511938 0.04298441 0.14658587 0.03321532 0.03869527
   ## AJ511891 0.13809967 0.05212052 0.08068128 0.06953674 0.04742001
   ## AJ511897 0.16008131 0.04765509 0.09558929 0.06925608 0.04360266 0.00000000
   ## AJ511904 0.16073464 0.04556049 0.09614380 0.06967621 0.04190798 0.16162972
   ## AJ511914 0.08519845 0.05065090 0.05103264 0.06922133 0.04518841 0.08539626
   ## AJ511897 AJ511904
   ## AJ511912
   ## AJ511930
   ## AJ511935
   ## AJ511928
   ## AJ511940
   ## AJ511906
   ## AJ511937
   ## AJ511879
   ## AJ511921
   ## AJ511938
   ## AJ511891
   ## AJ511897
   ## AJ511904 0.23067884
   ## AJ511914 0.09308166 0.09352930

四、序列比对
序列比对有多种软件可以选择，这里利用 muscle (Edgar, 2004) 进行多序列联配。用 ape 包的 muscle 函数进行序列比对之前，需要提前安装对应的软件。muscle 软件下载地址如下：http://www.drive5.com/muscle/ 。此外，使用muscle函数时需要指定软件的安装路径，或者将软件路径写入系统的环境变量。例如在windows系统下，muscle软件的路径为D:\Software，则把D:\Software写入系统的环境变量即可。在macOS 和Linux 系统可以直接将muscle软件复制到 /usr/local/bin 目录。在命令行中输入muscle 显示软件的使用帮助页面表示设置成功。
woodmouse_sel_align <- muscle(woodmouse_sel)
woodmouse_sel_align
## 15 DNA sequences in binary format stored in a matrix.
##
## All sequences of same length: 965
##
## Labels:
## AJ511910
## AJ511912
## AJ511930
## AJ511935
## AJ511928
## AJ511940
## ...
##
## Base composition:
## a c g t
## 0.309 0.259 0.124 0.307
## (Total: 14.47 kb)
序列比对完成后可以使用 trans 函数将核苷酸序列翻译为氨基酸序列。trans函数现支持6种不同的遗传密码类型。具体遗传密码的使用可以参考 https://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=cgencodes。 翻译的时候如果序列长度不是3的倍数，会出现警告信息。woodmouse_sel 数据属于脊椎动物线粒体 Cytb 基因，code 的数值应选择 2。
AA <- trans(woodmouse_sel_align, code = 2)
## Warning in trans(woodmouse_sel_align, code = 2): sequence length not a multiple
## of 3: 2 nucleotides dropped
AA
## 15 amino acid sequences in a matrix
## All sequences of the same length: 321
image(AA)


图 1: woodmouse_sel 翻译为氨基酸序列后的可视化结果

五、序列可视化
序列比对完成之后，使用image.DNAbin 可以对DNA序列可视化，使用 image.AAbin 函数可以对氨基酸序列可视化，对于核苷酸和氨基酸序列也可以直接用 image 函数。
### 显示缺失值为红色
image(woodmouse_sel_align, "n", "red")


图 2: 序列 woodmouse_sel_align 的可视化，缺失位点用红色表示

### 显示 G+C 的值为蓝色
image(woodmouse_sel_align, c("g", "c"), "blue")


图 3: 序列 woodmouse_sel_align 的可视化, G+C 用蓝色表示

### 碱基 a 为 黑色
image(woodmouse_sel_align, "a", "black")


图 4: 序列 woodmouse_sel_align 的可视化, 碱基 A 用黑色表示

此外，用 all.equal 函数可以比较两条序列的差异。若选择 plot = FALSE 参数，则直接返回序列差异的比较结果（文字描述）；若选择 plot = TRUE，则会用图形显示序列的差异。
### 定义 woodmouse_sel_align1 数据等于 woodmouse_sel_align 数据
woodmouse_sel_align1 <- woodmouse_sel_align
### 修改第一条序列1-10列的碱基为 “g”
woodmouse_sel_align1[1, c(1:10)] <- as.DNAbin("g")
all.equal(woodmouse_sel_align, woodmouse_sel_align1, plot = TRUE)


图 5: 序列差异的可视化结果

## $messages
## [1] "Labels in both objects identical."
##
## $different.sites
## [,1] [,2]
## [1,] 1 1
## [2,] 1 2
## [3,] 1 3
## [4,] 1 4
## [5,] 1 6
## [6,] 1 7
## [7,] 1 8
## [8,] 1 9
## [9,] 1 10
通过上图可以看出，all.equal 函数可以把序列有差异的列单独提取进行可视化，其他列则不再显示。

六、计算遗传距离
序列比对之后可以使用 dist.dna 函数计算成对比较遗传距离，模型选择会影响 dist.dna 的计算结果。本函数默认选择的模型是 “K80” （Kimura, 1980），此模型在DNA 条形码研究中经常用到。此外，对于存在缺失值的数据，pairwise.deletion 参数也会对结果有所影响。
models <- c("RAW", "JC69", "K80", "F81", "K81", "F84", "T92", "TN93", "GG95")
dist <- lapply(models, function(x) dist.dna(woodmouse_sel_align, model = x))
boxplot(dist, col = rainbow(9))


图 6: 不同模型下的遗传距离箱线图比较

七、构建进化树

1. NJ 法
   序列比对完成后，使用nj 函数可以基于遗传距离构建 NJ (neighbor-joining) 树. 此方法是 Saitou 和 Nei 于 1987年提出的，在DNA条形码研究中经常用到。
   tree <- nj(dist.dna(woodmouse_sel_align))
   plot(tree, edge.width = 1.5)
   
   
   图 7: 展示进化树 tree
   
   
   
   
   
2. bootstrap 分析
   上面步骤计算得到的NJ树如果想要得到节点支持率可以使用 boot.phylo 函数，重抽样次数一般设置为1000。计算完节点支持率后可以使用 drawSupportOnEdges 函数将支持率添加到树的分支上。
   ### 定义 fun 函数
   fun <- function(x) nj(dist.dna(x))
   tree <- fun(woodmouse_sel_align)
   set.seed(12345)
   bp <- boot.phylo(tree, woodmouse_sel_align, fun, quiet = F, B = 1000)
   ##
   Running bootstraps: 100 / 1000
   Running bootstraps: 200 / 1000
   Running bootstraps: 300 / 1000
   Running bootstraps: 400 / 1000
   Running bootstraps: 500 / 1000
   Running bootstraps: 600 / 1000
   Running bootstraps: 700 / 1000
   Running bootstraps: 800 / 1000
   Running bootstraps: 900 / 1000
   Running bootstraps: 1000 / 1000
   ## Calculating bootstrap values... done.
   plot(tree, edge.width = 1.5)
   ### 在分支上添加节点支持率
   drawSupportOnEdges(round(bp/10, 0), cex = 0.8)
   
   
   图 8: 展示进化树 tree，添加节点支持率
   

1. 树的写入
   树的写入可以使用的函数有 write.tree 和 write.nexus。write.tree 函数可以生成 Newick 格式树文件；write.nexus 则可生成 NEXUS 格式树文件。 在此我们可以把上面得到的进化树保存到本地。
   ### Newick 格式
   write.tree(tree, file = "woodmouse_sel.tre")
   ### NEXUS 格式
   write.nexus(tree, file = "woodmouse_sel.nex.tre")
   
   
2. 树的读取
   ape 有两个函数可以读取树文件，read.tree 和 read.nexus。 两个函数可以分别读取 Newick 和 NEXUS格式的树文件。我们可以用这两个函数读取刚才写入的两个树文件。
   tree.nwk <- read.tree("woodmouse_sel.tre")
   tree.nex <- read.nexus("woodmouse_sel.nex.tre")
   plot(tree.nwk, edge.width = 1.5)
   
   
   图 9: 展示进化树 tree.nwk
   
   

1. 置根
   在 ape 包中可以使用 root 函数置根系统发育树。置根需要选择外群，在这里为了便于演示，暂用 “AJ511938” 进行置根。实际操作时，选择指定的外群进行置根即可。
   tree_root <- root(tree, "AJ511938")
   plot(tree_root, edge.width = 1.5)
   
   
   图 10: 展示置根后的进化树 tree_root
   
   
2. 比较两个系统发育树的差异
   对于两个系统发育树，尤其是类群比较多的情况下，直接观察不容易发现两者的区别。使用 comparePhylo 函数可以比较两个系统发育树的差异。若参数选择 plot = TRUE，则通过图形显示两者差异，反之则只使用文字进行描述。
   a <- tree
   ### 去除两个末端节点
   b <- drop.tip(tree, 1:2)
   ### 比较两棵树的差异
   comparePhylo(a, b, plot = TRUE, force.rooted = TRUE)
   
   
   图 11: 两棵树的差异比较
   
   ## => Comparing a with b.
   ## Trees have different numbers of tips: 15 and 13.
   ## Tips in a not in b : AJ511910, AJ511912.
   ## Trees have different numbers of nodes: 13 and 12.
   ## a is unrooted, b is rooted.
   ## Both trees are not ultrametric.
   ## 2 clades in a not in b.
   ## 1 clade in b not in a.
   
   
3. 旋转一个节点的两个分支
   ape 包还有一个实用的函数是 rotate 函数。使用 rotate 函数可以旋转一个节点的两个分支。
   op <- par(no.readonly = TRUE)
   par(mfrow = c(1, 2))
   ### 旋转第24个节点位置的末端标签
   tree_rotate <- rotate(tree, 24)
   plot(ladderize(tree), sub = "A")
   nodelabels("", 24, frame = "c", bg = "red")
   plot(ladderize(tree_rotate), sub = "B")
   nodelabels("", 24, frame = "c", bg = "red")
   par(op)
   
   
   图 12: 比较进化树 tree 和 tree_rotate
   
   通过上图可以看出红色节点位置（node = 24）的末端节点标签发生了互换。
   
   
4. 系统发育树局部放大
   对于类群数目比较多的系统发育树，使用 zoom 函数可以放大系统发育树的指定部分。若参数选择 subtree = FALSE，则不展示其他末端节点，反之则会将其他末端节点依次进行合并，只显示合并后的末端节点数目。
   zoom(tree, 1:3, subtree = FALSE)
   
   
   图 13: 树的局部放大，只显示放大部分
   
   zoom(tree, 1:3, subtree = TRUE)
   
   
   图 14: 树的局部放大，其他部分只显示末端节点数目
   
   zoom(tree, grep("8", tree$tip.label), subtree = FALSE)
   
   
   图 15: 树的局部放大，显示末端标签含有 “8” 的分支
   
   
5. 去除系统发育树的末端节点
   使用 drop.tip 函数可以去除系统发育树的末端节点。这里可以使用末端节点的编号，也可以直接输入末端节点的名称。
   tree1 <- drop.tip(tree, 1:5, subtree = FALSE)
   comparePhylo(tree, tree1, plot = FALSE)
   ## => Comparing tree with tree1.
   ## Trees have different numbers of tips: 15 and 10.
   ## Tips in tree not in tree1 : AJ511910, AJ511912, AJ511930, AJ511935, AJ511928.
   ## Trees have different numbers of nodes: 13 and 9.
   ## tree is unrooted, tree1 is rooted.
   ## Both trees are not ultrametric.
   par(mfrow = c(1, 2))
   plot(tree, main = "tree", sub = "A", edge.width = 1.5)
   plot(tree1, main = "tree1", sub = "B", edge.width = 1.5)
   par(op)
   
   
   图 16: 展示去除特定末端节点后树形的变化
   
   
6. 生成超度量树
   为了推断分歧时间，需要把系统发育树转换为超度量树。在这里可以使用 chronos 函数，此函数可以使用 Penalised Likelihood （Sanderson, 2002; Kim and Sanderson, 2008) 和 Maximum Likelihood（Paradis, 2013）方法生成超度量树。为了便于演示，此处选择NJ树，实际使用时需要输入ML/BI树。
   tree <- nj(dist.dna(woodmouse_sel, pairwise.deletion = TRUE))
   str(tree)
   ## List of 4
   ## $ edge : int [1:27, 1:2] 16 16 19 21 23 24 24 23 25 26 ...
   ## $ edge.length: num [1:27] 0.004118 0.000615 0.002347 0.001033 0.003489 ...
   ## $ tip.label : chr [1:15] "AJ511910" "AJ511912" "AJ511930" "AJ511935" ...
   ## $ Nnode : int 13
   ## - attr(*, "class")= chr "phylo"
   ## - attr(*, "order")= chr "cladewise"
   ### 默认使用 correlated rate 模型
   tree_chronos <- chronos(tree)
   ##
   ## Setting initial dates...
   ## Fitting in progress... get a first set of estimates
   ## (Penalised) log-lik = -0.5692622
   ## Optimising rates... dates... -0.5692622
   ## Warning: false convergence (8)
   ##
   ## log-Lik = -0.56827
   ## PHIIC = 79.14
   tree$edge.length
   ## [1] 4.117731e-03 6.154783e-04 2.346552e-03 1.032984e-03 3.489051e-03
   ## [6] 6.447749e-03 1.939005e-03 2.990860e-02 2.562387e-03 1.265861e-03
   ## [11] 4.295423e-03 4.146896e-03 2.193855e-03 4.619535e-03 3.865063e-03
   ## [16] 4.689343e-03 1.173483e-02 2.515568e-03 7.214804e-05 1.991459e-04
   ## [21] 2.800554e-03 9.306119e-04 1.155711e-03 4.013051e-03 2.152766e-03
   ## [26] 1.038423e-03 1.038423e-03
   plot(tree_chronos, edge.width = 1.5)
   
   
   图 17: 展示超度量树 tree_chronos
   
   
7. 多样性分析
   生成超度量树之后，使用 ltt.plot 函数可以绘制 LTT 图形。LTT 图形可以描述支系随时间变化的趋势。
   ### 绘制LTT图
   ltt.plot(tree_chronos)
   title("Lineages Through Time Plot")
   
   
   图 18: 使用ltt.plot函数绘制的LTT图
   
   ### 将树和 LTT 图形绘制在一起
   par(mfrow = c(1, 2))
   plot(tree_chronos, show.tip.label = TRUE, edge.width = 1.5)
   ltt.plot(tree_chronos)
   par(op)
   
   
   图 19: 将树和 LTT 图共同展示
   
   ### 也可使用 ltt.coplot 函数
   ltt.coplot(tree_chronos, show.tip.label = TRUE, x.lim = 1.1)
   
   
   图 20: 使用ltt.plot函数绘制的LTT图
   
   ### 同时绘制多个 LTT 图形
   par(mfrow = c(1, 1))
   tree_chronos1 <- drop.tip(tree_chronos, 1:2)
   tree_chronos2 <- drop.tip(tree_chronos, 3:4)
   mltt.plot(tree_chronos, tree_chronos1, tree_chronos2)
   
   
   图 21: 同时展示多个 LTT 图形
   
   ### 模拟数据，rcoal 函数可以随机生成超度量树
   tree_rcoalx <- lapply(1:1000, function(x) rcoal(10))
   mltt.plot(tree_rcoalx, legend = FALSE)
   
   
   图 22: 同时展示1000个 LTT 图形
   
   
   

十、树的可视化

1. 改变树的展示方式
   ape 包的 plot.phylo 函数提供了系统发育树多种不同的展示方式，如下图所示：
   par(mfrow = c(3, 2))
   plot(tree, type = "p", main = "phylogram with branch lengths", sub = "A", use.edge.length = TRUE) # 默认
   plot(tree, type = "p", main = "phylogram without branch lengths", sub = "B", use.edge.length = FALSE)
   plot(tree, type = "c", main = "cladogram", sub = "C")
   plot(tree, type = "f", main = "fan", sub = "D")
   plot(tree, type = "u", main = "unrooted", sub = "E")
   plot(tree, type = "r", main = "radial", sub = "F")
   par(op)
   
   
   图 23: 系统发育树六种不同的展示方式
   
   
2. 末端标签对齐
   使用 align.tip.label = TRUE，可以将系统发育树的末端标签进行对齐操作。
   plot(tree, type = "p", align.tip.label = TRUE, edge.width = 1.5)
   
   
   图 24: 系统发育树末端标签对齐
   
   
3. 系统发育树的细节修改
   系统发育树的细节修改包括线型、分支颜色和粗细等，可以根据具体要求进行调节。
   ### 改变分支的线型
   par(mfrow = c(3, 2))
   for (i in 1:6) {
   plot(tree, edge.lty = i, sub = LETTERS[i], cex = 0.8)
   }
   par(op)
   
   
   图 25: 比较系统发育树的六种线型
   
   ### 改变分支的颜色
   plot(tree, edge.color = rainbow(27), edge.width = 1.5)
   
   
   图 26: 改变系统发育树分支的颜色
   
   ### 改变分支的粗细
   plot(tree, edge.width = 1:27/2)
   
   
   图 27: 改变系统发育树分支的粗细
   
   ### 改变末端标签的颜色
   plot(tree, tip.color = rainbow(15))
   
   
   图 28: 改变系统发育树末端标签的颜色
   
   ### 改变字体 (默认斜体)
   ### font 的数值可选择1-4，分别代表常规字体，粗体，斜体，粗斜体。
   par(mfrow = c(2, 2))
   for (i in 1:4) {
   plot(tree, font = i, sub = LETTERS[i], cex = 1.2)
   }
   par(op)
   
   
   图 29: 改变系统发育树末端标签的字体
   
   
4. 增加标签及注释
   nodelabels 函数可以为系统发育树增加内部节点标签。
   plot(tree, edge.width = 1.5)
   nodelabels(bg = "lightgray", frame = "c")
   
   
   图 30: 展示系统发育树的内部节点编号
   
   ### 加入节点标签
   plot(tree, edge.width = 1.5)
   nodelabels("circle", 17, frame = "c", bg = "red", adj = 0.5, cex = 0.6)
   nodelabels("rect", 19, frame = "r", bg = "yellow", font = 1, adj = 0)
   
   
   图 31: 展示系统发育树的特定节点标签
   
   ### 更改字体，为防止重叠可以使用 adj 参数进行调整
   plot(tree, use.edge.length = FALSE, font = 1, edge.width = 1.5)
   nodelabels(1:13, adj = -0.2, frame = "n", font = 1)
   nodelabels(14:26, adj = c(1.2, 1), frame = "n", font = 2)
   nodelabels(27:39, adj = c(1.2, -0.2), frame = "n", font = 3)
   
   
   图 32: 树的每个内部节点同时显示3个数字
   
   ### 节点标签使用饼图
   plot(tree, "c", use.edge.length = FALSE, edge.width = 1.5)
   nodelabels(pie = runif(13), piecol = rainbow(3))
   
   
   图 33: 树的节点标签使用饼图
   
   ### 节点标签使用柱状图
   plot(tree, "c", use.edge.length = FALSE, edge.width = 1.5)
   nodelabels(thermo = runif(13), piecol = rainbow(3))
   
   
   图 34: 树的节点标签使用柱状图
   
   ### 节点标签使用水平柱状图
   plot(tree, "c", use.edge.length = FALSE, edge.width = 1.5)
   nodelabels(thermo = runif(13), piecol = rainbow(3), horiz = TRUE)
   
   
   图 35: 树的节点标签使用水平柱状图
   
   此外，祖先状态重建结果的可视化也需要用到 nodelabels 函数。
   ### 随机生成15个超度量树
   tree_rcoal <- rcoal(15)
   ### rnorm(15) 可以生成符合正态分布的15个随机数 (mean = 0, sd = 1)
   x <- rnorm(15)
   ### 生成离散型变量
   x <- factor(c(rep(0, 5), rep(1, 10)))
   ans <- ace(x, tree_rcoal, type = "d")
   plot(tree_rcoal, type = "p", TRUE, edge.width = 1.5)
   co <- c("red", "blue")
   tiplabels(pch = 22, bg = co[as.numeric(x)], cex = 2, adj = 1)
   nodelabels(thermo = ans$lik.anc, piecol = co, cex = 0.6)
   
   
   图 36: 祖先特征重建结果
   
   ### 使用饼图
   plot(tree_rcoal, edge.width = 1.5)
   ### 定义颜色
   co <- c("red", "blue")
   tiplabels(pch = 21, bg = "yellow", cex = 1.5, adj = 1)
   nodelabels(pie = ans$lik.anc, piecol = co, cex = 0.6)
   
   
   图 37: 祖先特征重建结果，使用饼图
   
   使用 tiplabels 函数可以为系统发育树添加末端节点标签。
   ### 显示末端节点编号
   plot(tree, edge.width = 1.5)
   tiplabels()
   
   
   图 38: 展示系统发育树的末端节点编号
   
   plot(tree, "p", use.edge.length = FALSE, edge.width = 1.5, label.offset = 3, x.lim = 31)
   tiplabels(pch = 18, col = gray.colors(13), cex = 3, adj = 1.4)
   tiplabels(pch = 19, col = rainbow(13), adj = 2.5, cex = 3)
   
   
   图 39:树的末端同时展示两种不同的特征
   
   使用 edgelabels 函数可以在系统发育树的分支上添加标签。
   plot(tree, edge.width = 1.5)
   edgelabels()
   
   
   图 40: 展示树的分支编号
   
   plot(tree, edge.width = 1.5)
   edgelabels(round(tree$edge.length, 3), srt = 0)
   
   
   图 41: 在树的分支中展示支长
   
   plot(tree, "p", FALSE, edge.width = 1.5)
   edgelabels("a", adj = c(0.5, -0.25), bg = "grey")
   edgelabels("b", adj = c(0.5, 1.25), bg = "white")
   
   
   图 42: 在树的分支上下添加文字
   
   
5. 添加坐标轴
   使用 axisPhylo 函数可以为系统发育树添加坐标轴。
   par(mfrow = c(1, 2))
   plot(tree, edge.width = 1.5, sub = "A")
   axisPhylo()
   plot(tree, direction = "u", edge.width = 1.5, sub = "B")
   ### 改变方向
   axisPhylo(2, las = 1)
   par(op)
   
   
   图 43: 添加坐标轴并改变方向
   
   
6. 添加比例尺
   使用 add.scale.bar 函数可以添加比例尺。add.scale.bar 函数也可以定义比例尺的字号、字体和颜色。
   ### 自定义比例尺的字体和颜色
   plot(tree, edge.width = 1.5)
   add.scale.bar(cex = 1, font = 2, col = "red")
   
   
   图 44: 为系统发育树添加比例尺
   
   
7. 组合图形
   cophyloplot 将两个系统发育树绘制在一起，同时将相同的末端标签用线段连接。
   ### 生成两个随机树
   tree1 <- tree
   tree2 <- drop.tip(tree, 1:5, subtree = FALSE)
   
   ### 生成相关矩阵
   association <- cbind(tree2$tip.label, tree2$tip.label)
   
   cophyloplot(tree1, tree2, assoc = association, length.line = 1, space = 30, gap = 2)
   
   
   图 45: 同时展示 tree1和 tree2
   
   另外，通过改变 direction 的参数，可以改变系统发育树的展示方向。比如可以使系统发育树进行左右或上下组合。
   par(mfrow = c(2, 2))
   plot(tree, "p", FALSE, direction = "r", sub = "A", edge.width = 1.5)
   plot(tree, "p", FALSE, direction = "l", sub = "B", edge.width = 1.5)
   plot(tree, "c", FALSE, direction = "u", sub = "C", edge.width = 1.5)
   plot(tree, "c", FALSE, direction = "d", sub = "D", edge.width = 1.5)
   par(op)
   
   
   图 46: 改变树的方向
   
   
8. 将树和图形共同展示
   phydataplot 函数可以将树和图形进行组合，可供选择的类型有7种。
   ### 生成随机树，使展示的变量和树的标签对应
   x <- 1:15
   tree <- nj(dist.dna(woodmouse_sel_align, pairwise.deletion = TRUE))
   tree_chronos <- chronos(tree)
   ##
   ## Setting initial dates...
   ## Fitting in progress... get a first set of estimates
   ## (Penalised) log-lik = -0.5692622
   ## Optimising rates... dates... -0.5692622
   ## Warning: false convergence (8)
   ##
   ## log-Lik = -0.56827
   ## PHIIC = 79.14
   names(x) <- tree_chronos$tip.label
   ### 类型选择 'bars'
   plot(tree_chronos, x.lim = 6, edge.width = 1.5)
   phydataplot(x, tree_chronos, "bars", scaling = 0.25)
   
   
   图 47: 将树和柱状图共同展示
   
   ### 类型选择 “segments”
   plot(tree_chronos, x.lim = 6, edge.width = 1.5)
   phydataplot(x, tree_chronos, "s", lwd = 3, lty = 3, scaling = 0.25)
   
   
   图 48: 将树和线段共同展示
   
   ### 类型选择 “arrows”
   plot(tree_chronos, x.lim = 6, edge.width = 1.5)
   phydataplot(x, tree_chronos, "a", scaling = 0.25)
   
   
   图 49: 将树和箭头共同展示
   
   ### 类型选择 “dotchart”
   plot(tree_chronos, x.lim = 6, edge.width = 1.5)
   phydataplot(x, tree_chronos, "d", scaling = 0.25)
   
   
   图 50: 将树和点图共同展示
   
   ### 同时展示柱状图和点图
   plot(tree_chronos, x.lim = 20, edge.width = 1.5)
   phydataplot(x, tree_chronos, col = "yellow", scaling = 0.5, offset = 3)
   phydataplot(x, tree_chronos, "d", scaling = 0.5, offset = 12)
   
   
   图 51: 同时展示柱状图和点图
   
   ### runif 函数可以生成均匀分布的随机数
   x <- runif(15, 0, 0.5)
   y <- runif(15, 0, 0.5)
   z <- 1 - x - y
   m <- cbind(x, y, z)
   rownames(m) <- tree_chronos$tip.label
   plot(tree_chronos, x.lim = 3, edge.width = 1.5)
   col <- c("red", "green", "blue")
   phydataplot(m, tree_chronos, col = col, offset = 0.5)
   
   
   图 52: 将树和堆积柱状图共同展示
   
   ### 生成矩阵
   m <- matrix(rnorm(15 * 5), 15)
   rownames(m) <- tree_chronos$tip.label
   plot(tree_chronos, x.lim = 3, edge.width = 1.5)
   ### 展示箱线图
   phydataplot(m, tree_chronos, "bo", scaling = 0.25, offset = 0.5, boxfill = c("red",
   "skyblue"))
   
   
   图 53: 将树和箱线图共同展示
   
   ### 在树的旁边显示序列联配结果
   plot(tree_chronos, x.lim = 6, align.tip = FALSE, font = 1, edge.width = 1.5)
   phydataplot(woodmouse_sel_align[, 1:100], tree_chronos, "m", scaling = 1, border = NA)
   
   
   图 54: 将树和序列联配结果共同展示
   
   dist <- dist.dna(woodmouse_sel_align)
   par(mar = rep(1, 4))
   plot(tree_chronos, x.lim = 5, cex = 0.8, font = 1)
   phydataplot(as.matrix(dist), tree_chronos, "i", offset = 0.8)
   par(op)
   
   
   图 55: 将树和热图共同展示
   
   此外，对于系统发育树的可视化，推荐使用 R 包 ggtree（Yu et al., 2017, Methods in Ecology and Evolution）, 其使用方法详见 https://guangchuangyu.github.io/ggtree-book/chapter-ggtree.html。
   

致谢

感谢bio-protocol 的编辑工作以及两位匿名审稿专家提出建设性的修改意见。

参考文献

1. Edgar, R. C. (2004) MUSCLE: Multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res 32: 1792-1797.
2. Kimura, M. (1980). A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J Mol Evol 16: 111-120.
3. Kim, J. and Sanderson, M. J. (2008) Penalized likelihood phylogenetic inference: bridging the parsimony-likelihood gap. Syst Biol 57: 665-674.
4. Michaux, J. R., Magnanou, E., Paradis, E., Nieberding, C. and Libois, R. (2003) Mitochondrial phylogeography of the Woodmouse (Apodemus sylvaticus) in the Western Palearctic region. Mol Ecol 12: 685-697.
5. Paradis, E. and Schliep, K. (2019). ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35(3): 526-528.
6. Paradis, E. (2012) Developing and Implementing Phylogenetic Methods in R. In: Analysis of Phylogenetics and Evolution with R. Use R!. Springer, New York, NY.
7. Paradis, E. (2013) Molecular dating of phylogenies by likelihood methods: a comparison of models and a new information criterion. Mol Phylogenet Evol 67: 436-444.
8. Paradis, E., Claude, J. and Strimmer, K. (2004). APE: Analyses of Phylogenetics and Evolution in R language. Bioinformatics 20(2): 289-290.
9. Saitou, N. and Nei, M. (1987) The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol 4: 406-425.
10. Sanderson, M. J. (2002) Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach. Mol Biol Evol 19: 101-109.
11. Yu, G., Smith, D. K., Zhu, H., Guan, Y. and Lam, T.T.‐Y. (2017), GGTREE: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol Evol. 8(1): 28-36.