anvi'o泛基因组分析流程[官方版本(翻译)]

写在前面:鉴于本人想用这个软件分析泛基因组,但是并没有找到合适的中文说明,那就只好翻英文说明,并顺手翻译了一下。本文档主要根据官方原泛基因组分析流程直接翻译过来的。方便各位快速了解如何使用anvio进行泛基因组分析,并了解一些基本的原理。在重要的地方我进行一定的加粗处理。官方文档中有提到一些测试数据大家可以自行下载测试。文档内容较多,可能需要花费一些时间进行阅读,作者对各个功能讲解的非常详细,也给出了一些详细的例子,属于手把手教的那种了。但是有些地方可能会稍微有点啰嗦,大家可以选择性跳过一些内容。其实有些关于原理性的东西大家看不懂也不是很要紧。知道怎么用就行了!
本人英语水平一般,本文档翻译主要还是依赖翻译软件来的,有明显错误的地方我已经纠正了。如果有翻译不对的地方还请见谅。

这边有一个实例,本人已经使用多组数据顺利跑通程序anvio泛基因组分析 - 简书 (jianshu.com)

使用anvio进行泛基因组分析可以完成以下工作:

  1. 若干fasta文件,识别基因组中感兴趣的基因类群
  2. 将来自 anvi’o 项目基因组的宏基因组组装基因组(MAGs)与其他来源的基因组相结合; 进行同步分析;
  3. 互动式可视化泛基因组结果
  4. 基于基因聚类结果评估基因组的关系
  5. 以交互方式(或以编程方式)将基因簇划分为箱(bins)
  6. 针对指定的基因簇进系统发育分析
  7. 注释基因,并检查基因簇内的氨基酸比对情况
  8. 利用有关基因组或基因簇的上下游信息扩展全基因组
    9.定量分析基因簇的结构和生化亲缘关系
  9. 对泛基因组中的基因组进行功能富集分析
  10. 计算并可视化基因组之间的平均核苷酸同一性得分等

说明:使用anvio进行泛基因组分析输入文件仅需要fastaa文件(笔者一般直接使用faa文件进行分析)

介绍

本流程旨在引导用户完成以下流程
· 通过 anvi-gen-genomes-storage 程序讲基因组转化为anvio可识别的 *genomes-storage.db 格式
· 通过anvi-pan-genome构建泛基因组专用的db格式 pan.db
· 如何通过anvi-display-pan 命令唤出交互界面

然后,您可以使用交互式界面将基因簇放入集合中,或者使用程序 anvi-import-collection 导入基因簇的 bin,最后您可以使用程序 anvi-summarize 创建您的基因簇的静态 HTML 摘要。泛基因组。十分简单!

以下部分将详细介绍每个步骤,最后将是一个示例运行,但我们首先确保您已安装所有必需的依赖项,并且安装顺利进行。

依赖

在分析开始前请进行程序的自检,确保安装无误:

anvi-self-test --suite pangenomics

生成anvio特定的基因组存储格式(*genomes-storage.db)

genomes-storage.db是anvio存储基因组信息的特殊格式。这种数据格式可以通过外部基因组(external-genomes)或者已有的基因组(internal-genomes)生成;又或者两者结合。

· 外部基因组(external-genomes)是 FASTA 文件格式的任何内容(即从 NCBI 下载的基因组,或通过任何其他方式获得的基因组)。这意味着,您需要首先使用程序 anvi-gen-contigs-database 将每个 FASTA 文件转换为 anvi’o contigs 数据库。请阅读 contigs-db 工件,以确保您使用最有用的信息填充您的 contigs 数据库(例如使用功能注释您的基因等)
· 内部基因组(internal-genomes)是您存储在 anvi’o 集合中的任何基因组箱(在对 anvi’o 中的宏基因组进行分箱和/或精炼后)。使用内部基因组非常简单,因为您已经有一个 anvi'o 重叠群和一个 anvi'o 配置文件数据库,但如果您正在阅读本教程并且这对您来说还没有意义,请不要担心。

简而言之,外部基因组是指直接从fasta文件导入开始(从零开始),内部基因组指的是使用anvio进行分析过程中产生的db文件中抽提出来(anvio除了可以进行宏基因组的组装分箱等功能,可以理解为抽提在这些过程产生的MAGs出来)。

可以使用 anvi-gen-genomes-storage 程序创建新的 anvi’o 基因组存储,这将要求您提供要包含在该存储中的基因组的描述。外部基因组和内部基因组描述的文件格式略有不同。例如,这是一个示例 --external-genomes 文件:

示例数据

这是 --internal-genomes 的示例文件:


对于第一列中的名称,请仅使用字母、数字和下划线字符。

借助这两个文件,anvi'o 集合中描述的基因组和通过其他来源获得的基因组可以无缝组合和分析。基因组存储中一致性的最重要需求是确保每个内部和外部基因组在基因如何调用、功能如何分配等方面都是相同的。Anvi'o 会检查大多数事情,但它可以不能阻止你犯错误。例如,如果识别开放阅读框的基因调用程序在所有重叠群数据库中都不相同,则基因组存储中描述的基因不一定具有可比性。

专门为教程读者提供的真实示例:原绿球藻全基因组
为了可重复性,本教程的其余部分将遵循真实的示例。

我们将简单地创建本研究中使用的 31 个原绿球藻分离基因组的泛基因组

如果您想要在本地复现,您可以下载原绿球藻数据包 (doi:10.6084/m9.figshare.6318833 ),其中包括这些分离基因组的 anvi’o 重叠群数据库:

wget https://ndownloader.figshare.com/files/28834476 -O Prochlorococcus_31_genomes.tar.gz
tar -zxvf Prochlorococcus_31_genomes.tar.gz
cd Prochlorococcus_31_genomes
anvi-migrate *.db --migrate-dbs-safely

该目录包含 anvi’o contigs 数据库、一个外部基因组文件和一个 TAB 分隔的数据文件,其中包含每个基因组的附加信息(这是可选的,但稍后您将看到为什么它非常有用)。您可以通过以下方式生成本节中所述的基因组存储:

anvi-gen-genomes-storage -e external-genomes.txt \
                         -o PROCHLORO-GENOMES.db

进行泛基因组分析

基因组格式转换后,就可以使用程序 anvi-pan-genome 来运行实际的泛基因组分析。这是该命令最简单的形式:

 anvi-pan-genome -g MY-GENOMES.db -n PROJECT_NAME

原绿球菌泛基因组分析 [续]
让我们使用我们使用 31 个原绿球藻分离株创建的基因组存储来运行泛基因组:

anvi-pan-genome -g PROCHLORO-GENOMES.db \
                --project-name "Prochlorococcus_Pan" \
                --output-dir PROCHLORO \
                --num-threads 6 \
                --minbit 0.5 \
                --mcl-inflation 10 \
                --use-ncbi-blast

--project-name 之后的每个参数都是可选的(这里我们保持与文章中运行泛基因组的参数一致)。

下载的目录还包含一个名为“layer-additional-data.txt”的文件,它总结了每个基因组所属的进化枝。计算出泛基因组后,我们可以使用 anvi’o 程序 anvi-import-misc-data 将其添加到泛数据库中:

anvi-import-misc-data layer-additional-data.txt \
                      -p PROCHLORO/Prochlorococcus_Pan-PAN.db \
                      --target-data-table layers

输出

New layers additional data...
===============================================
Data key "clade" .............................: Predicted type: str
Data key "light" .............................: Predicted type: str

New data added to the db for your layers .....: clade, light.

我们都在寻找丰富泛基因组展示的方法,而 anvi’o 的附加数据表是实现这一目标的绝佳方法。请花一点时间在这里了解有关它们的更多信息。

当执行anvi-pan-genome后,程序执行以下步骤

  1. 使用基因组存储中的所有基因组。如果想关注特定的子集,可以使用参数 --genome-names

  2. 默认情况下仅使用单个核心。根据您正在分析的基因组数量,此过程可能非常耗时,因此您应该考虑通过参数 --num-threads 增加要使用的线程数量

  3. 默认情况下,在“快速”模式下使用 DIAMOND (Buchnfink et al., 2015)(或者您可以通过使用参数 --sensitive 来要求 DIAMOND 为“敏感”)来计算每个基因组中每个氨基酸序列与每个基因组的相似性。所有基因组中的其他氨基酸序列(这显然需要您安装 DIAMOND)。或者,您可以使用标志 --use-ncbi-blast 来使用 NCBI 的blastp 进行氨基酸序列相似性搜索

  4. 调用所有基因,无论它们是否完整。尽管这对于完整基因组来说不是一个大问题,但宏基因组组装基因组 (MAG) 在重叠群的末尾和开头会有许多不完整的基因调用。到目前为止,我们的实验表明它们不会引起重大问题,但如果您想排除它们,可以使用 --exclude-partial-gene-calls 参数。

  5. 使用最初在 ITEP 中实施的 minbit heuristic(Benedict 等,2014)来消除两个氨基酸序列之间的弱匹配。您会看到,泛基因组工作流程首先通过相似性搜索来识别有些相似的氨基酸序列,然后根据这些相似性解析基因簇。在这种情况下,弱相似性可以连接不应该连接的基因簇。尽管网络分区算法可以从这些弱连接中恢复,但最好在每一步中消除尽可能多的噪声。因此,minbit 启发式提供了一种设置 a 的方法,以消除两个氨基酸序列之间的弱匹配。我们从 ITEP(Benedict MN 等人,doi:10.1186/1471-2164-15-8)那里学到了它,这是另一个泛基因组的综合分析工作流程,并决定使用它,因为它很有意义。
    简而言之,如果有两个氨基酸序列 A 和 B,则 minbit 定义为 BITSCORE(A, B) / MIN(BITSCORE(A, A), BITSCORE(B, B))。因此,如果两个序列在“较短”氨基酸序列的整个长度上非常相似,则它们之间的 minbit 得分将变为 1.0;如果 (1) 它们在非常短的一段时间内匹配,甚至与“较短”氨基酸序列的长度相比,则它们之间的 minbit 得分将变为 0.0氨基酸序列较短或(2)序列同一性之间的匹配度较低。默认 minbit 为 0.5,但您可以使用参数 --minbit 更改它。

  6. 使用 MCL 算法(van Dongen 和 Abreu-Goodger,2012)识别氨基酸序列相似性搜索结果中的簇。我们默认使用 2 作为 MCL 膨胀参数。该参数定义了算法在识别基因簇过程中的灵敏度。更高的敏感性意味着更多的聚类,但当然更多的聚类并不意味着更好的进化关系推断。有关此参数及其对簇粒度的影响的更多信息,请参见 http://micans.org/mcl/man/mclfaq.html#faq7.2,但显然,我们宏基因组学人员需要更多地讨论这一点。到目前为止,在 Meren 实验室,如果我们要比较许多相关性较远的基因组(即,基因组属于不同的科或更远的科),则使用 2;如果我们要比较非常密切相关的基因组(即,同一基因组的“品系”),则使用 10。 “物种”(基于您喜欢的这些术语的任何定义)。您可以使用参数 --mcl-inflation 更改它。请自行实验,并考虑后续解释问题!

  7. 利用每个基因簇,即使它们在分析中仅出现在一个基因组中。当然,单例或双例的重要性取决于您分析中的基因组数量,或者具体的科学问题。但是,如果您设定义截止值,可以使用参数 --min-occurrence,默认情况下为 1。增加此截止值将提高聚类速度并使可视化更易于管理,但同样,应在每个研究的背景下考虑此参数

  8. 使用欧氏距离(euclidean distance)和ward linkage来组织基因簇和基因组。您可以使用 --distance 和 --linkage 参数更改它们。

  9. 如果已有目录,请尝试使用以前的搜索结果。这样您就可以使用 --minbit、--mcl-inflation 或 --min-occurrence 参数,而无需重新进行氨基酸序列搜索来进而减少时间。如果您困惑 anvi’o 使用现有的输出文件,请参阅本期的说明。如果您对使用以前的搜索结果不感兴趣,可以删除输出目录,或使用 --overwrite-output-destinations 标志从头开始重新搜索。

完成后,将出现一个包含分析结果的新目录。您可以使用 anvi’o 附加数据表子系统在泛配置文件数据库中添加或删除(附加数据项)[https://merenlab.org/2017/12/11/additional-data-tables/]

结果可视化

分析完成后,您将使用程序 anvi-display-pan 来显示结果。

这是该命令最简单的形式:

anvi-display-pan -p PROJECT-PAN.db -g PROJECT-PAN-GENOMES.db

anvi-display-pan 程序与 anvi-interactive 程序非常相似,交互界面只是标准的 anvi’o 交互式界面针对泛基因组分析进行了细微调整。当然,anvi-display-pan 将允许您设置要提供服务的 IP 地址和端口号、添加额外的视图数据、额外的图层和/或额外的树等等。请在终端中运行 anvi-display-pan -h 来熟悉它。

以下是我们在本教程前面部分中创建的 31 个原绿球藻分离株基因组的泛基因组:

anvi-display-pan -g PROCHLORO-GENOMES.db \
                 -p PROCHLORO/Prochlorococcus_Pan-PAN.db
看起来很丑。但没关系(目前)。例如,为了稍微改进一下,您可能希望通过从“样本”选项卡 >“样本”中选择“gene_cluster rates”树,根据基因聚类结果来组织基因组

排序菜单(order menu)


order menu

当您再次绘制它时会发生以下情况(请注意右侧出现的树):
更改后效果

尽管相对原始的图显得有一定的规律,但是整体绘图还是不好看!
因此,这正是您需要开始更有效地使用界面的地方。例如,使用交互界面的设置面板中的附加设置项进行一些更改后,这是相同的泛基因组:


修饰后

如果您下载了原绿球藻数据包,您还会发现一个 anvi’o 状态文件,您可以使用程序 anvi-import-state 将其导入到您的 pan 数据库中:

anvi-import-state -p PROCHLORO/Prochlorococcus_Pan-PAN.db \
                  --state pan-state.json \
                  --name default

No excuses for bad looking pangenomes!!!(anvi'o 作者原话)

有关如何改进泛基因组的分步说明的示例,请参阅此可重复工作流程中的完善泛基因组部分。

检查基因簇

分析中的每个基因簇都将包含源自一个或多个基因组的一个或多个氨基酸序列。虽然可能存在一个“核心”部分,其中所有基因簇都会出现在每个基因组中,但也很常见地发现包含来自单个基因组的多个基因调用的基因簇(即,所有多拷贝基因)给定的基因组最终将出现在同一个基因簇中)。如果对某些基因簇感到好奇,并想要了解更多有关它们的信息。可以右键单击任何基因簇,然后您会看到此菜单(或者甚至可能更多,具体取决于您阅读本文的时间):


示例

例如,如果单击“检查基因簇”,您将看到进入该基因簇的每个基因组的所有氨基酸序列(与主显示中排列的基因组的顺序和背景颜色相同):


核酸序列

如您所知,anvi’o 交互式界面允许您从树中进行选择。因此,您可以将基因簇组选择到分类簇中(不要介意左侧面板上的数字,这显然是一个错误,并将在您的版本中修复):


原图为动图,转载过来没法搬gif就这样了
原图为动图,转载过来没法搬gif就这样了

您可以创建具有多个选择的多个分类簇,如果您愿意,甚至可以为它们指定有意义的名称:


添加分类簇名称

虽然可以选择使用树状图手动选择基因簇,但也可以使用搜索界面来识别它们,该界面允许您定义非常具体的搜索条件:


通过面板选择

在界面中突出显示这些选择,也可以将它们添加到集合中以便稍后进行总结。
此外,还可以根据功能搜索基因簇:


根据功能搜索基因簇

同样,可以将这些基因簇添加到具有任何喜欢的名称的集合中,并稍后汇总这些集合。

推断基因簇的同源性

注:感谢 Mahmoud Yousef 的努力,此功能自 anvi’o v5.2 起可用。
基因簇是好的,但并非所有基因簇都是一样的。通过简单地检查几个基因簇内的比对,您可以发现不同基因组的氨基酸序列之间不同程度的不一致。

同源性概念

基因簇可能包含来自不同基因组的几乎相同的氨基酸序列,这将是高度同质的基因簇。另一个基因簇可能包含来自不同基因组的高度差异的氨基酸序列,这将是高度非同质的,依此类推。

人们可以通过关注序列比对的两个主要属性来推断基因簇内序列同质性的性质:功能同质性(即跨基因的氨基酸残基如何保守对齐)和几何同质性(即间隙/残基分布如何)无论氨基酸的身份如何,看起来都像在基因簇内)。虽然通过手动检查基因簇内的比对序列来了解基因簇的同质性是相当简单的,但还不可能自动量化这些信息。但 anvi'o 现在可以满足您的需求。

事实上,了解基因簇内的同质性可以产生关于在密切相关的分类群中塑造基因组背景的力量的详细生态或进化见解,或者帮助进一步审查基因簇以进行手动或编程的下游分析。本节的目的是向您展示我们如何在 anvi’o 中解决这个问题并展示其功效。

anvi’o 中的功能和几何均匀性估计

Anvi’o pangenomes 包含两层,总结每个基因簇的同质性指数(以及结合两者的附加层)。这是我们的原绿球菌全基因组上下文中的一个示例(请参阅最外面的两个附加层):

几何均匀性估计

均匀性指数是根据比对计算的。如果基因簇中给定基因组的比对由于某种原因失败,则其同质性指数将显示 -1。

几何均匀性指数

几何均匀性指数表示基因簇中基因之间的几何构型程度。具体来说,我们寻找由比对确定的基因簇的缺口/残基模式。根据定义,比对过程将基因的相似部分相互比对,并通过间隙字符表示缺失的残基(或基因上下文的不同结构配置)。如果缺口/残基分布模式在整个基因簇中大部分是均匀的,则该基因簇将具有较高的几何同质性,并且最大值1.0表示比对中没有缺口。

Anvi'o 通过结合两个级别的基因簇内容分析来计算几何同质性指数:位点级别分析(即垂直对齐位置)和基因级别分析(即水平对齐位置,因为它们位于相同的位置)。我们将基因簇中的信息转换为二进制矩阵,其中间隙和残基简单地由 1 和 0 表示,并且我们利用逻辑运算符 xor 来识别和枚举具有不同模式的所有比较。在站点级同质性检查中,算法从左到右扫描并识别对齐列之间的这些差异。在基因级同质性检查中,算法从一个基因扫描到下一个基因,并识别基因中模式不相同的所有差异,从而有效地检查基因簇中间隙的分布。可以使用标志 --quick-homeity 跳过基因级几何同质性计算。虽然这不如默认方法那么准确或全面,但它可能会节省您一些时间,具体取决于您正在处理的基因组数量。

功能同质性指数

相反,功能同质性指数考虑对齐的残基(通过忽略间隙),并尝试通过考虑不同残基的生化特性(这可能影响蛋白质水平上基因的功能保守性)来量化位点中残基之间的差异。为此,我们根据氨基酸的生化特性将氨基酸分为七个不同的“保守组”。这些基团是:非极性、芳香族、极性不带电、极性和非极性特征(这些氨基酸也属于除此之外的基团)、酸性、碱性和大部分非极性(但包含一些极性特征)有关更多信息,请参阅此文章。然后,该算法遍历整个基因簇,并根据氨基酸残基的生化特性彼此的接近程度,为基因中同一位置的每对氨基酸分配 0 到 3 之间的相似性分数。所有分配的相似性得分的总和表示基因簇的功能同质性指数,如果所有残基都相同,则将达到最大值1.0。

两个指数的范围都是 0 到 1,其中 1 是完全同质的,0 是完全异质的。如果算法因运行时错误而中断(由于意外问题,例如出于某种原因并非所有基因的长度都相同等),则索引的错误值将默认为 -1。因此,如果您在摘要输出中看到 -1,则意味着我们由于某种原因未能理解该基因簇中的比对

事实上,影响蛋白质折叠的复杂过程和氨基酸残基之间错综复杂的化学相互作用应该提醒我们,这些相似性评估仅仅是数字建议,并不一定反映准确的生化见解,这不是这些同质性的目标指数。

在全基因组中使用同质性指数

考虑到所有这些,让我们再次看看我们的全基因组。它们有什么用?好吧,人们可以利用这些同质性指数做很多事情,期望本教程涵盖所有这些是不太可能的。尽管如此,这里还是试图强调它的某些方面。

为了快速了解一些最不均质的基因簇,可以根据同质性的增加或减少对基因簇进行排序。假设我们想通过增加功能同质性来对其进行排序(您可能已经知道,这可以通过主设置面板中的“项目顺序”组合框来完成):

基因簇基于同质性进行排序

如果您打开显示屏右侧的“鼠标”面板,您实际上可以将鼠标悬停在该图层上并查看所有基因簇索引的确切值。上图中光标下方显示的基因簇具有相对较低的功能同质性估计,但具有较高的几何同质性。我们可以检查基因簇来亲自看看:
检查结果

您会明白为什么相对较高的几何均匀性分数是有意义的。其中三个基因具有相同的缺口/残基模式,而另一个基因末端的缺口稍微偏离了一些,使得几何分数接近 0.75。另一方面,排列的氨基酸的颜色编码也给我们暗示它们之间缺乏功能同质性。我们可以对几何均匀性指数做同样的事情。自己尝试一下:根据几何同质性指数对泛基因组进行排序,并检查得分相对较低的基因簇。

那么,我们可以做些什么来更深入地分析我们的基因簇呢?
Anvi’o 提供了一种非常强大的方法来过滤基因簇,既可以通过命令行程序

 anvi-get-sequences-for-gene-clusters

也可以通过交互式数据探索界面:


通过交互界面进行基因过滤

探索性分析

假设您希望找到一个代表单拷贝核心基因的基因簇,该基因的几何同质性与功能同质性之间存在很大差异。例如,您想要在所有基因组中高度保守的东西,它受到结构限制以保持其比对同质,但它有很大的多样化空间,从而影响其功能异质性。你想要很多。但 anvi'o 能做到吗?
那么,对于这组非常具体的约束,您可以首先根据递减的几何同质性指数对所有基因簇进行排序,然后输入以下值来设置过滤器,然后应用它并突出显示匹配的基因簇:


挑取单拷贝核心基因组

按逆时针顺序排列的第一个基因簇是最符合列出的标准的基因簇。经过仔细检查,


检查基因簇

人们可以看到,尽管基因簇中的序列具有几何同质性,但跨基因组的比对残基之间存在巨大的功能变异性(此处仅显示了比对的一部分):
检查碱基差异

那些读过我们关于该主题(示例数据的文章)的研究的人可能不会惊讶地发现这个特定基因簇的 COG 功能注释解析为与细胞壁/膜/包膜生物发生相关的酶。

构建系统发育树

这是同质性指数的另一个使用示例。我们经常使用单拷贝核心基因簇进行系统发育分析,以估计基因组之间的进化关系。通过高级过滤器或通过基因簇的手动分箱可以轻松识别单拷贝核心基因簇:


手动分箱选择单拷贝核心基因

但是,单拷贝核心基因簇将包含许多对齐不良的基因簇,您可能不希望将其用于系统发育分析,以便最大限度地减少因生物信息学决定在何处放置缺口而产生的噪音的影响。另一方面,这个集合中会有许多几乎相同的基因簇,这对于推断系统发育关系是相当无用的。幸运的是,您可以使用同质性指数和高级搜索选项来识别那些几何上完美但功能多样的基因簇:


使用同质性指数和高级搜索选项来识别那些几何上完美但功能多样的基因簇

并即时执行快速而繁琐的分析,看看他们如何直接在界面上组织您的基因组:
构建系统发育树

或者,您可以获得包含这些基因簇的对齐和连接氨基酸序列的 FASTA 文件,以进行更合适的系统发育分析:

anvi-get-sequences-for-gene-clusters -p PROCHLORO/Prochlorococcus_Pan-PAN.db \
                                       -g PROCHLORO-GENOMES.db \
                                       -C default -b Better_core \
                                       --concatenate-gene-clusters \
                                       -o better_core.fa

不用说,每个基因簇同质性指数的估计也会出现在 anvi-summarize 的摘要文件中,以满足统计学的需求。

分割泛基因组

在某些情况下,我们可能希望将给定的泛基因组拆分为多个独立的泛基因组,例如仅包含核心基因簇的泛基因组,或仅包含单基因组的泛基因组等。当人们希望进行非常高分辨率的分析时,这种高级功能最有用。
通过精心设计使用功能和几何同质性指数来仔细检查全基因组中的核心基因来进行系统基因组分析(这也可以在不分割全基因组的情况下完成,但我们发现有时它会让事情变得更容易)。

为此,我们有 anvi-split,一个非常高效的程序来分割泛基因组结果,它也适用于泛数据库中的集合和容器。使用此程序,您可以关注给定泛基因组中的 bin 中定义的任何基因簇组,并将它们拆分为独立更小的泛基因组。

如果没看懂后续通过实际的例子进行说明。以下步骤将使用原绿球藻泛基因组进行演示。假设在原绿球菌泛基因组中,您碰巧有三个 bin:CoreHL CoreSingleton,所有这些都存储在名为 default 的集合中:

分割泛基因组

这些 bin 可以使用 anvi-split 这种方式从这个泛基因组中分割成它们自己的小泛数据库:

anvi-split -p Prochlorococcus-PAN.db \
           -g Prochlorococcus-GENOMES.db \
           -C default \
           -o SPLIT_PANs

生成的目录将包含以下文件:

ls -lR SPLIT_PANs
total 0
drwxr-xr-x  3 meren  staff  96 Apr 26 16:37 Core
drwxr-xr-x  3 meren  staff  96 Apr 26 16:37 Hl_Core
drwxr-xr-x  3 meren  staff  96 Apr 26 16:37 Singletons

SPLIT_PANs/Core:
total 3240
-rw-r--r--  1 meren  staff  1658880 Apr 26 16:37 PAN.db

SPLIT_PANs/Hl_Core:
total 1560
-rw-r--r--  1 meren  staff  798720 Apr 26 16:37 PAN.db

SPLIT_PANs/Singletons:
total 4672
-rw-r--r--  1 meren  staff  1384448 Apr 26 16:37 PAN.db

正如您所看到的,这些是单独的泛基因组,确实可以使用 anvi-display-pan 进行可视化。例如运行以下命令,

anvi-display-pan -p SPLIT_PANs/Core/PAN.db \
                 -g Prochlorococcus-GENOMES.db
分割后结果

或者运行下面的命令

anvi-display-pan -p SPLIT_PANs/Singletons/PAN.db \
                 -g Prochlorococcus-GENOMES.db
相应的结果

量化泛基因组的功能富集

一旦我们有了全基因组,我们通常想做的关键事情之一就是研究与我们的基因簇相关的功能,特别是那些附属于我们全基因组中基因组子集的功能。 Anvi’o 可以通过程序帮助您识别泛基因组中包含的某些进化枝或子进化枝丰富的功能。这是通过我们新的和改进的程序 anvi-compute-function-enrichment-in-pan 完成的。

该程序通过将基因组分成多个组来利用您提供的信息,并找到任何这些组的特征功能,如仅在个别基因组中富集的功能,但在该组之外的基因组中主要缺失的功能

要使用此功能,您必须至少拥有一个分类附加层信息(可以通过 anvi-import-misc-data 轻松完成),以及至少一个用于基因组存储的功能注释源(如果每个运行 anvi-gen-genomes-storage 时使用的 contigs 数据库使用相同的功能源进行注释)。如果您不确定基因组存储中是否有任何可用的功能,您可以对其运行 anvi-db-info 并查看输出。

anvi'o 功能丰富分析的原理简短介绍(不关心的可以直接跳过)

为了使该分析与与全基因组相关的所有其他内容兼容,我们决定应该将重点放在基因簇上,因为它们是我们全基因组的核心。这意味着功能分析的第一步是尝试将每个基因簇与一个功能关联起来。这是必要的,因为默认情况下基因簇没有功能注释,因为功能注释是在 anvi’o 中在单个基因水平上完成的。

在理想情况下,单个基因簇中的所有基因都将被注释为具有相同的相同功能。尽管实际上情况通常如此,但并不总是如此。在存在与基因簇相关的多个功能的情况下,我们选择最常见的功能,即基因簇中最多数量的基因匹配的功能(如果存在多个功能的联系,那么我们只需选择一个任意)。如果基因簇中没有任何基因被注释,则该基因簇没有与其相关的功能。

当然,当我们将每个基因簇与单个功能相关联时,我们最终可能会在泛基因组中得到具有相同功能的多个基因簇。根据我们的经验,大多数功能都与单个基因簇相关,但仍然有很多功能与多个基因簇相关,这在包含远缘相关基因组的全基因组中更为常见。在这些情况下,为了找到基因组中给定功能的出现,我们只需“合并”与同一功能相关的所有基因簇的出现(对于计算读者来说,我们只需将出现不同基因组的基因簇频率向量求和即可)。

细心的读者会注意到,我们在下文中区分了“功能注释”和“功能关联”。当我们提到“功能注释”时,我们指的是功能注释源(即COG、EggNOG等)对单个基因具有功能的注释,而基因簇的“功能关联”是基因簇的关联具有如上所述的单一功能。

好的,现在我们有了基因组中的功能频率表,我们将其用作功能富集测试的输入。该测试由 Amy Willis 在 R 中实现(您可以在此处找到脚本),并使用带有 logit 链接函数的广义线性模型(Generalized Linear Model)来计算每个函数的富集分数和 p 值。使用 qvalue 包对 p 值进行错误检测率校正以考虑多个测试。

除了丰富测试之外,我们还使用简单的启发式方法来查找与每个函数关联的组。这种关联仅对真正丰富的功能才有意义,否则应被忽略。我们简单地确定,对于每个功能,相关组是基因组功能的出现次数大于均匀分布下的预期出现次数的组(如:如果该功能在所有群体的基因组中出现的可能性相同)。从数学上来说(如果你喜欢这类东西),如果我们表示E_ij 作为在第i组的j功能下基因组数量。在零分布下,我们认为零分布是均匀分布。

估算公式

将实际出现的次数计作 O_ij ;然后我们认为 O_ij>E_ij
现在我们可以使用所有这些信息来探索我们的数据,请参阅下面的详细信息!

让我们用原绿球藻的例子来了解我们可以用它做什么。

首先,我们将比较低光基因组与高光基因组,以查看使用“光”分类附加层数据是否有任何一组所独有的功能(如果这对您来说没有意义,请返回到上面的泛基因组图之一,查看附加数据层“light”):

anvi-compute-functional-enrichment-in-pan -p PROCHLORO/Prochlorococcus_Pan-PAN.db \
                                          -g PROCHLORO-GENOMES.db \
                                          --category light \
                                          --annotation-source COG14_FUNCTION \
                                          -o enriched-functions.txt

以下是输出文件 riched-functions.txt 的结构(还有更多列,向右滚动即可看到它们):


结果(部分)

结果(剩余部分)

以下描述了此输出中的每一列:
COG14_FUNCTION 此列具有计算富集的特定函数的名称。在此示例中,我们选择使用 COG14_FUNCTION 进行功能注释,因此列标题为 COG14_FUNCTION。您可以使用 --annotation-source 指定 PAN 数据库中的功能注释源,然后根据该注释源进行分析。即使您的基因组存储中有多个功能注释源,该程序的单次运行也只能使用一个源。如果您愿意,您可以多次运行它,并且每次都使用不同的注释源。如果您不记得基因组存储中有哪些注释源可用,您可以使用标志 --list-annotation-sources。

richment_score 是一个分数,用于衡量该函数对于属于特定群体的基因组与泛基因组中的所有其他基因组的独特性有多少。该乐谱由艾米·威利斯 (Amy Willis) 开发。有关我们如何生成此分数的更多详细信息,请参阅下面艾米的注释。

unadjusted_p_value 是富集测试的 p 值(多次测试未调整)。

adjustment_q_value 是调整后的 q 值,用于控制由于多次测试而导致的错误检测率 (FDR)(这是必要的,因为我们分别运行每个函数的富集测试)。

Associated_groups 是分类数据中与该函数关联的组(或标签)的列表。请注意,如果富集分数较低,则这是没有意义的(如果功能未富集,则它实际上与任何特定组没有关联)。请参阅上面 Alon 的解释以了解它们是如何计算的。

function_accession 是函数登录号。

gene_clusters_ids 是与此功能相关的基因簇。请注意,每个基因簇将与单个功能相关联,但一个功能可以与多个基因簇相关联。

每个组的 p_LL、p_HL(在此示例中,有两个组:LL 和 LH)将有一列,其中包含我们检测到函数的组成员部分。

每个组的 N_LL、N_HL 这些列指定该组中基因组的总数。
注:我们此处的示例仅包含两个类别(LL 和 HL),但您可以根据需要拥有任意多个不同的类别。请记住,如果您的某些群体的基因组非常少,那么统计测试将不会非常可靠。使测试可靠的组中基因组的最小数量取决于许多因素,但如果您的任何组的基因组少于 8 个,我们建议谨慎选择。

艾米·威利斯 (Amy Willis) 关于功能丰富评分的注释

这里提出并在 anvi’o 中实现的功能丰富评分是比例相等的 Rao 检验统计量。本质上,它将每个类别(此处为高光和低光)视为逻辑回归(二项式 GLM)中的解释变量,并测试分类变量在解释函数出现时的显着性。

该测试解释了这样一个事实:从一类观察到的基因组可能多于另一类。像往常一样,拥有更多的基因组使测试更加可靠。有很多不同的方法可以进行此测试,但我们做了一些调查,发现 Rao 测试在控制 1 类错误率的所有测试中具有最高的功效。耶!

由于许多用户会考虑测试多种功能的丰富性,因此默认情况下,我们通过控制错误发现率来调整多重测试。因此,如果您使用此统计数据,请在论文中报告 q 值而不是 p 值。

现在让我们搜索表中最重要的功能之一“Ser/Thr 蛋白激酶 RdoA 参与 Cpx 应激反应,MazF 拮抗剂”,该功能针对 LL 组的成员进行了富集,我们可以在表中看到它匹配四个基因簇。
基因富集

事实上,如果我们仔细观察,我们会发现这个功能与四个低光进化枝中每一个都是唯一的基因簇相匹配,并且是我们泛基因组中低光基因组的单拷贝核心基因。

让我们看看另一个丰富的功能:核酸外切酶 VII,大亚基。当我们搜索这个函数时,我们应该小心,因为函数的名称中包含逗号,而交互界面中的函数搜索选项将逗号视为分隔多个要同时搜索的函数。因此我只是搜索了核酸外切酶 VII,结果如下:


核酸外切酶Ⅶ

我们可以看到,搜索匹配了核酸外切酶 VII、大亚基和小亚基的匹配结果,总共有 22 个匹配结果。


核酸外切酶

大亚基与 CORE LL 中的单个基因簇相匹配,小亚基与每个分支特异性核心中的一个基因簇相匹配(类似于我们上面看到的 Ser/Thr 蛋白激酶)。这两个基因也是低光成员特有的单拷贝核心的一部分,而高光成员则不存在。

计算基因组的平均核苷酸同一性(ANI)(以及其他基因组相似性指标!)

Anvi’o 还包含一个程序 anvi-compute-genome-similarity,它使用各种相似性指标(例如 PyANI 来计算基因组中的平均核苷酸同一性),并使用 sourmash 来计算基因组中的混搭距离。它需要外部基因组文件、内部基因组文件或指向 FASTA 文件路径的 fasta 文本文件的任意组合(假设每个 FASTA 为 1 个基因组)。此外,anvi-compute-genome-similarity 可选择接受 pan-db 以将所有结果添加到其中作为附加层数据。

以下是我们的原绿球藻全基因组的示例:

anvi-compute-genome-similarity --external-genomes external-genomes.txt \
                               --program pyANI \
                               --output-dir ANI \
                               --num-threads 6 \
                               --pan-db PROCHLORO/Prochlorococcus_Pan-PAN.db

完成后,我们可以再次使用 anvi-display-pan 可视化泛基因组,看看有什么新内容:

anvi-display-pan -g PROCHLORO-GENOMES.db \
                 -p PROCHLORO/Prochlorococcus_Pan-PAN.db

当你第一眼看到它时,你不会发现任何异常。但是,如果您转到“图层”选项卡,您将看到以下添加内容:

ANI结果

如果单击 ANI_percentage_identity 等复选框,您将看到一组新条目将添加到图层数据条目列表中。然后,如果单击“绘制”或“重绘图层数据”按钮,您应该会看到 ANI 添加到显示中:
添加ANI信息

您可能需要更改界面中的最小值,以便更好地表示您自己的泛基因组中整个基因组的 ANI。

关于如何排序基因组说明

如何根据基因簇对基因组进行排序一直是我的一个疑问。是否应该使用基因簇存在/不存在模式来组织基因组(其中忽略旁系同源物并使用二进制表来对基因组进行聚类)?或者应该依赖基因簇频率数据来排序(其中确实考虑了旁系同源物,因此它们的表不再是二进制的,而是包含每个基因簇中每个基因组的基因数量的频率)?

感谢 Özcan 对代码库的新添加,当我处理本教程的这一部分有关泛基因组中“如何对基因组进行排序”的问题时,我得到了意想不到的观察结果。

当我根据基因簇频率数据排序原绿球藻基因组时,ANI 矩阵如下所示:


据基因簇频率数据排序原绿球藻基因组

相比之下,当我根据基因簇存在/不存在数据对基因组进行排序时,ANI 矩阵如下所示:


依据ANI对基因组进行排列

如果您正在使用分离的基因组,也许您应该考虑根据基因簇频率对它们进行排序,因为您期望基于基因簇的基因组的正确组织通常应与其整体相似性相匹配。以下是一些小注意事项:

在这种情况下,很容易看出,通过基因簇存在/不存在数据进行组织的效果很差,因为混合了高光组 II(紫色)和高光组 I(便便颜色)基因组。我们可以清楚地识别问题,因为我们知道这些基因组所属的进化枝,在您不知道这些关系的情况下,基因簇频率似乎可以更准确地组织您的基因组。

令人遗憾的是,MAG 通常缺乏多拷贝基因,并且它们的旁系同源物较少(因为短读段的组装会创建非常复杂的 de Bruijn 图,从而消除给定基因组中的所有冗余),因此它们的正确组织不会真正实现从基因簇频率数据中受益,至少不如分离基因组那么多,并且我们在基因簇存在/不存在数据的组织中看到的错误将对我们的泛基因组产生更大的影响。

看到高光组 II 中那些子进化枝的出现了吗?如果您有兴趣,请看一下我们关于原绿球菌宏基因组的论文(特别是“宏基因组揭示了具有不同适应性水平的密切相关的分离株”部分),其中我们将 HL-II 分为多个分支,并表明尽管差异很小在这些基因组之间,这些分支表现出显着不同的环境分布模式。我们可以在该研究中根据基因簇频率和环境分布模式定义的子分支与 ANI 揭示的群体非常匹配。这可以追溯到微生物学中一些最基本的问题,即如何定义微生物生命的生态相关单位。我不知道我们是否接近做出任何定义,但我可以告诉你,我们正在使用的那些“系统发育标记”……我不确定它们是否真的那么有效。

生成泛基因组分析报告

当您将选择的基因簇存储为集合时,anvi'o 将允许您总结这些结果。(即使您想简单地总结泛基因组中的所有内容而不在界面中进行任何选择,您仍然需要泛数据库中的集合。但幸运的是,您可以使用程序 anvi-script-add-default-collection 添加包含每个基因簇的默认集合。)

摘要步骤为您提供了两个重要的东西:一个静态 HTML 网页,您可以压缩并与同事共享或作为补充数据文件添加到您的文章中,以及输出目录中描述每个基因簇的全面制表符分隔文件。

您可以使用 anvi-summarize 程序对集合进行汇总,该命令的通用形式如下所示:

anvi-summarize -p PROJECT-PAN.db \
                 -g PROJECT-PAN-GENOMES.db \
                 -C COLLECTION_NAME \
                 -o PROJECT-SUMMARY

如果打开摘要目录中的 index.html 文件,您将看到一个输出,其中包含有关分析的一些基本信息:
报告的基本信息

以及基因簇的制表符分隔文件:
table数据

该文件的结构如下所示,并且使您能够以更详细的方式研究您的基因簇(您可能需要向右滚动才能查看更多表格):


详细数据

该文件会将每个基因组中的每个基因与您通过界面或程序 anvi-import-collection 所做和命名的每个选择联系起来,并且还允许您访问每个基因的氨基酸序列和功能。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容