R语言进行PCoA分析_vlambda技术博客

vlambda
2020-05-21

R语言进行PCoA分析

相信大家在做微生物多样性研究时经常听到PCA分析、PCoA分析，NMDS分析，CCA分析，RDA分析。它们对物种（或基因、功能）的分析具有重要作用，因而频频出现在16S测序及宏基因组测序中。那么你知道这些分析之前到底有什么区别吗？在什么情况下应该用什么分析呢？首先，以上分析本质上都属于排序分析（Ordination analysis）。排序(ordination)的过程就是在一个可视化的低维空间(通常是二维)重新排列这些样方,使得样方之间的距离最大程度地反映出平面散点图内样方之间的关系信息。常用的排序方法如下：

1：只使用物种组成数据的排序称作非限制性排序(unconstrained ordination)

(1)主成分分析(principal components analysis,PCA)

(2)对应分析(correspondence analysis, CA)

(3)去趋势对应分析(Detrended correspondence analysis, DCA)

(3)主坐标分析(principal coordinate analysis, PCoA)

(4)非度量多维尺度分析(non-metric multi-dimensional scaling, NMDS)

2：同时使用物种和环境因子组成数据的排序叫作限制性排序(constrained ordination)

(1)冗余分析(redundancy analysis,RDA)

(2)典范对应分析(canonical correspondence analysis, CCA)

我们先来仔细看看PCA与PCoA分析：

在非限制性排序中，16S和宏基因组数据分析通常用到的是PCA分析和PCoA分析。两者的区别在于：PCA分析是基于原始的物种组成矩阵所做的排序分析，而PCoA分析则是基于由物种组成计算得到的距离矩阵得出的。在PCoA分析中，计算距离矩阵的方法有很多种，一般基于Bray-Curtis使用比较广泛。

另外，我们还知道无论是PCA还是PCoA，一般都需要降维处理(一般物种数目都超过3个，样本数目都超过4个)，而降维就会产生数据损失。多数情况下，我们在做降维处理的时候，期望维数越低越好，这样我们就可以最大程度地保真原始数据。比如一维和二维数据一般不需要降维处理(直接呈现);再比如要想把三维坐标系的数据降维到一维坐标系上，我们首先要把三维空间的数据降维到二维空间上(此时损失一部分)，再将二维空间的数据降维到一维空间上(再损失一部分)。

那这样就好办了。如果样本数目比较多，而物种数目比较少，那肯定首选PCA;如果样本数目比较少，而物种数目比较多，那肯定首选PCoA。

PCoA主坐标分析

　　注：PC1 和PC2 是两个主坐标成分，PC1 表示尽可能最大解释数据变化的主坐标成分，PC2 为解释余下的变化度中占比例最大的主坐标成分，PC3 等依次类推。

如今计算能力如此发达，做PCA和PCoA基本都是分分钟的事情，不妨2个都可进行分析，这里不再赘述。(那之前说的这些还有什么用?古人云：知其然，知其所以然)

下面重点介绍一下PCoA

PCoA（principal co-ordinates analysis）是一种研究数据相似性或差异性的可视化方法，通过一系列的特征值和特征向量进行排序后，选择主要排在前几位的特征值，PCoA 可以找到距离矩阵中最主要的坐标，结果是数据矩阵的一个旋转，它没有改变样品点之间的相互位置关系，只是改变了坐标系统。通过PCoA 可以观察个体或群体间的差异。分析软件： R 语言PCoA 分析和作PCoA 图。

1.在Excel准备好原始数据，组别为行，检测变量为列，注意一定要遵循此规则

R语言进行PCoA分析

2.同样我们需要准备好原始数据的分组信息：

R语言进行PCoA分析

3.读取数据，先读取原始数据，分组信息暂时不读入：

>library(openxlsx)

>mydata<- read.xlsx("pcoa_genera.xlsx",sheet = 1)

>str(mydata)#因为后面我们会用到我们数据集的行与列的数据，因此我们通常利用str（）函数来观察数据集的信息。

4.依次加载以下三个R包

>library(permute)

>library(lattice)

> library(vegan)# vegan 包是进行群落数据分析最常用的R包

5.计算距离

>distance<-vegdist(mydata[1:14,2:21],method='bray') # 注意我们的数据是从第一行第二列开始，因此我们输入的数据为(mydata[1:14,2:21]，另外我们计算距离的方法为Bray-Curtis，特此说明。

>distance

6.保存结果

>write.csv(as.matrix(distance),'distance.csv')

7.随后我们对导出的结果再次导入到控制台

>mydata<- read.xlsx("Bray-curtis_HBR.xlsx",sheet = 1)

或者不需要导入，直接进行如下操作更加便捷

8.加载ape包

>library(ape)

>res=pcoa(distance, correction="none")

>head(res$value)

> head(res$value)

R语言进行PCoA分析

>biplot(res)

R语言进行PCoA分析

初步的图我们就做出来了，但是特别的粗糙，不美观，因此我们需要利用ggplot2进行优化。我们还可以利用 vegan 中的 ANOSIM 相似性分析是一种非参数检验，用来检验组间（两组或多组）差异是否显著大于组内差异，从而判断分组是否有意义。接下来我们导入分组。

>group<- read.xlsx("pca_genera_name.xlsx",sheet = 1)

R语言进行PCoA分析

>dune.ano <- with(mydata[1:14,2:21], anosim(distance, group$site))

>summary(dune.ano)#distance是我们的计算的矩阵，group$site是我们的分组名字，mydata是计算矩阵的出去表头的数据。

R语言进行PCoA分析

>plot(dune.ano)

>library(ggplot2)

>P=dune.ano$signif

>R=round(dune.ano$statistic,3)

>data=res$vectors[,1:2]

>data=data.frame(data)

>colnames(data)=c('x','y')

>eig=as.numeric(res$value[,1])

>type=group$site #我们分组信息

>data=cbind(data,type)

>p =ggplot(data, aes(x=x, y=y, color=type))+geom_point(alpha=.7, size=2) +labs(x=paste("PCoA 1 (",format(100* eig[1] / sum(eig), digits=4), "%)", sep=""),y=paste("PCoA 2 (", format(100*eig[2] / sum(eig), digits=4), "%)", sep=""),title="PCoA")

>p

经过ggplot2美化后的PCoA图就出来了，不得不感叹，ggplot2是真的强大，在R语言的学习中有专门针对ggplot2包的学习资料，笔者就购买了一本有关ggplot2的书，在这里一并推荐给大家：R数据可视化分析

后面如果需要调整图的颜色，大小，字体，可以利用AI实现。

END

vlambda博客
学习文章列表