SMURF流程之q2-sidle（二）– 序列重建

继续前面的文档学习，地址在这里啦！官方文档‎ SMURF 算法的核心是基于基于 kmer 的短区域重建到全长框架中。有两个步骤，首先是ASV在单个区域基于kmer进行比对，然后完整的序列集组装成重建的计数表。

区域比对

第一步是每个区域把序列比对到数据库，使用 align-regional-kmers 命令，我们前面使用--kmer-db-fp选项设置了数据库，使用  --rep-seq-fp选项传递ASV序列，最后是区域定义，来自前面你给区域起的别名，要完全一致。比对是一个开心的可并行任务，我们可以通过多多线程提升性能(--p-n-workers参数)。

# 继续循环解决啦
for i in {1..6}
 do
qiime sidle align-regional-kmers 
--i-kmers ../database/../database/sidle-db-P${i}-100nt-kmers.qza 
--i-rep-seq P${i}_rep-seq-100nt.qza 
--p-region P${i} 
--p-n-workers 4 
--o-regional-alignment alignment/P${i}_align-map.qza
done

‎这将输出一个对齐文件和任何不会与数据库对齐的 ASV 序列，以保留您自己的记录。‎另外，你可以设置与参考序列不同的长度，使用--p-max-mismatch参数（2-130）ps.我有点怀疑这样做有没有问题呀！‎您可能会发现，如果您有更长的 kmers，您可能需要相应地增加此参数。较低的（更严格的）值将增加丢弃的序列的数量，而较高的数字可能意味着您的匹配项质量较低。您可能会发现，如果您有更长的 kmers，您可能需要相应地增加此参数。较低的（更严格的）值将增加丢弃的序列的数量，数字越高可能意味着匹配项的质量越低。‎

表重建

这步里，首先，是区域片段重新比对到完整的数据库序列，然后相对丰度使用一个优化的步骤进行计算，最后生成重建的计数表。per-nucleotide-error参数和maximum-mismatch参数，决定了比对和错配的比例，Ps.我的理解是过程中会产生错误累积。min-abundance决定了要排队的数据库序列的最小丰度，这个参数在某种意义上取决于测序深度和你的偏好。最后，还是可以通过多多线程提升性能(--p-n-workers参数)。下面就是重建命令了，有点长呀：

qiime sidle reconstruct-counts 
 --p-region P1 
  --i-kmer-map ../database/sidle-db-P1-100nt-map.qza 
  --i-regional-alignment P1_align-map.qza 
  --i-regional-table P1_table-100nt.qza 
 --p-region P2 
  --i-kmer-map ../database/sidle-db-P2-100nt-map.qza 
  --i-regional-alignment P2_align-map.qza
  --i-regional-table P2_table-100nt.qza 
   --p-region P3 
  --i-kmer-map ../database/sidle-db-P3-100nt-map.qza 
  --i-regional-alignment P3_align-map.qza 
  --i-regional-table P3_table-100nt.qza 
   --p-region P4 
  --i-kmer-map ../database/sidle-db-P4-100nt-map.qza 
  --i-regional-alignment P4_align-map.qza 
  --i-regional-table P4_table-100nt.qza 
   --p-region P5 
  --i-kmer-map ../database/sidle-db-P5-100nt-map.qza 
  --i-regional-alignment P5_align-map.qza 
  --i-regional-table P5_table-100nt.qza 
     --p-n-workers 8 
 --o-reconstructed-table league_table.qza 
 --o-reconstruction-summary league_summary.qza 
 --o-reconstruction-map league_map.qza

好像这步是超级消耗资源的，我设置了个8核心，资源消耗竟然是＊5的，用上了全部40核心的资源。。。不过几分钟就结束了，那么如果设置1个核心，可能是需要8倍时间的。该命令将生成一个计数表、一个包含有关映射到某个区域的数据库Kmer数量以及ASV ID的详细信息的文件，以及一个需要进行分类重建的映射。

＃可视化看下
qiime feature-table summarize 
 --i-table reconstruction/league_table.qza 
 --o-visualization reconstruction/league_table.qzv

这里我是用了一个样本，133个汇总的feature（也就是以前的OTU），40598条等效数据库序列。您会注意到一些功能ID包含|字符，例如，1764594|195532|4471854。这意味着这两个数据库的序列在重建过程中不能被解析，因此我们将该序列分配给这两个区域。重建中使用的区域越多，就越有可能准确地重建数据库序列。第二个输出是summary。总结可以用来评估重建的质量；更多细节见SMURF文章原稿。默认情况下，摘要会将退化kmer视为唯一序列；您可以使用count-degenerates参数更改行为；如果为false，则仅当kmer属于唯一参考序列时才会对其进行计数。您可以通过将元数据制表来查看摘要。

＃ 汇总摘要
qiime metadata tabulate 
 --m-input-file league_summary.qza 
 --o-visualization league_summary.qzv

几个V区的计数明细和汇总

物种重建

重建特征表之后，需要重建物种，特别的，如果多个数据库序列不能得到处理时。该功能获取重建过程中生成的数据库映射和与数据库相关联的分类，并返回重建的分类。有三种可能，第一，有相同的全部物种分类信息；第二，在一些点上有区别；第三，前面一直相同，直到一个缺失。相同的情况最简单，二取一即可。有些情况下，可能物种分类在高水平上（这里是种）会不一样，例如下面这个，这种情况下会取到属，而不是种。

1236    k__Bacteria; p__Firmictues; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__Blautia; s__obeum
1237    k__Bacteria; p__Firmictues; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__Roseburia; s__

--database 参数允许选择一个数据库进行，如 greengenes, silva 或不提供。如果是greengenes 或 silva，一些专门的数据库清理将会自动进行，特别是缺失和不明处理参数。例如：

k__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Entrobacteriales; f__Enterobacteriaceae; g__; s__

将会整理成：

k__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Entrobacteriales; f__Enterobacteriaceae; g__unsp. f. Enterobacteriaceae; s__unsp. f. Enterobacteriaceae

这里我们用的是greengenes，

qiime sidle reconstruct-taxonomy 
 --i-reconstruction-map league_map.qza 
 --i-taxonomy 。../database/sidle-db-taxonomy.qza 
 --p-database 'greengenes' 
 --p-define-missing 'inherit' 
 --o-reconstructed-taxonomy ./league_taxonomy.qza

进化树重建

最后一步便是重建进化树了，如果不能解析参考序列，则不能简单地从数据库继承系统发育树。因此，我们需要重建一个新的系统树。我们以两种方式处理序列：

• 任何可以完全解析的数据库序列都可以保持其在参考树中的位置。
• 无法解析的序列需要以某种方式进行处理。

我们可以随机选择一个序列来映射重建的区域。然而，当有几个序列组合在一起时，这可能不起作用。因此，如果我们不能解析数据库序列，我们可以从组合的数据中计算一个融合序列，在我们能够绘制的区域中提取它们，然后这些一致的序列可以被插入到使用SEPP或类似的系统发育参考主干中(必须使用相同的数据库版本)。

因此，第一步，对不能解析的序列重建共识序列，这一步官方教程的结果应该有问题，看报错信息进行了更正。

qiime sidle reconstruct-fragment-rep-seqs 
 --p-region P1 
  --i-kmer-map ../database/sidle-db-P1-100nt-map.qza 
 --p-region P2 
  --i-kmer-map ../database/sidle-db-P2-100nt-map.qza 
   --p-region P3 
  --i-kmer-map ../database/sidle-db-P3-100nt-map.qza 
   --p-region P4 
  --i-kmer-map ../database/sidle-db-P4-100nt-map.qza 
   --p-region P5 
  --i-kmer-map ../database/sidle-db-P5-100nt-map.qza 
  --i-reconstruction-map league_map.qza 
  --i-reconstruction-summary league_summary.qza 
  --i-aligned-sequences ../database/sidle-db-aligned-sequences.qza 
  --o-representative-fragments league-rep-seq-fragments.qza

现在，我们可以将序列放入参考树中了，首先下载参考树：

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 
 --i-representative-sequences league-rep-seq-fragments.qza 
 --i-reference-database sepp-refs-gg-13-8.qza 
 --o-tree league-tree.qza 
 --o-placements league-placements.qza

最后，就可以继续常规的qiime2分析流程了，多样性，统计等等，能不能用picrust2分析呢？或许也可以的。

好了，一波三折，终于完成了这个教程的学习，不得不说这个技术好像还不够成熟，特别是qiime2的这个插件。有什么其他好的方法，欢迎交流，我的微信：

前面的16S流程：