SMURF流程之q2-sidle（一）

用户1075469

发布于 2021-03-11 14:54:19

5260

发布于 2021-03-11 14:54:19

前面说到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

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2021-02-06，如有侵权请联系 cloudcommunity@tencent.com 删除

html

php

本文分享自科技记者微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

html

php

登录后参与评论

0 条评论

热度