宏病毒组分析 ：病毒基因组的质量评估、纯化与vOTUs构建

在宏基因组数据分析中，从复杂的序列混合物中识别出病毒Contigs是第一步。然而，这些初步鉴定出的序列并非都是完整的病毒基因组，它们可能仅仅是病毒基因组的片段、存在宿主基因污染，甚至是被整合到宿主染色体中的前病毒（provirus）。此外，由于测序深度和样本多样性，我们还可能拥有大量高度相似的病毒序列，它们实际上代表了同一“物种”的不同株系。

因此，对这些初步鉴定的病毒序列进行质量评估、纯化以及去冗余聚类，是构建高质量病毒操作分类单元（viral Operational Taxonomic Units, vOTUs）的关键步骤。本篇文章将深入介绍两个核心工具：CheckV 用于病毒基因组的质量控制和纯化，以及 dRep (或 CD-HIT) 用于将相似病毒序列聚类成vOTUs，为后续的病毒分类、功能和宿主预测奠定坚实的基础。

本章的核心作用是将初步筛选的病毒Contigs转化为一套高质量、非冗余的vOTUs集合，确保每个vOTU都具有可靠的生物学意义，并最大程度地排除非病毒序列和低质量序列的干扰。

病毒基因组质量评估与纯化 (Viral Genome Quality Assessment & Refinement)

此步骤的目的是评估初步筛选出的病毒Contigs的完整性、识别潜在的宿主基因污染，并区分游离病毒和整合性前病毒。

软件工具：CheckV

CheckV 是一款专门为病毒宏基因组学设计的工具，能够评估病毒Contigs的完整度、预测其是游离病毒还是前病毒，并识别和剪切与宿主序列的嵌合片段。

• 核心作用与优势：

完整性评估：基于保守病毒蛋白或基因组结构推断病毒Contig的完整度（completeness），包括“Complete”、“High-quality”、“Medium-quality”、“Low-quality”等分级。

宿主污染识别：通过与参考数据库比对，识别病毒Contig两端可能存在的宿主基因片段，并将其剪除以获得纯化的病毒序列。

前病毒识别：能够鉴定出整合在宿主基因组中的前病毒序列。

环状基因组推断：对一些完整性很高的Contigs，可以推断其是否为潜在的环状基因组。

• 环境部署

CheckV 的安装相对简单，但其数据库较大，需要提前下载。

# 1. 创建并激活Conda环境
mamba create -n checkv_env -c conda-forge -c bioconda checkv
conda activate checkv_env

# 2. 下载`CheckV`数据库（首次使用需要，约16GB，耗时较长，请确保网络稳定和磁盘空间充足）
# 将 /path/to/your/checkv_db_directory 替换为您希望存放数据库的真实路径
checkv download_database /path/to/your/checkv_db_directory

• 软件核心用法

CheckV 的核心命令是 end_to_end，它将病毒Contigs作为输入，并自动完成所有评估和纯化步骤。

• 示例数据准备：

假设我们上一篇教程中，通过VirSorter2和VIBRANT合并去重后，得到了一个名为 merged_viral_contigs.fasta 的文件，其中包含了所有初步鉴定的病毒Contigs。

# 运行CheckV全流程
checkv end_to_end \
    merged_viral_contigs.fasta \           # 输入文件：初步鉴定的病毒Contigs
    /path/to/checkv_output \               # 输出目录
    -d /path/to/your/checkv_db_directory \ # CheckV数据库路径
    -t 16                                  # 使用的线程数

• 输出数据与结果解析

CheckV 会在指定的输出目录 /path/to/checkv_output 中生成多个文件。其中最关键的文件包括：

quality_summary.tsv：这是最重要的结果文件，以表格形式汇总了每个输入Contig的质量评估信息。

• 示例数据（quality_summary.tsv部分内容）：

contig_id contig_length provirus provirus_integrated provirus_prophage provirus_nuc_perc provirus_genes_perc completeness contamination warnings checkv_quality miuvig_quality
contig_1001 25432 Yes No No 90.5 92.1 Partial 0.0 NA Low-quality NA
contig_1002 45899 No NA NA NA NA Complete 0.0 NA Complete NA
contig_1003 15678 Yes Yes Yes 98.1 99.0 Complete 5.2 NA Medium-quality NA
contig_1004 8976 No NA NA NA NA Partial 0.0 NA Low-quality NA
contig_1005 120050 No NA NA NA NA High-quality 0.0 NA High-quality NA

• 结果解析：

contig_id：原始Contig的ID。

contig_length：Contig长度。

provirus：是否被识别为前病毒（Yes/No）。

completeness：估计的基因组完整度（例如：Partial, Complete, High-quality等）。

contamination：宿主序列污染百分比。

checkv_quality：CheckV基于完整度和污染度给出的最终质量分级。

miuvig_quality：根据MIUViG（Minimum Information about an Uncultivated Virus Genome）标准的分级。

• 如何筛选正常结果：

通常，我们会根据checkv_quality或miuvig_quality字段进行筛选。例如：

高质量vOTUs：选择checkv_quality为 "Complete", "High-quality", "Medium-quality" 的Contigs。或者更严格地，只选择completeness在90%以上且contamination低于5%的。

注意前病毒：provirus字段为“Yes”的Contigs可能是整合到宿主基因组中的病毒。根据研究目的，可以选择保留或单独分析它们。

viruses.fna：经过纯化处理后的游离病毒Contigs序列（通常移除了两端的宿主序列）。这是我们后续聚类的主要输入。

proviruses.fna：鉴定为前病毒的序列。

contigs_cut.fna：如果检测到宿主污染并进行了剪切，这个文件会包含被剪切后的病毒Contigs。

通过CheckV的评估和筛选，我们得到了一个更纯净、质量更高的病毒Contigs集合，为下一步的去冗余操作做好了准备。

病毒基因组聚类与去冗余 (Viral Genome Clustering & Dereplication)

此步骤的目的是将通过CheckV筛选出的高质量病毒Contigs进行聚类，从而获得代表病毒“物种”水平的非冗余操作分类单元（vOTUs）。

软件工具：dRep (主流推荐) / CD-HIT (备选方案)

dRep (De-replication)

是目前处理基因组去冗余和聚类的黄金标准工具。它不仅仅是基于简单的序列相似性，而是通过计算平均核苷酸一致性 (Average Nucleotide Identity, ANI) 来对基因组进行聚类。对于细菌，95% ANI通常被认为是物种划分的阈值；对于病毒，这个阈值可能需要根据具体病毒类群进行调整，但95% ANI也是常用的物种水平阈值。

• 核心作用与优势： ○ 基于ANI聚类：比简单的序列相似性（如BLAST）更能准确反映基因组间的进化距离，从而更好地划分物种。 ○ 自动选取代表基因组：在每个聚类簇中，dRep能根据基因组质量（如完整度、长度等）自动选择一个最佳的代表基因组。 ○ 灵活的阈值设置：允许用户自定义ANI阈值。
• 环境部署

# 1. 创建并激活Conda环境
mamba create -n drep_env -c conda-forge -c bioconda drep
conda activate drep_env

• 软件核心用法

示例数据准备：

我们将使用从CheckV结果中筛选出的高质量病毒Contigs文件作为输入。假设我们已经提取了所有 checkv_quality 为 "Complete", "High-quality", "Medium-quality" 的Contigs，并合并为一个文件 checkv_hq_viruses.fasta。

# 运行dRep进行病毒基因组去冗余
dRep dereplicate /path/to/drep_output \ # 输出目录
    -g /path/to/checkv_hq_viruses.fasta \ # 输入文件：CheckV筛选出的高质量病毒Contigs
    -p 16 \                               # 使用的线程数
    -sa 0.95 \                            # (重要) 设置ANI聚类阈值为95%（病毒物种常用）
    --ignoreGenomeQuality                 # (重要) 忽略dRep内置的基因组质量评估，我们已用CheckV完成

参数实践建议：

-sa 0.95：代表ANI相似度达到95%的基因组被认为属于同一物种。这个阈值在病毒学中是一个常用的经验值，但对于某些病毒类群，可能需要根据具体研究调整（例如，对于RNA病毒或进化速率极快的病毒，可能需要更低的阈值）。

--ignoreGenomeQuality：dRep本身也包含基因组质量评估功能（基于单拷贝核心基因），但这主要针对细菌。对于病毒，CheckV的评估更专业和全面。因此，我们通过此参数告知dRep直接使用我们输入的基因组进行聚类，无需再进行内部的质量过滤。

• 输出数据与结果解析

dRep会在指定的输出目录 /path/to/drep_output 中生成多个文件和子目录。最核心的结果是：

1. dereplicated_genomes/ 目录：

这个目录下包含了每个vOTU的代表序列文件（FASTA格式）。每个文件对应一个最终的vOTU。这些就是我们从大量原始Contigs中提炼出来的、非冗余的、高质量的病毒“物种”基因组。

1. Cdb.csv：

这是一个重要的表格文件，记录了聚类的详细信息，即每个原始输入基因组（genome列）被分到了哪个vOTU簇（cluster列），以及该簇的代表基因组是哪个（primary_cluster_represent列）。

示例数据（Cdb.csv部分内容）：

genome,cluster,primary_cluster_represent,completeness,contamination,length,N50,mean_coverage
contig_1001,cluster_001,contig_1001,Complete,0.0,45899,25000,10x
contig_1005,cluster_001,contig_1001,High-quality,0.0,42000,20000,8x
contig_2003,cluster_002,contig_2003,Complete,0.0,60000,30000,15x
contig_3010,cluster_003,contig_3010,Medium-quality,0.0,18000,10000,5x

• 结果解析：

genome：原始输入checkv_hq_viruses.fasta中的Contig ID。

cluster：该Contig所属的聚类簇ID。

primary_cluster_represent：该聚类簇的代表基因组ID。

其他列（completeness, contamination等）是dRep在运行时从输入基因组名称或从辅助文件（如CheckV的摘要文件）中提取的质量信息，用于辅助选择代表序列。

• 如何检查结果是否正常： ○ 查看聚类大小：检查Cdb.csv中每个cluster有多少个genome。如果有很多非常大的簇（数百甚至上千个基因组），可能ANI阈值设置过低；如果大部分簇只有一个基因组，可能ANI阈值设置过高，或者数据本身多样性非常高。 ○ 代表序列质量：检查dereplicated_genomes目录中的代表序列，它们应该都是CheckV筛选过的高质量Contigs。

log_ani.csv：记录了簇内所有基因组与代表基因组之间的ANI值，可以用来进一步验证聚类的合理性。

CD-HIT (Cluster Database at High Identity with Tolerance)

CD-HIT 是一款广泛使用的序列聚类工具，基于局部序列相似性，适用于DNA和蛋白质序列。虽然不如dRep在基因组层面精确，但其速度快，对于快速去冗余或初步筛选仍有价值。

• 核心作用与优势：速度快、内存占用低，适合处理超大规模数据集。
• 环境部署：mamba install -c bioconda cd-hit
• 软件核心用法：

cd-hit-est \
    -i /path/to/checkv_hq_viruses.fasta \ # 输入文件
    -o /path/to/cdhit_output.fasta \      # 输出去冗余后的代表序列
    -c 0.95 \                             # 设置序列相似性阈值（95%）
    -aS 0.9 \                             # (重要) 较短序列与较长序列的比对长度覆盖率阈值
    -d 0 \                                # 序列描述信息不截断
    -T 16                                 # 线程数

• 参数实践建议：

-c 0.95: 表示95%的核苷酸相似性。

-aS 0.9: 表示短序列必须至少有90%的长度与长序列比对上才被认为是同源的。这对于避免短片段与长基因组的局部高相似性造成的误聚类很重要。

• 输出数据与结果解析：

CD-HIT 的结果相对简单，主要用于快速获得一个非冗余的序列集。但由于它不是基于ANI，对于基因组级别的“物种”定义不如dRep精确，且无法自动评估和选择最佳代表序列。因此，**对于病毒vOTUs的构建，强烈推荐使用dRep**。

cdhit_output.fasta：包含每个聚类簇的代表序列。

cdhit_output.fasta.clstr：详细描述了每个聚类簇的组成，包括哪些原始序列被聚到了一起。

本章小结与展望

通过本篇文章的学习，我们已经掌握了从初步鉴定的病毒Contigs到构建高质量vOTUs的完整流程：

1. 利用 CheckV 对病毒Contigs进行严格的质量评估，识别并移除了宿主污染，区分了游离病毒和前病毒，并对病毒基因组的完整度进行了量化分级。
2. 接着，我们使用 dRep（或者备选的 CD-HIT）对通过质量控制的病毒Contigs进行聚类，基于平均核苷酸一致性（ANI）构建了代表病毒“物种”水平的非冗余vOTUs集合。

现在，我们手中拥有了一份可靠、高质量且具有生物学意义的病毒基因组列表。在下一篇文章中，我们将聚焦于这些vOTUs的身份识别和功能探索，包括病毒物种分类注释、病毒基因组功能注释，以及备受关注的病毒宿主预测，以进一步揭示它们在宏基因组中的生态位和作用。