Ubuntu+qiime2 16s扩增子变异分析与rrnDB使用全过程

所需环境

Windows系统，商城下载Ubuntu即可，不过默认装在C盘，我之前没有注意，然后在使用的过程中，又把Ubuntu迁到了D盘。所以最好在安装的时候注意一下，免得后面麻烦。qiime2的安装有很多方法，最友好的方法是conda安装，由于conda比较大，当时没有换源下载，比较慢，最后选择了miniconda3。以及为了批量从NCBI上下载数据，还需要安装sratoolkit插件。以上过程其他博客都有介绍，不再赘述。

批量下载文件

首先文献中会给出accession number，搜索后如下：

选中所有的sample，点击右上角send to，选择run selector，go之后得到如下界面：

选择Total这一行中的accession list，metadata（后面要用）。然后输入：

prefetch --option-file SRR_Acc_List.txt

默认下载的数据在sratoolkits文件夹下sra文件夹中，为了条理更清晰，可以自己创建文件夹，把数据剪切过去。

批量读入文件

自己创建一个excel表格，需要包含sample-id、绝对位置和方向这三列。sample-id就是accession list中的id，绝对位置是数据存放的位置，方向通常包括前向（forword）和后向（reverse），也就是测序的方向。需要说明如果是成对数据，需要先进行拆分，拆分出来后缀为1的文件则为前向，2的文件为后向。将excel表格转换成tsv文件，网上有很多在线转换网页。我的数据是single的，不需要拆分。

表格格式如下：

cd到表格和数据存放的文件夹下，单端输入：

qiime tools import --type 'SampleData[SequencesWithQuality]'  --input-path manifest.tsv --output-path single-end-demux.qza  --input-format SingleEndFastqManifestPhred33 

把输出的single-end-demux.qza可视化一下，看看质量怎么样，觉得接下来切除多少。

qiime demux summarize  --i-data single-end-demux.qza  --o-visualization demux.qzv

在qiime2 view上看结果：

可以看到序列几乎都高于最高值的一半，所以我不需要切割序列。

（Ubuntu中$表示当前所处位置的绝对路径；cd ++：表示进入 ++文件夹；cd …表示退回上一层目录；ls表示当前目录下所包含的文件，也就是list的缩写）

DADA2 去噪

qiime dada2 denoise-single --i-demultiplexed-seqs single-end-demux.qza --p-trim-left 0 --p-trunc-len 100 --o-representative-sequences rep-seqs-dada2.qza --o-table table.qza --o-denoising-stats dada2-stats.qza 

--p-trim-left 0 #表示从左侧0开始
--p-trunc-len 100 #表示保留的长度为100，也就是没有切割，因为我的序列长度就是100

rrnDB的使用：

rep-seqs-dada2.qza —— 为代表序列，也就是按照这个序列划分了ASV表格
table.qza —— 也就是ASV表格
dada2-stats.qza —— 为统计状态

可视化代表序列：
qiime feature-table tabulate-seqs
–i-data rep-seqs-dada2.qza
–o-visualization rep-seqs-dada2.qzv
即可得到feature ID对应的序列是什么：

点击download your sequence as a raw FASTA file，就得到：

将这个fasta文件upload到rrnDB数据库比对，可以设置置信度，得到：

第一个文件是根据RDP分类器进行的物种注释结果，第二个文件是没有矫正的物种分度结果，第三个数据是矫正之后的物种分度结果。二三表格只有sequences.fasta不同，用未校正的sequences.fasta除以矫正后的sequences.fasta得到的就是其所对应的rrna拷贝数了。（原因是sequences.fasta表示序列的丰度，但是未校正时存在多拷贝的情况，所以矫正就是除掉了多拷贝数。）

训练RDP分类器

如果自己训练RDP分类器得到物种注释等下游工作分析，需要下载RDP classifier的参考OTU表和参考分类表，由于Ubuntu网速较慢，我自己在Windows上下载拷贝到了Ubuntu。关于Ubuntu的文件存放位置，输入\wsl即可：

在此网页下载，下翻找到qiimeformat的zip，也就是下图第二个：

将它们读入，转换成qza文件：

qiime tools import \
--type 'FeatureData[Sequence]' \
--input-path RefOTUs.fasta \
--output-path RefOTUs.qza

qiime tools import \
--type 'FeatureData[Taxonomy]' \
--input-format HeaderlessTSVTaxonomyFormat \
--input-path Ref_taxonomy.txt \
--output-path Ref_taxonomy.qza

引物通常会给出，注意查找：

• -p-f:前向引物（forward）；- -p-r:后向（reverse）

qiime feature-classifier extract-reads \
  --i-sequences RefOTUs.qza \
  --p-f-primer GTGCCAGCMGCCGCGGTAA \
  --p-r-primer GGACTACHVGGGTWTCTAAT \
  --p-trunc-len 100 \
  --o-reads ref-seqs.qza

qiime feature-classifier fit-classifier-naive-bayes \
  --i-reference-reads ref-seqs.qza \
  --i-reference-taxonomy Ref_taxonomy.qza \
  --o-classifier classifier.qza

输入之前得到的代表序列，得到分类结果：

qiime feature-classifier classify-sklearn --i-classifier classifier.qza --i-reads rep-seqs-dada2.qza --o-classification taxonomy.qza 

将分类结果可视化：

qiime metadata tabulate --m-input-file taxonomy.qza --o-visualization taxonomy.qzv

如下：

绘制每个sample的物种组成bar图：

qiime taxa barplot --i-table table.qza --i-taxonomy taxonomy.qza --m-metadata-file sample-metadata.tsv --o-visualization taxa-bar-plots.qzv 

结果如下：