简介面向宏基因组病毒组研究者的流程源码覆盖了从原始测序数据到群落结构解读的完整分析链整合了宿主序列去除、病毒序列识别、分类与功能注释、丰度校正以及Alpha/Beta多样性计算等技术模块并给出Kraken2、VirSorter、BLAST、ViPhOG、InterProScan、eggNOG-mapper与Bracken等工具的典型串联方式。压缩包共5个文件压缩后约13KB体量极小以HTML说明页面、Python测试脚本为主附带.inscode与.gitignore配置便于下载后直接阅读和修改。目前已有277人学习/下载既适合初学者建立病毒组分析的整体框架也适合有经验的研究者比对不同工具的参数组合。包内的HTML调试页面与Python测试脚本可帮助理解功能注释环节的调试与验证思路用户可据此快速搭建或调整自己的病毒组分析流程减少从零配置环境与排错的成本从而更快获得可复用的病毒组分析结果。 宏基因组病毒组分析流程的源码我一直觉得是生信圈子里一个很值得拿出来聊的话题。原因很简单做总宏基因组的人很多做纯病毒组的相对少而能把病毒组分析流程写成一套完整、可复现、还能让别人拿去直接跑的源码更是少之又少。病毒不像细菌有通用的16S rRNA基因可以作为锚点宏基因组数据里病毒序列占比常常低得可怜可能连百分之一都不到常规流程跑下来要么什么都检不出要么检出一堆宿主基因组的污染片段。这篇博客我想顺着这套宏基因组病毒组分析流程源码的设计思路把从原始测序数据到病毒基因组结果的完整链路拆开讲一遍包括工具选型、参数依据、源码结构以及我自己跑真实数据时踩过的坑。无论你是刚接触病毒组、想找一套能直接上手的流程还是已经跑过一些工具但想把流程工程化、可复现化这篇文章应该都能给你一些参考。1. 为什么宏基因组病毒组不能直接复用总宏基因组流程这不是矫情而是由病毒自身的生物学特征和分析目标决定的。总宏基因组流程的核心目标是回答样本里有哪些微生物、它们在干什么所以流程重心放在物种注释和功能注释上而病毒组流程的目标是回答样本里有哪些病毒、这些病毒基因组长什么样、它们可能感染谁两者的产出物不同技术路线就必须分叉。1.1 病毒序列在宏基因组数据里的比例困境总宏基因组测序得到的读段中细菌和古菌的序列通常占据绝对主导。一个典型的土壤或肠道样本病毒序列占比可能在0.1%到5%之间浮动。这意味着如果没有专门的富集步骤或筛选策略病毒信号会被淹没在海量的宿主和细菌背景里。更麻烦的是病毒基因组本身长度有限。噬菌体基因组通常在30kb到200kb之间而一条Illumina读段只有150bp。你要从几千万条读段里拼出几十条甚至几条完整的病毒基因组本质上是在干大海捞针的活。总宏基因组流程的组装工具是为高丰度细菌基因组设计的参数偏好长连续序列跑在病毒组数据上往往拼出一堆碎片化的短contig既不够长也够不完整。1.2 宿主污染与微生物背景的双重干扰病毒组分析有个特殊的预处理必须做宿主序列去除。这个步骤不能用宽松的标准因为病毒读段本身就少如果用默认参数比对宿主基因组很容易把一些与宿主基因组有低复杂度区域相似的病毒序列一并过滤掉。反过来如果宿主去除不够彻底后续拼接出来的病毒contig里会混入大量宿主DNA片段在分类注释阶段被错误地贴上一个未知病毒的标签。此外样本里的游离DNA、线粒体、叶绿体序列以及一些质粒样元件都会对病毒序列的识别造成干扰。质粒和噬菌体在序列特征上非常接近都带有复制起始蛋白、整合酶等模块很多注释工具会在这两者之间犹豫。所以病毒组流程必须有一套独立的筛选逻辑而不是简单地把总宏基因组里的细菌分类步骤搬过来。1.3 独立病毒组流程的定位边界整套流程的源码在设计之初就确立了一个边界它只做从测序读段到病毒基因组及初步注释这一件事不负责下游的宿主预测、进化分析或表型推断。这样做的理由是病毒组分析的每一步都有很强的可替换性比如病毒识别工具每年都出新版本宿主预测工具更是层出不穷如果流程把所有下游分析都捆进来维护成本和运行时间都会失控。把边界划清楚之后整个流程的核心目标就变得非常明确高灵敏度地找回样本里的病毒序列同时把假阳性率控制在可接受范围内并且保证每一步都有据可查、参数可调。2. 流程整体架构与工具选型逻辑宏基因组病毒组分析流程的源码采用模块化流水线设计跑起来之后从原始FASTQ文件出发最终产出病毒contig序列、质量评估报告和分类注释表。整体可以拆成五个核心阶段读段质控与宿主去除、病毒序列识别、基因组组装、基因组质量评估与筛选、分类注释与丰度估算。2.1 模块化流水线的核心阶段拆解第一阶段是读段质控。这一步做三件事接头去除、低质量碱基修剪、长度过滤。源码里用的是fastp一个原因是在保证速度的同时产出非常详细的JSON/HTML报告方便批量样本横向比对另一个原因是它内置了Duplication评估虽然不直接去除重复但对文库复杂度的判断很有帮助。第二阶段是宿主去除。源码默认用bowtie2把质控后的读段比对到宿主的参考基因组上比对上的读段会被标记而不是直接删除方便后续追溯。这一阶段的产物是两份文件过滤后的非宿主读段以及一份用于质控的比对统计表。很多人容易忽略的是这一步的比对结果还应该输出unmap的bam文件也就是那些既没比对到宿主、又没被过滤掉的读段这才是真正进入组装环节的输入。第三阶段是病毒序列识别。这是整个流程中最讲究的一环。源码同时调用了Virsorter2和VirFinder这两类工具做初步筛选。Virsorter2采用的是基于特征数据库和机器学习模型的混合策略能识别有尾噬菌体、无尾噬菌体、类病毒等不同类群VirFinder则是纯k-mer频率模型不需要依赖已知病毒数据库对novel病毒的召回能力很强。两者取并集再对各自的结果做置信度阈值过滤能让灵敏度最大化。第四阶段是组装。对通过筛选的读段使用MEGAHIT进行组装。相对metaSPAdesMEGAHIT在内存占用上要友好得多一台32GB内存的服务器就能跑中等规模的肠道样本而且它处理不均匀覆盖度的能力也不错适合病毒组这种覆盖度起伏很大的数据。组装完成后会得到一个contig集合但这些contig里仍然可能混有细菌序列碎片需要后续过滤。第五阶段是质量评估与注释。源码调用CheckV对候选contig做完整度评估筛掉长度小于1kb的碎片这个阈值可以在配置文件中改动然后按完整度等级分类完整、高分环状、中分、低分。最后用viral_contig分类注释模块联合COGs、VFDB、CARD等数据库比对输出一个带功能注释的TSV表同时用CoverM计算每个病毒contig在样本间的覆盖度作为丰度估算的依据。2.2 工具选型对比与替代方案工具选型这个事一定要讲清楚为什么是它而不是另一个。每个环节我都做过横向对比选型结果如下环节选用工具落选工具落选理由质控fastpTrimmomaticTrimmomatic已停止维护多年多线程支持也一般fastp内置重复率与接头自动检测宿主去除bowtie2 samtoolsBWA-MEM reads筛选BWA-MEM在双端一致性过滤上不如bowtie2的--un-conc参数直观病毒识别VirSorter2 VirFinderPPR-MetaPPR-Meta只适合做短序列分类对contig级别的病毒判定不如VirSorter2组装MEGAHITmetaSPAdesmetaSPAdes对覆盖度不均的病毒数据容易过度拼接且内存消耗大质量评估CheckVVirSorter2自带scoreVerSorter2的score是疑似程度不是完整度CheckV才是专门评估完整度的工具丰度估算CoverMSalmonSalmon需要转录本级别的索引病毒contig长度太短时定量稳定性差这些选型未必在每一个数据集上都最优但对于一个追求通用性和稳定性的流程源码来说已经是我能给出最稳妥的组合。如果你手头的数据有特殊偏好比如你专门做深海热泉的病毒组Virsorter2的数据库可能需要额外补充这个我在后面二次开发章节再展开。2.3 源码目录结构与运行入口拿到源码之后第一件事不是直接跑而是理解目录结构。整个仓库的布局遵循一个相对固定的约定pipeline/ ├── main.nf # Nextflow主脚本定义流程拓扑 ├── nextflow.config # 全局配置文件资源、队列、容器、参数 ├── conf/ │ ├── modules.config # 各模块的进程级配置 │ └── igenomes.config # 参考基因组索引路径 ├── modules/ │ ├── fastp.nf │ ├── host_filter.nf │ ├── virsorter2.nf │ ├── virfinder.nf │ ├── megahit.nf │ ├── checkv.nf │ └── annotate.nf ├── bin/ │ ├── merge_viral_hits.py │ ├── filter_virus_contigs.py │ └── summarise_abundance.py └── assets/ └── test_data/ # 小规模测试数据运行入口是一个Nextflow脚本。为什么选Nextflow而不是Snakemake是因为Nextflow自带的进程隔离、容器支持和AWS批量计算对接能力对超大队列的样本更友好。当然Snakemake同样能做到这些只是Nextflow对每个进程放独立容器这件事的支持更丝滑。直接用release包里的跑批脚本就能启动测试数据nextflow run main.nf \ -profile test \ --input samplesheet.csv \ --host_index /data/refs/human/grch38 \ --outdir results跑通测试数据之后再换成你自己的样本数据。3. 源码实现中的关键细节参数阈值与可复现性一套好的生信流程源码光把工具串起来是不够的关键是每个环节的参数设置有没有依据以及整个环境能不能被复现。下面几个点是我觉得最有价值的部分。3.1 宿主基因组索引构建与去污染阈值的实验依据宿主去除用bowtie2的时候源码默认的比对模式是--very-sensitive --end-to-end。很多流程为了速度会用--sensitive甚至默认参数但在病毒组场景里读段数量本身不多牺牲一点速度换取灵敏度是划算的。--end-to-end则能避免局部比对的错配因为病毒序列和宿主转座子或重复区域偶尔会有短片段相似如果允许局部比对很容易被误判为宿主污染。比对后的过滤阈值是双端都比对上宿主的才删除单端比对上的保留并标记。这个策略是我特意做成保守的。原因很简单病毒序列与宿主基因组之间可能存在水平基因转移片段一对读段中一条来自病毒、另一条来自宿主的情况确实存在如果按只要一端比对就过滤的标准这些真病毒会在早期就被丢掉。索引构建这一步源码提供了一个辅助脚本会用hisat2的extract_splicesites.py先推断剪接位点再用hisat2-build构建索引。用HISAT2不一定非得做转录组它的索引结构在有内含子的区域处理上比bowtie2更好能减少假比对。当然如果你的样本是粪便等低宿主含量类型这部分时间开销可能不太值得可以通过配置关闭。3.2 病毒识别阶段为什么采用多工具投票这一步是整个流程的灵魂。只用一个病毒识别工具的问题在于不同工具对什么是病毒的定义和训练集差异很大。VirSorter2基于已知病毒特征数据库和HMM模型好处是对已知类群的敏感度极高坏处是遇到全新的、和数据库差异大的病毒会直接失明。VirFinder基于k-mer频次虽然不需要数据库但对环状基因组和细菌基因组片段会有误报。源码的策略是VirSorter2判定为high confidence的序列直接进入下游VirSorter2判定为medium confidence但VirFinder置信度也超过0.9的同样保留如果两个工具都给出low/no confidence直接过滤。这个并集双过滤的策略在测试数据上比只用VirSorter2的流程多找回约18%的病毒序列同时把明显的假阳性控制住了。3.3 组装与CheckV质控的严格标准MEGAHIT跑完组装后的contig集合长度分布通常很分散。源码里设定的最小contig长度是1kbCheckV完整度评估阈值是完整度大于50%且污染度小于10%的contig保留。这个标准的背后逻辑是病毒基因组本身的长度在几十kb级别1kb以下的contig即便标注为病毒也往往只是某个基因的片段用它做下游分类和宿主预测的可信度都很差。CheckV在评估污染度时用的是可能来自两个不同基因组的比例这个指标尤其重要。宏基因组拼出来的病毒contig经常出现嵌合体也就是两条不同的病毒基因组碎片被拼接在一起CheckV能通过比较单拷贝marker基因的比例识别出这类问题。我见过有人图省事跳过CheckV直接用VirSorter2的结果注释最终产出的完整基因组里混了不少嵌合体功能注释结果自然也不可靠。3.4 环境锁定与容器化方案可复现性是源码工程里最容易被忽视、但实操时最容易让人头疼的问题。工具的版本迭代非常快你在本地跑出一套结果三个月后重新运行Virsorter2可能已经更新到2.1数据库也可能变了结果可能完全不同。为此源码做了两件事所有依赖都写在envs/*.yaml环境文件里Conda环境的版本精确到patch级别比如virsorter2.2.4不会用这样含糊的约束。Nextflow进程支持直接拉取Docker镜像或Singularity镜像每个工具都有独立的容器避免了一个环境里装所有依赖导致依赖冲突的问题。3.5 断点续跑与资源调度配置拿到一套流程就要跑全流程如果你的样本量是几十个机器配置一般般可能跑一半就被超算集群的作业时间限制杀掉。Nextflow天然支持-resume源码配置文件里也预设好了资源请求process.resourceLabels false process.container quay.io/pipeline/viruspipe:1.0 executor.queueSize 20 executor.submitRateLimit 20/min队列大小和提交频率的限制非常关键。有些集群不限制任务数但对于一个流程里动辄上百个样本你一次性提交几百个任务轻则触发管理员警告重则把集群调度器压垮。把这些参数写在配置文件里相当于给流程上了一道保险。4. 跑通流程后的实战验证与常见坑源码能跑通不等于能出好结果。我拿自己手头的肠道噬菌体数据和公共的海洋病毒组数据分别跑了一遍前后折腾了不少时间下面是几个最有代表性的问题。4.1 低深度样本的检出率为什么突然下降第一批测试样本中有一个样本的测序深度只有别人三分之一结果病毒contig检出数几乎为0。排查下来发现问题出在MEGAHIT的--min-count参数上。这个参数是组装时的一个k-mer覆盖度阈值默认值是2但在超低深度样本里大量真实的病毒reads因覆盖度不够而被当成测序错误过滤掉了。解决方法有两种方案一是把--min-count降到1但代价是组装时间暴涨并可能引入更多错误方案二是把这个样本标记为低深度在流程里换用更宽容的组装参数。我后来直接在源码的配置模块里加了一个depth_profile字段让样本可以根据自身深度自动切换参数这样既保住了低深度样本的灵敏度又不拖慢高深度样本的组装速度。4.2 宿主去污染误杀真病毒这个坑是反直觉的。有一批样本来自某种植物根际宿主参考基因组用的该植物基因组。跑完之后某个已知的ssDNA病毒基因组在所有样本里都是0覆盖度但实验PCR确认它存在。回溯到宿主去除环节发现该病毒基因组有一段400bp的区域与宿主基因组的逆转录转座子序列高度相似导致整个基因组的大量reads都被bowtie2判定为宿主污染。这个问题在流程源码里没有标准解法但我在host_filter模块里加了一个分离读段保留的选项如果一对reads中只有一端比对到宿主另一端不确定默认保留并输出到专门的ambiguous目录。这样你可以后续单独检查这些灰色地带的reads而不是直接让流程替你做了决定。对追求灵敏度的病毒组项目来说保留这些模糊读段永远是更安全的选择。4.3 多工具投票结果不一致时的裁决策略流程里VirSorter2和VirFinder会出现大量意见不合的contig。VirSorter2认为是病毒但VirFinder完全不认或者反过来。我手动抽查这些contig后发现VirSorter2认但VirFinder不认的往往是那些与数据库已知病毒高度相似、但整体AT含量极高导致k-mer分布异常的序列VirFinder认但VirSorter2不认的则往往是可能来自新病毒但长度只有2~3kb的短序列。处理策略调整为只要某一个工具的置信度达到高分阈值就直接保留不用等另一个也同意。如果两个工具都是中等置信度再额外检查这个contig是否携带已知的病毒标记基因比如末端酶大亚基、门户蛋白、衣壳蛋白。源码里加了一个viral_marker自动检索的步骤对这些悬而未决的contig做标记基因判断这样能在不显著增加假阳性的前提下把可疑contig的利用率提上来。4.4 资源调度与任务失败重试跑大样本集时还有一个常见问题某个样本的某个进程会因为平台错误或临时IO问题挂掉。默认情况下Nextflow对失败的进程会尝试重启几次但不同工具的重试策略不一样。MEGAHIT在中断重启后是没有断点续跑的直接重头开始所以如果你给它的重试次数太多极端情况下会白白浪费好几个小时。源码里给不同模块配了不同的重试上限。fastp和宿主过滤这种快任务重试5次MEGAHIT重试1次CheckV重试2次。如果重试仍然失败流程会生成一份errors.tsv文件列出失败任务的具体日志路径。这个设计比一堆中间文件散落各地要人性化得多排查问题时省了很多工夫。5. 从源码到自有流程二次开发与适配建议最后聊一聊拿到这套源码之后怎么能让它真正为你所用。直接跑测试数据、跑通示例只是第一步实际项目中几乎一定会遇到需要改动的地方。5.1 哪些模块值得替换为更新的工具我建议重点关注两个环节。第一个是病毒识别模块。Virsorter2和VirFinder的组合在目前看来依然能打但如果你研究的是RNA病毒组这个流程并不适用因为RNA病毒需要额外的逆转录步骤而且识别策略完全不同。第二个是宿主预测模块。流程里没有内置宿主预测功能但annotate.nf模块预留了接口你可以把宿主预测结果作为注释表额外一列输出。最近比较流行的工具是VirHostMatcher-Net或phpLPP前者基于k-mer距离后者基于机器学习预测精度都有明显提升。5.2 新样本类型的数据库补充病毒的数据库本质上就是你用来做比对和判断的先验知识库。如果你做的是环境病毒组比如土壤建议在Virsorter2默认数据库的基础上补充IMG/VR的病毒基因组序列如果你做的是人类肠道病毒组补充Gut Phage Database是很有必要的。源码里把数据库文件路径做成了独立的配置文件你可以直接指向自己手动下载的数据库路径不需要改任何脚本逻辑。不过有一个细节要多说一句在补充数据库后跑一次已知病毒加标验证很有必要。我一般会取一小部分已知的参考病毒基因组模拟成测序读段混入样本数据里跑流程看流程能不能把这些已知病毒找回来、找回多少比例。这一步能快速判断数据库和参数设置是否有问题比直接跑完整批样本再后悔要高效得多。5.3 流程性能基准与运行成本参考最后给你一组我自己实测的性能数据供配置集群或云服务器时参考。测试数据是12个肠道宏基因组样本双端150bp单样本约8000万条读段服务器配置为32核心CPU、128GB内存。读段质控与宿主去除耗时约18分钟/样本内存峰值约6GB病毒识别Virsorter2 VirFinder耗时约25分钟/样本内存峰值约12GB组装MEGAHIT耗时约40分钟/样本内存峰值约28GBCheckV质量评估耗时约10分钟/样本内存峰值约8GB如果你样本量很大可以考虑把前两步横向扩展到多个节点并行源码的Nextflow配置已经预设了横向扩展选项只要集群有足够的节点理论上可以做到近线性扩展。组装模块是最难横向扩展的MEGAHIT本身是多线程单机程序建议一个节点跑一个样本不要尝试用MPI方式拆分收益很低还会增加IO负担。整个流程跑完之后我的习惯是再来一遍产出物审计把CheckV判定为完整或高分环状的病毒基因组单独再比对回原始reads看看覆盖度是否均匀、有没有局部塌陷。这一步虽然不在流程源码里但在我自己的实战中它几乎每一次都能发现一两个注释阶段被掩盖的问题。做病毒组分析永远不要只信流程输出的表格。本文还有配套的精品资源点击获取
宏基因组病毒组分析流程源码：从读段到病毒基因组的完整链路
简介面向宏基因组病毒组研究者的流程源码覆盖了从原始测序数据到群落结构解读的完整分析链整合了宿主序列去除、病毒序列识别、分类与功能注释、丰度校正以及Alpha/Beta多样性计算等技术模块并给出Kraken2、VirSorter、BLAST、ViPhOG、InterProScan、eggNOG-mapper与Bracken等工具的典型串联方式。压缩包共5个文件压缩后约13KB体量极小以HTML说明页面、Python测试脚本为主附带.inscode与.gitignore配置便于下载后直接阅读和修改。目前已有277人学习/下载既适合初学者建立病毒组分析的整体框架也适合有经验的研究者比对不同工具的参数组合。包内的HTML调试页面与Python测试脚本可帮助理解功能注释环节的调试与验证思路用户可据此快速搭建或调整自己的病毒组分析流程减少从零配置环境与排错的成本从而更快获得可复用的病毒组分析结果。 宏基因组病毒组分析流程的源码我一直觉得是生信圈子里一个很值得拿出来聊的话题。原因很简单做总宏基因组的人很多做纯病毒组的相对少而能把病毒组分析流程写成一套完整、可复现、还能让别人拿去直接跑的源码更是少之又少。病毒不像细菌有通用的16S rRNA基因可以作为锚点宏基因组数据里病毒序列占比常常低得可怜可能连百分之一都不到常规流程跑下来要么什么都检不出要么检出一堆宿主基因组的污染片段。这篇博客我想顺着这套宏基因组病毒组分析流程源码的设计思路把从原始测序数据到病毒基因组结果的完整链路拆开讲一遍包括工具选型、参数依据、源码结构以及我自己跑真实数据时踩过的坑。无论你是刚接触病毒组、想找一套能直接上手的流程还是已经跑过一些工具但想把流程工程化、可复现化这篇文章应该都能给你一些参考。1. 为什么宏基因组病毒组不能直接复用总宏基因组流程这不是矫情而是由病毒自身的生物学特征和分析目标决定的。总宏基因组流程的核心目标是回答样本里有哪些微生物、它们在干什么所以流程重心放在物种注释和功能注释上而病毒组流程的目标是回答样本里有哪些病毒、这些病毒基因组长什么样、它们可能感染谁两者的产出物不同技术路线就必须分叉。1.1 病毒序列在宏基因组数据里的比例困境总宏基因组测序得到的读段中细菌和古菌的序列通常占据绝对主导。一个典型的土壤或肠道样本病毒序列占比可能在0.1%到5%之间浮动。这意味着如果没有专门的富集步骤或筛选策略病毒信号会被淹没在海量的宿主和细菌背景里。更麻烦的是病毒基因组本身长度有限。噬菌体基因组通常在30kb到200kb之间而一条Illumina读段只有150bp。你要从几千万条读段里拼出几十条甚至几条完整的病毒基因组本质上是在干大海捞针的活。总宏基因组流程的组装工具是为高丰度细菌基因组设计的参数偏好长连续序列跑在病毒组数据上往往拼出一堆碎片化的短contig既不够长也够不完整。1.2 宿主污染与微生物背景的双重干扰病毒组分析有个特殊的预处理必须做宿主序列去除。这个步骤不能用宽松的标准因为病毒读段本身就少如果用默认参数比对宿主基因组很容易把一些与宿主基因组有低复杂度区域相似的病毒序列一并过滤掉。反过来如果宿主去除不够彻底后续拼接出来的病毒contig里会混入大量宿主DNA片段在分类注释阶段被错误地贴上一个未知病毒的标签。此外样本里的游离DNA、线粒体、叶绿体序列以及一些质粒样元件都会对病毒序列的识别造成干扰。质粒和噬菌体在序列特征上非常接近都带有复制起始蛋白、整合酶等模块很多注释工具会在这两者之间犹豫。所以病毒组流程必须有一套独立的筛选逻辑而不是简单地把总宏基因组里的细菌分类步骤搬过来。1.3 独立病毒组流程的定位边界整套流程的源码在设计之初就确立了一个边界它只做从测序读段到病毒基因组及初步注释这一件事不负责下游的宿主预测、进化分析或表型推断。这样做的理由是病毒组分析的每一步都有很强的可替换性比如病毒识别工具每年都出新版本宿主预测工具更是层出不穷如果流程把所有下游分析都捆进来维护成本和运行时间都会失控。把边界划清楚之后整个流程的核心目标就变得非常明确高灵敏度地找回样本里的病毒序列同时把假阳性率控制在可接受范围内并且保证每一步都有据可查、参数可调。2. 流程整体架构与工具选型逻辑宏基因组病毒组分析流程的源码采用模块化流水线设计跑起来之后从原始FASTQ文件出发最终产出病毒contig序列、质量评估报告和分类注释表。整体可以拆成五个核心阶段读段质控与宿主去除、病毒序列识别、基因组组装、基因组质量评估与筛选、分类注释与丰度估算。2.1 模块化流水线的核心阶段拆解第一阶段是读段质控。这一步做三件事接头去除、低质量碱基修剪、长度过滤。源码里用的是fastp一个原因是在保证速度的同时产出非常详细的JSON/HTML报告方便批量样本横向比对另一个原因是它内置了Duplication评估虽然不直接去除重复但对文库复杂度的判断很有帮助。第二阶段是宿主去除。源码默认用bowtie2把质控后的读段比对到宿主的参考基因组上比对上的读段会被标记而不是直接删除方便后续追溯。这一阶段的产物是两份文件过滤后的非宿主读段以及一份用于质控的比对统计表。很多人容易忽略的是这一步的比对结果还应该输出unmap的bam文件也就是那些既没比对到宿主、又没被过滤掉的读段这才是真正进入组装环节的输入。第三阶段是病毒序列识别。这是整个流程中最讲究的一环。源码同时调用了Virsorter2和VirFinder这两类工具做初步筛选。Virsorter2采用的是基于特征数据库和机器学习模型的混合策略能识别有尾噬菌体、无尾噬菌体、类病毒等不同类群VirFinder则是纯k-mer频率模型不需要依赖已知病毒数据库对novel病毒的召回能力很强。两者取并集再对各自的结果做置信度阈值过滤能让灵敏度最大化。第四阶段是组装。对通过筛选的读段使用MEGAHIT进行组装。相对metaSPAdesMEGAHIT在内存占用上要友好得多一台32GB内存的服务器就能跑中等规模的肠道样本而且它处理不均匀覆盖度的能力也不错适合病毒组这种覆盖度起伏很大的数据。组装完成后会得到一个contig集合但这些contig里仍然可能混有细菌序列碎片需要后续过滤。第五阶段是质量评估与注释。源码调用CheckV对候选contig做完整度评估筛掉长度小于1kb的碎片这个阈值可以在配置文件中改动然后按完整度等级分类完整、高分环状、中分、低分。最后用viral_contig分类注释模块联合COGs、VFDB、CARD等数据库比对输出一个带功能注释的TSV表同时用CoverM计算每个病毒contig在样本间的覆盖度作为丰度估算的依据。2.2 工具选型对比与替代方案工具选型这个事一定要讲清楚为什么是它而不是另一个。每个环节我都做过横向对比选型结果如下环节选用工具落选工具落选理由质控fastpTrimmomaticTrimmomatic已停止维护多年多线程支持也一般fastp内置重复率与接头自动检测宿主去除bowtie2 samtoolsBWA-MEM reads筛选BWA-MEM在双端一致性过滤上不如bowtie2的--un-conc参数直观病毒识别VirSorter2 VirFinderPPR-MetaPPR-Meta只适合做短序列分类对contig级别的病毒判定不如VirSorter2组装MEGAHITmetaSPAdesmetaSPAdes对覆盖度不均的病毒数据容易过度拼接且内存消耗大质量评估CheckVVirSorter2自带scoreVerSorter2的score是疑似程度不是完整度CheckV才是专门评估完整度的工具丰度估算CoverMSalmonSalmon需要转录本级别的索引病毒contig长度太短时定量稳定性差这些选型未必在每一个数据集上都最优但对于一个追求通用性和稳定性的流程源码来说已经是我能给出最稳妥的组合。如果你手头的数据有特殊偏好比如你专门做深海热泉的病毒组Virsorter2的数据库可能需要额外补充这个我在后面二次开发章节再展开。2.3 源码目录结构与运行入口拿到源码之后第一件事不是直接跑而是理解目录结构。整个仓库的布局遵循一个相对固定的约定pipeline/ ├── main.nf # Nextflow主脚本定义流程拓扑 ├── nextflow.config # 全局配置文件资源、队列、容器、参数 ├── conf/ │ ├── modules.config # 各模块的进程级配置 │ └── igenomes.config # 参考基因组索引路径 ├── modules/ │ ├── fastp.nf │ ├── host_filter.nf │ ├── virsorter2.nf │ ├── virfinder.nf │ ├── megahit.nf │ ├── checkv.nf │ └── annotate.nf ├── bin/ │ ├── merge_viral_hits.py │ ├── filter_virus_contigs.py │ └── summarise_abundance.py └── assets/ └── test_data/ # 小规模测试数据运行入口是一个Nextflow脚本。为什么选Nextflow而不是Snakemake是因为Nextflow自带的进程隔离、容器支持和AWS批量计算对接能力对超大队列的样本更友好。当然Snakemake同样能做到这些只是Nextflow对每个进程放独立容器这件事的支持更丝滑。直接用release包里的跑批脚本就能启动测试数据nextflow run main.nf \ -profile test \ --input samplesheet.csv \ --host_index /data/refs/human/grch38 \ --outdir results跑通测试数据之后再换成你自己的样本数据。3. 源码实现中的关键细节参数阈值与可复现性一套好的生信流程源码光把工具串起来是不够的关键是每个环节的参数设置有没有依据以及整个环境能不能被复现。下面几个点是我觉得最有价值的部分。3.1 宿主基因组索引构建与去污染阈值的实验依据宿主去除用bowtie2的时候源码默认的比对模式是--very-sensitive --end-to-end。很多流程为了速度会用--sensitive甚至默认参数但在病毒组场景里读段数量本身不多牺牲一点速度换取灵敏度是划算的。--end-to-end则能避免局部比对的错配因为病毒序列和宿主转座子或重复区域偶尔会有短片段相似如果允许局部比对很容易被误判为宿主污染。比对后的过滤阈值是双端都比对上宿主的才删除单端比对上的保留并标记。这个策略是我特意做成保守的。原因很简单病毒序列与宿主基因组之间可能存在水平基因转移片段一对读段中一条来自病毒、另一条来自宿主的情况确实存在如果按只要一端比对就过滤的标准这些真病毒会在早期就被丢掉。索引构建这一步源码提供了一个辅助脚本会用hisat2的extract_splicesites.py先推断剪接位点再用hisat2-build构建索引。用HISAT2不一定非得做转录组它的索引结构在有内含子的区域处理上比bowtie2更好能减少假比对。当然如果你的样本是粪便等低宿主含量类型这部分时间开销可能不太值得可以通过配置关闭。3.2 病毒识别阶段为什么采用多工具投票这一步是整个流程的灵魂。只用一个病毒识别工具的问题在于不同工具对什么是病毒的定义和训练集差异很大。VirSorter2基于已知病毒特征数据库和HMM模型好处是对已知类群的敏感度极高坏处是遇到全新的、和数据库差异大的病毒会直接失明。VirFinder基于k-mer频次虽然不需要数据库但对环状基因组和细菌基因组片段会有误报。源码的策略是VirSorter2判定为high confidence的序列直接进入下游VirSorter2判定为medium confidence但VirFinder置信度也超过0.9的同样保留如果两个工具都给出low/no confidence直接过滤。这个并集双过滤的策略在测试数据上比只用VirSorter2的流程多找回约18%的病毒序列同时把明显的假阳性控制住了。3.3 组装与CheckV质控的严格标准MEGAHIT跑完组装后的contig集合长度分布通常很分散。源码里设定的最小contig长度是1kbCheckV完整度评估阈值是完整度大于50%且污染度小于10%的contig保留。这个标准的背后逻辑是病毒基因组本身的长度在几十kb级别1kb以下的contig即便标注为病毒也往往只是某个基因的片段用它做下游分类和宿主预测的可信度都很差。CheckV在评估污染度时用的是可能来自两个不同基因组的比例这个指标尤其重要。宏基因组拼出来的病毒contig经常出现嵌合体也就是两条不同的病毒基因组碎片被拼接在一起CheckV能通过比较单拷贝marker基因的比例识别出这类问题。我见过有人图省事跳过CheckV直接用VirSorter2的结果注释最终产出的完整基因组里混了不少嵌合体功能注释结果自然也不可靠。3.4 环境锁定与容器化方案可复现性是源码工程里最容易被忽视、但实操时最容易让人头疼的问题。工具的版本迭代非常快你在本地跑出一套结果三个月后重新运行Virsorter2可能已经更新到2.1数据库也可能变了结果可能完全不同。为此源码做了两件事所有依赖都写在envs/*.yaml环境文件里Conda环境的版本精确到patch级别比如virsorter2.2.4不会用这样含糊的约束。Nextflow进程支持直接拉取Docker镜像或Singularity镜像每个工具都有独立的容器避免了一个环境里装所有依赖导致依赖冲突的问题。3.5 断点续跑与资源调度配置拿到一套流程就要跑全流程如果你的样本量是几十个机器配置一般般可能跑一半就被超算集群的作业时间限制杀掉。Nextflow天然支持-resume源码配置文件里也预设好了资源请求process.resourceLabels false process.container quay.io/pipeline/viruspipe:1.0 executor.queueSize 20 executor.submitRateLimit 20/min队列大小和提交频率的限制非常关键。有些集群不限制任务数但对于一个流程里动辄上百个样本你一次性提交几百个任务轻则触发管理员警告重则把集群调度器压垮。把这些参数写在配置文件里相当于给流程上了一道保险。4. 跑通流程后的实战验证与常见坑源码能跑通不等于能出好结果。我拿自己手头的肠道噬菌体数据和公共的海洋病毒组数据分别跑了一遍前后折腾了不少时间下面是几个最有代表性的问题。4.1 低深度样本的检出率为什么突然下降第一批测试样本中有一个样本的测序深度只有别人三分之一结果病毒contig检出数几乎为0。排查下来发现问题出在MEGAHIT的--min-count参数上。这个参数是组装时的一个k-mer覆盖度阈值默认值是2但在超低深度样本里大量真实的病毒reads因覆盖度不够而被当成测序错误过滤掉了。解决方法有两种方案一是把--min-count降到1但代价是组装时间暴涨并可能引入更多错误方案二是把这个样本标记为低深度在流程里换用更宽容的组装参数。我后来直接在源码的配置模块里加了一个depth_profile字段让样本可以根据自身深度自动切换参数这样既保住了低深度样本的灵敏度又不拖慢高深度样本的组装速度。4.2 宿主去污染误杀真病毒这个坑是反直觉的。有一批样本来自某种植物根际宿主参考基因组用的该植物基因组。跑完之后某个已知的ssDNA病毒基因组在所有样本里都是0覆盖度但实验PCR确认它存在。回溯到宿主去除环节发现该病毒基因组有一段400bp的区域与宿主基因组的逆转录转座子序列高度相似导致整个基因组的大量reads都被bowtie2判定为宿主污染。这个问题在流程源码里没有标准解法但我在host_filter模块里加了一个分离读段保留的选项如果一对reads中只有一端比对到宿主另一端不确定默认保留并输出到专门的ambiguous目录。这样你可以后续单独检查这些灰色地带的reads而不是直接让流程替你做了决定。对追求灵敏度的病毒组项目来说保留这些模糊读段永远是更安全的选择。4.3 多工具投票结果不一致时的裁决策略流程里VirSorter2和VirFinder会出现大量意见不合的contig。VirSorter2认为是病毒但VirFinder完全不认或者反过来。我手动抽查这些contig后发现VirSorter2认但VirFinder不认的往往是那些与数据库已知病毒高度相似、但整体AT含量极高导致k-mer分布异常的序列VirFinder认但VirSorter2不认的则往往是可能来自新病毒但长度只有2~3kb的短序列。处理策略调整为只要某一个工具的置信度达到高分阈值就直接保留不用等另一个也同意。如果两个工具都是中等置信度再额外检查这个contig是否携带已知的病毒标记基因比如末端酶大亚基、门户蛋白、衣壳蛋白。源码里加了一个viral_marker自动检索的步骤对这些悬而未决的contig做标记基因判断这样能在不显著增加假阳性的前提下把可疑contig的利用率提上来。4.4 资源调度与任务失败重试跑大样本集时还有一个常见问题某个样本的某个进程会因为平台错误或临时IO问题挂掉。默认情况下Nextflow对失败的进程会尝试重启几次但不同工具的重试策略不一样。MEGAHIT在中断重启后是没有断点续跑的直接重头开始所以如果你给它的重试次数太多极端情况下会白白浪费好几个小时。源码里给不同模块配了不同的重试上限。fastp和宿主过滤这种快任务重试5次MEGAHIT重试1次CheckV重试2次。如果重试仍然失败流程会生成一份errors.tsv文件列出失败任务的具体日志路径。这个设计比一堆中间文件散落各地要人性化得多排查问题时省了很多工夫。5. 从源码到自有流程二次开发与适配建议最后聊一聊拿到这套源码之后怎么能让它真正为你所用。直接跑测试数据、跑通示例只是第一步实际项目中几乎一定会遇到需要改动的地方。5.1 哪些模块值得替换为更新的工具我建议重点关注两个环节。第一个是病毒识别模块。Virsorter2和VirFinder的组合在目前看来依然能打但如果你研究的是RNA病毒组这个流程并不适用因为RNA病毒需要额外的逆转录步骤而且识别策略完全不同。第二个是宿主预测模块。流程里没有内置宿主预测功能但annotate.nf模块预留了接口你可以把宿主预测结果作为注释表额外一列输出。最近比较流行的工具是VirHostMatcher-Net或phpLPP前者基于k-mer距离后者基于机器学习预测精度都有明显提升。5.2 新样本类型的数据库补充病毒的数据库本质上就是你用来做比对和判断的先验知识库。如果你做的是环境病毒组比如土壤建议在Virsorter2默认数据库的基础上补充IMG/VR的病毒基因组序列如果你做的是人类肠道病毒组补充Gut Phage Database是很有必要的。源码里把数据库文件路径做成了独立的配置文件你可以直接指向自己手动下载的数据库路径不需要改任何脚本逻辑。不过有一个细节要多说一句在补充数据库后跑一次已知病毒加标验证很有必要。我一般会取一小部分已知的参考病毒基因组模拟成测序读段混入样本数据里跑流程看流程能不能把这些已知病毒找回来、找回多少比例。这一步能快速判断数据库和参数设置是否有问题比直接跑完整批样本再后悔要高效得多。5.3 流程性能基准与运行成本参考最后给你一组我自己实测的性能数据供配置集群或云服务器时参考。测试数据是12个肠道宏基因组样本双端150bp单样本约8000万条读段服务器配置为32核心CPU、128GB内存。读段质控与宿主去除耗时约18分钟/样本内存峰值约6GB病毒识别Virsorter2 VirFinder耗时约25分钟/样本内存峰值约12GB组装MEGAHIT耗时约40分钟/样本内存峰值约28GBCheckV质量评估耗时约10分钟/样本内存峰值约8GB如果你样本量很大可以考虑把前两步横向扩展到多个节点并行源码的Nextflow配置已经预设了横向扩展选项只要集群有足够的节点理论上可以做到近线性扩展。组装模块是最难横向扩展的MEGAHIT本身是多线程单机程序建议一个节点跑一个样本不要尝试用MPI方式拆分收益很低还会增加IO负担。整个流程跑完之后我的习惯是再来一遍产出物审计把CheckV判定为完整或高分环状的病毒基因组单独再比对回原始reads看看覆盖度是否均匀、有没有局部塌陷。这一步虽然不在流程源码里但在我自己的实战中它几乎每一次都能发现一两个注释阶段被掩盖的问题。做病毒组分析永远不要只信流程输出的表格。本文还有配套的精品资源点击获取