install.packages("ggplot2")

install.packages("ggpubr")

install.packages("vegan")

计算Shannon-香农指数和Simpson-辛普森指数的命令在vegan包中，计算各组显著性的命令在ggpubr包中；画图使用ggplot命令，在行使每个命令之前一定要加载相应的包，如下:

library(ggplot2)

library(ggpubr)

library(vegan)

拿到一个otu表格，要先计算香农指数和辛普森指数，操作如下:

otu=read.table('D:/r-working/feature-table.taxonomy.txt',row.names = 1,skip=1,header=T,comment.char ='',sep='\t')

#读取out表格

#'D:/feature table.taxonomy.txt'为文件路径，注意斜线方向

#row.names = 1指定第一列为行名

#skip=1跳过第一行不读

#header=T指定第一个有效行为列名

#sep='\t'表示指定制表符为分隔符

#comment.char=''表示设置注释符号为空字符‘’，这样#后面的内容就不会被省略

otu=otu[,-ncol(otu)]

#去除表格的最后一列，无用信息

otu=t(otu)

#表格转置，必须将样品名作为行名

shannon=diversity(otu,"shannon")

#计算香农指数,先加载vegan包

shannon

#查看香农指数

simpson=diversity(otu,"simpson")

#计算辛普森指数,先加载vegan包

simpson

#查看辛普森指数

alpha=data.frame(shannon,simpson,check.names=T)

#合并两个指数

write.table(alpha,"D:/r-working/alpha-summary.xls",sep='\t',quote=F)

#存储数据，注意路径使用反斜杠

将各样本进行分组，并进行画图，操作如下:

map<-read.table('D:/r-working/mapping_file.txt',row.names = 1,header=T,comment.char ='',sep='\t',check.names=F)

#读取分组表格

group<-map["Group1"]

#提取需要的分组,'Group1'是表中的分组列名，包括A,B,C三组

alpha<-alpha[match(rownames(group),rownames(alpha)),]

#重排alpha的行的顺序,使其与group的样本id(行名)一致

data<-data.frame(group,alpha,check.rows=T)

#合并两个表格.'<-'与'='同属赋值的含义.

p=ggplot(data=data,aes(x=Group1,y=shannon))+geom_boxplot(fill=rainbow(7)[2])

#data = data指定数据表格

#x=Group1指定作为x轴的数据列名

#y=shannon指定作为y轴的数据列名

#geom_boxplot()表示画箱线图

#fill=rainbow(7)[2]指定填充色

此处用到ggplot2包画箱线图，将画图函数赋值给p后，可以用‘+’不断进行图层叠加，给图片p增加新的特性

p

#查看p

mycompare=list(c('A','B'),c('A','C'),c('B','C'))

#指定多重比较的分组对

mycompare

p<-p+stat_compare_means(comparisons=mycompare,label = "p.signif",method = 'wilcox')

#添加显著性标记的第一种方法，在此之前先加载ggpubr包

p<-p+ylim(2,5.5)

#调整图像的外观

有两种方法可以安装

1.先安装BiocManager

install.packages("BiocManager")

library("BiocManager")

BiocManager::install("phyloseq")

library("phyloseq")

2.source("https://bioconductor.org/biocLite.R")

biocLite("phyloseq")

#安装phyloseq

library("phyloseq")

安装并加载了phyloseq包后，开始读取数据，前面计算α多样性，用到的是read.table……

qiimedata <- import_qiime(otufilename = "feature-table.taxonomy.txt", mapfilename = "mapping_file.txt", treefilename = "tree.rooted.nwk", refseqfilename = "dna-sequences.fasta")

#读取数据，参数都是文件名，注意加后缀

#otufilename指定out表格，mapfilename指定map文件（分组数据)

#treefilename指定有根进化树文件

#refseqfilename指定代表序列文件

otu<-qiimedata@[email protected]

#从qiimedata里面提取otu

sum_of_otus<-colSums(t(otu))

#t_转置,colsums计算列的和,即计算各个otu检测到的总序列数，为了筛掉一些总序列数过低的otu（可能是测序错误）

sum_of_otus

#查看otu总序列数

selected_otu<-names(sum_of_otus)[sum_of_otus>10]

#获取总序列数大于10的otu id

sub_qiimedata <- prune_taxa(selected_otu, qiimedata)

#筛选总序列数大于10的otu的phyloseq数据

weighted_unifrac<-distance(sub_qiimedata,method = 'wunifrac')

#计算样本间加权unifrac

unweighted_unifrac<-distance(sub_qiimedata,method = 'unifrac')

#计算样本间非加权unifrac

bray_curtis <- distance(sub_qiimedata, method='bray')

write.table(as.matrix(bray_curtis),"bray_curtis.txt",sep = '\t',quote = FALSE,col.names = NA)

#保存距离矩阵

#计算样本间Bray-Curtis距离矩阵，method 可选" wunifrac ", " unifrac " ，"jaccard"等

pcoa_of_bray_curtis<-ordinate(physeq=sub_qiimedata,distance = 'bray',method = "PCoA")

#基于Bray-Curtis距离矩阵的PCoA排序分析

p<-plot_ordination(sub_qiimedata, pcoa_of_bray_curtis, type="samples", color="Group1",shape = "Group1")

#将PCoA排序分析结果可视化

library("ggplot2")

p<-p+ scale_colour_manual(values=c("#DC143C","#808000","#00CED1")) + geom_point(size=2) +ggtitle("PCoA of Bray-Curtis distance")+theme(text = element_text(size = 15))

#修改图形大小,ggtitle加标题,stat_ellipse加椭圆

#用scale_colour_manual(values=c())自定义颜色，可查颜色的16进制对照表

p

nmds_of_bray_curtis<-ordinate(physeq=sub_qiimedata,distance = 'bray',method = "NMDS")

#基于Bray-Curtis距离矩阵的NMDS排序分析

p1<-plot_ordination(qiimedata, nmds_of_bray_curtis, type="samples", color="Group1")

#将NMDS排序分析结果可视化

# color=“Group1”指定不同分组的点染不同颜色

p1

p1<-p1+ geom_point(size=3) +ggtitle("NMDS of Bray-Curtis distance") + stat_ellipse()+theme(text = element_text(size = 15))

#对图片进行适当修饰， stat_ellipse()加椭圆， ggtitle()加标题

ggsave(plot = p1,“nmds_of_bary_curtis.pdf",dpi = 300,width

PCoA中的两个点距离，接近β多样性指数

PCA(Principal Components Analysis)即主成分分析，也称主分量分析或主成分回归分析法，首先利用线性变换，将数据变换到一个新的坐标系统中然后再利用降维的思想，使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上，第二大方差在第二个坐标(第二主成分)上。这种降维的思想首先减少数据集的维数，同时还保持数据集的对方差贡献最大的特征，最终使数据直观呈现在二维坐标系。

PCoA(Principal Co-ordinates Analysis)分析即主坐标分析，可呈现研究数据相似性或差异性的可视化坐标，是一种非约束性的数据降维分析方法，可用来研究样本群落组成的相似性或相异性。它与PCA类似，通过一系列的特征值和特征向量进行排序后，选择主要排在前几位的特征值，找到距离矩阵中最主要的坐标，结果是数据矩阵的一个旋转，它没有改变样本点之间的相互位置关系，只是改变了坐标系统。两者的区别为PCA是基于样本的相似系数矩阵(如欧式距离)来寻找主成分，而PCoA是基于距离矩阵(欧式距离以外的其他距离)来寻找主坐标。

NMDS图中两个点的距离的排序，接近β多样性指数的排序