lyao222lll的个人博客分享 http://blog.sciencenet.cn/u/lyao222lll

博文

R包vegan的群落PCA及tb-PCA分析

已有 1758 次阅读 2019-11-16 21:04 |个人分类:R统计分析|系统分类:科研笔记| _R语言, _多元统计, _排序分析

Rvegan的群落PCAtb-PCA分析

 

主成分分析(Principal Component AnalysisPCA)是基于特征向量的线性无约束排序方法。在群落分析中,分析对象是原始的定量数据,排序图展示样方之间的欧几里得距离。

PCA的基本概念及其在群落分析中的应用,参考前文

 

本篇以某16S扩增子测序所得细菌群落数据,以及测量的环境因子数据为例,在生态学分析的角度,简介在R中运行PCA的过程。

示例文件、R脚本等,已上传至百度盘:

https://pan.baidu.com/s/1RKOSji5gPHqSey_AlJSwEA

1.png

1)文件“phylum_table.txt”为通过16S测序所得的细菌类群丰度表格,此处统计至细菌“门水平”(即界门纲目科属种的“门”这一水平)。(为了方便演示,我没使用OTU水平的)

2)文件“env_table.txt”记录了多种土壤理化指标信息,即环境数据。

接下来期望通过PCA分析这些数据。

 


环境数据的PCA分析


数据集中存在9种定量环境变量,它们有什么相关性?各样方的环境特征是否存在相似性?

 

PCA执行


vegan包中,可通过rda()执行PCA,该函数的输入格式大致如下。

library(vegan)
 
#详情 ?rda
#带协变量的 RDA
rda(Y, X, W)
#不带协变量的 RDA,即常规 RDA
rda(Y, X)
#PCA
rda(Y)

其中,Y为响应变量矩阵,X为解释变量矩阵,W为偏RDA分析需要输入的协变量矩阵。XW是在RDA分析中才需要提供的,对于PCA而言,只需提供响应变量矩阵Y即可,即由对象(行)和变量(列)组成的矩阵。

因此,对环境数据的PCA分析,方法如下。

#读取环境数据
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
 
##PCA 排序
#环境变量需要标准化,详情 ?rda
pca_env <- rda(env, scale = TRUE)
 
#I 型标尺
summary(pca_env, scaling = 1)
#II 型标尺
#summary(pca_env, scaling = 2)

rda()scale是对数据进行标准化。“scale=F”,不对变量标准化,排序对象是离差矩阵(包含方差和协方差的变量之间的关联矩阵);“scale=T”,对变量标准化,排序对象是相关矩阵(不同量纲的变量之间的相关系数矩阵)。由于各环境变量具有不同的量纲,为了使不同物理单位的变量趋向一致,且降低极端值的影响,因此需要标准化。

本篇末尾处,也举例比较了标准化前后PCA的差异,并简述了环境数据标准化的必要性。

 

PCA结果解读


summary()展示结果时,根据所重要要关注的内容,可选根据I型标尺或II型标尺缩放排序坐标后展示出。如果对排序样方之间的距离更感兴趣,则考虑I型标尺;如果对变量之间的相关关系更感兴趣,则考虑II型标尺。关于I型标尺和II标尺,可参考前文

这里展示了按I型标尺缩放后排序结果(参scaling = 1),部分内容如下所示。

2.png

 

关键内容提取


可以将重要的信息提取出来,例如:

#各主成分(PCA轴)特征根,特征根是固定的,和标尺选择无关
pca_eig <- pca_env$CA$eig
#除以特征根总和,即可得各主成分(PCA轴)解释量
pca_exp <- pca_env$CA$eig / sum(pca_env$CA$eig)
 
#提取对象排序坐标,以 I 型标尺为例,以前两轴为例
#scores() 提取样方坐标
site.scaling1 <- scores(pca_env, choices = 1:2, scaling = 1, display = 'site')
#或者在 summary 中提取
site.scaling1 <- summary(pca_env, scaling = 1)$sites[ ,1:2]
site.scaling1
#write.table(site.scaling1, 'site.scaling1.txt', col.names = NA, sep = '\t', quote = FALSE)
 
#提取变量排序坐标,以 II 型标尺为例,以前两轴为例
#scores() 提取环境变量坐标
env.scaling2 <- scores(pca_env, choices = 1:2, scaling = 2, display = 'sp')
#或者 summary 中提取
env.scaling2 <- summary(pca_env, scaling = 2)$species[ ,1:2]
env.scaling2
#write.table(env.scaling2, 'env.scaling2.txt', col.names = NA, sep = '\t', quote = FALSE)

3.png

 

PCA双序图


最后根据PCA所得的样方排坐标以及环境变量的排序坐标,绘图就可以了。

通常,排序对象直接在对应坐标处绘制为点;变量则呈现为向量,由原点(0,0)起始,指向变量得分的对应坐标处,向量的方向表示了该变量数值增加的方向。虽然也可以将变量在顶点坐标处直接展示为点(而非箭头形式),但考虑到它们可能会与对象坐标点混淆,因此一般不推荐这样表示。

#简洁版双序图,详情 ?biplot
pc1 <- paste('PC1:', round(pca_exp[1]*100, 2), '%')
pc2 <- paste('PC2:',round(pca_exp[2]*100, 2), '%')
 
par(mfrow = c(1, 2))
biplot(pca_env, choices = c(1, 2), scaling = 1, main = 'I 型标尺', col = c('red', 'blue'), xlab = pc1, ylab = pc2)
biplot(pca_env, choices = c(1, 2), scaling = 2, main = 'II 型标尺', col = c('red', 'blue'), xlab = pc1, ylab = pc2)

4.png

env.scaling1 <- summary(pca_env, scaling = 1)$species[ ,1:2]

#ordiplot() 作图,I 型标尺为例,详情 ?ordiplot
png('pca_env.png', width = 600, height = 600, res = 100, units = 'px')

ordiplot(pca_env, type = 'n', choices = c(1, 2), scaling = 1, main = 'I 型标尺', xlab = pc1, ylab = pc2)
points(pca_env, dis = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 4), rep('orange', 4), rep('green3', 4)), col = NA, cex = 1.2)
arrows(0, 0, env.scaling1[ ,1], env.scaling1[ ,2], length = 0.1, lty = 1, col = 'blue')
text(pca_env, dis = 'sp', choices = c(1, 2), scaling = 1, col = 'blue', cex = 0.8)

dev.off()

5.png

#ggplot2 作图,I 型标尺为例
#提取样方和环境变量排序坐标,前两轴
site.scaling1 <- data.frame(summary(pca_env, scaling = 1)$sites[ ,1:2])
env.scaling1 <- data.frame(summary(pca_env, scaling = 1)$species[ ,1:2])
 
#添加分组信息
site.scaling1$sample <- rownames(site.scaling1)
site.scaling1$group <- c(rep('A', 4), rep('B', 4), rep('C', 4))
env.scaling1$group <- rownames(env.scaling1)
 
#ggplot2 作图
library(ggplot2)
 
p <- ggplot(site.scaling1, aes(PC1, PC2)) +
geom_point(aes(color = group)) +
scale_color_manual(values = c('red', 'orange', 'green3')) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.title = element_blank(), legend.key = element_rect(fill = 'transparent')) + 
labs(x = pc1, y = pc2) +
geom_vline(xintercept = 0, color = 'gray', size = 0.5) + 
geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
geom_segment(data = env.scaling1, aes(x = 0, y = 0, xend = PC1, yend = PC2), arrow = arrow(length = unit(0.1, 'cm')), size = 0.3, color = 'blue') +
geom_text(data = env.scaling1, aes(PC1 * 1.1, PC2 * 1.1, label = group), color = 'blue', size = 3)
 
#ggsave('pca_env.pdf', p, width = 4.5, height = 4)
ggsave('pca_env.png', p, width = 4.5, height = 4)

6.png

 


物种组成数据的PCA分析


物种组成数据比较特殊。由于PCA本质上是一种展示对象之间欧几里得距离的线性方法,所以可能不适合原始的物种多度数据分析(特别是含有许多0值的情形,欧几里得距离对双零问题很敏感)。但是,可以事先对物种多度数据执行Hellinger转化,将转化后的物种多度数据用于PCA。原因可参考前文

这种通过首先对原始物种组成数据作转化后,再执行PCA的方法,也称为基于转化的PCATransformation-based Principal Component Analysistb-PCA)。

更推荐的方法是使用主坐标分析(Principal Coordinate AnalysisPCoA)代替PCA,通过计算样方间的定量非对称距离测度(如Bray-curtis距离),实现多维物种空间中的样方排序。本篇不讨论PCoA,继续看tb-PCA

 

我们通过tb-PCA分析物种组成数据,细菌门水平丰度表。

此时,rda()函数中的scale参数(标准化)未再使用“scale=TRUE”。由于物种数据的量纲单位都一致,且Hellinger转化本身即可降低高丰度物种的权重,所以我使用了“scale=FALSE”。

#读取物种数据
species <- read.delim('phylum_table.txt', sep = '\t', row.names = 1, stringsAsFactors = FALSE, check.names = FALSE)
 
#物种数据 Hellinger 预转化,详情 ?decostand
species_hel <- decostand(species, method = 'hellinger')
 
#使用转化后的物种数据执行 PCA,无需标准化
pca_sp <- rda(species_hel, scale = FALSE)
 
#summary(pca_sp, scaling = 1)
#summary(pca_sp, scaling = 2)
 
#这里的排序图就只简单展示样方,不展示物种了
par(mfrow = c(1, 2))
ordiplot(pca_sp, scaling = 1, display = 'site', type = 'text')
ordiplot(pca_sp, scaling = 1, display = 'site', type = 'points')
 
#其它细节可参考上文

7.png

 

其它细节信息,如结果解读、关键信息提取、排序图可视化等,参考上文内容即可,方法一致,就不再演示了。

有一点需要注意,由于物种数据在这里作了转化,所以排序图中的物种相关性并不完全等于原始数据的物种相关性。

 


其它相关内容


以下是一些常见的问题,在此列出来帮助理解PCA

 

评估有解读价值的 PCA


PCA只是一种探索性分析方法,主要目的是在低维空间尽可能多地展示数据的主要趋势特征。最终选择多少个观测轴、选择哪些轴并没有统一的标准,自己看着来。

通常的做法就是选取前2-4个排序轴,然后提取对象或变量在这些轴中的坐标,绘制排序图查看对象或变量沿这些轴的分布是否出现了某种规律。借此评估数据的趋势特征。

排位偏后的轴(通常就是第5轴之后)一般就不考虑了,因为它们的特征值太低。

 

可以通过柱形图由高到低展示所有的特征值,然后根据特征值做个评判,多少个排序轴值得保留和解读。(这也只是一种参考,并不是强制的)

#Kaiser-Guttman 准则
#计算特征值的平均值,可据此保留高于平均值的主成分轴
pca_eig[pca_eig > mean(pca_eig)]
 
#断棍模型
#将单位长度的“棍子”随机分成与 PCA 轴数一样多的几段,由高往低排序
#可据此选取特征根大于所对应断棍长度的轴,或者总和大于所对应断棍长度总和的前几轴
n <- length(pca_eig)
bsm <- data.frame(j=seq(1:n), p = 0)
bsm$p[1] <- 1/n
for (i in 2:n) bsm$p[i] <- bsm$p[i-1] + (1/(n + 1 - i))
bsm$p <- 100*bsm$p/n
bsm
 
#绘制每轴的特征根和方差百分比 
par(mfrow = c(2, 1))
barplot(pca_eig, main = '特征根', col = 'bisque', las = 2)
abline(h = mean(pca_eig), col = 'red')
legend('topright', '平均特征根', lwd = 1, col = 2, bty = 'n')
barplot(t(cbind(100 * pca_eig/sum(pca_eig), bsm$p[n:1])), beside = TRUE, main = '% 变差', col = c('bisque', 2), las = 2)
legend('topright', c('% 特征根', '断棍模型'), pch = 15, col = c('bisque', 2), bty = 'n')

8.png

 

评估高贡献的变量


哪些变量对数据特征具有较高的贡献度?

可以将变量向量垂直投影至两侧的坐标轴上,对于主要的PCA轴,根据变量投影到其上的相对长度,选择高贡献的变量。

 

还有就是通过平衡贡献圈判断,仅可在I型标尺中适用,平衡贡献圈概念可参考前文向量长度超过平衡贡献圈半径的变量,其重要性相对较高。

II型标尺可帮助评估高度共线性、高度相关的变量,据此可以手动对变量去冗余。

#某大佬的函数
source('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/scripts/NumEcolR2/cleanplot.pca.R')
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
pca_env <- rda(env, scale = TRUE)
 
par(mfrow = c(1, 2))
cleanplot.pca(pca_env, scaling = 1)        #I 型标尺中加入平衡贡献圈,评估重要的变量
cleanplot.pca(pca_env, scaling = 2)

9.png

 

关于环境变量的标准化


PCA分析必须从相同量纲的变量矩阵开始,原因是需要将变量总方差分配给特征根。因此,变量必须有相同的物理单位,方差和才有意义(方差的单位是变量单位的平方)。或者,变量是无量纲的数据,例如标准化后的数据。

然后我们不妨比较一下,标准化前后的环境数据,执行PCA后的差异。

#对比环境变量标准化前后的差异
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
 
pca_env1 <- rda(env, scale = FALSE)    #不执行标准化
pca_env2 <- rda(env, scale = TRUE)      #执行标准化
 
#I 型标尺
par(mfrow = c(1, 2))
cleanplot.pca(pca_env1, scaling = 1)
cleanplot.pca(pca_env2, scaling = 1)
 
round(pca_env1$CA$eig/sum(pca_env1$CA$eig)*100, 2)
round(pca_env2$CA$eig/sum(pca_env2$CA$eig)*100, 2)

结果显而易见,标准化前(左)相对于标准化后(右)的环境数据,方差很大,受高权重数值(由于各变量物理单位不同,变量间数值大小差异巨大)的影响甚重,特别是第一主成分轴的特征值极高。所以环境数据标准化是很有必要的。

10.png

11.png

 

关于物种组成数据的Hellinger转化


 一般而言,对于物种组成数据较少使用PCA这种线性模型去做,PCoANMDS等是比较推荐的。

如果仍要使用PCA的话,事实上,并不是说物种多度数据一定要作转化才行。如果物种多度数据比较均匀,且含有很少的0值(这是最关键的,避开了欧几里得距离对“双零”敏感这一问题),则不转化也是可以的。

比方说示例数据“phylum_table.txt”,由于是门水平(一个很大的类别)的类群统计,本身0值较少,其实也可以不执行Hellinger转化的。尽管如此,可能高丰度和低丰度种群的丰度差别还是非常大,数据并不“均匀”,此时可以通过在运行rda()时,通过“scale=TRUE”消除高丰度物种与稀有物种在重要性方面的差异,降低欧几里得距离对高数值的敏感性。


但如果0值过多,使用PCA时还是推荐Hellinger转化,我们换个数据说明下吧。

网盘文件“otu_table.txt”是OTU水平的物种丰度表,其中含有非常多的0值,然后我们比较下Hellinger转化前后的PCA差异。

#对比物种数据 Hellinger 转化前后的差异
otu <- read.delim('otu_table.txt', sep = '\t', row.names = 1, stringsAsFactors = FALSE, check.names = FALSE)
otu_hel <- decostand(otu, method = 'hellinger')
 
pca_sp1 <- rda(otu, scale = FALSE)      #使用 Hellinger 转化前的数据
pca_sp2 <- rda(otu_hel, scale = FALSE)        #使用 Hellinger 转化后的数据
pca_sp3 <- rda(otu, scale = TRUE)        #使用 Hellinger 转化前的数据,但执行标准化
 
#特征值提取
pca_exp1 <- pca_sp1$CA$eig / sum(pca_sp1$CA$eig)
pc1_sp1 <- paste('PC1:', round(pca_exp1[1]*100, 2), '%')
pc2_sp1 <- paste('PC2:',round(pca_exp1[2]*100, 2), '%')
 
pca_exp2 <- pca_sp2$CA$eig / sum(pca_sp2$CA$eig)
pc1_sp2 <- paste('PC1:', round(pca_exp2[1]*100, 2), '%')
pc2_sp2 <- paste('PC2:',round(pca_exp2[2]*100, 2), '%')
 
pca_exp3 <- pca_sp3$CA$eig / sum(pca_sp3$CA$eig)
pc1_sp3 <- paste('PC1:', round(pca_exp3[1]*100, 2), '%')
pc2_sp3 <- paste('PC2:',round(pca_exp3[2]*100, 2), '%')
 
#I 型标尺
par(mfrow = c(2, 2))
#ordiplot(pca_sp1, scaling = 1, display = 'site', type = 'text')
#ordiplot(pca_sp2, scaling = 1, display = 'site', type = 'text')
#ordiplot(pca_sp3, scaling = 1, display = 'site', type = 'text')
ordiplot(pca_sp1, dis = 'site', type = 'n', choices = c(1, 2), scaling = 1, main = 'Hellinger 前,不标准化', xlab = pc1_sp1, ylab = pc2_sp1)
points(pca_sp1, dis = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 12), rep('orange', 12), rep('green3', 12)), col = NA, cex = 1.2)
ordiplot(pca_sp2, dis = 'site', type = 'n', choices = c(1, 2), scaling = 1, main = 'Hellinger 后,不标准化', xlab = pc1_sp2, ylab = pc2_sp2)
points(pca_sp2, dis = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 12), rep('orange', 12), rep('green3', 12)), col = NA, cex = 1.2)
ordiplot(pca_sp3, dis = 'site', type = 'n', choices = c(1, 2), scaling = 1, main = 'Hellinger 前,标准化', xlab = pc1_sp3, ylab = pc2_sp3)
points(pca_sp3, dis = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 12), rep('orange', 12), rep('green3', 12)), col = NA, cex = 1.2)

Hellinger转化前,且未标准化的数据,存在很大的方差、组内离散度高、趋势不明显且略微呈现了“马蹄形效应”;Hellinger转化后或者未Hellinger转化但标准化后的数据,样方分布明显,但相较之下,Hellinger转化后的区分程度更高,主成分轴的解释量具代表性。

12.png

 

PCA图中添加聚类树


还可以对样方执行聚类,将聚类树与PCA排序图相结合展示,也不失为一种可视化好方法。

#排序图中添加聚类树
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
pca_env <- rda(env, scale = TRUE)
 
#I 型标尺
p <- ordiplot(pca_env, dis = 'site', type = 'n', choices = c(1, 2), scaling = 1)
points(pca_env, dis = 'site', choices = c(1, 2), scaling = 1, pch = 21, bg = c(rep('red', 4), rep('orange', 4), rep('green3', 4)), col = NA, cex = 1.2)
 
#例如 UPGMA 聚类
env_upgma <- hclust(dist(scale(env)), method = 'average')
ordicluster(p, env_upgma, col = 'gray')

13.png

 

排序图中的对象点或变量向量可以同比例缩放


有时发现vegan输出的排序图中,所展示的对象或变量坐标,与计算结果中输出的数值对应不起来?但是仔细看的话,尽管横纵坐标的数值不对应,但各对象之间或各变量之间的相对位置并没有变,似乎是同比例缩放了。

在绘制排序图时,有时会考虑将排序对象或变量的得分(即坐标)乘以一个常量,即对排序图中的点或向量坐标同比例放大或缩小一定的数值后展示,以产生易于解读的生态排序结果(可参见LegendreLegendre1998404页)。要么同时缩放所有对象,要么同时缩放所有变量,不能只改变部分对象或变量的位置。这种做法对于解答生态学问题通常是没有影响的。

比方说A图使用原始样方及物种得分(排序坐标)绘图展示,但物种向量的长度过短影响观测。因此考虑将物种得分同比例放大4倍,即对物种向量顶点的原始横纵坐标均乘以4后表示,即得到B图的样式更加易于辨认。

14.png

 

那么,变量向量缩放了,怎么再通过平衡贡献圈评判重要的变量呢?很简单,平衡贡献圈的半径也同比例缩放一下就可以了。

比方说上述使用的函数“cleanplot.pca”(细心的同学可能也注意到了),就考虑到使用一个内部常数对结果重新标定,这时变量向量的箭头长度及平衡贡献圈的半径不完全等于原始值,而是等比例原始值。

 


参考资料


DanielBorcard, FranoisGillet, PierreLegendre, et al. 数量生态学:R语言的应用(赖江山 译). 高等教育出版社, 2014.

Legendre P, Legendre L. Numerical Ecology. Second English edition. Developments in Environmental Modelling, 1998, 20, Elsevier

 



http://blog.sciencenet.cn/blog-3406804-1206413.html

上一篇:Illumina、PacBio和Nanopore的基因组polish
下一篇:对应分析(CA)和去趋势对应分析(DCA)简述

1 李学友

该博文允许注册用户评论 请点击登录 评论 (1 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-2-29 07:54

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部