SMURF流程之q2-sidle（三）--reads准备

zd200572

于 2021-02-06 16:32:19 发布

阅读量717

点赞数 1

分类专栏：生物信息文章标签： 16S

本文链接：https://blog.csdn.net/zd200572/article/details/113457234

版权

生物信息专栏收录该内容

58 篇文章 23 订阅

订阅专栏

完成了前面的数据库准备，下面就是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