前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >SMURF流程之q2-sidle(二)-- 序列重建

SMURF流程之q2-sidle(二)-- 序列重建

作者头像
用户1075469
发布2021-03-11 14:54:09
3980
发布2021-03-11 14:54:09
举报
文章被收录于专栏:科技记者

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

区域比对

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

代码语言:javascript
复制
# 继续循环解决啦
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参数)。下面就是重建命令了,有点长呀:

代码语言:javascript
复制
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的详细信息的文件,以及一个需要进行分类重建的映射。

代码语言:javascript
复制
#可视化看下
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属于唯一参考序列时才会对其进行计数。您可以通过将元数据制表来查看摘要。

代码语言:javascript
复制
# 汇总摘要
qiime metadata tabulate \
 --m-input-file league_summary.qza \
 --o-visualization league_summary.qzv

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

物种重建

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

代码语言:javascript
复制
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,一些专门的数据库清理将会自动进行,特别是缺失和不明处理参数。例如:

代码语言:javascript
复制
k__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Entrobacteriales; f__Enterobacteriaceae; g__; s__

将会整理成:

代码语言:javascript
复制
k__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Entrobacteriales; f__Enterobacteriaceae; g__unsp. f. Enterobacteriaceae; s__unsp. f. Enterobacteriaceae

这里我们用的是greengenes,

代码语言:javascript
复制
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或类似的系统发育参考主干中(必须使用相同的数据库版本)。

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

代码语言:javascript
复制
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

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

代码语言:javascript
复制
wget  -O "sepp-refs-gg-13-8.qza" \
 "https://data.qiime2.org/2020.11/common/sepp-refs-gg-13-8.qza"

然后,插入序列:

代码语言:javascript
复制
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的这个插件。有什么其他好的方法,欢迎交流,我的微信:

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2021-02-06,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 科技记者 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 区域比对
  • 表重建
  • 物种重建
  • 进化树重建
相关产品与服务
数据库
云数据库为企业提供了完善的关系型数据库、非关系型数据库、分析型数据库和数据库生态工具。您可以通过产品选择和组合搭建,轻松实现高可靠、高可用性、高性能等数据库需求。云数据库服务也可大幅减少您的运维工作量,更专注于业务发展,让企业一站式享受数据上云及分布式架构的技术红利!
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档