16s分析之PCoA分析学习笔记 / 四六文摘

今天我们来一起学习一下PCoA分析：PCoA可以使用很多种距离的相异或者相似矩阵；

基于关联矩阵，所以它也会有表征其重要性的特征根；

######################################################################

##########################################################

#导入我们需要的R包

library("ggplot2")

library("vegan")

#找到我们需要的文件，mapping文件和beta多样性矩阵文件，现在我们先就bray_curtis做一个PCoA

setwd("E:/Shared_Folder/HG-QIIME-OTU")

design =read.table("map_lxdjhg.txt", header=T, row.names= 1,sep="\t")

setwd("E:/Shared_Folder/HG-QIIME-OTU/Beta")

bray_curtis =read.table("bray_curtis_CSS_otu_table12072.txt",sep="\t",row.names= 1,header=T, check.names=F)

#这里使用我们提前弄好的mapping文件和距离矩阵做一个匹配，主要就是两个文件处理编号检查一下，防止出现错误

idx =rownames(design) %in% colnames(bray_curtis)

#输出逻辑值

idx

#匹配上的筛选出来mapping

sub_design =design[idx,]

#初次作图，避免意外，展示一下，内心安稳

sub_design

#匹配上的矩阵筛出来，所以，做聚类矩阵的时候我们无所谓多少个处理，到时候，我们编一个mapping文件就可以了

bray_curtis =bray_curtis[rownames(sub_design), rownames(sub_design)]

#这里请记住pcoa函数，这里我们简单了解一下这个函数，第一个参数表示距离矩阵，第二个k代表选择多少个维度，在我们当然选择二维就够了，eig就是特征值，图片上常见的主成分得分就是使用这个值计算的

pcoa =cmdscale(bray_curtis, k=2, eig=T) # k is dimension, 3 is recommended; eig iseigenvalues

此时PCoA分析就算是做完了

下一步就是出图的了：

#提取我们作图需要的数据

points =as.data.frame(pcoa$points) # 获得坐标点

colnames(points)= c("x", "y") #坐标点命名

eig = pcoa$eig#提取特征值为了计算坐标轴得分

#将分组信息和我们的坐标轴合并

points =cbind(points, sub_design[match(rownames(points), rownames(sub_design)), ])

#我们是在做科研，显著性分析必然少不了

#Adonis又称置换多元方差分析或非参数多元方差分析。它利用半度量（如Bray-Curtis）或度量的距离矩阵（如Euclidean）对总方差进行分解，分析不同分组因素对样品差异的解释度，并使用置换检验对其统计学意义进行显著性分析

adonis(bray_curtis~design$SampleType,permutations = 999,method="bray")

#也可以使用anosim，这里都可以，关于统计，之后我会退出微生态中常用的统计及其运用这里就不详细说了，这里我们选择这种检验

anosim.result<-anosim(bray_curtis,design$SampleType,permutations = 999)

#要通过summary提取出来我们需要的东西

summary(anosim.result)

下一步就是出图了，出图我们使用ggplot包

#还是首先将保存命令加上

tiff(file="beta_bray_curtis.tif",res = 300, compression = "none", width=180,height=140,units="mm")

#开始作图

p =ggplot(points, aes(x=x, y=y, color=SampleType)) +

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="bray_curtis PCoA")

#这里我是通过手动将显著性分析值输出到图片上

p +theme_classic()+stat_ellipse(type = "t", linetype = 2)+

annotate("text",x=0.004,y=-0.15,parse=TRUE,size=4,label="'R:'*0.876",family="serif",fontface="italic",colour="darkred")+

annotate("text",x=0,y=-0.17,parse=TRUE,size=4,label="'p:'*0.001",family="serif",fontface="italic",colour="darkred")

#呼应保存图片命令，结束画板

dev.off()

成图结果展示：

到这里pcoa就大功告成了

我们作图的矩阵文件是如何求得的：

方法一：

#R包vegan十分强大，这里使用vegdist函数x代表标准化后的otu表格

vegdist(x,method="bray")

备注：我们应该使用什么样的otu表格：

当我们聚类得到otu表格直接来做vegdist吗？当然不行，基于序列数的otu表格并不是代表真正的丰度序列信息，一个简单的只能表征相对丰度的表格；如果所有的样品序列总数一样，倒也没事，但是我们测序得到的序列深度往往跨度在几千到几万条之间，

所以我们对otu表格进行标准化，那么我为什么不使用重抽样将测序深度抽成相同的呢，主要就是损失信息太多，所以之前我使用CSS标准化；

方法二：

#使用Qiime

#准备我们的CSS标准化

normalize_table.py-i otu_table.biom -a CSS -o CSS_otu.biom

#用对齐文件成功标准化

#这里给出转化为txt的文件方便使用R语言计算

biom convert -iCSS_otu_table.biom -o CSS_otu_table.txt --table-type="OTU table"--to-tsv

#当然我们可以使用Qiime直接求出

beta_diversity.py-i CSS_otu_table.biom -o ./Beta_d -t rep_set.tree -m bray_curtis,weighted_unifrac,unweighted_unifrac

学习永无止境，分享永不停歇！

16s分析之PCoA分析学习笔记

相关推荐