前面已经完成了qiime2-sidle插件的安装，测试方法就是输入qiime：
sidle Plugin for kmer-based marker gene reconstruction 出现了上面的选项，应该就说明已经安装成功了。
数据库准备 数据库准备是一劳永逸的，前面我们已经完成了数据库过滤的步骤
准备一个区域数据库 这一步是提取一个区域的数据库，基于K-mer，为了提升内存效率，把简并碱基和重复kmer作为一条序列。
# 首先，使用feature-classifier插件的extract-reads提取，那么对于6V区的SMURF应该是：# 2.2 序列提取# 74f 315R V1-V2qiime feature-classifier extract-reads \ 上sidle插件，文档介绍这里和原方法略有区别，原算法是允许2个碱基不一致，而这里用的80%,作者说可以用–p-identity 参数，但是没找到。。。
# 可以给region起个别名，看喜好起即可 这步会输出含简并的序列,移除重复序列和一个kmer映射文件，可以qiime可视化，给出原始数据库序列名称(db-seq)、序列名称(seq-name)、区域标识(Kmer)、引物(fwd-primer和rev-primer)、区域(Region)和序列长度(Trim-Length)之间的关系。
# 如果参考序列不能覆盖正反向引物，可以把--p-trim-length设置为负值，增加个参数--p-reverse-complement-result
完成了前面的数据库准备，下面就是reads的准备，基本过程就是把reads拆成对应不同引物的几个部分，后面再重建合并在一起啦。首先声明，这个方法还在开发和完善之路，最近一次更新在这个月，可能结果会有变动，应该说还处于beta版本中，不建议在生产环境中使用。
Reads准备 尽管SMURF依赖于质控过滤，还是推荐继续使用去噪的方法。推荐使用现有的方法进行预处理，当然可以用别的方法替代，这只是一个例子。
拆数据 我们这里用的是PE150的示例数据，没有进行样本拆分，那么使用多q2-cutadapt的trim-paired切去引物，–p-discard-untrimmed移除没有引物的序列。如果正反向序列不能拼接，当做单独的序列进行处理。
首先，导入数据，导入前重命名一下，以适应导入格式要求，个人觉得这样导入最简单。
# 重命名 拆分区域数据
# P1 从得到的数据看，各个区域的数据分布不是太一致，是不是多重PCR做的，引物偏好性有一定影响，所以这里面的影响也值得思考哈。
-rw-r--r-- 1 zjd zjd 2.6M Jan 31 11:00 P1_demux.qza 去噪 这里发现官方推荐的两去噪方法均报错呢，于是用了qiime2 vsearch解决。
# ######## 以下是我走过的弯路（可以略过）：
＃dada2for i in {1..6} deblur
##############选项2，deblur，耗时也较长# 序列质控6152.75s 继续前面的文档学习，地址地这里啦！ 官方文档
区域比对 第一步是每个区域把序列比对到数据库，使用 align-regional-kmers 命令，我们前面使用 --kmer-db-fp 选项设置了数据库，使用 --rep-seq-fp 选项传递ASV序列，最后是区域定义，来自前面你给区域起的别名，要完全一致。比对是一个开心的可并行任务，我们可以通过多多线程提升性能( --p-n-workers 参数)。
# 继续循环解决啦for i in {1..6} 这将输出一个对齐文件和任何不会与数据库对齐的 ASV 序列，以保留您自己的记录。另外，你可以设置与参考序列不同的长度，使用 --p-max-mismatch 参数（2-130）ps.我有点怀疑这样做有没有问题呀！
表重建 这步里，首先，是区域片段重新比对到完整的数据库序列，然后相对丰度使用一个优化的步骤进行计算，最后生成重建的计数表。 per-nucleotide-error 参数和 maximum-mismatch 参数，决定了比对和错配的比例，Ps.我的理解是过程中会产生错误累积。 min-abundance 决定了要排队的数据库序列的最小丰度，这个参数在某种意义上取决于测序深度和你的偏好。 --p-n-workers 参数)。
qiime sidle reconstruct-counts \ 好像这步是超级消耗资源的，我设置了个8核心，资源消耗竟然是＊5的，用上了全部40核心的资源。。。不过几分钟就结束了，那么如果设置1个核心，可能是需要8倍时间的。映射到某个区域的数据库Kmer数量以及ASV ID 的详细信息的文件，以及一个需要进行分类重建的映射。
＃可视化看下 这里我是用了一个样本，133个汇总的feature（也就是以前的OTU），40598条等效数据库序列。不能被解析 ，因此我们将该序列分配给这两个区域。重建中使用的区域越多，就越有可能准确地重建数据库序列。summary 。总结可以用来评估重建的质量；更多细节见SMURF文章 原稿 。默认情况下，摘要会将退化kmer视为唯一序列；您可以使用 count-degenerates 参数更改行为；如果为false，则仅当kmer属于唯一参考序列时才会对其进行计数。您可以通过将元数据制表来查看摘要。
＃ 汇总摘要 几个V区的计数明细和汇总
物种重建 重建特征表之后，需要重建物种，特别的，如果多个数据库序列不能得到处理时。该功能获取重建过程中生成的数据库映射和与数据库相关联的分类，并返回重建的分类。
1236 k__Bacteria; p__Firmictues; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__Blautia; s__obeum --database 参数允许选择一个数据库进行，如 greengenes, silva 或不提供。如果是greengenes 或 silva，一些专门的数据库清理将会自动进行，特别是缺失和不明处理参数。例如：
k__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Entrobacteriales; f__Enterobacteriaceae; g__; s__1 将会整理成：
k__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Entrobacteriales; f__Enterobacteriaceae; g__unsp. f. Enterobacteriaceae; s__unsp. f. Enterobacteriaceae1 这里我们用的是greengenes，
qiime sidle reconstruct-taxonomy \ 进化树重建 最后一步便是重建进化树了，如果不能解析参考序列，则不能简单地从数据库继承系统发育树。因此，我们需要重建一个新的系统树。我们以两种方式处理序列：
我们可以随机选择一个序列来映射重建的区域。然而，当有几个序列组合在一起时，这可能不起作用。因此，如果我们不能解析数据库序列，我们可以从组合的数据中计算一个融合序列 ，在我们能够绘制的区域中提取它们，然后这些一致的序列可以被插入到使用SEPP 或类似的系统发育参考主干中(必须使用相同的数据库版本)。
因此，第一步，对不能解析的序列重建共识序列，这一步官方教程的结果应该有问题，看报错信息进行了更正。
qiime sidle reconstruct-fragment-rep-seqs \ 现在，我们可以将序列放入参考树中了，首先下载参考树：
wget -O "sepp-refs-gg-13-8.qza" \ "https://data.qiime2.org/2020.11/common/sepp-refs-gg-13-8.qza" 然后，插入序列：
qiime fragment-insertion sepp \ 最后，就可以继续常规的qiime2分析流程了，多样性，统计等等，能不能用picrust2分析呢？或许也可以的。
好了，一波三折，终于完成了这个教程的学习，不得不说这个技术好像还不够成熟，特别是qiime2的这个插件。
转载本文请联系原作者获取授权，同时请注明本文来自赵加栋科学网博客。 https://blog.sciencenet.cn/blog-623545-1270944.html
上一篇：
SMURF(5R)-Science封面文章使用的16S新流程(二) 下一篇：
Kraken2 Vs qiime2 16S物种注释
