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

qiime2会直接给出各个分类水平的丰度表，用R语言也可以从otu表中抽提出各个水平的丰度表。

读取otu丰富度

代码思路：

1 strsplit后，phylum不是na也不是""则记下，否则记行数到delete

2 添加phylum列，删除注释列和记到delete的phylum

3 按Phylum排序，获取Phylum unique列表

4提取第一个phylum所有行到新表，apply列求和

5遍历剩下的phylum列表，提取每个phylum，apply按列求和，添加到新表

6 数据归一化，保存

拓展:

R语言读取Excel文件的各种方法

最近初学R语言，在R语言读入EXCEL数据格式文件的问题上遇到了困难，经过在网上搜索解决了这一问题，下面归纳几种方法，供大家分享：

第一： R中读取excel文件中的数据的路径：

假定在您的电脑有一个 excel 文件，原始的文件路径是： D:workdata1

如果直接把这个路径拷贝到R中，就会出现错误，原因是：

是escape character（转义符），\才是真正的字符，或者用/

因此，在R中有两种方法读取该路径：

1：在R中输入一下路径： D:\work\data\1 2：在R中输入一下路径： D:\work\data\1 第二： R中读取excel文件中的数据的方法 ：

read.table(),read.csv(),read.delim()直接读取EXCEl文件时，都会遇到一下问题：“在读取‘.xls’的TableHeader时遇到不完全的最后一行”。解决的方法有以下几种：假如文件1.1中是一个6乘以2的矩阵，元素为：

方法1： xls另存为csv格式然后用read.csv ：

具体过程如下：

>data<-read.csv("D:\work\data\1.csv") >data X1 X23 1 2 24 2 3 25 3 4 26 4 5 27 5 6 28 >data<-read.csv("D:\work\data\1.csv",header = F) >data V1 V2 1 1 23333 2 2 24 3 3 25 4 4 26 5 5 27 6 6 28 >data<-read.csv("D:\work\data\1.csv",header = T) >data X1 X23333 1 2 24 2 3 25 3 4 26 4 5 27 5 6 28 也就是说 header = T（TURE）是默认的状态 ，在这默认状态下，输出的data矩阵是一个5乘以2的矩阵，第一行作为了data的名字，如果 header = F（FALSE）， 则会现实原始的矩阵结果。

方法2： xls另存为txt格式然后用read.table ： 如例子所示：

>data<-read.table("D:\work\data\1.txt",header = T) >data X1 X23 1 2 24 2 3 25 3 4 26 4 5 27 5 6 28

>data<-read.table("D:\work\data\1.txt",header = F) >data V1 V2 1 1 23 2 2 24 3 3 25 4 4 26 5 5 27 6 6 28 方法3：打开EXCEL，全选里面的内容，点击复制，然后在R中输入一下命令：数据分析培训

data <- read.table("clipboard", header = T, sep = 't') 结果如下所示：

>data <- read.table("clipboard", header = T, sep = 't') >data X1 X23 1 2 24 2 3 25 3 4 26 4 5 27 5 6 28 >data <- read.table("clipboard", header = F, sep = 't') >data V1 V2 1 1 23 2 2 24 3 3 25 4 4 26 5 5 27 6 6 28 使用这种方法的时候一定要注意复制！剪切板里面没有内容是无法运行的！

主成分分析(Principal Components Analysis，PCA) ，也称主分量分析或主成分回归分析法，是一种无监督的数据降维方法。PCA通过线性变换将原始数据变换为一组各维度线性无关的表示，可用于提取数据的主要特征分量，常用于高维数据的 降维 。这种降维的思想首先减少数据集的维数，同时还保持数据集的对方差贡献最大的特征，最终使数据直观呈现在二维坐标系。

直观上，第一主成分轴 优于 第二主成分轴，具有最大可分性。

主坐标分析（Principal Coordinates Analysis，PCoA），即经典多维标度（Classical multidimensional scaling），用于研究数据间的相似性。

主成分分析（Principal components analysis，PCA）是一种统计分析、简化数据集的方法。它利用正交变换来对一系列可能相关的变量的观测值进行线性变换，从而投影为一系列线性不相关变量的值，这些不相关变量称为主成分（Principal Components）。具体地，主成分可以看做一个线性方程，其包含一系列线性系数来指示投影方向（如图）。PCA对原始数据的正则化或预处理敏感（相对缩放）。PCA是最简单的以特征量分析多元统计分布的方法。通常情况下，这种运算可以被看作是揭露数据的内部结构，从而更好的解释数据的变量的方法。

主坐标分析（Principal Coordinates Analysis，PCoA），即经典多维标度（Classical multidimensional scaling），用于研究数据间的相似性。PCoA与PCA都是降低数据维度的方法，但是差异在在于PCA是基于原始矩阵，而PCoA是基于通过原始矩阵计算出的距离矩阵。因此，PCA是尽力保留数据中的变异让点的位置不改动，而PCoA是尽力保证原本的距离关系不发生改变，也就是使得原始数据间点的距离与投影中即结果中各点之间的距离尽可能相关（如图）。

R中有很多包都提供了PCA和PCoA,比如常用的ade4包。本文将基于该包进行PCA和PCoA的分析，数据是自带的deug，该数据提供了104个学生9门课程的成绩（见截图）和综合评定。综合评定有以下几个等级：A+,A,B,B-,C-,D。

让我们通过PCA和PCoA来看一看这样的综合评定是否合理，是否确实依据这9门课把这104个学生合理分配到不同组（每个等级一个组）。

前文已经介绍了PCA是基于原始数据，所以直接进行PCA分析即可。相信大家都比较熟悉散点图的绘制方法，这里不再细讲，PCA分析完毕后我们直接作图展示结果。

整体看起来还不错，就是B-和C-的学生似乎难以区分。

有时候PCA和PCoA的结果差不多，有时候某种方法能够把样本有效分开而另一种可能效果不佳，这些都要看样本数据的特性。

除转录组研究以外，在16S微生物的研究中我们会根据物种丰度的文件对数据进行PCA或者PCoA分析，也是我们所说的β多样性分析。根据PCA或者PCoA的结果看感染组和对照组能否分开，以了解微生物组的总体变化情况。

β多样性分析的概念

Beta多样性指的是样本间多样性。在肠道菌群分析中，Beta多样性是衡量个体间微生物组成相似性的一个指标。通过计算样本间距离可以获得β多样性计算矩阵，后续一般会利用PCoA、进化树聚类等分析对此数值关系进行图形展示。主要基于OTU的群落比较方法，有欧式距离、bray curtis距离、Jaccard 距离，这些方法优势在于算法简单，考虑物种丰度（有无）和均度（相对丰度），但其没有考虑OTUs之间的进化关系，认为OTU之间不存在进化上的联系，每个OTU间的关系平等。另一种算法Unifrac距离法，是根据系统发生树进行比较，并根据16s的序列信息对OTU进行进化树分类， 一般有加权和非加权分析。

QIIME2中重要的Beta多样性指数：

Jaccard距离：群落差异的定性度量，即只考虑种类，不考虑丰度。

Bray-Curtis距离：群落差异的定量度量，较常用。

Unweighted UniFrac距离：包含特征之间的系统发育关系的群落差异定性度量。

Weighted UniFrac距离：包含特征之间的系统发育关系的群落差异定量度量。

解压缩通过qiime2输出的 .qza文件，获得绘图的matrix和pcoa结果文件

将pcoa结果整理成下表，保存为 ***_site.txt

注意没有legend，需要AI加入。

后期需要继续摸索，其实可以加legend的，只是目前自己的技术做不到。。。

PCA思想解析：

https://www.jianshu.com/p/09bae5cbdc53

---