SMURF流程之q2-sidle(一)

前面说到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
# 去除序列中有过多简并的,这里是3个,99%的标准
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
# 去除没有门或者界信息的物种 剩余5400条左右
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



本篇文章来源于微信公众号: 微因

发表评论