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


学习永无止境,分享永不停歇!

(0)

相关推荐