用Phylomatic和PhyloCom进行群落系统进化分析(20180820)

Phylomatic (http://phylodiversity.net/phylomatic/) 和PhyloCom (http://phylodiversity.net/phylocom/) 是Cam Webb、Steve Kembel和David Ackerly编写的软件，用于群落系统发育分析。

Phylomatic可基于植物名录以及相应进化树骨架生成进化树。不过，通过Phylomatic建立的进化树，一般只能准确反映科或者属以上分类等级的进化关系（取决于该属是否在骨架进化树上是否有信息）。由于不能准确反映科或属以下的系统发育关系，Phylomatic生成的进化树，一般只用于群落系统发育，或者空间尺度较大，类群较多，且系统发育关系较远的系统发育多样性分析中。如果分类单元数比较少，在Genbank中又有足够的序列， 则应该优先考虑自行建立分子系统树，再进行相应分析，具体请参考 （Roquet et al. 2013) 。

通过植物名录建立进化树， 目前也有两个程序包可以选择， 一个是浙江大学金毅老师编写的 S.PhyloMaker R 脚本 (https://github.com/jinyizju/S.PhyloMaker), 另一个是 Scott Chamberlain 编写的 brranching R 程序包 (https://github.com/ropensci/brranching)。 S.PhyloMaker 的数据是基于Zanne 2014进化树生成的，可建立包含几万个分类单元的进化树，但速度较慢。 brranching 主要是基于Phylomatic的，可以直接输入植物名录获得进化树，速度较快，但是在分类单元数量上有一定限制。以上两个程序包的用法，用户可自行参考指南，本文只介绍通过Phylomatic建立进化树。

Phylocom是用来进行群落系统发育分析的软件。该软件可以为Phylomatic软件得到的无枝长进化树，用BLADJ方法拟合出枝长（这是因为Phylomatic软件建立的进化树， 有科属种的标签，而另外一个提供的ages文件， 有科属种的分化时间），不过BLADJ算法在2014年以后已经很少使用。此外，Phylocom还可以：计算系统发育多样性（PD），系统发育结构(community structure)，群落的系统发育距离(community phylogenetic distance)，分析群落的性状进化(AOT)等。下面就介绍如何使用Phylomatic和Phylocom进行相应的分析。

假设有如下植物名录要建立进化树：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
Castanopsis eyrei
Castanopsis carlesii
Castanopsis fargesii
Castanopsis tibetana
Cyclobalanopsis gracilis
Cyclobalanopsis glauca
Lithocarpus glaber
Myrica rubra
Elaeocarpus decipiens
Syzygium buxifolium
Daphniphyllum oldhamii
Loropetalum chinense
Rhododendron latoucheae
Rhododendron ovatum
Vaccinium carlesii
Schima superba
Adinandra millettii
Eurya muricata
Camellia fraterna
Ternstroemia gymnanthera
Machilus thunbergii
Neolitsea aurata var. chekiangensis
Chimonanthus salicifolius

一般流程如下

1. 准备 科/属/种名录

请注意: 《中国植物志》 以及Flora of China等植物志, 由于成书较早，所以采用的分类系统都不是APG。而APG分类系统是目前为止能够较为客观体现科水平系统发育关系的分类系统（最新版本APG IV），因而是推荐采用的。不过目前能获得的准确科属列表为APGIII的版本，与APG IV差别不大。为了批量查询植物所在APG的科名，可使用本人编写的R程序包 plantlist (https://github.com/helixcn/plantlist)。该程序包的简要介绍参见 http://blog.sciencenet.cn/blog-255662-846673.html 。

注意： 由于plantlist给出的是APGIII的科属列表，如果用户要使用APG IV，请自行修订。

1.1 plantlist程序包的安装

1
2
library(devtools)
install_github("helixcn/plantlist")

1.2 生成科属种列表

R 命令为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
sp <- c( "Castanopsis eyrei","Castanopsis carlesii",
"Castanopsis fargesii","Castanopsis tibetana",
"Cyclobalanopsis gracilis","Cyclobalanopsis glauca",
"Lithocarpus glaber","Myrica rubra",
"Elaeocarpus decipiens","Syzygium buxifolium",
"Daphniphyllum oldhamii","Loropetalum chinense",
"Rhododendron latoucheae","Rhododendron ovatum",
"Vaccinium carlesii","Schima superba",
"Adinandra millettii","Eurya muricata",
"Camellia fraterna","Ternstroemia gymnanthera",
"Machilus thunbergii","Neolitsea aurata var. chekiangensis",
"Chimonanthus salicifolius")

library(plantlist)
sp2 <- TPL(sp)
res <- taxa.table(sp2)
writeLines(res, "taxa.table.txt")
getwd()

运行以上R代码，科属种列表就保存在taxa.table.txt 文件中，保存的路径在R中有所显示。

所生成文件taxa.table.txt用Notepad++ (https://notepad-plus-plus.org/download/ ) 记事本打开， 内容如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
Fagaceae/Castanopsis/Castanopsis_eyrei
Fagaceae/Castanopsis/Castanopsis_carlesii
Fagaceae/Castanopsis/Castanopsis_fargesii
Fagaceae/Castanopsis/Castanopsis_tibetana
Fagaceae/Cyclobalanopsis/Cyclobalanopsis_gracilis
Fagaceae/Cyclobalanopsis/Cyclobalanopsis_glauca
Fagaceae/Lithocarpus/Lithocarpus_glaber
Myricaceae/Myrica/Myrica_rubra
Elaeocarpaceae/Elaeocarpus/Elaeocarpus_decipiens
Myrtaceae/Syzygium/Syzygium_buxifolium
Daphniphyllaceae/Daphniphyllum/Daphniphyllum_oldhamii
Hamamelidaceae/Loropetalum/Loropetalum_chinense
Ericaceae/Rhododendron/Rhododendron_latoucheae
Ericaceae/Rhododendron/Rhododendron_ovatum
Ericaceae/Vaccinium/Vaccinium_carlesii
Theaceae/Schima/Schima_superba
Theaceae/Camellia/Camellia_fraterna
Pentaphylacaceae/Adinandra/Adinandra_millettii
Pentaphylacaceae/Eurya/Eurya_muricata
Pentaphylacaceae/Ternstroemia/Ternstroemia_gymnanthera
Lauraceae/Machilus/Machilus_thunbergii
Lauraceae/Neolitsea/Neolitsea_aurata_var._chekiangensis
Calycanthaceae/Chimonanthus/Chimonanthus_salicifolius

注： 查询植物科名，也可以用 taxize程序包（https://ropensci.org/tutorials/taxize_tutorial/）

2 Phylomatic建树

在Zanne 2014进化树发表之前，用Phylomatic建立进化树一般只是按照APGIII的骨架，生成无枝长的newick进化树。然后再在Phylocom软件中，用BLADJ模块，将Wikstroem等人推断的各类群分化时间拟合到无枝长的进化树上。Wikstroem的分化时间由于发表年代久远， 被子植物各类群的分化时间有一些变动，通过Phylocom BLADJ的方法校正分化时间已经显得过时。从2015年开始, Phylomatic网站整合了 Zanne 2014的进化树骨架, 可以直接基于该进化树直接获得有枝长的进化树。

方法如下:

  1. 打开 http://phylodiversity.net/phylomatic/
  2. 将 科/属/种 列表 拷贝到 taxa = 框中
  3. storedtree = 选择 Zanne 2014
  4. 勾选 clean = TRUE， 其余选择默认.
  5. 点击Send按钮获得结果

将网页中的输出结果拷贝到文本文件中， 并另存为 phylo文件， 请注意，如果后续用Phylocom进行分析，则该phylo文件不能有任何扩展名。该进化树为Newick格式， 不含Singletons，因此可以用R的ape读取并进行数据处理，或者Figtree软件绘图。

获得的进化树如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
((Chimonanthus_salicifolius:108.8,(Neolitsea_aurata_var._chekiangensis:88.1284,
Machilus_thunbergii:88.1284)Lauraceae:20.672)Laurales:71.8279,
(((((((Cyclobalanopsis_gracilis:20.205,Cyclobalanopsis_glauca:20.205)
Cyclobalanopsis:20.205,(Lithocarpus_glaber:3.30372,
(Castanopsis_eyrei:0.94544,(Castanopsis_carlesii:0.320126,
(Castanopsis_tibetana:0.28527,Castanopsis_fargesii:0.28527):0.034856):
0.625313):2.35828):37.1063)Fagaceae:21.8954,Myrica_rubra:62.3054):
50.9532,Elaeocarpus_decipiens:113.259)Fabidae:4.27613,Syzygium_buxifolium:
117.535)Rosidae:0.138911,(Daphniphyllum_oldhamii:107.877,Loropetalum_chinense:
107.877):9.79685)Superrosidae:1.49186,((Rhododendron_latoucheae:62.0517,
Vaccinium_carlesii:62.0517,Rhododendron_ovatum:62.0517)Ericaceae:20.7165,
((Ternstroemia_gymnanthera:43.8733,(Adinandra_millettii:33.7061,
Eurya_muricata:33.7061):10.1672):27.9519,(Schima_superba:57.076,
Camellia_fraterna:57.076):14.7492):10.943):36.3974)
Gunneridae:61.4628):259.665;

注：用brranching R程序包可以直接输入物种名，获得进化树。

3. 用Phylocom计算PD, NRI 和NTI （以Example文件为例）

注意：这部分计算使用的是Example文件，与上文的植物名录已经没有关系。

3.1 下载phylocom软件

打开 https://github.com/phylocom/phylocom/releases 下载zip文件，解压缩, 可见以下几个文件夹：

  1. Example_data 中是练习数据
  2. Mac 中是苹果机的运行软件
  3. R，是驱动Phylocom的R脚本
  4. Src 是Phylocom的源程序，phylocom是C语言写的
  5. W32 是Windows平台下可以运行的exe程序
  6. Phylocom_manual.pdf 是软件说明书
  7. README 注意事项

3.2 创建工作路径

Phylocom必须在纯ASCII码的路径下运行，运行路径不能有中文。在C盘根目录下，创建一个名为phylocom的文件夹。将w32中的phylocom.exe文件，拷贝到文件夹 C:/phylocom下。 将example_data文件夹中的phylo和sample文件，拷贝到该文件夹C:/phylocom中。

3.3 利用phylocom软件计算PD、NRI、NTI等指数

sample文件第一列是样方的名称，第二列是个体数，第三列是物种名。物种名一定要与phylo文件中的物种名完全对应（包括大小写）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
clump1    1    sp1
clump1    1    sp2
clump1    1    sp3
clump1    1    sp4
clump1    1    sp5
clump1    1    sp6
clump1    1    sp7
clump1    1    sp8
clump2a    1    sp1
clump2a    1    sp2
.... (中间省略)
random    2    sp15
random    3    sp17
random    1    sp22
random    2    sp24

将phylo和sample文件夹拷贝到 C:\phylocom 之后，

点击 开始> cmd > 输入 cd C:\phylocom>

输入 phylocom comstruct > comstruct.output 计算群落的系统发育结构的NRI和NTI指数

输入 phylocom pd > pd.output.txt 计算系统发育多样性 Faith’s PD

输入 phylocom comdist > comdist.output.txt 计算群落两两之间系统发育距离 MPD的MNDT等

注意： 进行系统发育多样性相关的计算， 可以参考picante 程序包 (http://blog.sciencenet.cn/home.php?mod=space&uid=255662&do=blog&id=1116212)

参考文献

  1. Chamberlain, S. (2018). brranching: Fetch ‘Phylogenies’ from Many Sources. R package version 0.3.0. https://CRAN.R-project.org/package=brranching
  2. Qian, H. and Y. Jin. (2016) An updated megaphylogeny of plants, a tool for generating plant phylogenies and an analysis of phylogenetic community structure. Journal of Plant Ecology 9(2): 233–239.
  3. Roquet, C., Thuiller, W., & Lavergne, S. (2013). Building megaphylogenies for macroecology: taking up the challenge. Ecography, 36(1), 13-26.
  4. Webb, C. O., & Donoghue, M. J. (2005). Phylomatic: tree assembly for applied phylogenetics. Molecular Ecology Notes, 5(1), 181-183.
  5. Zhang, J.L. (2018). plantlist: Looking Up the Status of Plant Scientific Names based on The Plant List Database. R package version 0.3.7/r31. https://R-Forge.R-project.org/projects/plantlist/