这一篇是衔接上一篇的，就是要用ggplot2程序包对PCA和PCoA进行可视化。代码我直接照搬过来了，只是绘图的时候用ggplot函数。ggplot2包实现了一个在R中基于全面一致的语法创建图形时的系统。这提供了在R中画图时经常缺乏的图形创造的一致性并允许我们创建具有创新性和新颖性的图标类型，ggplot2中，图是采用‘+’号串联的，每个函数修改属于自己的部分。

#PCA分析

#读入物种数据

phylum

#行为物种，列为样方(原是此形式无需转置)

phylum

##加载 vegan 包，如果尚未安装vegan包，需要在此步骤之前执行：install.packages('vegan')

library(vegan)

#物种数据 Hellinger 预转化(处理包含很多 0 值的群落物种数据时，推荐使用)

phylum_hel

#PCA 排序(vegan 包中 rda() 执行)

tbpca

以上完成了pca排序，需要把数据整理成ggplot可以识别的形式，所以需要一下一些操作。需要自己用excel写一个分组分件。像这样，我把每个样本都分到了对应的分组里。

#提取样方和环境因子排序坐标，前两轴，I 型标尺,还记得上篇的scaling=1型和scaling=2型吗，这里又出现

#了summary(tbpca,scaling=1)为了把这几步都干了啥解释清楚，必须展示一下#summary(tbpca,scaling=1)，执行之后会出现一段很长的展示，有这几项重点关注哈，一个是species #scores 一个是sites scores。显然我们是是想比较不同样品的差异，所以我们选者了sites scores，下面这句##代码的意思是将的tbpca在1型标尺下的site score的前两列创建为数据框并赋值给tbpca(tbpca$表示取#tbpca的某一部分，如果我写tbpca$species那么我取出来的就是他的species scores了，后面的[]是说取行和#列的，逗号规定行，逗号后定义列，举个例子把[1,2]就是说1行2列取出来就是一个数字，也可以说是一行一##列的矩阵，可以拿自己的数据试一下，看见了吧就是这个效果—>

#PCA 前两轴解释量(%)

tbpca_eig

#提取样方和环境因子排序坐标，前两轴，I 型标尺

tbpca

#读取样本分组数据"group1.csv"，

group

#合并样本分组信息，构建 ggplot2 作图数据集

tbpca$sample

tbpca

species scores

site scores

#ggplot2作图 ggplot2中图是采用串联起来(+)号函数创建的，aes()函数的功能是指定每个变量扮演的角##色(说人话就是变量和图像上的点、颜色、形状等的对应关系)。在这里本写作aes(x=PC1,y=PC2)略写为##下面的形式，ggplot()函数使用一个或多个几何函数向图中添加几何对象(简写geom)，ggplot()函数指##定要绘制的数据源的变量，几何函数则指定这些变量如何在视觉上进行表示(使用点、条、线和阴影区)，##scale_color_manual()函数来设定四组样本的点的颜色。theme(主题)，theme()函数中的选项可以让我们调整字体，背景、颜色和网格线等。

library(ggplot2)

ggplot(tbpca_site.scaling1, aes(PC1,PC2,color = group)) +

geom_point() +

scale_color_manual(values = c('red', 'orange', 'green3','gold')) +

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 = paste('PCA1: ', tbpca_eig[1], '%'), y = paste('PCA2: ', tbpca_eig[2], '%')) +

geom_vline(xintercept = 0, color = 'gray', size = 0.5) +

geom_hline(yintercept = 0, color = 'gray', size = 0.5)

然后你就得到了这样一个图了

下面是PCoA部分

#自己数据PCoA(ggplot2)

##加载 vegan 包

library(vegan)

#物种数据 Hellinger 预转化(处理包含很多 0 值的群落物种数据时，推荐使用)

phylum

#使用Bray-curtis方法计算距离矩阵

bray

#PCoA 排序(vegan 包中 cmdscale() 执行)

pcoa

ggplot麻烦的地方要来了，要把你的数据整理成ggplot函数可以识别的形式。这部分就前文的PCA的部分相似了。

#读取样本分组数据

group

pcoa_eig

pcoa_eig$sample

pcoa_eig

pcoa_eig1

#接着是作图，PCA部分也都细说过各部分了。

ggplot(pcoa_eig, aes(Dim1, Dim2)) +

geom_point(aes(color = group)) +

scale_color_manual(values = c('red', 'orange', 'green3','gold'))+

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 = paste('PCoA1: ', pcoa_eig1[1], '%'), y = paste('PCoA2: ',  pcoa_eig1[2], '%')) +

geom_vline(xintercept = 0, color = 'gray', size = 0.5) +

geom_hline(yintercept = 0, color = 'gray', size = 0.5)

然后你就得到了这样一个图了

好的，就这样吧~

转载本文请联系原作者获取授权，同时请注明本文来自靳泽星科学网博客。

链接地址：http://wap.sciencenet.cn/blog-3419243-1231602.html

上一篇：R 语言 PCA PCoA 区别

下一篇：生信 python 练习题 把每条FASTA序列分割成80个字母一行的序列