最近好几位读者在公众号后台留言PCoA。今天的推文就来简要的介绍一下PCoA在R语言里实现的代码。具体的原理我自己也没有看明白。目前的理解是：这个PCoA(主坐标分析)和主成分分析（PCA）作用是一样的，都是将多变量的数据进行降维，方便在二维平面进行展示，使我们能够直观的观察到样本之间可能存在的某种关系。PCoA和PCA一样，在R语言里都是一条命令就能够完成分析。但是PCoA有一个不同点是：PCoA的输入数据要求是距离矩阵。所以如果要做PCoA分析第一步是用我们自己的数据计算距离矩阵。
示例数据
今天推文的示例数据使用R语言vegan包里的varespec数据集
首先加载vegan这个包
library(vegan)
然后通过data()函数来获得这个数据
data("varespec")
这个数据集是24行44列，每行是一个样地，每列是样地里的植物种类。数值具体代表什么我还没有看明白，现在的理解就是度量这种植物在这个样地丰富度的一个标准吧。
下面开始实际操作
- 首先是将数据集赋值给新的变量，并以site1-24对新的数据集的行进行命名
df<-varespec
rownames(df)<-paste0("site",1:24)
- 计算距离
bray_dist<-vegdist(df,method = "bray")
- 使用ape这个包中的pcoa()函数做PCoA分析
library(ape)
df.pcoa<-pcoa(bray_dist,correction = "cailliez")
df.pcoa$vectors能够获得用于画图的数据
image.png
df.pcoa$values可以获得坐标轴上显示的百分比
image.png
最后用ggplot2来画这个图
df.plot<-data.frame(df.pcoa$vectors)
head(df.plot)
library(ggplot2)
x_label<-round(df.pcoa$values$Rel_corr_eig[1]*100,2)
y_label<-round(df.pcoa$values$Rel_corr_eig[2]*100,2)
x_label
y_label
ggplot(data=df.plot,aes(x=Axis.1,y=Axis.2))+
geom_point()+
theme_bw()+
theme(panel.grid = element_blank())+
geom_vline(xintercept = 0,lty="dashed")+
geom_hline(yintercept = 0,lty="dashed")+
labs(x=paste0("PCoA1 ",x_label,"%"),
y=paste0("PCoA2 ",y_label,"%"))
image.png
通过上图我们可以看到这些样地大体上可以分为两组，如果自己手头有样地的分组数据就可以看看这个结果是不是和自己的分组数据一致。
下面人为的给他分个组，然后添加一个表示分组的椭圆
df.plot$group<-ifelse(df.plot$Axis.1<0,"AAA","BBB")
ggplot(data=df.plot,aes(x=Axis.1,y=Axis.2,
color=group,shape=group))+
geom_point(size=5)+
theme_bw()+
theme(panel.grid = element_blank())+
geom_vline(xintercept = 0,lty="dashed")+
geom_hline(yintercept = 0,lty="dashed")+
labs(x=paste0("PCoA1 ",x_label,"%"),
y=paste0("PCoA2 ",y_label,"%"))+
stat_ellipse(data=df.plot,
geom = "polygon",
aes(fill=group),
alpha=0.3)+
scale_fill_manual(values = c("#e31a1c","#1f78b4"))
image.png
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！
网友评论