1. 项目概述从宏基因组中揪出病毒的“侦探工具”在宏基因组学的研究里我们常常面对的是一个由无数微生物DNA/RNA片段混合而成的“超级大杂烩”。这里面有细菌、古菌、真核微生物当然还有我们这次要寻找的主角——病毒。病毒尤其是那些感染细菌的噬菌体它们个头小、基因组成分复杂多变而且很多没有通用的保守标记基因比如细菌的16S rRNA想从海量的测序数据里把它们准确识别出来就像在沙滩上找几颗特定形状的沙子。VirSorter2就是目前学界公认的、用来干这活儿的一把“尖刀”。它不是一个简单的序列比对工具而是一个集成了多种机器学习模型和大量病毒特征数据库的综合性检测流程。简单说它通过“以毒攻毒”的方式用已知的病毒特征去“嗅探”未知数据中的病毒序列无论是完整的病毒基因组还是整合在宿主基因组里的前病毒都难逃它的法眼。对于从事微生物生态、病毒组学、环境科学甚至医学相关研究的同行来说熟练掌握VirSorter2就等于拿到了打开病毒世界大门的钥匙。接下来我就以一个过来人的身份带你从零开始搞定它的安装并分享一些实战中真正好用的技巧和避坑指南。2. 环境准备与安装打好地基避免后续“楼塌了”安装VirSorter2本身不复杂但它的运行依赖于一个正确配置的Python环境和一系列生物信息学软件。很多新手卡在这一步不是因为命令难而是环境没弄干净。2.1 Conda环境为VirSorter2打造专属“单间”强烈建议使用CondaAnaconda或Miniconda来管理环境。这能完美解决软件依赖冲突这个老大难问题。别直接在系统环境或你的基础环境里安装否则将来升级其他工具时很可能“打架”。# 1. 创建一个新的conda环境命名为virsorter2名字可自定并指定python版本为3.8-3.10之间3.7太老3.11可能有不兼容风险 conda create -n virsorter2 python3.9 -y # 2. 激活这个环境 conda activate virsorter2激活后你的命令行提示符前面通常会显示(virsorter2)这表明你已经进入了这个专属环境。之后的所有操作都请确保在这个环境下进行。2.2 安装VirSorter2两种主流方式官方推荐通过pip安装这是最直接的方法。# 使用pip安装最新版的VirSorter2 pip install virsorter安装完成后可以验证一下virsorter --version如果成功显示版本号如2.2.4恭喜你核心工具安装成功了。但先别急着跑它只是个“指挥官”还需要“士兵”数据库和依赖工具。注意有时网络问题会导致pip安装缓慢或失败。可以尝试使用国内镜像源例如pip install virsorter -i https://pypi.tuna.tsinghua.edu.cn/simple2.3 数据库下载装备你的“病毒特征库”VirSorter2的识别能力很大程度上取决于它的数据库是否完整。首次运行任何命令时它会自动检查并尝试下载数据库但自动下载经常因为网络权限问题失败。因此手动下载并配置是更稳妥的做法。# 1. 设置数据库存放目录选择一个空间充足的路径比如~/db/virsorter2 export VS2DBDIR~/db/virsorter2 mkdir -p $VS2DBDIR # 2. 使用VirSorter2自带的脚本下载数据库 # 这一步会下载几个G的数据需要耐心等待并确保网络通畅。 virsorter setup --db-dir $VS2DBDIR --force--force参数会强制重新下载。如果中途断网可以重新执行此命令它会自动续传。下载的内容主要包括病毒蛋白家族pVOGs数据库用于基于隐马尔可夫模型HMM的搜索。病毒基因组的聚类信息用于序列相似性比较。其他分类和模型文件。下载完成后务必将VS2DBDIR这个环境变量永久添加到你的shell配置文件中如~/.bashrc或~/.zshrc并执行source命令使其生效这样每次打开终端都能找到数据库。echo export VS2DBDIR~/db/virsorter2 ~/.bashrc source ~/.bashrc2.4 依赖工具检查确保“士兵”到位VirSorter2在运行过程中会调用一些外部工具最主要的是MMseqs2用于超快速的序列搜索和聚类和Prodigal用于基因预测。通过pip安装VirSorter2时通常会自动安装这些Python包的 wrapper但有时系统级的二进制文件更稳定。# 在virsorter2的conda环境中直接安装这些工具 conda install -c bioconda mmseqs2 prodigal -y安装后检查一下mmseqs version prodigal -v确保都能正确输出版本信息。至此软件和数据库的安装才算真正完成。3. 核心原理与参数解析理解它在做什么在动手跑数据之前花几分钟理解VirSorter2的工作流程能让你在结果解读和问题排查时更加得心应手。它不是个黑箱。3.1 工作流程概览VirSorter2对一个输入文件通常是宏基因组组装的contig序列fasta格式的处理可以简化为以下几步基因预测使用Prodigal对每条contig进行编码基因CDS的预测。特征提取对预测出的基因进行多维度特征提取包括同源性特征使用MMseqs2将基因比对到病毒蛋白数据库如pVOGs看它像不像已知的病毒基因。基因密度与长度特征病毒基因组通常基因密度高基因间区短。基因方向与排列特征某些病毒有特定的基因排列模式。“宿主样”基因特征检测是否含有细菌等宿主特有的基因如核糖体蛋白基因这是判断是否为污染或前病毒的重要依据。机器学习分类将提取到的特征向量输入到预训练好的随机森林Random Forest分类模型中。这个模型已经用大量已知的病毒和非病毒序列训练过。打分与分类模型对每条contig输出一个得分0-1之间并根据得分和额外规则如是否位于片段末端将其分类到不同的类别中。输出结果给出每条contig的详细分类信息、得分、以及被识别出的病毒基因标记。3.2 关键运行参数详解了解核心参数才能定制化你的分析。以下是virsorter run命令最常用的一些参数virsorter run \ -i input.fasta \ # 输入文件组装的contigs -o output_dir \ # 输出目录会自动创建 --min-length 1500 \ # 最小序列长度短于该值的contig被忽略默认1500bp可调 --include-groups dsDNAphage,ssDNA,lavidaviridae,RNA \ # 指定检测的病毒组默认是dsDNAphage双链DNA噬菌体 --seqname-suffix-off \ # 关闭在输出序列名后添加后缀保持原名方便后续追踪 --prep-for-dramv \ # 为下游工具DRAM-v优化输出格式强烈推荐开启 -w dir_for_intermediate_files \ # 指定中间文件目录默认在输出目录下可单独设置以管理空间 -j number_of_threads \ # 使用的CPU线程数加快分析速度 --use-conda-off \ # 如果已手动安装所有依赖关闭conda自动管理依赖避免冲突 --db-dir $VS2DBDIR # 指定数据库路径如果环境变量已设置可省略参数选择心得--min-length病毒基因组通常大于1.5kb设置太低会引入大量噪音太高可能漏掉小病毒。对于高质量组装保持1500或提高到3000都是常见选择。--include-groups这是最容易出错的地方之一。如果你的样本可能包含真核病毒比如从人肠道、海洋浮游生物中提取的必须加上dsDNAphage,ssDNA,lavidaviridae,RNA,dsDNA,ssDNA,retro等更全的组。默认只找噬菌体会漏掉很多真核病毒--prep-for-dramv这个选项务必开启。它会生成一个*_for-dramv.fasta文件和一个*_for-dramv.tsv注释表这两个文件是下游进行病毒基因组功能注释使用DRAM-v的完美输入能省去大量格式转换的麻烦。4. 完整实操流程从数据到结果假设我们有一个名为metagenome_contigs.fasta的组装结果文件现在我们来完整跑一遍流程。4.1 数据准备与质控VirSorter2的输入是组装的contigs/scaffolds。在运行前建议先对组装文件进行简单的质控去除过短序列虽然VirSorter2有--min-length参数但事先用seqkit或bbmap过滤掉极短的序列如500bp可以减少不必要的计算量。确保序列ID简单避免在序列ID中使用空格、|、.等特殊字符最好只用字母、数字和下划线。复杂的ID有时会导致后续脚本解析出错。检查文件格式确保是标准的FASTA格式。# 使用seqkit过滤并简化序列名示例 conda install -c bioconda seqkit -y seqkit seq -m 1500 metagenome_contigs.fasta | seqkit replace -p .* -r metagenome_contigs_filtered.fasta4.2 运行VirSorter2在终端中执行以下命令。这里我们使用12个线程并检测所有主要病毒类型。# 激活环境并设置数据库路径如果已永久设置则无需重复 conda activate virsorter2 export VS2DBDIR~/db/virsorter2 # 运行VirSorter2 virsorter run \ -i metagenome_contigs_filtered.fasta \ -o ./virsorter2_results \ --min-length 1500 \ --include-groups dsDNAphage,ssDNA,lavidaviridae,RNA,dsDNA,ssDNA,retro \ --seqname-suffix-off \ --prep-for-dramv \ -w ./virsorter2_workdir \ -j 12 \ --use-conda-off运行时间取决于数据量大小和线程数。一个包含10万条contig的数据集在12核服务器上可能需要数小时。屏幕上会滚动显示当前进度。4.3 结果文件解读运行结束后进入./virsorter2_results目录你会看到一系列文件其中最重要的有以下几个final-viral-score.tsv核心结果文件。TSV格式包含每条contig的详细分类信息。seqname: 序列ID。length: 序列长度。hallmark: 检测到的病毒标志性基因数量。viral: 预测为病毒基因的数量。cellular: 预测为宿主细胞基因的数量。dsDNAphage,ssDNA, etc.: 在各病毒组中的得分。max_score: 所有组中的最高得分。max_score_group: 取得最高得分的病毒组。category:分类类别最关键1: 确定的病毒序列有病毒标志性基因。2: 可能的病毒序列无标志性基因但特征高度疑似。3: 整合的前病毒区域位于细菌contig内部的部分。45: 分别是类别1和2的片段位于contig末端。6: 可疑序列得分低可能是假阳性。final-viral-combined.fa所有被预测为病毒类别1, 2, 4, 5的contig序列集合。这是你初步得到的“病毒基因组”库。*_for-dramv.fasta和*_for-dramv.tsv为DRAM-v工具准备好的文件和注释表。4.4 结果筛选与后续分析通常我们会根据category和max_score进行筛选。一个常见的保守策略是只保留类别1和2即确定的/可能的完整病毒序列。# 使用awk从结果文件中提取类别1和2的序列名 awk -F\t $11 ~ /^(1|2)$/ {print $1} final-viral-score.tsv high_confidence_viral_contigs.list # 使用seqkit根据名单提取序列 seqkit grep -f high_confidence_viral_contigs.list final-viral-combined.fa high_confidence_viral.fasta得到的high_confidence_viral.fasta就可以用于下游分析了比如去冗余与聚类使用cd-hit-est或MMseqs2对病毒contigs进行聚类以代表物种水平。功能注释使用DRAM-v专门为病毒设计的注释流程对病毒基因组进行功能预测包括辅助代谢基因AMGs的鉴定。分类学注释使用vContact2或CAT/BAT等工具进行病毒的分类学归属。宿主预测使用iPHoP、VirHostMatcher等工具预测这些病毒的潜在宿主。5. 高级技巧与避坑指南这部分是文档里不会写但实战中能让你效率倍增、避免翻车的关键。5.1 处理大型宏基因组数据集当你的contig数量达到数十万甚至百万级别时直接运行可能会内存不足或耗时极长。分批次处理将大的fasta文件拆分成多个小文件例如每个文件5万条contig分别运行VirSorter2最后合并结果。注意合并时需要重新统一处理序列ID的重复问题。# 使用seqkit拆分示例 seqkit split -s 50000 metagenome_contigs.fasta充分利用-w参数将中间文件目录-w指向一个高速、大容量的存储位置如SSD或本地硬盘避免因网络存储NFS的I/O延迟拖慢速度。内存监控VirSorter2在比对阶段MMseqs2比较吃内存。确保服务器有足够物理内存如64G以上。可以使用top或htop命令监控进程。5.2 分类结果的可视化与解读直接看TSV文件不直观。可以用R或Python进行快速可视化。# R语言示例绘制病毒序列长度分布和得分分布 library(ggplot2) library(dplyr) data - read.delim(final-viral-score.tsv) # 筛选高置信度病毒 high_conf - data %% filter(category %in% c(1,2)) ggplot(high_conf, aes(xlog10(length), fillmax_score_group)) geom_histogram(bins50, alpha0.7) theme_minimal() labs(titleViral Contig Length Distribution, xlog10(Length bp), yCount) ggplot(high_conf, aes(xmax_score)) geom_density(fillsteelblue, alpha0.5) theme_minimal() labs(titleDistribution of Max VirSorter2 Scores, xScore, yDensity)通过可视化你可以快速判断病毒序列的长度集中范围、不同病毒组的比例以及得分分布是否合理通常希望高置信度序列得分集中在0.7-1.0。5.3 常见报错与解决方案ERROR: database folder is empty or does not exist原因数据库路径VS2DBDIR未设置或设置错误或者数据库没有成功下载。解决检查echo $VS2DBDIR输出是否正确。手动进入该目录查看是否有pfam、viral等子文件夹。如果没有重新运行virsorter setup。[Errno 28] No space left on device原因中间文件目录-w指定或默认在输出目录下所在的磁盘空间已满。解决清理磁盘空间或将-w指向一个空间充足的磁盘分区。VirSorter2运行产生的中间文件可能比原始输入大10倍以上。Command ‘[‘conda‘, ‘run‘, ‘-n‘, ‘virsorter‘ ...]’ returned non-zero exit status 127.原因Conda环境问题。可能是名为virsorter的conda环境不存在或者--use-conda-off参数未使用但系统conda配置有问题。解决最稳妥的方法是确保已手动安装所有依赖mmseqs2, prodigal并在运行命令中明确加上--use-conda-off参数。运行极其缓慢卡在某个步骤可能原因通常是MMseqs2比对步骤。检查是否使用了网络存储NFS。MMseqs2的临时I/O非常频繁网络延迟是性能杀手。解决使用-w参数将工作目录指向本地硬盘。同时确保$TMPDIR环境变量也指向本地高速存储。5.4 与CheckV联用评估病毒基因组完整性VirSorter2预测出的病毒contig可能是完整的也可能是片段。CheckV是评估病毒基因组完整性和纯度的黄金标准工具。强烈建议将VirSorter2的输出作为CheckV的输入。# 假设已安装CheckV (conda install -c bioconda checkv) checkv end_to_end high_confidence_viral.fasta ./checkv_output -d /path/to/checkv_database -t 12CheckV会给出每条序列的“完整性”百分比并识别出序列两端的宿主污染直接终端重复序列Direct Terminal Repeats。最终你可以得到一个质量分级的病毒基因组集合用于发表级分析。6. 性能优化与集群部署对于超大规模项目在服务器集群上运行是常态。使用Snakemake或Nextflow编写流程将VirSorter2、CheckV、去冗余等步骤封装成可复现的流程方便并行和重跑。Slurm作业提交示例#!/bin/bash #SBATCH --job-namevirsorter2 #SBATCH --outputlogs/virsorter_%j.out #SBATCH --errorlogs/virsorter_%j.err #SBATCH --time48:00:00 #SBATCH --mem100G #SBATCH --cpus-per-task24 conda activate virsorter2 export VS2DBDIR/shared/db/virsorter2 virsorter run \ -i $INPUT_FASTA \ -o $OUTPUT_DIR \ --min-length 1500 \ --include-groups dsDNAphage,ssDNA,lavidaviridae,RNA,dsDNA,ssDNA,retro \ --seqname-suffix-off \ --prep-for-dramv \ -w $TMPDIR \ -j $SLURM_CPUS_PER_TASK \ --use-conda-off这里的关键是将-w指向$TMPDIR通常是节点本地的高速临时存储并申请足够的内存和CPU。数据库共享在集群上将庞大的数据库$VS2DBDIR放在共享存储如Lustre, NFS上所有计算节点都能访问避免每个任务重复下载。7. 下游分析串联实例从VirSorter2到生态解读单独看病毒序列列表意义有限必须与上下游分析结合。这里分享一个简单的串联分析思路。VirSorter2筛选得到high_confidence_viral.fasta。CheckV评估得到完整和高质量的病毒基因组checkv_quality_summary.tsv筛选完整性90%且污染5%的序列。去冗余使用cd-hit-est在95%相似度、85%覆盖度下聚类取最长序列为代表。cd-hit-est -i viral_complete.fasta -o viral_reps.fasta -c 0.95 -aS 0.85 -G 0 -M 0 -T 24DRAM-v注释对代表序列进行功能注释挖掘AMGs。dram-v.py annotate -i viral_reps.fasta -o dramv_annotation --threads 24 dram-v.py distill -i dramv_annotation/annotations.tsv -o dramv_distillvContact2网络分析基于基因共享网络对病毒进行聚类近似属/种水平并可能与已知的病毒参考基因组关联。与宿主关联如果同时有宏基因组组装基因组MAGs可以使用CRISPR spacer匹配或序列相似性如通过blastn将病毒与潜在宿主关联起来。统计分析将病毒在不同样本中的丰度通过read mapping计算、与宿主的关联信息、携带的AMG功能等与环境因子进行关联分析如RDA、Mantel检验最终阐释病毒在生态系统中的潜在作用。整个流程看似步骤繁多但一旦用流程管理工具如Snakemake串联起来就能实现自动化、可复现的分析极大提升研究效率。VirSorter2作为这个流程的起点其输出的准确性和完整性直接决定了下游所有分析的可靠性。因此花时间理解它的原理、优化它的参数、妥善处理它的结果是每一个病毒组学研究者值得投入的必修课。