宏基因组中的病毒组分析（1）病毒序列的鉴定geNomad

从宏基因组数据中鉴定病毒序列的主要方法可以分为：序列比对法，特征分析方法和机器学习方法

1、geNomad简介

这里介绍的geNomad，一种结合序列比对和机器学习用于识别和分类病毒序列的工具。geNomad 擅长处理复杂和高噪声的环境数据，适用于多样化的宏基因组样本。geNomad运行速度快，可对大型数据集进行分析。https://github.com/apcamargo/genomad

geNomad 的主要目标是鉴定测序数据（分离株、宏基因组和宏转录组）中的病毒和质粒。包括：

• 病毒基因组的分类分配。
• 鉴定整合到宿主基因组中的病毒（前病毒）。
• 蛋白质的功能注释。

2、安装和使用

a、创建环境 geNomad
b、安装数据库 genomad db_genomad
c、一步执行

2.1.安装

这里介绍使用mamba进行安装，这个方法比较方便，可以一起安装所需要的所有依赖项

# 为geNomad创建环境
mamba create -n genomad -c conda-forge -c bioconda genomad
# 打开geNomad环境
mamba activate genomad
#利用mamba进行下载
mamba install geNomad

geNomad适用于分离基因组、宏基因组和宏转录组中的病毒序列，需要下载数据库（大概1.36G）：

方法一：直接下载

genomad download-database .

如果想要更改数据库的位置，只需要在数据库前添加路径即可。

genomad  <file_path> download-database .

方法二：手动下载

从Zenodo（geNomad database (zenodo.org)）下载数据库并手动提取，推荐手动下载，更快一些

2.2.运行

geNomad可以一步执行，命令如下：

genomad end-to-end [OPTIONS] INPUT OUTPUT DATABASE

解释如下：

end-to-end:一步执行
[OPTIONS]：一些可以定义的选项
         --cleanup：强制 geNomad 删除执行过程中生成的中间文件，节省空间
         --splits 8：geNomad 搜索一个占用大量内存空间的蛋白质谱的大型数据库。为了防止由于内存不足 
                     导致执行失败，我们可以使用 parameter 将搜索拆分为 chuncks。如果您在大型服务 
                     器中运行 geNomad，则可能不需要拆分搜索，从而提高执行速度
INPUT：输入文件，一般为核苷酸序列的FASTA文件，geNomad支持.gz.bz2.xz
OUTPUT：输出文件，可以规定输出文件的路径
DATABASE：数据库genomad_db，如果数据库位置改变的话，这里也要相应改变，添加改变后的路径

我选择的运行命令是：

genomad end-to-end --cleanup --splits 8  <input> <output>  /mnt/f/database/genomad_db

input和output是输入文件和输出文件的位置，我的电脑配置是64G的运行内存和CPU12代i7-12700H，平均每个跑2小时，内存消耗16-18G，每个序列的数据量不同，时间也不同，快得时候半小时就跑完了，可以从两个方面提高运行速度：

1、不带--splits 8参数

2、使用合适的 --threads，默认是64，但一般可能没有这么多。

另外，使用 --disable-nn-classification 选项禁用基于神经网络的分类，这也将禁用分数聚合并强制 geNomad 仅依赖基于标记的分类器。可能对 geNomad 的分类性能产生负面影响。不推荐
使用 --sensitivity 参数降低 MMseqs2 搜索的灵敏度，该搜索将标记分配给基因。这将使注释模块更快，但也会减少分配给标记的基因数量。可能对 geNomad 的分类性能产生负面影响。不推荐

3、输出文件

geNomad 的分析结果将写入设定的输出文件目录：

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

每个模块都自带一个日志文件，.log结尾的文件。我们一般需要查看的是summary文件，里面包含了鉴定出的质粒（plasmid）和病毒（virus）的列表。

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_2885510_2934610   49101    Provirus              2885510-2934610   69        11             0.9776        NA    14            76.0892             Viruses;Duplodnaviria;Heunggongvirae;Uroviricota;Caudoviricetes;;
NZ_CP045015.1|provirus_3855947_3906705   50759    Provirus              3855947-3906705   79        11             0.9774        NA    16            75.1552             Viruses;Duplodnaviria;Heunggongvirae;Uroviricota;Caudoviricetes;;
NZ_CP045018.1                            51887    No terminal repeats   NA                57        11             0.9774        NA    14            67.7749             Viruses;Duplodnaviria;Heunggongvirae;Uroviricota;Caudoviricetes;;
…

这个表格文件列出了 geNomad 发现的所有病毒。以下是每列包含的内容：

seq_name：输入 FASTA 文件中序列的标识符。
length：序列的长度（或前病毒，在整合病毒的情况下）。
topology：病毒序列的拓扑结构。
coordinates：宿主序列中前病毒区域的索引坐标。将用于未预测会集成的病毒。
n_genes：序列中编码的基因数。
genetic_code： 预测的遗传密码。可能的值为：11（细菌和古细菌的标准代码）、4（重新编码的 TGA 终止密码子）或 15（重新编码的 TAG 终止密码子）。
virus_score： 衡量 geNomad 对序列是病毒的置信度的指标。分数接近 1.0 的序列比分数较低的序列更有可能是病毒。
fdr：分类的估计错误发现率 （FDR）（即直到此行的序列中假阳性的预期比例）。默认情况下，校准处于关闭状态。
n_hallmarks：与标志性 geNomad 标记匹配的基因数量。标志是以前与病毒功能相关的基因，它们的存在强烈表明该序列确实是一种病毒。
marker_enrichment：表示序列中病毒标志物总富集的分数。该值随着序列中病毒标记物数量的增加而变化，因此具有多个标记物的序列将具有更高的分数。染色体和质粒标记物会降低分数。
taxonomy： 病毒基因组的分类分配。病毒可以分类学上分配到科水平，但不能分配给该科中的特定属或种。分类以固定数量的字段（对应于分类等级）显示，用分号分隔，空字段留空。

4、不足

geNomad标记的主要目的是分类。主要是染色体、质粒或病毒。并不是每个病毒基因都会用geNomad标记进行注释。如果想尽可能地注释序列中的基因，可以使用Pfam或COG等数据库。

参考来源：apcamargo/genomad: geNomad: Identification of mobile genetic elements (github.com)