宏基因组分析流程 | 全新的 bioBakery 3

本文整理自 Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with bioBakery 3:https://www.biorxiv.org/content/10.1101/2020.11.19.388223v1。本人水平有限,很多时候不能准确表达文章的意思,同时对侧重点做了删减,建议感兴趣的同学阅读原文。

为了进一步扩大微生物群落研究的范围,新版本的 bioBakery 平台中引入了一套全新的计算方法。BioBakery 3 流程包括新版本的数据质控软件 KneadData,分类学注释分析软件 MetaPhlAn 3,功能分析软件 HUMAnN 3,菌株分析软件 StrainPhlAn3 和 PanPhlAn 3 以及系统发育分析软件 PhyloPhlAn 3。同时这些工具几乎都使用了新版本的 ChocoPhlAn 3 泛基因组数据库,该数据库包含系统组织和注释的微生物基因组和基因家族簇。基准测试显示,每个单独的工具都比其先前版本和其他类似算法更准确且有效,有的甚至在灵敏度和特异性上提高了两倍以上。所有组件都已开源,具体请访问 http://segatalab.cibio.unitn.it/tools/biobakery 和 http://huttenhower.sph.harvard.edu/biobakery

bioBakery 提供了完整的宏组学工具套件和分析环境,包括用于宏组学的上游分析方法,下游统计分析,可复现的集成工作流程,储存于开源库(GitHub,Conda,PyPI, 和 R/Bioconductor)的标准化软件包和文档,可网格和云部署的镜像(AWS,GCP 和 Docker),在线学习教程和演示数据以及一个公共讨论社区。对于任何样本集,都可以通过该流程生成质量控制,分类分析,功能分析,菌株分析以及所得的输出数据和报告文本,同时保持版本控制和日志记录。新版本中已更新了所有方法,并进行了相应优化。例如,优化 Docker 镜像大小以方便云部署,并将流程分别移植到 AWS(Amazon Web Services)和 Terra/Cromwell(Google Compute Engine)。所有基本镜像和依赖也已更新,包括最新的 Python(v3.7+)和 R(v4.0+)。

高质量的参考数据集可改善宏基因组分析

bioBakery 3 套件中的大多数软件都利用了全新的参考数据集 ChocoPhlAn 3。ChocoPhlAn 使用公开可用的基因组和标准化的基因和基因家族信息来生成 marker。这些标记将用于 MetaPhlAn 3,StrainPhlAn 3 和 PanPhlAn 3 进行宏基因组的分类和菌株水平分析,PhyloPhlAn 3 进行基因组和 MAGs 的系统发育分析以及 HUMAn 3 进行宏基因组的功能分析。

ChocoPhlAn 3 是一个包含了来自 16.8k 个物种的 99.2k 个高质量,全注释的参考微生物基因组以及对应的 87.3M UniRef90 基因家族注释数据库。通过这些数据,ChocoPhlAn 会为每个物种生成带功能注释的泛基因组。这些泛基因组为随后的整个 bioBakery 3 流程提供了统一的参考注释资源,HUMAnN 3 和 PanPhlAn 3 直接基于完整的泛基因组进行整体功能和菌株谱分析,而其他工具则只使用一些附加信息进行注释,比如,PhyloPhlAn 3 着眼于保守的核心基因家族(即该基因家族存在于一个物种的所有菌株中)以推断准确的系统发育,MetaPhlAn 3 和 StrainPhlAn 3 则将核心基因家族细化为物种特有的独特基因家族,以产生用于宏基因组物种鉴定和菌株水平遗传表征的生物标记。

MetaPhlAn 3 提高了定量分类学分析的准确性

MetaPhlAn 使用进化枝特异性标记基因的覆盖度估算宏基因组中微生物类群的相对丰度。MetaPhlAn 3 整合了 13.5k 个 species(比 MetaPhlAn 2 多出两倍多),并整合了全新的 ChocoPhlAn 3 数据库,其包含从 16.8k 物种泛基因组中选择的一组全新的 1.1M 标记基因(每个 species 平均 84±47)。采用 UniRef90 基因家族可进一步扩展核心基因的鉴定程序,然后将潜在核心基因与所有可用的整个微生物基因组进行比对,以确保这些标记的独特性。此外,核心算法也进行了多项改进和升级,比如标记物比对的质量控制过程和可探究由未知微生物组成的宏基因组含量。

通过 OPAL 基准框架,使用来自 2nd CAMI Challenge 的 113 个合成样本的 118 个合成宏基因组评估了 MetaPhlAn 3 的分类学性能。这些代表了来自五个人体相关身体部位和鼠类肠道的典型微生物群,并为它们补充了另外 五 个新生成的合成的非人类相关的宏基因组。除了 MetaPhlAn 3 外,还比较了 MetaPhlAn 2.7,mOTU 2.51 和 Bracken 2.5。多个评价指标均显示 MetaPhlAn 3 胜过其他算法。

对于 F1 score,MetaPhlAn 3 的表现无论在哪种样本中都优于其他方法。相对于其他工具,MetaPhlAn 3 检出的假阳性菌种数量极少。对于 recall,Bracken 和 mOTU 在某些情况下要优于 MetaPhlAn 3,但有很高的假阳性。MetaPhlAn 通过设置更高比例的阳性物种(--stat_q 参数)可降低假阳率,但总体而言,默认设置下的 F1 值仍然高于其他算法。

宏基因组分析流程 | 全新的 bioBakery 3

除了更准确的物种检测外,在量化物种丰度上 MetaPhlAn 3 也比 MetaPhlAn 2,mOTUs2 和 Bracken 更准确。

宏基因组分析流程 | 全新的 bioBakery 3

在计算效率方面,MetaPhlAn 3 比 MetaPhlAn 2 快 3 倍以上,并且几乎与 Bracken 的速度(每秒11k reads)差不多。MetaPhlAn 3 的内存使用量比 MetaPhlAn 2(2.1Gb)略高(2.6Gb),但优于其他方法(mOTU 为 4.3 Gb,Bracken 为 32.5 Gb)。

宏基因组分析流程 | 全新的 bioBakery 3

HUMAnN 3 准确量化物种对群落功能的贡献

HUMAnN 3 使用来自 ChocoPhlAn 泛基因组的 UniRef90 功能注释数据以进行物种功能分析。同样用上文介绍的 30 个 CAMI 和 5 个其他合成宏基因组比较了 HUMAnN 2 和最近发表的 Carnelian 算法的性能。之所以选择 Carnelian 是因为它是在 HUMAnN 2 之后发表的,同时其算法策略也类似 HUMAnN ,即直接从宏基因组测序数据而不是从组装的重叠群估计分子功能的相对丰度。虽然 HUMAnN 2 和 3 都可以通过宏基因组估计多种功能注释的相对丰度,但选择了4级酶促(EC)注释作为与 Carnelian 比较的基础,因为 Carnelian 为 EC 量化提供了基准指标。

宏基因组分析流程 | 全新的 bioBakery 3

通过 HUMAnN 3 流程,产生了 30 个 CAMI 宏基因组中社区水平 EC 丰度的高准确率估算值(Bray-Curtis 相似性平均值 = 0.93±0.03)。HUMAnN 2 的注释准确率紧随其后为 0.70±0.04,而 Carnelian 的准确率仅 0.49±0.04。虽然 HUMAnN 3 受益于更新的序列数据库,但 HUMAnN 2 的数据库(2014年)也比 Carnelian 早了数年。

与其他功能分析流程相比(比如 Carnelian),HUMAnN 3(和2)的主要优势之一是它们能够根据贡献物种对群落功能概况进行分层。此功能在 HUMAnN 3 中更为准确和有用,因为他有更广泛的泛基因组注释。在 CAMI 数据中,对至少具有 1x 平均深度的物种,在 HUMAnN 3 中 EC 的准确率为 0.81±0.16,而 HUMAnN 2 为 0.51 ± 0.15 。与 HUMAnN 2 相比,HUMAnN 3 倾向于在此覆盖度内检测到更多预期物种。HUMAnN 的物种间功能灵敏度对于样本中 1x 深度以下的物种自然较低,但对每个物种的准确率仍然很高,并且新版本中,v3 与 v2 相比略有提高 (0.95 ± 0.08 vs. 0.91 ± 0.07)。

Carnelian 是这三种方法中计算效率最高的方法,分析 CAMI 的时间为 26.4±2.7 CPU 小时,而 HUMAnN 2 为 38.1±12.8 CPU小时,HUMAnN 3 为52.5±19.2 CPU小时。峰值内存使用量(MaxRSS)的趋势与此类似,其中 Carnelian 要求 11.9±0.0 GB,而 HUMAnN 2 为 17.0±0.3 GB,HUMAnN 3 为 21.5±1.9 GB。这些差异主要归因于算法检索的数据库大小:尽管 Carnelian 仅关注 EC 可注释的序列子集,但 HUMAnN 的目标是首先量化数千万个唯一的 UniRef90,其中最终只有 12.5% 是由 EC 注释。与 HUMAnN 2 相比,HUMAnN 3 运行时间的增加同样归因于其更大的翻译检索数据库(UniRef90 87.3M v.s. 23.9M),因为即使在前面的核苷酸水平检索层中已解释了大多数数据,翻译后的检索层也是 HUMAnN 算法的速率限制步骤。这种现象也解释了 HUMAnN 的运行时间差异,因为运行时间与翻译后的检索层数据成反比。值得注意的是,当翻译检索步骤,HUMAnN 3 可以在 5.8±0.8 CPU小时中解释大多数 CAMI 宏基因组读数(每个样本70.9±9.6%),当然这通常仅适用于已知被相关参考序列充分覆盖的菌群。

使用一组富含非人类相关物种的合成数据进行评估,可以看到这三种方法之间的相对准确率和运行效率趋势相似。因此,HUMAnN 3 的强大性能不仅限于由宿主相关物种组装而成的微生物群落。此外,相对于 HUMAnN 2,MetaPhlAn 3 对非宿主相关物种的灵敏度提高了,从而提高了 HUMAnN 3 的准确性和性能(通过更快更准确的泛基因组检索步骤注释了绝大多数 reads)。

宏基因组分析流程 | 全新的 bioBakery 3

算法原理

bioBakery workflow

GitHub: https://github.com/biobakery/biobakery

bioBakery3 是一套从宏组学数据中分析微生物群落的计算方法,可产生分类,功能,系统发育和菌株水平概况,以直接解释或包括在下游统计分析中。

宏基因组分析流程 | 全新的 bioBakery 3

用 KneadData 进行 reads 质量控制后,MetaPhlAn3 将估计样品中存在的微生物种类及其相对丰度。进一步,StrainPhlAn3 通过完善由 MetaPhlAn3 鉴定的菌种的菌株水平基因型来加深遗传学表征。而 HUMAnN3 则更侧重于鉴定和定量在宏基因组中编码或在宏转录组中表达的分子功能,PanPhlAn3 可以解析其菌株的基因型。PhyloPhlAn3 提供了一种全面的方法来解释基于组装算法的基因组草图。这些 bioBakery3 模块都基于 ChocoPhlAn3 数据集,目前共包括 99,227 个基因组和 87.3M 个基因家族,这几乎是 bioBakery 第一个发行版所包含数据的 100 倍。

用 AnADAMA 2 串联流程

GitHub: https://github.com/biobakery/anadama2

bioBakery3 分析流程使用 AnADAMA v2 进行串联。简单来说,AnADAMA 打包了 Python 的依赖管理器 doit( http://pydoit.org ),该模块为分析任务定义,版本和输出跟踪,变更管理,文档,网格,云部署和自动报告提供了一种简单但可扩展的方式。如果修改了工作流程或更改了输入文件,则只会重新运行那些受影响的任务。使用默认的记录器和命令行报告器记录所有任务的基本信息,确保流程的可重复性。记录的信息包括流程的各项命令和参数,可执行文件的版本以及每个任务的输出文件。同时也可用于将其他 bioBakery3 模块连接在一起。

用 KneadData 进行 reads 质控

GitHub: https://github.com/biobakery/kneaddata

bioBakery 3 包括一个原始数据质控模块 KneadData,该模块可自动化清洗原始宏基因组和宏转录组数据。主要内容包括:

使用 Trimmomatic 去除1)低质量的碱基(默认值:4-mer 窗口,平均 Phred 质量 <20)2)截短的 reads(默认值:<50% 的预修剪长度)3)adapter 和 barcode使用 FastQC 去除过表达序列(默认为 > 0.1% 频率),使用 TRF 去除低复杂度的序列。但其实软件默认情况下并没有加上这个参数 --run-trim-repetitive,因为扩增子序列通常具有大量重复 reads,会干扰结果。但强烈建议在宏基因组流程中加上该参数。使用 bowtie2 比对并去除参考基因组序列。软件自带一个扩展的人类参考基因组(包含 hg37,已知的污染序列 和 decoy 序列)。

有研究(Human contamination in bacterial genomes has created thousands of spurious proteins)对 NCBI RefSeq 数据库中的细菌和古细菌基因组进行了大规模扫描,结果发现有 2250 个基因组被人类序列污染。污染物序列主要来自高拷贝的重复区域。在某些情况下,污染的重叠群被错误地注释为包含蛋白质编码序列,该序列随着时间的流逝已传播到多个原核和真核基因组中,从而形成了伪蛋白“家族”。

在转录组数据中使用 SILVA 数据库去除微生物核糖体和结构 RNA 序列。

建议在进行后续分析前使用 KneadData 进行质控,bioBakery workflows 工作流程会默认对所有类型的数据执行此操作。

构建 ChocoPhlAn 3 泛基因组数据库

构建 ChocoPhlAn 数据库的目的就是根据物种分类水平来组织微生物参考基因组,并整合相关序列和注释数据。在检索 UniProt 基因组和基因注释后,使用所有经过初始质控的微生物参考基因组来生成特定物种的泛基因组(比如,物种的基因家族需存在于其至少一个基因组中)。然后从整个泛基因组库中鉴定出核心基因组(即,一个物种的所有基因组中都存在的基因家族),并将其用作 PhyloPhlAn 3 中的标记。进一步处理核心基因组,以提取独特的标记基因(即,与一个物种唯一相关的核心基因家族 ),构成 MetaPhlAn 3 和 StrainPhlAn 3 的标记数据库。最后,对带有功能注释的泛基因组进行整理,构成 PanPhlAn 3 和 HUMAnN 3 的参考数据库。

泛基因组:某一物种全部基因的总称,包括核心基因组(core genome),在所有菌株中都存在的基因;以及非必须基因组(dispensable genome), 在部分菌株中存在的基因。在实际研究中,有时候,泛基因组也可以分成核心基因组(在所有菌株中都存在的基因)、非必须基因组(在2个以及2个以上的菌株中存在的基因),以及菌株特有基因(strains-specific gene,仅在某一个菌株中存在的基因)。

数据检索

ChocoPhlAn 依赖于 UniProt 核心数据资源(release January 2019)和 NCBI 分类学和基因组存储库(release January 2019)。ChocoPhlAn 共纳入两种基本序列数据类型,即微生物的原始基因组以及这些基因组上鉴定的微生物蛋白质/基因数据。ChocoPhlAn 用微生物分类整合基因组数据,用蛋白家族整合微生物蛋白质数据。

ChocoPhlAn 采用 NCBI 分类数据库进行分类注释。完整数据可从 NCBI FTP 服务器 ftp.ncbi.nlm.nih.gov/pub/taxonomy/ 下载。ChocoPhlAn 使用 unidentifiedsp.Candidatusbacterium 来标记 species。也有一些关键词作为低质量的菌种的标志,具体来说,用于过滤低质量分类注释的正则表达式如下所示:

(C|c)andidat(e|us)| _sp(_.*|$)|(.*_|^)(b|B)acterium(_.*|)|.*(eury|)archaeo(n_|te|n$).*|.*(endo|)symbiont.*|.*genomosp_.*|.*unidentified.*|.*_bacteria_.*|.*_taxon_.*|.*_et_al_.*|.*_and_.*|.*(cyano|proteo|actino)bacterium_.*)

根据 UniProt 蛋白质组资源中的 GCA 编号从 ftp.ncbi.nlm.nih.gov/genomes/genbank/assembly_summary_genbank.txt 文件中获取对应 URL,再从 NCBI GenBank FTP 服务器 ftp.ncbi.nlm.nih.gov/genomes/all/GCA 下载每个蛋白质组的参考基因组( fasta 格式,后缀 .fna)和相关的基因组注释(.gff)。丢弃了 12,598 个缺少 GenBank 登录号的蛋白质组,最终得到 99,227 个基因组(997 古细菌,97,941 细菌,339 真核生物)。

从 UniProt Knowledgebase(UniProtKB)和 UniProt Archive(UniParc)数据库中检索与至少一个 UniProt 蛋白质组相关并包含在 ChocoPhlAn 数据库中的微生物蛋白质(和基因)。UniProtKB 中包含的蛋白质已从 UniProt 蛋白质组中所有可用参考基因组的 CDS 的翻译衍生而来。ChocoPhlAn 3 还检索并包含了 UniProtKB 条目中存在的相关数据(从 ftp.uniprot.org/pub/databases/uniprot/ 下载 uniprot_sprot.xml.gzuniprot_trembl.xml.gz 和 uniparc_all.xml.gz),比如,功能,系统发育和蛋白质结构域注释(KEGG,KO,EggNOG,GO,EC,Pfam),与外部数据库(GenBank,ENA和BioCyc)进行交叉引用条目的登录名和该基因的名称。

共处理了 UniProtKB 和 UniParc 中 203.9M 蛋白质,其中 126.9M 的蛋白质与 UniProt 蛋白质组条目相关。Bacteria domain 的蛋白质数量最多(194.8M),而古细菌和真核生物分别占 5.0M 和 4.0M 蛋白质。

为了减少数据库冗余,使用 UniProt 提供的 UniProtKB UniRef90。UniRef90 是通过将 UniRef100(将相同的 UniProtKB 蛋白结合在一个簇中)用 CD-HIT 聚类得到的(2019 年 8月后换用 MMseqs2 进行聚类)。UniRef90 中至少有 90% 的序列同一性。UniRef50 是通过将 UniRef90 种子序列进行聚类而生成的,每个簇都包含具有至少 50% 同一性的蛋白质。UniRef90 和 UniRef50 都要求每种蛋白质与簇的最长序列重叠至少80%。在 ChocoPhlAn 3 中考虑的 UniRef 条目包含一个代表性蛋白质的序列,该簇中包括的所有条目的登录 ID,UniProtKB 和 UniParc 记录的登录以及其他关联的 UniRef 簇的登录号都包含在 UniProt 条目中。

共处理了 292.1M 个 UniRef cluster(UniRef100,UniRef90和UniRef50分别为172.3M,87.3M和32.5M),并与 ChocoPhlAn 3 中的每种蛋白质和每个基因组相关联。

生成泛蛋白质组

接着为至少由一个 UniProt 蛋白质组表示的物种生成他们的泛蛋白质组,即该物种具有蛋白编码潜力的非冗余序列。通过整合基因组中独特的 UniRef90 和 UniRef50 蛋白家族,可以得到每个物种的蛋白数据。

对于每种泛蛋白,将计算几个分数:

核心评分 coreness:该物种的泛蛋白质组中包含的基因组数目唯一性评分 uniqueness:其他具有相同泛蛋白质的物种的泛蛋白质组数目uniqueness_sp 得分:唯一性评分的一种变换形式,即不包括那些以前被标记为低质量物种的物种external_genomes:拥有相同泛蛋白的其他物种的“泛蛋白质组”的基因组数目(而不是物种或物种的泛蛋白质组)。

这些分数是针对 UniRef50 和 UniRef90 蛋白家族计算的。

在 ChocoPhlAn 3 中,共纳入了 22,096 种泛蛋白质组和总共 87.3M UniRef90 核心蛋白。

生成 MetaPhlAn 3 标记

MetaPhlAn 依赖于一套独特的且特定于物种的核苷酸标记。首先使用 species-level genome bin(SGB)过滤掉了先前被标记为低质量分类的物种。归入同一 SGB 的低质量物种被合并,仅考虑代表 SGB,这里共合并了 1,328 个物种(6%),因为它们可能导致假阳性的分类学注释。对于 NCBI 分类中包含在“species-group”中的多个物种显示出大量具有较高唯一性评分(>30)的标记,主要包括下面这些菌种群:Streptococcus anginosus group,Lactobacillus casei group,Bacillus subtilis group,Enterobacter cloacaecomplex,Pseudomonassyringae group,Pseudomonas stutzeri group,Pseudomonas putida group,Pseudomonasfluorescens group,Pseudomonas aeruginosa group,Streptococcus dysgalactiae group 以及 Bacillus cereus group。在这些情况下,泛基因组都是通过合并所有物种级别的泛基因组并将它们视为一个物种来构建的。

第一步,使用 UniRef90 中长度在 150 至 1500 个氨基酸之间的蛋白质构建泛蛋白组数据库。从核心评分和唯一性评分开始,应用迭代的方法来查找标记,最多保留 150 个标记,并仅保留至少有 10 种标记的物种。根据“唯一性”评分将候选标记分为唯一标记和准标记:“唯一性”为零的标记被报告为“唯一标记”。如果无法识别出唯一的标记,则将使用更宽松阈值得到所谓的“准标记”,即“唯一性”评分大于零。

迭代的方法主要是结合“核心”,“唯一性”和“ external_genomes”评分来定义四个标记质量等级。

等级 A:核心评分超过 80%,不超过与 2 个其他泛蛋白质组共享(Uniqueness_NR90 和Uniqueness_NR50),并且出现在不超过 10 个基因组中(UniRef90,在 UniRef50 中为 5)(即 External_genomes_NR90External_genomes_NR50等级 B:核心评分介于 70% 和 80% 之间,Uniqueness_NR90 Uniqueness_NR50值为 5 , External_genomes_NR90 External_genomes_NR50值分别低于 15 和 10 个基因组等级 C:不符合先前条件的标记将包含在 C 层中,其中包括核心评分在 50% 和 70% 之间,Uniqueness_NR90小于10,Uniqueness_NR50小于15, External_genomes_NR90小于25,External_genomes_NR50小于20。等级 U:对于泛蛋白质组中仅包含一个基因组的物种,其核心定义微不足道的标记被分类为“ U”级,它们的唯一性评分为零。

对迭代过程使用如下定义的评分对满足阈值的候选标记进行排序:

宏基因组分析流程 | 全新的 bioBakery 3

该评分将优先选择进化枝中高度保守(核心评分高)但与其他物种的数量尽可能少(唯一性评分低)的候选标记。为每个候选标记注释其等级,如果识别出 50 个以上的候选标记,那么将从排名列表中选择至多 150 个标记。如果没有找到足够的标记(少于50个),则使用下一层的阈值重复该过程。如果使用 C 级阈值仍没有找到符合的标记,则将物种丢弃。

然后,将通过此过程选择的每个标记的核苷酸序列作为 MetaPhlAn 数据库的条目。为了完善通过“唯一性”参数估算的物种数量,将标记序列拆分为150bp 的非重叠块,并使用 bowtie2( version 2.3.4.3,参数 -a --very-sensitive --no-unal --no-hq --no-sq)进行比对。如果在识别出的目标参考基因组中发现了至少 150 个标记序列的连续核苷酸,将基于“唯一性”参数注释一个新识别的物种。

对由 Co-Abundance gene group(CAG)获得的基因组物种的标记物进行了额外的质控步骤。为了减少假阳性的数量,如果与提供 CAG 基因组分类的物种共享超过 50% 的标记,就删除该 CAG 物种。

每个标记都在 MetaPhlAn 数据库中关联了一个条目,该数据库包括标记的物种,共享该标记的物种列表,序列长度和物种分类。病毒标记继承自 v20_m200 MetaPhlAn2 数据库。

共确定了 13,475 个物种的 1.1M 标记。

用 MetaPhlAn 3 进行物种注释

GitHub: https://github.com/biobakery/MetaPhlAn

质控后的宏基因组 reads 由 MetaPhlAn 3 调用 bowtie2 比对到 1.1M 标记基因数据库中。默认的 bowtie2 参数是 very-sensitive ,但可通过 MetaPhlAn 3 进行自定义。在 MetaPhlAn 3 中,输入文件可以是单个 FASTQ 文件,也可以是压缩文件亦或事先比对好的数据。程序内部,MetaPhlAn 3 会估算每个标记的覆盖度,并将进化枝的覆盖度计算为同一进化枝的各个标记覆盖度的平均值。然后对进化枝的覆盖度进行标准化,以获得每个分类单元的相对丰度。

在 MetaPhlAn 3 中优化了 robust average 的参数,这将剔除位于顶部和底部分位数的标记丰度(stat_q 参数),默认为 0.2(即不包括丰度最高的 20% 的标记以及丰度最低的 20% 的标记)。为了提高比对质量,MetaPhlAn 3 在比对之前和之后丢弃了低质量的序列和比对结果( reads 短于 70bp 且 MAPQ 值小于 5 的比对)。

MetaPhlAn 3 还引入了一项新功能,可以估计当前数据库中不存在的未知分类单元。这是通过从 reads 总数中减去根据分类单元平均基因组长度归一化的每个分类单元的平均 reads 深度来计算的。此外,默认情况下,MetaPhlAn 3 的新输出格式包括每个进化枝的 NCBI 分类标准 ID,从而重新注释的情况下,可以更好地比较工具和跟踪物种名称。

最后,除了默认的 MetaPhlAn 输出格式外,现在还可以使用 CAMI 输出格式报告配置文件,该格式可用于执行 OPAL 框架的基准测试。为了支持后期分析,可从软件仓库(metaphlan/utils/calculate_unifrac.R)中获得 R 脚本,用于计算 MetaPhlAn 结果中的 weighted 和 unweighted UniFrac 距离,以及系统发育分析。

用 StrainPhlAn 3 进行菌种分析

StrainPhlAn 通过重建 MetaPhlAn markers 的样本特异性序列并将其用于多序列比对和系统发育建模从而在菌株水平上进行基因分型。StrainPhlAn 3 也在多个方面改进了之前的版本:

1.集成了经过改进和验证的共识序列生成流程2.PhyloPhlAn 3 的集成,从而提高了系统发育建模的质量和分析的灵活性3.改进了算法,可以过滤没有足够标记支持的样本以及在菌株和样本之间不够保守的标记。

StrainPhlAn 3 将 MetaPhlAn 3 的结果和 MetaPhlAn 3 数据库作为输入。对于每个样品,StrainPhlAn 3 通过在标记的每个位置观察比对到标记并覆盖该位置的 reads 中频率最高的核苷酸,来重建物种特异性标记的高质量共有序列。默认会丢弃少于 8 个 reads 或覆盖范围(即 reads 所覆盖的标记部分)小于 80% 的共有标记( --breadth_threshold 参数)。Ambiguous base 定义为比对中质量低于 30 或多态性较高(major allele dominance 低于 80%)的位置,并结合覆盖范围的阈值来判断为 unmapped positions。

重建标记后,过滤算法将丢弃少于 20 个标记的样本以及在少于 80% 的样本中存在的标记(分别为 --sample_with_n_markers 和 --marker_in_n_samples 参数)。然后,通过删除前导和尾随的 50 个碱基(--trim_sequences 参数)来修剪标记,由于比对过程中的边界效应,这些标记通常覆盖度较低。最后,通过 PhyloPhlAn 3 处理过滤的样品和标记,以进行系统发育分析。默认情况下,使用 BLASTn 将重建的序列比对到标记数据库上,通过 MAFFT 执行多序列比对,并通过 RAxML 生成系统发育树。为了重建菌株水平的系统发育,将 PhyloPhlAn 设置为使用 --diversity low 参数运行。

StrainPhlan 3 产生的系统发育树也可用于鉴定样品中的相同菌株,例如可以在菌株传播分析中加以利用,新添加的 strain_transmission.py 脚本对此提供了支持。该脚本处理由 StrainPhAn 生成的系统发育树以及描述样本之间的元数据以推断菌株传播事件(例如纵向样本或具有相关关系的样本(例如母婴)的样本之间的关系)。首先,使用系统发育树,生成成对距离矩阵,并通过树的总分支长度对其进行归一化。使用距离矩阵和相关的元数据,可以推断出定义相同菌株的阈值,从而选择非相关样本距离的分布的第一个百分位数(即将理论错误发现率的上限设为 1%)。如果提供了纵向样本,则每个对象仅考虑一个样本,并且元数据中未包括的样本被视为不相关。最后,距离小于推断阈值的相关样本配对将被报告为潜在的传播事件。

用 HUMAnN 3 进行功能注释

GitHub: https://github.com/biobakery/humann

HUMAnN 使用包含 UniRef90 注释的泛基因组数据库对 MetaPhlAn 中的每个样品检测到的所有物种进行微生物群落的功能分析。为了获得每个泛蛋白质组的核苷酸代表序列,通过选择在分类上分配给所需物种的 UniProtKB 或 UniParc 条目,为每种泛蛋白确定了代表簇。每个代表簇用于从来源参考基因组中提取核苷酸序列以及来自不同系统的一些功能注释,包括 GO,KEGG,Pfam,EC 和 eggNOG。除了功能注释外,数据库还将每个 UniRef90 簇与其对应的 UniRef50 簇相关联,以提供多个级别的注释分辨率。

HUMAnN 3 实现了许多新选项以微调分层检索步骤(例如将自定义搜索参数分别传递给 bowtie2 和 DIAMOND)。在将 reads 分配给物种泛基因组时,首先加入了两个新选项:

1.覆盖度过滤。HUMAnN 2 中基于覆盖度过滤的选项已用于对翻译后的搜索结果进行质控,结果表明这样可以显著提高 UniRef90 水平的特异性,并对灵敏度的影响很小。在泛基因组检索的背景下,也观察到了相似的结果。所以 HUMAnN 3 在其基因组和翻译的检索步骤中都设置了(可单独调整)数据库序列覆盖度过滤器(默认均为 50%)。2.允许一个 read 与多个泛基因比对,而不是 bowtie2 的默认设置(在 HUMAnN 2 中使用)仅比对到单个目标。但实验表明,允许 reads 最多比对到 5个泛基因(通过 bowtie2 的 -k 5 设置实现的新选项2)对准确性的影响很小,所以这在 HUMAnN 3 中未默认启用。

其实也考虑了在翻译检索过程中调整 DIAMOND 0.9 的严格性和内存使用情况的新选项。其中影响最大的是将 UniRef90 的每个 read 比对的识别阈值从90%(HUMAnN 2 默认值)降低到 80%(HUMAnN 3 新默认值)。虽然选择前一个值是为了尊重 UniRef90 家族成员之间的平均一致性,但 80%的阈值更能包容滑动窗口内的变化。80% 的阈值与 HUMAnN 的数据库序列覆盖度过滤器配合使用,可以在翻译搜索过程中正确比对到更多的 reads,而不会影响特异性。

虽然 HUMAnN 2 接受 DIAMOND(默认)每次查询的前 20 个数据库目标,但重新评估前 1 个和前 5 个目标,以及在最佳匹配分数的 1%,2% 或 10% 之内的目标。根据 UniRef90 特异性的显着提高和最小的灵敏度损失,选择“最佳匹配得分低于 1%”(DIAMOND 的 --top 1 选项)作为 HUMAnN 3 的新默认设置。最后,软件探讨了通过 --block-size(-b) 和 --index-chunks(-c) 参数来调整 DIAMOND 的内存。结果表明,大内存对速度的提升并不明显,因此 HUMAnN 3 将继续支持 DIAMOND 的默认低内存配置。

用 PanPhlAn 3 进行功能注释和菌株谱分析

GitHub: https://github.com/SegataLab/panphlan

PanPhlAn 通过识别单个宏基因组样本中特定物种的基因库组成来进行菌株水平的宏基因组分析。它使用 bowtie2 将宏基因组与目标物种的泛基因组相对应。覆盖度归一化后(通过将一个基因家族中所有基因的基因覆盖度相加,再除以该家族的平均基因长度),PanPhlAn 会在每个样本中建立一个基因家族的覆盖度曲线,并评估这些基因家族中存在与否。这将产生一个基因家族存在/缺失的二元矩阵。

与以前的版本相比,在 PanPhlAn 3 中采用了一个新的 ChocoPhlAn 3 泛基因组数据库,其中包含 MetaPhlAn 3 中的 2,298 个物种,他们均至少有 2 个参考基因组可用。对于具有超过 200 个参考基因组的物种,使用 200 个基因组的代表性子集来建立泛基因组,以最大化它们之间的 Mash 距离。PanPhlAn 数据库的基因组由所有重叠群的 FASTA 文件,预先建好的 bowtie2 索引和制表符分隔的文件组成,该文件包含基因家族的 UniRef90 ID 以及重叠群上的基因名称,基因组坐标,以及功能和结构注释。

此外,新功能还包括一个脚本,该脚本可快速显示存在/缺失矩阵,并具有将样本中基因家族的概况聚类的功能。可以基于一组基因的长度之和与沿重叠群的总跨度之比,为每个聚类计算经验 p 值。因此,可以确定一个显着的“紧密”基因组,并计算经验 p 值,以评估这些家族沿着重叠群的遗传邻近性是否被认为是重要的。这简化了宏基因组样本中移动元素的检测和鉴定。

用 PhyloPhlAn 3 进行系统发育分析

GitHub: https://github.com/biobakery/phylophlan

PhyloPhlAn 3是一种易于使用的方法,可以对微生物基因组和宏基因组组装的基因组(MAG)进行分类分析和系统发育分析。PhyloPhlAn 数据库利用了ChocoPhlAn 3 鉴定的核心基因集和参考基因组集,并从 111,825 种 UniProt 蛋白质组学中提取出每个分类学物种。简而言之,PhyloPhlAn 3 数据库中包含的核心基因用于鉴定输入基因组和 MAG 中的序列同源物,然后进行比对,连接并用于系统发育重建。还可以包括先前分析的一组 MAGs 以提供新组装的MAGs的系统进化语境。因此,PhyloPhlAn 3提供了将基于装配的方法和系统发育分析集成到bioBakery 3分析框架中的方法。

生物医学科研方法

【生信文献200篇】06 高通量药物筛选 + Rhabdoid tumors

2021-1-18 9:14:44

生物医学科研方法

这本4分+老牌2区SCI一审最快仅12天!

2020-12-19 17:25:15