geNomad
geNomad:移动遗传元件的鉴定
功能
geNomad 的主要目标是在测序数据(分离株、宏基因组和宏转录组)中鉴定病毒和质粒。它还提供了几个额外的功能,可以帮助您进行分析:
- 病毒基因组的分类学归属。
- 鉴定整合在宿主基因组中的病毒(前病毒)。
- 蛋白质的功能注释。
文档
有关安装说明、geNomad 工作原理的信息以及如何执行它的详细说明,请查阅完整文档:https://portal.nersc.gov/genomad/
Web 应用
geNomad 可作为 Web 应用程序在 Galaxy 和 NMDC EDGE 平台上使用。在那里,您可以上传序列数据,直接在浏览器中查看结果,并下载输出以进行进一步分析。
引用 geNomad
如果您在工作中使用了 geNomad,请考虑引用其论文:
Identification of mobile genetic elements with geNomad
Camargo, A. P., Roux, S., Schulz, F., Babinski, M., Xu, Y., Hu, B., Chain, P. S. G., Nayfach, S., & Kyrpides, N. C. — Nature Biotechnology (2023), DOI: 10.1038/s41587-023-01953-y.
快速入门
我们建议用户在开始使用 geNomad 之前阅读 documentation。不过,如果您时间紧迫,可以按照此快速分步示例进行操作。
安装
首先,您需要安装 geNomad。有几种方法可以做到这一点(有关更多信息,请参阅文档),但两个方便的选择是使用 Pixi 或 Mamba。它们都会为您处理所有依赖项的安装。
Pixi 允许您将 geNomad 安装为全局可用的命令,以便轻松执行。
pixi global install -c conda-forge -c bioconda genomad
使用 Mamba,你需要为 geNomad 创建一个环境,并在能够使用它之前激活该环境。
# Create an environment for geNomad
mamba create -n genomad -c conda-forge -c bioconda genomad
# Activate the geNomad environment
mamba activate genomad
另一个选项是通过 Docker 使用 geNomad。
# Pull the image
docker pull antoniopcamargo/genomad
# Run the image
docker run --rm -ti -v "$(pwd):/app" antoniopcamargo/genomad
下载数据库
geNomad 依赖于一个数据库,该数据库包含用于分类序列的标记物的特征信息、其分类学信息、其功能注释等。因此,您应首先将数据库下载到当前目录:
genomad download-database .
数据库将包含在 genomad_db 目录中。
如果您愿意,也可以从 Zenodo 下载数据库并手动解压。
运行 geNomad
现在您可以开始了!geNomad 通过按顺序执行一系列模块来工作(您可以在 pipeline documentation 中了解有关此的更多信息),但我们提供了一个方便的 end-to-end 命令,可以一次性为您执行整个流程。
在本示例中,我们将使用 Klebsiella pneumoniae 基因组(GCF_009025895.1)作为输入。您可以使用任何包含核苷酸序列的 FASTA 文件作为输入。geNomad 适用于分离株基因组、宏基因组和宏转录组。
执行 geNomad 的命令结构如下:
genomad end-to-end [OPTIONS] INPUT OUTPUT DATABASE
因此,要运行完整的 geNomad 流程(end-to-end 命令),以核苷酸 FASTA 文件(GCF_009025895.1.fna.gz)和数据库(genomad_db)作为输入,我们将执行以下命令:
genomad end-to-end --cleanup --splits 8 GCF_009025895.1.fna.gz genomad_output genomad_db
结果将写入 genomad_output 目录中。
关于上述命令的三个重要细节:
- 使用了
--cleanup选项来强制 geNomad 删除执行过程中生成的中间文件。这将为您节省一些存储空间。 - 此处使用了
--splits 8参数,以便在 notebook 中运行此示例。geNomad 会搜索一个占用大量内存的蛋白质谱大型数据库。为防止因内存不足导致执行失败,我们可以使用--splits参数将搜索拆分为块。如果您在大型服务器上运行 geNomad,可能无需拆分搜索,从而提高执行速度。 - 请注意,我用作输入的 FASTA 文件是经过压缩的。这是可能的,因为 geNomad 支持以
.gz、.bz2或.xz格式压缩的输入文件。
[!NOTE] 默认情况下,geNomad 会应用一系列分类后过滤器以移除可能的假阳性。例如,序列的质粒或病毒得分必须至少为 0.7,且长度短于 2,500 bp 的序列必须编码至少一个标志性基因。如果您想禁用分类后过滤器,请在您的命令中添加
--relaxed标志。另一方面,如果您希望分类非常保守,可以使用--conservative标志。这将使分类后过滤器更加严格,防止缺乏强支持的序列被分类为质粒或病毒。您可以在此处查看默认、宽松和保守的分类后过滤器 here。
理解输出
在此示例中,geNomad 分析的结果将写入 genomad_output 目录,其结构如下:
genomad_output
├── GCF_009025895.1_aggregated_classification
├── GCF_009025895.1_aggregated_classification.log
├── GCF_009025895.1_annotate
├── GCF_009025895.1_annotate.log
├── GCF_009025895.1_find_proviruses
├── GCF_009025895.1_find_proviruses.log
├── GCF_009025895.1_marker_classification
├── GCF_009025895.1_marker_classification.log
├── GCF_009025895.1_nn_classification
├── GCF_009025895.1_nn_classification.log
├── GCF_009025895.1_summary
╰── GCF_009025895.1_summary.log
如上所述,geNomad 通过顺序执行多个模块来工作。这些模块中的每一个都会生成一个日志文件(<prefix>_<module>.log)和一个子目录(<prefix>_<module>)。
在本示例中,我们仅查看 GCF_009025895.1_summary 内的文件。<prefix>_summary 目录包含汇总整个流程中生成结果的文件。如果您只需要输入中识别出的质粒和病毒列表,这就是您正在寻找的内容。
genomad_output
╰── GCF_009025895.1_summary
├── GCF_009025895.1_plasmid.fna
├── GCF_009025895.1_plasmid_genes.tsv
├── GCF_009025895.1_plasmid_proteins.faa
├── GCF_009025895.1_plasmid_summary.tsv
├── GCF_009025895.1_summary.json
├── GCF_009025895.1_virus.fna
├── GCF_009025895.1_virus_genes.tsv
├── GCF_009025895.1_virus_proteins.faa
╰── GCF_009025895.1_virus_summary.tsv
首先,让我们看看 GCF_009025895.1_virus_summary.tsv:
seq_name length topology coordinates n_genes genetic_code virus_score fdr n_hallmarks marker_enrichment taxonomy
-------------------------------------- ------ ------------------- --------------- ------- ------------ ----------- --- ----------- ----------------- -----------------------------------------------------------------
NZ_CP045015.1|provirus_3855947_3906705 50759 Provirus 3855947-3906705 79 11 0.9772 NA 16 73.7974 Viruses;Duplodnaviria;Heunggongvirae;Uroviricota;Caudoviricetes;;
NZ_CP045015.1|provirus_2885031_2934610 49580 Provirus 2885031-2934610 70 11 0.9769 NA 13 73.2757 Viruses;Duplodnaviria;Heunggongvirae;Uroviricota;Caudoviricetes;;
NZ_CP045018.1 51887 No terminal repeats NA 57 11 0.9760 NA 14 65.4720 Viruses;Duplodnaviria;Heunggongvirae;Uroviricota;Caudoviricetes;;
…
此表格文件列出了 geNomad 在您的输入中发现的所有病毒,并为您提供了一些关于它们的便捷信息。以下是各列的内容:
seq_name: 输入 FASTA 文件中序列的标识符。前病毒将采用以下命名方案:<sequence_identifier>|provirus_<start_coordinate>_<end_coordinate>。length: 序列的长度(对于整合病毒,则为前病毒的长度)。topology: 病毒序列的拓扑结构。可能的值为:No terminal repeats、DTR(直接末端重复)、ITR(反向末端重复)或Provirus(整合在宿主基因组中的病毒)。coordinates: 前病毒区域在宿主序列中的 1 起始坐标。对于未预测为整合的病毒,该值将为NA。n_genes: 序列中编码的基因数量。genetic_code: 预测的遗传密码。可能的值为:11(细菌和古菌的标准密码)、4(重编码的 TGA 终止密码子)或 15(重编码的 TAG 终止密码子)。virus_score: 衡量 geNomad 对序列是病毒的置信度的指标。得分接近 1.0 的序列比得分较低的序列更可能是病毒。fdr: 分类的估计错误发现率(FDR)(即,直到该行序列中假阳性的预期比例)。为了估计 FDR,geNomad 需要 score calibration,该功能默认关闭。因此,在此示例中,该列将仅包含NA值。n_hallmarks: 匹配到标志性 geNomad 标记的基因数量。标志性基因是先前与病毒功能相关的基因,它们的存在是序列确实是病毒的强烈指示。marker_enrichment: 一个代表序列中病毒标记物总富集程度的分数。随着序列中病毒标记物数量的增加,该值会增大,因此具有多个标记物的序列将获得更高的分数。染色体和质粒标记物会降低该分数。taxonomy: 病毒基因组的分类学归属。谱系遵循 ICTV 分类学发布 MSL39 中包含的分类系统。病毒可以分类到科级别,但不能分类到该科内的特定属或种。分类系统以固定数量的字段(对应分类等级)表示,字段之间用分号分隔,空字段保持空白。
在我们的示例中,geNomad 识别出多个整合到 K. pneumoniae 基因组中的前病毒以及一个染色体外噬菌体。由于它们都具有高分和标记富集,我们可以确信这些确实是病毒。它们均被预测使用遗传密码 11,并被归类到 Caudoviricetes 纲,该纲包含所有带尾细菌噬菌体。对于这些病毒的 taxonomy 字段,在 Caudoviricetes 之后有两个连续的分号,因为 geNomad 只能将它们分配到纲级别,导致目和科级别为空。
另一个重要文件是 GCF_009025895.1_virus_genes.tsv。在执行过程中,geNomad 使用包含染色体、质粒和病毒特异性标记的数据库,对输入序列编码的基因进行注释。<prefix>_virus_genes.tsv 文件总结了已识别病毒所编码基因的注释。
gene start end length strand gc_content genetic_code rbs_motif marker evalue bitscore uscg plasmid_hallmark virus_hallmark taxid taxname annotation_conjscan annotation_amr annotation_accessions annotation_description
--------------- ----- ----- ------ ------ ---------- ------------ ----------- ----------------- ---------- -------- ---- ---------------- -------------- ----- -------------- ------------------- -------------- -------------------------------- --------------------------------------------------------------------------------------
NZ_CP045018.1_1 1 399 399 1 0.536 11 None GENOMAD.108715.VP 2.675e-31 120 0 0 1 1246 Caudoviricetes NA NA PF05100;COG4672;TIGR01600 Phage minor tail protein L
NZ_CP045018.1_2 401 1111 711 1 0.568 11 AGGAG GENOMAD.168265.VP 1.523e-39 149 0 0 0 1246 Caudoviricetes NA NA PF14464;COG1310;K21140;TIGR02256 Proteasome lid subunit RPN8/RPN11, contains Jab1/MPN domain metalloenzyme (JAMM) motif
NZ_CP045018.1_3 1143 1493 351 1 0.382 11 AGGAG GENOMAD.147875.VV 7.818e-13 66 0 0 0 1246 Caudoviricetes NA NA COG5633;TIGR03066 NA
NZ_CP045018.1_4 1509 2120 612 1 0.477 11 GGA/GAG/AGG GENOMAD.143103.VP 2.238e-48 173 0 0 1 1246 Caudoviricetes NA NA PF06805;COG4723;TIGR01687 Phage-related protein, tail component
NZ_CP045018.1_5 2183 13516 11334 1 0.566 11 None GENOMAD.159864.VP 7.104e-262 901 0 0 0 1246 Caudoviricetes NA NA PF12421;PF09327 Fibronectin type III protein
NZ_CP045018.1_6 13585 15084 1500 1 0.550 11 AGGAG GENOMAD.195756.VP 1.695e-13 76 0 0 0 1246 Caudoviricetes NA NA NA NA
NZ_CP045018.1_7 15163 16128 966 -1 0.469 11 GGAGG NA NA NA 0 0 0 1 NA NA NA NA NA
…
此文件中的列如下:
gene: 基因的标识符(<sequence_name>_<gene_number>)。通常,基因编号从 1 开始(序列中的第一个基因)。然而,整合在宿主染色体中部的前噬菌体编码的基因可能从不同的编号开始,具体取决于其在染色体中的位置。start: 基因的 1 起始坐标。end: 基因的 1 结束坐标。length: 基因位点的长度(以碱基对为单位)。strand: 编码该基因的单链。可以是 1(正向链)或 -1(反向链)。gc_content: 基因位点的 GC 含量。genetic_code: 预测的遗传密码(详见摘要文件的说明)。rbs_motif: 检测到的核糖体结合位点基序。marker: 最佳匹配的 geNomad 标记。如果该基因不匹配任何标记,值将为NA。evalue: 该基因编码的蛋白质与最佳匹配的 geNomad 标记之间的比对 E 值。bitscore: 该基因编码的蛋白质与最佳匹配的 geNomad 标记之间的比对 Bitscore。uscg: 分配给该基因的标记是否对应于通用单拷贝基因(UCSG,定义见 BUSCO v5)。这些基因预期存在于染色体中,在质粒和病毒中较为罕见。可以是 1(基因是 USCG)或 0(基因不是 USCG)。plasmid_hallmark: 分配给该基因的标记是否代表质粒特征。virus_hallmark: 分配给该基因的标记是否代表病毒特征。taxid: 分配给该基因的标记的分类学标识符(您可以忽略此项,因为它旨在供 geNomad 内部使用)。taxname: 与分配的 geNomad 标记相关的分类单元名称。在此示例中,我们可以看到注释的蛋白质均具有 Caudoviricetes 的特征(这就是为什么前病毒被分配到该纲的原因)。annotation_conjscan: 如果与基因匹配的标记是接合相关基因(如 CONJscan 中定义),此字段将显示分配给该标记的 CONJscan 登录号。annotation_amr: 如果与基因匹配的标记被注释了抗菌素耐药性(AMR)功能(如 NCBIfam-AMRFinder 中定义),此字段将显示分配给该标记的 NCBIfam 登录号。annotation_accessions: 部分 geNomad 标记具有功能注释。此列告诉您分配给该标记的 Pfam、TIGRFAM、COG 和 KEGG 条目。annotation_description: 描述分配给标记的功能的文本。
在上面的示例中,我们可以看到由 NZ_CP045018.1 编码的前七个基因的信息。最后一个条目未匹配到任何 geNomad 标记。前六个均被分配至蛋白质家族,其中一些是带尾噬菌体的典型特征(例如小尾蛋白),这让我们确信这些确实是 Caudoviricetes。
这里的一个重要细节是,geNomad 标记的主要用途是分类。它们被设计为特异性针对染色体、质粒或病毒,从而能够区分属于这些类别的序列。因此,你不应该期望每一个病毒基因都会被标注一个 geNomad 标记。如果你希望尽可能彻底地注释序列中的基因,应使用诸如 Pfam 或 COG 之类的数据库。
摘要目录中另外两个与病毒相关的文件是 GCF_009025895.1_virus.fna 和 GCF_009025895.1_virus_proteins.faa。它们分别是已识别病毒序列及其蛋白质的 FASTA 文件。前病毒会自动从宿主序列中切除。
接下来是质粒,与其识别相关的数据可以在 <prefix>_plasmid_summary.tsv、<prefix>_genes.tsv、<prefix>_plasmid.fna 和 <prefix>_plasmid_proteins.faa 文件中找到。这些文件与它们的病毒对应文件大多非常相似。<prefix>_plasmid_summary.tsv(如下所示)中的差异如下:
<prefix>_virus_summary.tsv(coordinates和taxonomy)中存在的病毒特异性列不存在。conjugation_genes列列出了可能参与接合的基因。需要注意的是,此类基因的存在不足以判断给定质粒是接合性质粒还是可移动性质粒。如果您有兴趣鉴定接合性质粒,我们建议您使用 geNomad 识别的质粒,并通过 CONJscan 进行分析。amr_genes列列出了注释为具有抗菌耐药功能的基因。您可以在 AMRFinderPlus 网站上查看与每个登录号相关的具体功能。
seq_name length topology n_genes genetic_code plasmid_score fdr n_hallmarks marker_enrichment conjugation_genes amr_genes
------------- ------ ------------------- ------- ------------ ------------- --- ----------- ----------------- ----------------------------------------------------------------------------------------------------- -----------------
NZ_CP045020.1 28729 No terminal repeats 36 11 0.9954 NA 6 26.4290 F_traE NA
NZ_CP045022.1 50635 No terminal repeats 61 11 0.9946 NA 9 44.8458 T_virB1;T_virB3;virb4;T_virB5;T_virB6;T_virB8;T_virB9 NA
NZ_CP045019.1 44850 No terminal repeats 52 11 0.9943 NA 4 27.8509 F_traE NA
NZ_CP045016.1 82240 No terminal repeats 110 11 0.9935 NA 11 36.6048 T_virB8;T_virB9;F_traF;F_traH;F_traG;T_virB1 NF000270;NF012171
NZ_CP045021.1 5251 No terminal repeats 7 11 0.9932 NA 2 3.1408 MOBP1 NA
NZ_CP045017.1 61331 No terminal repeats 76 11 0.9929 NA 16 35.2570 I_trbB;I_trbA;MOBP1;I_traI;I_traK;I_traL;I_traN;I_traO;I_traP;I_traQ;I_traR;traU;I_traW;I_traY;F_traE NA