【R语言编程】---利用三代测序绘制菌群聚类热图与物种丰度图

Python017

【R语言编程】---利用三代测序绘制菌群聚类热图与物种丰度图,第1张

前言： 仍然是三代测序数据的分析，宏基因组的文章中经常出现聚类热图和物种丰度图，用来直观地识别与某些疾病或者表型相关的菌群构成。

1.读取数据

一共有11个样本，每一个样本的测序reads都经过Nanopore官方的Epi2Me程序鉴定了物种，下表中第一列是被鉴定的菌种，第二列是该样本中每个物种产生的reads数目。

首先导入到R语言中，合并所有样本到一个数据框：

2.绘制热图

经过上一步，我们得到了列名为样本，行名为菌种的reads数据框，然后就可以绘制热图，进行聚类分析了：

绘制结果：

3.绘制物种丰度图

丰度图，其实就是堆积图，把每个样本的reads数目转换为百分数，然后作图就可以了：

绘制结果：

基因测序分析微生物菌群结构 NA是什么意思

微生物群落测序是指对微生物群体进行高通量测序，通过分析测序序列的构成分析特定环境中微生物群体的构成情况或基因的组成以及功能。借助不同环境下微生物群落的构成差异分析我们可以分析微生物与环境因素或宿主之间的关系，寻找标志性菌群或特定功能的基因。对微生物群落进行测序包括两类，一类是通过16s rDNA，18s rDNA，ITS区域进行扩增测序分析微生物的群体构成和多样性；还有一类是宏基因组测序，是不经过分离培养微生物，而对所有微生物DNA进行测序，从而分析微生物群落构成，基因构成，挖掘有应用价值的基因资源。以16s rDNA扩增进行测序分析主要用于微生物群落多样性和构成的分析，目前的生物信息学分析也可以基于16s rDNA的测序对微生物群落的基因构成和代谢途径进行预测分析，大大拓展了我们对于环境微生物的微生态认知。目前我们根据16s的测序数据可以将微生物群落分类到种（species）（一般只能对部分菌进行种的鉴定），甚至对亚种级别进行分析，几个概念：16S rDNA（或16S rRNA）：16S rRNA 基因是编码原核生物核糖体小亚基的基因，长度约为1542bp，其分子大小适中，突变率小，是细菌系统分类学研究中最常用和最有用的标志。16S rRNA基因序列包括9个可变区和10个保守区，保守区序列反映了物种间的亲缘关系，而可变区序列则能体现物种间的差异。16S rRNA基因测序以细菌16S rRNA基因测序为主,核心是研究样品中的物种分类、物种丰度以及系统进化。OTU：operational taxonomic units (OTUs)在微生物的免培养分析中经常用到，通过提取样品的总基因组DNA，利用16S rRNA或ITS的通用引物进行PCR扩增，通过测序以后就可以分析样品中的微生物多样性，那怎么区分这些不同的序列呢，这个时候就需要引入operational taxonomic units，一般情况下，如果序列之间，比如不同的 16S rRNA序列的相似性高于97%就可以把它定义为一个OTU，每个OTU对应于一个不同的16S rRNA序列，也就是每个OTU对应于一个不同的细菌（微生物）种。通过OTU分析，就可以知道样品中的微生物多样性和不同微生物的丰度。测序区段：由于16s rDNA较长（1.5kb），我们只能对其中经常变化的区域也就是可变区进行测序。16s rDNA包含有9个可变区，分别是v1-v9。一般我们对v3-v4双可变区域进行扩增和测序，也有对v1-v3区进行扩增测序。

主成分分析(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