|
目录
摘要工具与方法流程与命令数据录入，去重与质控结果展示
生成features table（OTU）结果展示
发育树生成部分结果展示
注释分类结果展示
总结
摘要
前两个月项目特别多，最近终于有机会闲下来写点文章，把之前搭建的流程梳理一遍。前同事分析16S/ITS使用的qiime1，我接手后感觉不太适应，希望能够使用新版本来搭建，于是花了几天时间重新搭了这个流程，
工具与方法
使用工具：qiime2 使用版本：qiime2-2021.4 参考文档：https://docs.qiime2.org/2021.4/（最新版本2021.8）
流程与命令
梳理一下16S全套分析的流程，这里借用了参考模板,需要说明的是qiime2现在已经不是计算OTUs，而是ASVs,得到feature table。本流程分为上下两章，本章为上章，主要介绍质控，生成特征表，以及注释富集等步骤。 
数据录入，去重与质控
qiime tools import --type 'SampleData[PairedEndSequencesWithQuality]' --input-path all_sample.txt --output-path paired-demux.qza --input-format PairedEndFastqManifestPhred33V2 #数据录入
qiime cutadapt trim-paired --p-cores ${threads} --i-demultiplexed-sequences paired-demux.qza --p-front-f ${forward_adapt} --p-front-r ${reverse_adapt} --o-trimmed-sequences paired-end-demux.qza --verbose #质控，需要输入数据的adapter
qiime demux summarize \--i-data paired-end-demux.qza \--o-visualization demux.qzv #获取质控结果
第二步，调用cutadapt时，adapter可以找供应商提供，如果没有提供的话，可以使用常规接口（16S，ITS有通用adapter），或者用vim命令打开数据看看首尾端的重复序列，直接复制过来。 展示一下demux.qzv的结果，这个结果个人感觉意义不大，因此我没有整理给客户,但是对于后续分析还是比较重要，稍后会解释。
结果展示
 
生成features table（OTU）
qiime dada2 denoise-paired \--i-demultiplexed-seqs paired-end-demux.qza \--p-trunc-len-f 0 \--p-trunc-len-r 0 \--o-table table.qza \--o-representative-sequences rep-seqs.qza \--o-denoising-stats denoising-stats.qza ## 使用dada2软件包识别特征序列，创建特征表
qiime feature-table summarize --i-table table.qza --o-visualization table.qzv --m-sample-metadata-file sample-metadata.txt #获取特征表统计结果
qiime tools extract --input-path table.qzv --output-path ../02.OTU_table #输出统计表结果文件
qiime tools export --input-path table.qza --output-path phyloseq
biom convert -i phyloseq/feature-table.biom -o phyloseq/otu_table.tsv --to-tsv #输出统计表tsv格式
qiime tools export --input-path rep-seqs.qza --output-path phyloseq #输出特征序列
qiime2生成特征表的算法有两个，dada2和Deblur,常规使用的是前者，看参考文档说的是dada2准确性更高一些。这行命令输入的 --p-trunc-len-f 和--p-trunc-len-r两个参数分别代表双端序列的截取位置。一般质量较好的测序数据我们默认为0就行，代表所有序列的整条片段全部纳入分析。如何判定序列质量较好？这就需要我们回到上一步的结果来看。  根据上图可以看到，大概在200bp之后质量分数会降到10左右，大部分质量分数都较好或者在25，符合咱们的分析要求。当然，我之前也遇到过一个测序质量很差的数据。当时没有意识到这个问题，结果没有截取序列的情况下可能有太多噪音，就没有特征表生成。 另外要补充一点，qiime2生成的可视化文件是.qzv格式，想看结果必须到qiime2 view网站提交查看。而qiime tools extract命令提供了将.qzv结果文件里的可视化图表以文件形式进行输出。之后所有的分析结果都会用到这个命令。
结果展示
发育树生成
qiime phylogeny align-to-tree-mafft-fasttree --i-sequences rep-seqs.qza --o-alignment aligned-rep-seqs.qza --o-masked-alignment masked-aligned-rep-seqs.qza --o-tree unrooted-tree.qza --o-rooted-tree rooted-tree.qza
qiime tools export --input-path unrooted-tree.qza --output-path phyloseq/unrooted #unrooted-tree
qiime tools export --input-path rooted-tree.qza --output-path phyloseq/rooted #rooted-tree
这一步通过特征序列来生成发育树，包括无根发育树和有根发育树两个。结果是.nwk格式，可以用第三方网站打开，MEGA也可以打开，就是感觉意义不大。这个和后面beta分析中的UPGMA生成的样品同源发育树还有些区别，需要注意区分。另外网上有文档介绍修改特征序列的名称，把乱码改成统一名称，这个看大家个人需求。
部分结果展示
注释分类
qiime feature-classifier classify-sklearn --i-classifier ${taxonomy_dataset} --i-reads rep-seqs.qza --o-classification taxonomy.qza
qiime metadata tabulate --m-input-file taxonomy.qza --o-visualization taxonomy.qzv #转换为可视化结果
qiime tools export --input-path taxonomy.qza --output-path phyloseq #输出.tsv表格结果
qiime taxa barplot --i-table table.qza --i-taxonomy taxonomy.qza --m-metadata-file sample-metadata.txt --o-visualization taxa-bar-plots.qzv #taxa-barplot
qiime tools extract --input-path taxa-bar-plots.qzv --output-path ../03.Taxonomy #输出注释富集柱状图
cp ../script/heatmap.R ../03.Taxonomy/*/data/
cd ../03.Taxonomy/*/data/
Rscript heatmap.R
cd ../../../01.data
qiime gneiss correlation-clustering --i-table table.qza --o-clustering cluster_tree #生成聚类树
qiime tools export --input-path cluster_tree.qza --output-path ../03.Taxonomy/
这一步注释的时间会稍微长一些，需要耐心等待，需要说明的是，参数${taxonomy_dataset} 代表你比对的数据库，我这个流程中，16S和ITS的区别就在数据集的差别。一个来自silva，另一个来自Unite。当然网站提供了不同软件，不同处理形式的数据库，比如99，97，冗余，非冗余等等，大家可以好好了解一下。 dataset_ITS="Unite-ITS-97-classifier.qza" dataset_16S="silva-138-99-nb-classifier.qza"
结果展示
  这个注释富集柱状图只支持在线展示，因此输出结果是网页格式。好处是可以自行调整界门纲目科属种的富集结果，并且可以根据之前的分组来进行调整。至于热图的绘制，在另一篇文章中有介绍。
总结
本篇介绍了微生物的前期质控，特征表分析，以及注释分析几个步骤，下章将会主要介绍Alpha，Beta多样性分析，以及探讨差异分组的和稀释曲线的结果。有兴趣的小伙伴可以加群交流，或者添加作者VX：bbplayer2021,这边邀请进群。 
|