16s RNA ：16s rDNA 基因存在于所有细菌的基因组中，具有高度的保守性。 该序列包含9个高变区和10个保守区，通过对某一段高变区序列进行PCR扩增后进行测序，得到1500bp左右的序列。
宏基因组：是将微生物基因组DNA随机打断成500bp的小片段，然后在片段两端加入通用引物进行PCR扩增测序，再通过组装的方式，将小片段拼接成较长的序列。
针对16s RNA 数据主要使用 QIIME2 软件进行。
此分析笔记流程概览
——8 以后的分析待后续更新吧！
数据下载：ENA 浏览器 (ebi.ac.uk)
0.数据下载 PRJNA359624
点击获取下载脚本，到服务器上运行
nohup ./ena-file-download-read_run-PRJNA359624-fastq_ftp-20240111-0625.sh &
确定下载成功，一共52个样本
分析流程如下：
要自己安装fastqc哦
1.使用fastqc进行质控
mkdir fastqout
fastqc -c *.fastq.gz -o fastqout/
看见报告正常生成
multiqc 整合质控报告
multiqc fastqout/
查看质控报告，在当前目录下生成，
质量数据没啥问题，下一步
2.切除引物
- 送公司测序返还的数据一般都是拆分过并去除了引物的，可以自己再做一下质检，引物切除了么？？
发现此处我们下载的已经切除，不用再去。
3. 数据导入：在qiime2中分析测序数据
3.0 安装 qiime2
参考官网 本机安装 QIIME 2 — QIIME 2 2023.9.2 文档
wget https://data.qiime2.org/distro/amplicon/qiime2-amplicon-2023.9-py38-linux-conda.yml
conda env create -n qiime2-amplicon-2023.9 --file qiime2-amplicon-2023.9-py38-linux-conda.yml
conda activate qiime2-amplicon-2023.9
确定安装成功，我大概花了20分钟吧~
3.1 先查看下测序类型质量值体系
质量值体系分为 Phred33 和 Phred 64两种，如下图所示，一般看fastq文件的质量值那行包含！和？（对应ASCII值33和63）等，即为Phred33体系（一般都为Phred33）。
发现全为?
3.2 写参数配置文件
下面为了运行 qiime2 需要自己写个文件 manifest file, 按以下格式写，
开头的行是注释行会被自动忽略，例如以下命名为为se-33-manifest的文件，保存为txt等文件
第一列为sample-ID, 第二列是绝对路径，如果是conda环境里的路径需要调整一下。第三列是序列文件对应的是reverse还是forward。
我们这次示例的是个是双端合并后的文件，所以不需要区分正向还是反向，direction一列，每个都写forward。然后导出为UTF-8（逗号分隔）文件。格式写成后长这样。
sample-id,absolute-filepath,direction
SRR5142316,/xtdisk/PRJNA359624/SRR5142316.fastq,forward
SRR5142325,/xtdisk/PRJNA359624/SRR5142325.fastq.gz,forward
SRR5142334,/xtdisk/PRJNA359624/SRR5142334.fastq.gz,forward
如是双端测序，分了两个fastq文件，格式写成后大概长这样，需要定指定 forward和reverse
sample-id,absolute-filepath,direction
SRR5142316,/xtdisk/PRJNA359624/SRR5142316_R1.fastq.gz,forward
SRR5142316,/xtdisk/PRJNA359624/SRR5142316_R2.fastq.gz,reverse
SRR5142325,/xtdisk/PRJNA359624/SRR5142325_R1.fastq.gz,forward
SRR5142334,/xtdisk/PRJNA359624/SRR5142334_R2.fastq.gz,reverse
- 一定要有表头 sample-id,absolute-filepath,direction
3.3 运行 qiime tools，导入数据，生成数据摘要
qiime tools import \
--input-path se-33-manifest.txt \
--output-path qiime_out_manifest.qza \
--type SampleData[JoinedSequencesWithQuality] \
--input-format SingleEndFastqManifestPhred33
*双端序列 --input-format PairedEndFastqManifestPhred33
--type 'SampleData[PairedEndSequencesWithQuality]'
*双端序列合并 --input-format SingleEndFastqManifestPhred33
--type SampleData[JoinedSequencesWithQuality]
结果文件出来了，很快~
3.4 可视化数据质量
#可视化文件
qiime demux summarize \
--i-data qiime_out_manifest.qza \
--o-visualization qiime_out_manifest.qzv
qiime_out_manifest.qzv文件需要从Linux中下载后再拖拽到qiime2 view网页中才能打开。
此处可以得到质检矢量图，通过放大观察可以清楚的判断碱基质量明显下降的位置，从而辅助确定下一步中流程中的参数reads1_cutpoint和reads2_cutpoint。
- 生成的.qzv文 拖拽进网页查看(推荐) QIIME 2 查看
- 或是使用qiime tools view qiime_out_manifest.qzv查看
这里发现数据质量都是30，其实这个数据有点奇怪。
此处可以得到质检矢量图，通过放大观察可以清楚的判断碱基质量明显下降的位置，从而辅助确定下一步中的参数t。
4 过滤， 切割、去嵌合体、拼接等
这步可使用Deblur 或 DATA 2。
Deblur VS DADA2： 方法不同，使用相同数据集的序列和特征数量的也会有差异。
Deblur 要求所有序列的长度相同。另一方面，DADA2 的特征序列具有不同的长度。DADA2 纠正错误，Deblur 删除错误。这种差异对到达终点的总读取量有很大影响，但对频率的影响要小得多，而频率更重要。
如果使用 DADA2来合并和消除双端数据的噪声，请在用DADA2去噪之前不要合并您的序列；DADA2希望读长尚未合并的序列，并将在去噪过程中为您双端合并。deblur只接受single-end的序列，如果你把没有jion的序列传递给deblur，它会只处理forward seqs,eblur在denoising时需要输入整齐一样长度的序列，所以需要trim成相同的长度。
目前：无论选择使用哪种方法，生物学解释（至少从表面上看）似乎都非常相似
目前关于哪个更好，没有达成共识，根据自己的数据情况来选择。此教程我们的数据是合并后的数据，
那么选择Deblur，注意一些数据分析比较原则，如果是整合数据集，那么保证每个样本使用相同的分析流程。
4.1 质量过滤
qiime quality-filter q-score \
--i-demux qiime_out_manifest.qza \
--o-filtered-sequences demux-joined-filtered.qza \
--o-filter-stats demux-joined-filter-stats.qza
输出对象：
-
demux-joined-filter-stats.qza: 统计结果。 -
demux-joined-filtered.qza: 数据过滤后结果。
此方法的参数尚未在双端合并的数据上进行广泛的基准测试，因此建议尝试使用不同的参数设置。
4.2 去噪
现在已经准备好用Deblur去噪你的序列了。您应该从质量分数图中为--p-trim-length选择合适的序列长度值。这将把所有序列修剪到这个长度，并丢弃任何小于这个长度的序列。
大约 251-255（EMP 引物组）捕获了几乎所有的 V4 序列。因此，修剪到 251-255 是合理的。
Deblur的开发者们建议设置一个质量分数开始迅速下降的长度（recommend setting this value to a length where the median quality score begins to drop too low）。
qiime deblur denoise-16S \
--i-demultiplexed-seqs demux-joined-filtered.qza \
--p-trim-length 250 \
--p-sample-stats \
--o-representative-sequences rep-seqs.qza \
--o-table table.qza \
--o-stats deblur-stats.qza
输出对象：
-
rep-seqs.qza: 代表序列。 -
deblur-stats.qza: 统计过程。 -
table.qza: 特征表。
4.3 查看结果特征表
统计结果可视化
qiime feature-table summarize \
--i-table table.qza \
--o-visualization table.qzv
代表序列统计
qiime feature-table tabulate-seqs \
--i-data rep-seqs.qza \
--o-visualization rep-seqs.qzv
5 多样性分析
5.1 建树用于多样性分析
qiime phylogeny align-to-tree-mafft-fasttree \
--i-sequences rep-seqs.qza \
--o-alignment aligned-rep-seqs.qza \
--o-masked-alignment masked-aligned-rep-seqs.qza \
--o-tree unrooted-tree.qza \
--o-rooted-tree rooted-tree.qza
5.2 Alpha多样性分析
α多样性是度量单个样本内有多少种微生物物种，以及每个物种所占比例的指标。
计算多样性(包括所有常用的Alpha和Beta多样性方法)，输入有根树、Feature表、样本重采样深度
取样深度看table.qzv文件确定（一般为样本最小的sequence count，或覆盖绝大多数样品的sequence count）
还需要构建一个 sample-metadata.txt file ，这里我们的演示，从NCBI下载
Run Selector :: NCBI (nih.gov)
改写"RUN" 为 sampleID，这里我们添加了一个Group 列，根据Library Name 分了两组Gout_Train 和HC_Train，注意格式，保存为分隔符为空格的txt 文件。
然后运行以下流程：
qiime diversity core-metrics-phylogenetic \
--i-phylogeny rooted-tree.qza \
--i-table table.qza \
--p-sampling-depth 1775 \
--m-metadata-file sample-metadata.txt \
--output-dir core-metrics-results
运行很快，
输出结果包括多种多样性结果，文件列表和解释如下：
beta多样性bray_curtis距离矩阵 bray_curtis_distance_matrix.qza
alpha多样性evenness(均匀度，考虑物种和丰度)指数 evenness_vector.qza
alpha多样性faith_pd(考虑物种间进化关系)指数 faith_pd_vector.qza
beta多样性jaccard距离矩阵 jaccard_distance_matrix.qza
alpha多样性observed_otus(OTU数量)指数 observed_otus_vector.qza
alpha多样性香农熵(考虑物种和丰度)指数 shannon_vector.qza
beta多样性unweighted_unifrac距离矩阵，不考虑丰度 unweighted_unifrac_distance_matrix.qza
beta多样性unweighted_unifrac距离矩阵，考虑丰度 weighted_unifrac_distance_matrix.qza
测试分类元数据(sample-metadata)列和alpha多样性数据之间的关联，输入多样性值、sample-medata，输出统计结果，统计faith_pd算法Alpha多样性组间差异是否显著
qiime diversity alpha-group-significance \
--i-alpha-diversity core-metrics-results/faith_pd_vector.qza \
--m-metadata-file sample-metadata.txt \
--o-visualization core-metrics-results/faith-pd-group-significance.qzv
# 统计evenness组间差异是否显著
qiime diversity alpha-group-significance \
--i-alpha-diversity core-metrics-results/evenness_vector.qza \
--m-metadata-file sample-metadata.txt \
--o-visualization core-metrics-results/evenness-group-significance.qzv
以evenness-group-significance.qzv为例，图中可点Category选择分类方法，查看不同分组下箱线图间的分布与差别。图形下面的表格，详细详述组间比较的显著性和假阳性率统计。
发现 Alpha 多样性和 我们的分组没啥关系，P值不显著
5.3 Beta 多样性分析
β多样性是度量不同样本间菌群组成的相似度大小的指标，即关注各样本间的菌群组成差异。
# 统计group 的组间是否有显著差异,其他的分组分析类似。
qiime diversity beta-group-significance \
--i-distance-matrix core-metrics-results/unweighted_unifrac_distance_matrix.qza \
--m-metadata-file sample-metadata.txt \
--m-metadata-column Group \
--o-visualization core-metrics-results/unweighted-unifrac-subject-significance.qzv \
--p-pairwise
可视化结果：
发现有P值。
# 可视化三维展示其它类，定量值的 分析
qiime emperor plot \
--i-pcoa core-metrics-results/unweighted_unifrac_pcoa_results.qza \
--m-metadata-file sample-metadata.txt \
--p-custom-axes weight \
--o-visualization core-metrics-results/unweighted-unifrac-emperor-weight.qzv
我们这里演示的数据集没有，跳过。
5.4 Alpha rarefaction plotting
# --p-max-depth should be determined by reviewing the “Frequency per sample” information presented in the table.qzv file
# --p-max-depth一般取table.qzv文件Frequency per sample的中位数左右
qiime diversity alpha-rarefaction \
--i-table table.qza \
--i-phylogeny rooted-tree.qza \
--p-max-depth 31131 \
--m-metadata-file sample-metadata.txt \
--o-visualization alpha-rarefaction.qzv
可视化将有两个图。顶部图是α稀疏图，主要用于确定样品的丰富度是否已被完全观察或测序。如果图中的线在沿x轴的某个采样深度处看起来“平坦化”（即接近零斜率），则表明收集超出该采样深度的其他序列将不可能会有其他的OTU（feature）产生。如果图中的线条没有达到平衡，这可能是因为尚未完全观察到样品的丰富程度（因为收集的序列太少），或者它可能表明在数据中存在大量的测序错误（被误认为是新的多样性）。底部图表示当特征表稀疏到每个采样深度时每个组中保留的样本数。样本被分成两组,图中显示即两条线.
6 菌群比对数据库：引物特异性菌群比对
Greengenes细菌16S全长序列数据库
下载地址:Data resources — QIIME 2 2023.9.2 documentation
下载得到gg_13_8_otus.tar.gz后将其解压,具体使用gg_13_8_otus/rep_set 下99_otus.fasta（序列）和taxonomy/99_otu_taxonomy.txt（分类）两个文件
tar -zxvf gg_13_8_otus.tar.gz
6.1. 转换格式
#将99_otus.fasta（序列）和99_otu_taxonomy.txt（分类）两个文件的格式转换QIIME2能识别和利用的格式
qiime tools import \
--type 'FeatureData[Sequence]' \
--input-path gg_13_8_otus/rep_set/99_otus.fasta \
--output-path 99_otus.qza
qiime tools import \
--type 'FeatureData[Taxonomy]' \
--input-format HeaderlessTSVTaxonomyFormat \
--input-path gg_13_8_otus/taxonomy/99_otu_taxonomy.txt \
--output-path 99-taxonomy.qza
6.2. 抽提V3V4模板序列
用测序引物序列从Greengenes数据库中的16S全长序列99_otus.qza中抽提出引物特异性的细菌参考序列，就会得到本研究特异性的参考序列
qiime feature-classifier extract-reads \
--i-sequences 99_otus.qza \
--p-f-primer CCTACGGGNGGCWGCAG \
--p-r-primer GACTACHVGGGTATCTAATCC \
--o-reads 99-v3v4-seqs.qza
大概20分钟
6.3. 训练V3V4分类器
qiime feature-classifier fit-classifier-naive-bayes \
--i-reference-reads 99-v3v4-seqs.qza \
--i-reference-taxonomy 99-taxonomy.qza \
--o-classifier gg-13-8-99-v3v4-classifier.qza
大概20分钟 等待中~
7 菌群特征比对，分类学注释
7.1. 比对
把DADA2分析得到的细菌特征代表性序列rep-seqs.qza和训练好的分类器gg-13-8-99-v3v4-classifier.qza进行比对，获得具体的细菌分类信息taxonomy.qza
qiime feature-classifier classify-sklearn \
--i-classifier gg-13-8-99-v3v4-classifier.qza \
--i-reads rep-seqs.qza \
--o-classification taxonomy.qza
7.2. 丰度统计
将细菌特征丰度表table.qza和细菌分类信息taxonomy.qza进行整合获得完整的细菌分类丰度表，包含界、门、纲、目、科、属、种多水平的细菌丰度信息
qiime taxa barplot \
--i-table table.qza \
--i-taxonomy taxonomy.qza \
--m-metadata-file sample-metadata.txt \
--o-visualization taxa-bar-plots.qzv
8 数据处理: 获取菌属分类表
8.1. 获取属水平菌丰度表
QIIME2 view网页中打开taxa-bar-plots.qzv，下载level6的CSV文件，如下：
8.2. 标准化菌属丰度表
把CSV文件导入到Excel中进行标准化，即每个菌属的原始丰度除以该菌所在样本的总菌属丰度得到标准相对菌属相对丰度
