计算β多样性指数需要用到phyloseq包。它的安装方式不同于简单的install.packages（“phyloseq”）

有两种方法可以安装

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图中两个点的距离的排序，接近β多样性指数的排序

这是最后一篇讲解有关矩阵操作的博客，介绍有关矩阵的函数，主要有 rowSums() , colSums() , rowMeans() , colMeans() , apply() , rbind() , cbind() , row() , col() , rowsum() , aggregate() , sweep() , max.col() 。

下面通过例子来了解这些函数的用法：

我们知道，通过下标索引 [i, j] 可以访问矩阵的某一部分，索引如果没有提供意味着“所有行”或“所有列”。来看个例子，比如：

在R中，可以用一些特殊的函数来进行矩阵的行、列计算。来看些例子：

上述矩阵的行、列计算，还可以使用 apply() 函数来实现。 apply() 函数的原型为 apply(X, MARGIN, FUN, ...) ，其中： X 为矩阵或数组； MARGIN 用来指定是对行运算还是对列运算， MARGIN＝1 表示对行运算， MARGIN＝2 表示对列运算； FUN 用来指定运算函数； ... 用来指定 FUN 中需要的其它参数。来看些例子：

用 apply() 函数来实现上面的例子

apply() 函数功能很强大，我们可以对矩阵的行或列进行其它运算，例如：

如果矩阵存在 NA 值，可通过设置 na.rm=TRUE 来忽略 NA 值，然后再计算。比如：

甚至我们还可以自定义运算函数，来看个例子：

在R中， rbind() 和 cbind() 函数可分别为矩阵添加行和列，来看一个例子：

在R中， row() 和 col() 函数将分别返回元素的行和列下标矩阵，来看个例子：

通过这两个函数，可以获取矩阵的对角元素以及上下三角矩阵，例如：

有时，你可能需要对每行进行分组，然后组内每列求和。在R中可以用 rowsum() 函数来解决，而且效率也非常高。先看个例子：

你也可以用 aggregate() 函数获得类似结果：

有人就会问“为啥没有列分组求和的操作？”，其实你可以先将矩阵转置，然后行分组求和；这两步就等同于列分组求和。

sweep() 函数的原型为 sweep(x, MARGIN, STATS, FUN = "-", check.margin = TRUE, ...) ，其中： x 为矩阵或数组； MARGIN 用来指定是对行运算还是对列运算， MARGIN＝1 表示对行运算， MARGIN＝2 表示对列运算； STATS 表示想要清除的统计量； FUN 用来指定运算函数，默认为减法 - ； check.margin 用来核实 x 的维度是否与 STATS 的匹配，如果事先知道它们匹配的话，将其设为 FALSE 将提高运算速度； ... 用来指定 FUN 中需要的其它参数。来看些例子：

事实上，通过改变 FUN 参数的具体形式或自定义函数， sweep() 函数可以实现很多不同操作，这里就不细讲了。

max.col() 函数返回矩阵每行最大值所在的列位置（即列下标），其原型为 max.col(m, ties.method = c("random", "first", "last")) ，其中： m 为矩阵；当存在多个最大值时， ties.method 指定用哪种方式来处理这种情况，默认为"random"（随机），"first"指使用第一个最大值，"last"指使用最后一个最大值。来看个官网例子：

我们也可以结合 apply() 和 which.max() 函数来实现 max.col(mm, 'first') 。看个例子，

R矩阵的最后一部分内容就讲到这。

如若有遗漏，后期将会添加至本博客。

---