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流程:

16S流程知多少
Nanopore 16S测序数据分析之last和centrifuge进行物种注释
使用nf-core的ampliseq(qiime2)流程分析16S数据
纳米孔Nanopore-16S数据分析学习笔记


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

发表评论