前面说到Science封面文章用的16S数据分析流程有qiime2的插件版本，可以解决基于matlab MCR standalone版本的报错，于是实践一下！https://github.com/jwdebelius/q2-sidle。conda的安装就不表了，教程挺多的。
1、环境准备
安装qiime2-2020.11作者说只测试了兼容这个版本，于是就装这个啦！
# 激活环境
source ~/data_home/Miniconda3/bin/activate
# 下载配置文件
wget https://data.qiime2.org/distro/core/qiime2-2020.11-py36-linux-conda.yml
# 修改配置为镜像，加速下载和安装，这里用的是北外的,conda的配置也要同样的镜像，防止冲突
vi qiime2-2020.11-py36-linux-conda.yml
# channels:
- qiime2/label/r2020.11
- https://mirrors.bfsu.edu.cn/anaconda/cloud/conda-forge
- https://mirrors.bfsu.edu.cn/anaconda/cloud/bioconda
- defaults
# 安装
conda env create -n qiime2-2020.11 --file qiime2-2020.11-py36-linux-conda.yml
# 激活环境
conda activate qiime2-2020.11
# 安装插件和相关依赖
conda install dask regex
conda install -c conda-forge -c bioconda -c qiime2 -c defaults xmltodict
pip install git+https://gitee.com/zd200572/RESCRIPt.git
pip install git+https://gitee.com/zd200572/q2-sidle
qiime dev refresh-cache
# 安装完成
# 数据库准备
wget https://gitee.com/zd200572/q2-sidle/raw/main/docs/tutorial_data/database.zip
unzip database.zip
cd database
# 两个文件
# sidle-db-full-sequences.qza, 全长数据库
# sequences and sidle-db-taxonomy.qza 物种注释
# 可视化看下 5649条序列，这么少？
qiime feature-table tabulate-seqs
--i-data sidle-db-full-sequences.qza
--o-visualization sidle-db-full-sequences.qzv
# 去除序列中有过多简并的，这里是３个，９９％的标准
qiime sidle filter-degenerate-sequences
--i-sequences sidle-db-full-sequences.qza
--p-max-degen 3
--o-filtered-sequences sidle-db-full-degen-filtered-sequences.qza
# 去除没有门或者界信息的物种 剩余５４００条左右
qiime taxa filter-seqs
--i-sequences sidle-db-full-degen-filtered-sequences.qza
--i-taxonomy sidle-db-taxonomy.qza
--p-exclude "p__;,k__;"
--p-mode contains
--o-filtered-sequences sidle-db-full-degen-filtered-phylum-def-sequences.qza
完成了前面的数据库准备，下面就是reads的准备，基本过程就是把reads拆成对应不同引物的几个部分，后面再重建合并在一起啦。首先声明，这个方法还在开发和完善之路，最近一次更新在这个月，可能结果会有变动，应该说还处于beta版本中，不建议在生产环境中使用。这里就有几种情况啦，一种是已经每个样本每个V区拆好的数据，另一种是每个样本几个V区混在一起的数据，或者完全没拆的数据。这里根据SMURF的示例，按第二种情况进行，应该是最常见的情况。下面是具体步骤：
Reads准备
尽管SMURF依赖于质控过滤，还是推荐继续使用去噪的方法。推荐使用现有的方法进行预处理，当然可以用别的方法替代，这只是一个例子。简单来说步骤就是分而治之，最后合并，作者也说了这个方法其实是可以用来meta分析的，但是我还是对meta分析持怀疑态度的，毕竟每个实验室使用的方法区别那么大，样本保存条件不一样，提取方法有区别，再有就是PCR扩增区域、引物和酶也是有区别，还有建库方法的不一致，这样下来，区别差到天涯海角了，meta分析只能得出一个有问题的结论吧！16S/宏基因组，微生物、微生态的研究亟需一个标准化吧，否则每个结果只能对每个实验室的方法负责了，完全失去了可比性。对于这篇文章中使用的引物，我发现多数不是最常用的引物，所以结果是不是有偏差，也是一个值得思考的问题呀！好，回归正题！
拆数据
我们这里用的是PE150的示例数据，没有进行样本拆分，那么使用多q2-cutadapt的trim-paired切去引物，–p-discard-untrimmed移除没有引物的序列。如果正反向序列不能拼接，当做单独的序列进行处理。
首先，导入数据，导入前重命名一下，以适应导入格式要求，个人觉得这样导入最简单。
# 重命名
mv Example_L001_R1_001.fastq.gz 1_Example_L001_R1_001.fastq.gz
mv Example_L001_R2_001.fastq.gz 1_Example_L001_R2_001.fastq.gz
# 导入
qiime tools import
--type 'SampleData[PairedEndSequencesWithQuality]'
--input-path .
--input-format CasavaOneEightSingleLanePerSampleDirFmt
--output-path demux.qza
拆分区域数据
# P1
qiime cutadapt trim-paired
--i-demultiplexed-sequences demux.qza
--p-front-f TGGCGGACGGGTGAGTAA
--p-front-r CTGCTGCCTCCCGTAGGA
--p-discard-untrimmed
--p-error-rate 0.15
--o-trimmed-sequences ./P1_demux.qza
# P2
qiime cutadapt trim-paired
--i-demultiplexed-sequences demux.qza
--p-front-f TCCTACGGGAGGCAGCAG
--p-front-r TATTACCGCGGCTGCTGG
--p-discard-untrimmed
--p-error-rate 0.15
--o-trimmed-sequences ./P2_demux.qza
# P3
qiime cutadapt trim-paired
--i-demultiplexed-sequences demux.qza
--p-front-f CAGCAGCCGCGGTAATAC
--p-front-r CGCATTTCACCGCTACAC
--p-discard-untrimmed
--p-error-rate 0.15
--o-trimmed-sequences ./P3_demux.qza
# P4
qiime cutadapt trim-paired
--i-demultiplexed-sequences demux.qza
--p-front-f AGGATTAGATACCCTGGT
--p-front-r GAATTAAACCACATGCTC
--p-discard-untrimmed
--p-error-rate 0.15
--o-trimmed-sequences ./P4_demux.qza
# P5
qiime cutadapt trim-paired
--i-demultiplexed-sequences demux.qza
--p-front-f GCACAAGCGGTGGAGCAT
--p-front-r CGCTCGTTGCGGGACTTA
--p-discard-untrimmed
--p-error-rate 0.15
--o-trimmed-sequences ./P5_demux.qza
# P6
qiime cutadapt trim-paired
--i-demultiplexed-sequences demux.qza
--p-front-f AGGAAGGTGGGGATGACG
--p-front-r CCCGGGAACGTATTCACC
--p-discard-untrimmed
--p-error-rate 0.15
--o-trimmed-sequences ./P6_demux.qza
从得到的数据看，各个区域的数据分布不是太一致，是不是多重PCR做的，引物偏好性有一定影响，所以这里面的影响也值得思考哈。
-rw-r--r-- 1 zjd zjd 2.6M Jan 31 11:00 P1_demux.qza
-rw-r--r-- 1 zjd zjd 466K Jan 31 11:32 P2_demux.qza
-rw-r--r-- 1 zjd zjd 13M Jan 31 11:33 P3_demux.qza
-rw-r--r-- 1 zjd zjd 11M Jan 31 11:34 P4_demux.qza
-rw-r--r-- 1 zjd zjd 3.5M Jan 31 11:34 P5_demux.qza
-rw-r--r-- 1 zjd zjd 22M Jan 31 11:35 P6_demux.qza
去噪
这里发现官方推荐的两去噪方法均报错呢，于是用了qiime2 vsearch解决。
# ########
# Vsearch
###########
for i in {1..6}
do
qiime vsearch join-pairs
--i-demultiplexed-seqs P${i}_demux.qza
--o-joined-sequences P${i}_demux-joined.qza
time qiime quality-filter q-score
--i-demux P${i}_demux-joined.qza
--o-filtered-sequences P${i}_demux-joined-filtered.qza
--o-filter-stats P${i}_demux-joined-filter-stats.qza
qiime vsearch dereplicate-sequences
--i-sequences P${i}_demux-joined-filtered.qza
--o-dereplicated-table P${i}_table.qza
--o-dereplicated-sequences P${i}_rep-seqs.qza
echo vsearch denovo start
date
#denovo pick otu
qiime vsearch cluster-features-de-novo
--i-table P${i}_table.qza
--i-sequences P${i}_rep-seqs.qza
--p-perc-identity 0.99
--o-clustered-table P${i}_table-dn-99.qza
--o-clustered-sequences P${i}_rep-seqs-dn-99.qza
echo vsearch end
done
以下是我走过的弯路（可以略过）：dada2和deblur，两个算法，习惯于用第一个啦，遗憾报错，换了个2020.2月版本依然，上后者。
＃dada2
for i in {1..6}
do
qiime dada2 denoise-paired
--i-demultiplexed-seqs P${i}_demux.qza
--p-trim-left-f 0
--p-trim-left-r 0
--p-trunc-len-f 0
--p-trunc-len-r 0
--o-table ${i}_table.qza
--o-representative-sequences P${i}_rep-seqs.qza
--o-denoising-stats P${i}_denoising-stats.qza
done
# Plugin error from dada2:
# An error was encountered while running DADA2 in R (return code 1), please inspect stdout and stderr to learn more.
# Debug info has been saved to /tmp/qiime2-q2cli-err-wqromfhv.log
deblur
##############选项2，deblur，耗时也较长
# 序列质控6152.75s
for i in {1..6}
do
qiime vsearch join-pairs
--i-demultiplexed-seqs P${i}_demux.qza
--o-joined-sequences P${i}_demux-joined.qza
time qiime quality-filter q-score
--i-demux P${i}_demux-joined.qza
--o-filtered-sequences P${i}_demux-joined-filtered.qza
--o-filter-stats P${i}_demux-joined-filter-stats.qza
# # # deblur去噪生成特征表
time qiime deblur denoise-16S
--i-demultiplexed-seqs P${i}_demux-joined.qza
--p-sample-stats
--o-representative-sequences P${i}_rep-seqs-deblur.qza
--o-table P${i}_table-deblur.qza
--o-stats P${i}_deblur-stats.qza
done
time qiime deblur denoise-16S
--i-demultiplexed-seqs P1_demux-joined-filtered.qza
--p-trim-length 236
--p-sample-stats
--o-representative-sequences P1_rep-seqs-deblur.qza
--o-table P1_table-deblur.qza
--o-stats P1_deblur-stats.qza
qiime sidle trim-dada2-posthoc
--i-table table-dada2.qza
--i-representative-sequences rep-seqs-dada2.qza
--p-trim-length 100
--o-trimmed-table table-dada2-100nt.qza
--o-trimmed-representative-sequences rep-seq-dada2-100nt.qza
qiime demux summarize
--i-data P1_demux-joined.qza --o-visualization demux-subsample.qzv
本篇文章来源于微信公众号: 微因