|
在R中分析物种和系统发育多样性 加拿大 Steve Kembel 原著 张金龙 编译
本课程中， 我们将学习如何分析阿尔伯塔草原的一组植物群落数据。数据中包括阿尔伯塔几个地点草原植物的个体数以及功能性状和系统发育关系等， 数据详细介绍参见以下论文: S.W. Kembel and J.F. Cahill, Jr. 2011. Independent evolution of leaf and root traits within and among temperate grassland plant communities. PLoS ONE 6(6): e19992. (doi:10.1371/journal.pone.0019992). 怎样使用课程材料本课程覆盖了导入数据，在R中分析生物多样性的全部内容。如果要学习全部课程， 可以从R入门部分开始。如果已经有R的基础， 则可以直接学习本讲义。本课程需要用到picante程序包，以及数据映像（即.Rdata），后者包括练习用的全部数据，有了这些数据，运行后面的命令就可直接得到结果。 生物多样性数据导入R首先， 查看一下我们需要的程序包是否能加载。需要用的程序包有ape, picante, vegan. 由于picante是依赖另外两个程序包的，所以加载picante的时候，另外两个程序包也会自动加载。 设定工作路径，请注意Windows和MacOS以及Linux操作系统上，路径斜杠的方向 setwd("C:/Users/jlzhang/Dropbox/04 to review/kembel")
library(picante)## Loading required package: ape ## Loading required package: vegan ## Loading required package: permute ## Loading required package: lattice ## This is vegan 2.5-2 ## Loading required package: nlme 群落数据群落数据包括每个地点，或者每块样地物种的个体数，或者相对多度等。在本课程的数据中，多度数据是每个不同生境20m*20m样方中植物的相对盖度。 群落数据的格式为data.frame, 行表示样地，列表示物种。 练习数据已经包括comm数据了， 所以直接用就行了。 如果数据是保存在csv文件中， 则需要用read.csv函数读取，格式如下： # 每一行代表一块样地，因此行名就是样地编号，第一行作为各列名称，所以这里设定header = TRUE。
comm <- read.csv("grassland.community.csv", header = TRUE, row.names = 1)用以上方法将数据读取到R中，注意行名和列名并不是数据本身，而是各行列的名称。行列有了名称， 在后续就更容易操作和展示。群落数据的各行列有了名称，才方便和其他类型的数据关联。 # 先查看comm的数据类型 class(comm) ## [1] "data.frame" # 查看该data.frame是由多少行列组成的。(rows x columns) dim(comm) ## [1] 27 76 # 查看行名 rownames(comm) ## [1] "mix-O-1" "mix-O-2" "mix-O-3" "mix-O-4" "mix-O-5" "mix-O-6" ## [7] "mix-O-7" "fes-K-8" "fes-K-9" "fes-K-10" "fes-K-11" "fes-K-12" ## [13] "fes-K-13" "fes-K-14" "fes-K-15" "fes-K-16" "fes-K-17" "mix-H-18" ## [19] "mix-H-19" "mix-H-20" "mix-H-21" "mix-H-22" "mix-H-23" "mix-H-24" ## [25] "mix-H-25" "mix-H-26" "mix-H-27" # 查看前六列的名 head(colnames(comm)) ## [1] "Antennaria_parvifolia" ## [2] "Artemisia_cana" ## [3] "Artemisia_frigida" ## [4] "Symphyotrichum_ericoides_var._ericoides" ## [5] "Bouteloua_gracilis" ## [6] "Carex_filifolia" # 提取数据的子集，1到5行，1到5列数据 comm[1:5, 1:5] ## Antennaria_parvifolia Artemisia_cana Artemisia_frigida ## mix-O-1 10 10 50 ## mix-O-2 0 10 50 ## mix-O-3 20 20 30 ## mix-O-4 0 0 0 ## mix-O-5 0 10 0 ## Symphyotrichum_ericoides_var._ericoides Bouteloua_gracilis ## mix-O-1 10 70 ## mix-O-2 10 90 ## mix-O-3 10 60 ## mix-O-4 0 90 ## mix-O-5 0 100 本数据中， 多度其实是样地中某种植物的盖度。多元统计中的很多方法对于多度很敏感，因此可以将绝对多度转换为相对多度。用vegan中的函数即可完成。 # 查看每个样方中盖度之合， 如果数据为个体数， 则表示每个样方中所有种的总个体数 apply(comm, 1, sum) ## mix-O-1 mix-O-2 mix-O-3 mix-O-4 mix-O-5 mix-O-6 mix-O-7 fes-K-8 ## 640 630 710 350 400 650 560 960 ## fes-K-9 fes-K-10 fes-K-11 fes-K-12 fes-K-13 fes-K-14 fes-K-15 fes-K-16 ## 960 980 830 980 980 830 640 1080 ## fes-K-17 mix-H-18 mix-H-19 mix-H-20 mix-H-21 mix-H-22 mix-H-23 mix-H-24 ## 710 440 590 540 340 420 400 600 ## mix-H-25 mix-H-26 mix-H-27 ## 540 590 420 # 除以每个样方的盖度和， 将盖度转换为相对盖度 comm <- decostand(comm, method = "total") # 每个样方的总盖度 apply(comm, 1, sum) ## mix-O-1 mix-O-2 mix-O-3 mix-O-4 mix-O-5 mix-O-6 mix-O-7 fes-K-8 ## 1 1 1 1 1 1 1 1 ## fes-K-9 fes-K-10 fes-K-11 fes-K-12 fes-K-13 fes-K-14 fes-K-15 fes-K-16 ## 1 1 1 1 1 1 1 1 ## fes-K-17 mix-H-18 mix-H-19 mix-H-20 mix-H-21 mix-H-22 mix-H-23 mix-H-24 ## 1 1 1 1 1 1 1 1 ## mix-H-25 mix-H-26 mix-H-27 ## 1 1 1 # 查看转换后的数据 comm[1:5, 1:5] ## Antennaria_parvifolia Artemisia_cana Artemisia_frigida ## mix-O-1 0.01562500 0.01562500 0.07812500 ## mix-O-2 0.00000000 0.01587302 0.07936508 ## mix-O-3 0.02816901 0.02816901 0.04225352 ## mix-O-4 0.00000000 0.00000000 0.00000000 ## mix-O-5 0.00000000 0.02500000 0.00000000 ## Symphyotrichum_ericoides_var._ericoides Bouteloua_gracilis ## mix-O-1 0.01562500 0.10937500 ## mix-O-2 0.01587302 0.14285714 ## mix-O-3 0.01408451 0.08450704 ## mix-O-4 0.00000000 0.25714286 ## mix-O-5 0.00000000 0.25000000 性状数据trait包括每个种叶和根的数据。读取到R中的方法与群落数据一样，但是在性状数据中，每个种是一行，每个性状是一列。 traits <- read.csv("species.traits.csv", header = TRUE, row.names = 1)
# 查看性状数据的前6列
head(traits)## SLA LeafArea LeafThickness SLV ## Achillea_millefolium 140.2663 9.275390 0.4163333 59.56525 ## Allium_textile 137.7006 2.445361 0.9147222 125.69496 ## Amelanchier_alnifolia 156.1014 14.064856 0.2900000 45.45227 ## Androsace_occidentalis 257.2050 0.274745 0.2535000 84.22189 ## Antennaria_neglecta 171.0442 1.731990 0.2810000 48.14442 ## Antennaria_parvifolia 193.8718 0.317200 0.2466667 47.64151 ## LeafTissueDens SRL SRV RootTissueDens ## Achillea_millefolium 0.018085428 74.14570 5.038776 0.2553510 ## Allium_textile 0.008137136 187.85485 14.013757 0.1049986 ## Amelanchier_alnifolia 0.022841232 20.87560 2.518939 0.5039683 ## Androsace_occidentalis 0.017706402 207.45582 3.291592 0.4071078 ## Antennaria_neglecta 0.020920233 124.73397 6.710526 0.1593750 ## Antennaria_parvifolia 0.021048340 44.93859 4.003997 0.2504167 ## RootDiam ## Achillea_millefolium 0.3123600 ## Allium_textile 0.3107833 ## Amelanchier_alnifolia 0.3760667 ## Androsace_occidentalis 0.1148714 ## Antennaria_neglecta 0.2749500 ## Antennaria_parvifolia 0.3495500 # 绘制两两相关图 pairs(traits) # 由于一些变量显著偏向一侧，因此进行对数转换 traits <- log10(traits) # 查看转换后的数据 pairs(traits)
## 样地的总体情况Metadata 主要包括样地的生境，采集地点，环境数据，如坡度、干湿度等 metadata <- read.csv("plot.metadata.csv", header = TRUE, row.names = 1)
# 查看前6行
head(metadata)## habitat site slope aspect slope.position rel.moisture ## mix-O-1 Mixedgrass Onefour 0 270 3.0 1 ## mix-O-2 Mixedgrass Onefour 20 130 1.5 2 ## mix-O-3 Mixedgrass Onefour 5 90 1.0 2 ## mix-O-4 Mixedgrass Onefour 5 40 2.0 1 ## mix-O-5 Mixedgrass Onefour 5 130 2.0 1 ## mix-O-6 Mixedgrass Onefour 1 90 3.0 1 进化树（或称为系统树）newick格式的进化树，用read.tree读取. Nexus格式的进化树，用read.nexus读取。 phy <- read.tree("grassland.phylogeny.newick")
class(phy)## [1] "phylo" phy ## ## Phylogenetic tree with 76 tips and 68 internal nodes. ## ## Tip labels: ## Antennaria_neglecta, Antennaria_parvifolia, Erigeron_glabellus, ## Erigeron_pumilus, Heterotheca_villosa, Symphyotrichum_falcatum_var._falcatum, ... ## Node labels: ## , , , , , , ... ## ## Rooted; includes branch lengths. 读取到R中的进化树为phylo为格式。 phylo格式的详细说明参见ape的主页(http://ape.mpl./)。 phylo对象的本质是list，包括进化树末端分类单元名称（tip label），枝长（edge length）等等， ape内置的相关函数可以针对phylo这种数据格式打印， 汇总， 绘图等。 # 列出phy对象的所有名称， 这些名称是list内部对象的名 names(phy) ## [1] "edge" "edge.length" "Nnode" "node.label" "tip.label" ## [6] "root.edge" # 查看末端分类单元名称， 取前5个 phy$tip.label[1:5] ## [1] "Antennaria_neglecta" "Antennaria_parvifolia" "Erigeron_glabellus" ## [4] "Erigeron_pumilus" "Heterotheca_villosa" # 进化树有多少末端分类单元？ Ntip(phy) ## [1] 76 # 绘制进化树，用cex参数调整物种名字体大小 plot(phy, cex = 0.5)
## 数据整理（cleaning）及匹配 本分析用到的数据包括群落、性状、进化树以及metadata。 ls() ## [1] "comm" "metadata" "phy" "traits" 示例中的数据，因为已经经过了认真整理，样地和样品中的物种名能够完全对应。但是很多情况下，物种名并不能完全匹配。例如， 样地数据可能只有一部分种出现在进化树中，有些种可能有性状数据，但未出现在进化树中。一些分析中， R假设物种在群落数据和进化树中出现的顺序是一致的。但是有可能出现输入错误等。 以上问题都要找出来。picante程序包中有几个函数来检测不同数据中物种名是否相同。 # 检查物种是否匹配 combined <- match.phylo.comm(phy, comm) # 返回结果是一个包含 phy和com的list， # 未在两个数据集中都出现的物种都已经去除， # 并且顺序都排好了。 # 这里用combined中的数据，替换原始数据。 phy <- combined$phy comm <- combined$comm # 性状数据也需要进行同样的处理 combined <- match.phylo.data(phy, traits) # 用结果中检查过的数据，替换掉原来的数据。 # 由于之前的返回结果是list， 这里用$调取phy和data phy <- combined$phy traits <- combined$data #再检查一下metadata与群落数据出现的顺序是否一致？ all.equal(rownames(comm), rownames(metadata)) ## [1] TRUE # 若顺序不同， metadata数据中样方出现的顺序就应该按照群落顺序排序metadata <- metadata[rownames(comm), ]# 这样，数据预处理就完成了。 后面即可开始各种生物多样性分析。 生物多样性数据的可视化及汇总物种丰富度和多样性每种生境上的物种数是否相同？ # 比较 fescue 和 mixedgrass 生境中物种数的差异 boxplot(specnumber(comm) ~ metadata$habitat, ylab = "# of species") # 用t检验比较不同生境中物种数是否存在差异 t.test(specnumber(comm) ~ metadata$habitat) ## ## Welch Two Sample t-test ## ## data: specnumber(comm) by metadata$habitat ## t = 5.1371, df = 17.179, p-value = 7.972e-05 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 6.274287 15.008066 ## sample estimates: ## mean in group Fescue mean in group Mixedgrass ## 22.70000 12.05882 # 物种累计曲线，用于查看样方总面积是否足够大 plot(specaccum(comm), xlab = "# of samples", ylab = "# of species")
# 群落多元分析 不同样地植物群落物种组成是怎样变化的？生境类型和环境变量与群落植物组成有什么关系？ 要研究这些问题， 需要用到多元排序方法。 相关函数在vegan程序包中，相关的函数的帮助文件和使用指南也非常详尽。 也可以参考 Borcard等人编写的 Numerical Ecology in R， 这本书在2018年出的第二版。 等级聚类通过聚类，物种组成相近的群落能聚合到一起。 先计算样方之间的Bray-Curtis距离，Bray-Curtis不但考虑了物种组成的差异， 同时考虑了每个种的多度。之后，用聚合等级聚类（agglomerative hierachical clustering algorithm）算法即可实现聚类。 # 计算样地之间的 Bray-Curtis距离 comm.bc.dist <- vegdist(comm, method = "bray") # 用UPGMA方法聚类 comm.bc.clust <- hclust(comm.bc.dist, method = "average") # 绘制聚类图 plot(comm.bc.clust, ylab = "Bray-Curtis dissimilarity")
两种类型的草地， 由于有不同的物种组成，分别聚成两类 排序R中的排序方法比较多，这里我们采用NMDS （non-metric multidimensional scaling）方法。 NMDS的方法是基于距离矩阵的。 # metaMDS函数会自动进行数据转换，并检查结果是否已经收敛。 comm.bc.mds <- metaMDS(comm, dist = "bray") ## Run 0 stress 0.07174232 ## Run 1 stress 0.07577799 ## Run 2 stress 0.07481014 ## Run 3 stress 0.08008957 ## Run 4 stress 0.08432338 ## Run 5 stress 0.07886623 ## Run 6 stress 0.07886637 ## Run 7 stress 0.0794153 ## Run 8 stress 0.07344627 ## Run 9 stress 0.07344482 ## Run 10 stress 0.08455081 ## Run 11 stress 0.07886629 ## Run 12 stress 0.08390732 ## Run 13 stress 0.07886603 ## Run 14 stress 0.07935565 ## Run 15 stress 0.08447472 ## Run 16 stress 0.07941531 ## Run 17 stress 0.07480942 ## Run 18 stress 0.07909592 ## Run 19 stress 0.08269313 ## Run 20 stress 0.07481005 ## *** No convergence -- monoMDS stopping criteria: ## 20: stress ratio > sratmax # 绘制stressplot，检查排序结果 stressplot(comm.bc.mds)
排序结果可通过多种方式展示 # 在Biplot中显示样地 ordiplot(comm.bc.mds, display = "sites", type = "text") # 添加物种名 ordipointlabel(comm.bc.mds) # ordiplot可以只打开绘图设备和画布，不绘制任何内容 mds.fig <- ordiplot(comm.bc.mds, type = "none") # 只添加样地，颜色为样地类型，pch这里设定的是点的形状。试试 # plot(1:30, pch = 1:30) points(mds.fig, "sites", pch = 19, col = "green", select = metadata$habitat == "Fescue") points(mds.fig, "sites", pch = 19, col = "blue", select = metadata$habitat == "Mixedgrass") # 添加生境置信椭圆 ordiellipse(comm.bc.mds, metadata$habitat, conf = 0.95, label = TRUE) # 将聚类结果在NMDS排序图中显示 ordicluster(comm.bc.mds, comm.bc.clust, col = "gray") 也可以用ordisurf函数显示物种多度 # 显示球葵属Sphaeralcea 的多度. cex参数可以调整点的大小 ordisurf(comm.bc.mds, comm[, "Sphaeralcea_coccinea"], bubble = TRUE, main = "Sphaeralcea coccinea abundance", cex = 3) ## ## Family: gaussian ## Link function: identity ## ## Formula: ## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE) ## ## Estimated degrees of freedom: ## 4.32 total = 5.32 ## ## REML score: -43.63335 添加环境信息和性状数据到排序中环境变量与排序轴之间是怎样的关系？这里使用envfit函数拟合 ordiplot(comm.bc.mds) plot(envfit(comm.bc.mds, metadata[, 3:6]))
用vegan中的cca或者rda等排序方法也是可以的。cca或者rda直接考虑了物种组成和多度， 并不用先计算样方之间的距离。 性状与进化树系统发育信号最近，系统发育保守性受到很多人关注。群落系统发育结构分析就假设物种的性状是系统发育保守的。 系统发育信号是度量不同物种的性状在进化树上相似程度的指数。Blomberg K统计量就是比较性状的观察值与布朗运动模型预测值的统计量 （Blomberg et al. 2003）。K 值接近1表明性状进化过程接近布朗运动， 表明有一定程度的系统发育信号或者呈现一定的保守性。K接近于0表明性状进化倾向于随机，K>1 表明性状保守。 K的显著性检验可以用以下方法： 将进化树上的物种名随机打乱，计算每次打乱物种名后，性状的系统发育独立差（Phylogenetic independent contrast）的方差, 在多次打乱物种名生成零分布后，真实进化树的系统发育独立差方差与之进行比较。相关检验在Kcalc，phylosignal和multiPhylosignal 函数中可实现。 下面用练习数据计算系统发育信号 # 可以用Kcalc依次计算traits数据框的所有列 apply(traits, 2, Kcalc, phy) ## SLA LeafArea LeafThickness SLV LeafTissueDens ## 0.2563107 0.4231067 0.2418525 0.3310724 0.3298808 ## SRL SRV RootTissueDens RootDiam ## 0.2290093 0.2698531 0.2544904 0.3150881 # multiPhylosignal函数可以检验多列的P值，但是该函数需要的进化树为严格的二分叉树 multiPhylosignal(traits, multi2di(phy)) ## K PIC.variance.obs PIC.variance.rnd.mean ## SLA 0.2563107 0.0006068179 0.0008428445 ## LeafArea 0.4231067 0.0055507994 0.0133247586 ## LeafThickness 0.2418525 0.0006669402 0.0008816348 ## SLV 0.3310724 0.0005728618 0.0010752657 ## LeafTissueDens 0.3298808 0.0006688721 0.0012296789 ## SRL 0.2290093 0.0028408725 0.0036853452 ## SRV 0.2698531 0.0011064783 0.0016710753 ## RootTissueDens 0.2544904 0.0010609020 0.0015356573 ## RootDiam 0.3150881 0.0005510416 0.0009813974 ## PIC.variance.P PIC.variance.Z ## SLA 0.037 -1.569983 ## LeafArea 0.001 -4.068208 ## LeafThickness 0.081 -1.324220 ## SLV 0.001 -3.031784 ## LeafTissueDens 0.001 -2.983913 ## SRL 0.065 -1.365218 ## SRV 0.007 -2.291677 ## RootTissueDens 0.012 -2.064882 ## RootDiam 0.001 -2.990749 结果中包括K和PIC.variance.P, 后者表示性状的非随机程度。大部分变量都完全随机表现出更强的系统发育信号。 性状进化的可视化plot(phy, direction = "up", show.tip.label = FALSE, show.node.label = TRUE, cex = 0.7) tiplabels(pch = 19, col = "black", cex = 3 * (traits[, "LeafArea"]/max(traits[,"LeafArea"])))
## 性状相关性的系统发育分析 若功能性状表现出较强的系统发育信号，就违背了数据完全独立的假设。此时可以在gls中设定corBrownian 关联矩阵，考虑物种之间的系统发育关系。（注：在分析过程中， 可以使用caper的pgls回归。） 其中gls是一般最小二乘(Generalised least squares)法的简称，方法类似ANOVA 或线性模型。无序类别变量以及连续变量都可以通过这种方式检验。 下面检验检测一下， 考虑物种之间系统发育关系时， specific root length (SRL) 和 root tissue density 之间的关系。 # 普通的gls模型， 不考虑系统发育关系 root.gls <- gls(RootTissueDens ~ SRL, data = traits) anova(root.gls) ## Denom. DF: 74 ## numDF F-value p-value ## (Intercept) 1 611.0773 <.0001 ## SRL 1 3.0592 0.0844 # 考虑系统发育关系 root.pgls <- gls(RootTissueDens ~ SRL, correlation = corBrownian(value = 1, phy), data = traits) anova(root.pgls) ## Denom. DF: 74 ## numDF F-value p-value ## (Intercept) 1 13.42163 5e-04 ## SRL 1 20.07297 <.0001 # 绘图
plot(RootTissueDens ~ SRL, data = traits, xlab = "SRL (specific root length)",
ylab = "Root tissue density")
abline(coef(root.gls), lwd = 2, col = "black")
abline(coef(root.pgls), lwd = 2, col = "red")
legend("bottomleft", legend = c("GLS fit", "Phylogenetic GLS fit"), lwd = 2,
col = c("black", "red"))
当不考虑系统发育关系时，SRL和root tissue density的关系很弱。考虑系统发育信号后，两者的相关性就明显了。 系统发育与性状多样性系统发育多样性# 计算Faith's PDcomm.pd <- pd(comm, phy) head(comm.pd) ## PD SR ## mix-O-1 1072.3697 16 ## mix-O-2 1475.4767 22 ## mix-O-3 1406.1708 21 ## mix-O-4 564.5899 6 ## mix-O-5 783.4028 10 ## mix-O-6 1028.5796 13 # 不同生境的 Faith's PDboxplot(comm.pd$PD ~ metadata$habitat, xlab = "Habitat", ylab = "Faith's PD") # 不同生境PD是否相同 t.test(comm.pd$PD ~ metadata$habitat) ## ## Welch Two Sample t-test ## ## data: comm.pd$PD by metadata$habitat ## t = 5.7161, df = 17.627, p-value = 2.195e-05 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 451.1886 976.8551 ## sample estimates: ## mean in group Fescue mean in group Mixedgrass ## 1602.2371 888.2153 # 查看PD和物种丰富度的关系 plot(comm.pd$PD ~ comm.pd$SR, xlab = "Species richness", ylab = "Faith's PD")
picante可以用来计算 MPD, MNTD, SESMPD and SESMNTD 具体定义， 请参见Phylocom的说明书。 群落系统发育最重要指数之一就是标准效应值。 公式如下
一般包括： SESMPD 等于 -NRI， SESMNTD 等于 -NTI picante中提供了不同的零模型，用来进行随机化，以计算SES。 # 系统发育距离矩阵 phy.dist <- cophenetic(phy) # ses.mpd comm.sesmpd <- ses.mpd(comm, phy.dist, null.model = "richness", abundance.weighted = FALSE, runs = 999) head(comm.sesmpd) ## ntaxa mpd.obs mpd.rand.mean mpd.rand.sd mpd.obs.rank mpd.obs.z ## mix-O-1 16 231.3054 237.7399 11.979528 252 -0.53712739 ## mix-O-2 22 239.5479 237.8670 9.079322 519 0.18513371 ## mix-O-3 21 236.5260 237.5431 9.939899 390 -0.10231908 ## mix-O-4 6 222.5255 239.8113 23.919368 178 -0.72267054 ## mix-O-5 10 234.2013 238.9068 16.809032 315 -0.27993465 ## mix-O-6 13 239.4120 239.3863 13.018818 431 0.00197437 ## mpd.obs.p runs ## mix-O-1 0.252 999 ## mix-O-2 0.519 999 ## mix-O-3 0.390 999 ## mix-O-4 0.178 999 ## mix-O-5 0.315 999 ## mix-O-6 0.431 999 # ses.mntd comm.sesmntd <- ses.mntd(comm, phy.dist, null.model = "richness", abundance.weighted = FALSE, runs = 999) head(comm.sesmntd) ## ntaxa mntd.obs mntd.rand.mean mntd.rand.sd mntd.obs.rank ## mix-O-1 16 94.98812 103.95233 18.80944 322 ## mix-O-2 22 97.41972 94.29428 14.93484 593 ## mix-O-3 21 98.71519 95.13883 15.00211 591 ## mix-O-4 6 136.86094 155.07310 43.84286 347 ## mix-O-5 10 107.36711 123.61058 30.60802 310 ## mix-O-6 13 118.94915 111.77332 23.51786 634 ## mntd.obs.z mntd.obs.p runs ## mix-O-1 -0.4765805 0.322 999 ## mix-O-2 0.2092718 0.593 999 ## mix-O-3 0.2383905 0.591 999 ## mix-O-4 -0.4153961 0.347 999 ## mix-O-5 -0.5306932 0.310 999 ## mix-O-6 0.3051223 0.634 999 输出结果
SES.mpd>0，且mpd.obs.p>0.95时，表示群落系统发育均匀。 SES.mpd<0，且mpd.obs.p<0.05时，表明与零模型相比系统发育聚集。 一般认为， MPD对对系统发育整体的聚集性或均匀性较为敏感，MNTD对靠近进化树末端的均匀性和聚集性更为敏感。 # 不同生境的 ses.mpdplot(comm.sesmpd$mpd.obs.z ~ metadata$habitat, xlab = "Habitat", ylab = "SES(MPD)") abline(h = 0, col = "gray")
t.test(comm.sesmpd$mpd.obs.z ~ metadata$habitat) ## ## Welch Two Sample t-test ## ## data: comm.sesmpd$mpd.obs.z by metadata$habitat ## t = 5.9818, df = 20.313, p-value = 7.067e-06 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 1.202306 2.487850 ## sample estimates: ## mean in group Fescue mean in group Mixedgrass ## 0.8486423 -0.9964360 # 不同生境的 ses.mntdplot(comm.sesmntd$mntd.obs.z ~ metadata$habitat, xlab = "Habitat", ylab = "SES(MNTD)") abline(h = 0, col = "gray")
t.test(comm.sesmntd$mntd.obs.z ~ metadata$habitat) ## ## Welch Two Sample t-test ## ## data: comm.sesmntd$mntd.obs.z by metadata$habitat ## t = 2.7692, df = 14.489, p-value = 0.01469 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 0.2216608 1.7234088 ## sample estimates: ## mean in group Fescue mean in group Mixedgrass ## 0.5978721 -0.3746627 # 某fescue群落中的物种 plot(phy, show.tip.label = FALSE, main = "Fescue community fes-K-11") tiplabels(tip = which(phy$tip.label %in% colnames(comm)[comm["fes-K-11", ] > 0]), pch = 19)
‘mix-H-23’ Mixedgrass 群落中包含了一系列系统发育聚集的种。 # 某mixedgrass群落中出现的个体 plot(phy, show.tip.label = FALSE, main = "Fescue community mix-H-23") tiplabels(tip = which(phy$tip.label %in% colnames(comm)[comm["mix-H-23", ] > 0]), pch = 19)
## 性状多样性 性状的分析方法与进化树分析的方法类似，一般都是先进行标准化，再计算欧式距离，再计算功能性状相似性，计算MPD或者MNTD # 性状先用scale标准化， 再用dist计算欧式距离计算 trait.dist <- as.matrix(dist(scale(traits), method = "euclidean")) comm.sesmpd.traits <- ses.mpd(comm, trait.dist, null.model = "richness", abundance.weighted = FALSE, runs = 999) plot(comm.sesmpd.traits$mpd.obs.z ~ metadata$habitat, xlab = "Habitat", ylab = "Trait SES(MPD)") abline(h = 0, col = "gray")
vegan中的treedive函数与picante中的pd类似。 系统发育beta多样性unifrac函数和phylosor函数，都是计算样方之间的Faith’s PD。 comdist和comdistnt的计算类似于MPD和MNTD，但是是针对样方之间的距离。计算过程中可以考虑多度，由于之前我们计算的Bray-Curtis距离已经考虑了多度。为了方便比较，后续的计算也考虑多度， 即abundance.weighted = TRUE # 系统发育beta多样性 comm.mntd.dist <- comdistnt(comm, phy.dist, abundance.weighted = TRUE) # 功能性状beta多样性 comm.mntd.traits.dist <- comdistnt(comm, trait.dist, abundance.weighted = TRUE) # 计算两个距离矩阵的相关性，用mantel检验 # 此处计算Bray-Curtis距离和mntd距离的相关性 mantel(comm.bc.dist, comm.mntd.dist) ## ## Mantel statistic based on Pearson's product-moment correlation ## ## Call: ## mantel(xdis = comm.bc.dist, ydis = comm.mntd.dist) ## ## Mantel statistic r: 0.8597 ## Significance: 0.001 ## ## Upper quantiles of permutations (null model): ## 90% 95% 97.5% 99% ## 0.0759 0.1084 0.1525 0.1822 ## Permutation: free ## Number of permutations: 999 # 此处计算Bray-Curtis距离和功能性状距离的相关性 mantel(comm.bc.dist, comm.mntd.traits.dist) ## ## Mantel statistic based on Pearson's product-moment correlation ## ## Call: ## mantel(xdis = comm.bc.dist, ydis = comm.mntd.traits.dist) ## ## Mantel statistic r: 0.9524 ## Significance: 0.001 ## ## Upper quantiles of permutations (null model): ## 90% 95% 97.5% 99% ## 0.0793 0.1206 0.1726 0.2266 ## Permutation: free ## Number of permutations: 999 系统发育和功能性状排序之前，我们用Bray-Curtis进行了NDMS排序，实际上，系统发育MNTD以及功能性状MNTD都可以用类似的方法排序， 以展示样地之间的关系 # 由于只有样地之间的距离， 这里用monoMDS. 具体请认真阅读 monoMDS的说明书。 # 系统发育距离 comm.mntd.mds <- monoMDS(comm.mntd.dist) mds.fig <- ordiplot(comm.mntd.mds, type = "none") ## species scores not available points(mds.fig, "sites", pch = 19, col = "green", select = metadata$habitat == "Fescue") points(mds.fig, "sites", pch = 19, col = "blue", select = metadata$habitat == "Mixedgrass") ordiellipse(comm.mntd.mds, metadata$habitat, conf = 0.95, label = TRUE)
# 功能性状 comm.mntd.traits.mds <- monoMDS(comm.mntd.traits.dist) mds.fig <- ordiplot(comm.mntd.traits.mds, type = "none") ## species scores not available points(mds.fig, "sites", pch = 19, col = "green", select = metadata$habitat == "Fescue") points(mds.fig, "sites", pch = 19, col = "blue", select = metadata$habitat == "Mixedgrass") ordiellipse(comm.mntd.traits.mds, metadata$habitat, conf = 0.95, label = TRUE)
无论是对物种beta多样性， 系统发育beta多样性还是功能性状beta多样性进行排序方式，两种生境类型的差别都比较大。 ##解释beta多样性 为了解释哪些因子对样方之间的系统发育关系有影响， 可以用permutational MANOVA （adonis）计算，R代码如下。输出结果类似ANOVA方差分析。 以下分别检验生境对于Bray-Curtis距离、系统发育距离、功能性状距离的解释能力。 # Bray-Curtis距离 adonis(comm.bc.dist ~ habitat, data = metadata) ## ## Call: ## adonis(formula = comm.bc.dist ~ habitat, data = metadata) ## ## Permutation: free ## Number of permutations: 999 ## ## Terms added sequentially (first to last) ## ## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) ## habitat 1 3.3656 3.3656 25.124 0.50123 0.001 *** ## Residuals 25 3.3490 0.1340 0.49877 ## Total 26 6.7146 1.00000 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # 系统发育多样性距离 adonis(comm.mntd.dist ~ habitat, data = metadata) ## ## Call: ## adonis(formula = comm.mntd.dist ~ habitat, data = metadata) ## ## Permutation: free ## Number of permutations: 999 ## ## Terms added sequentially (first to last) ## ## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) ## habitat 1 22814 22813.7 34.904 0.58267 0.001 *** ## Residuals 25 16340 653.6 0.41733 ## Total 26 39154 1.00000 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # 功能性状距离 adonis(comm.mntd.traits.dist ~ habitat, data = metadata) ## ## Call: ## adonis(formula = comm.mntd.traits.dist ~ habitat, data = metadata) ## ## Permutation: free ## Number of permutations: 999 ## ## Terms added sequentially (first to last) ## ## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F) ## habitat 1 10.9154 10.9154 63.959 0.71897 0.001 *** ## Residuals 25 4.2666 0.1707 0.28103 ## Total 26 15.1820 1.00000 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 |
|
|
