Varscan1文献解读：读段计数框架在体细胞突变检测中的原理与实战

1. 从一次肿瘤样本复核说起：为什么还要读Varscan的文献

去年帮一个做肿瘤方向的朋友复核一批配对样本的体细胞突变结果，他用的是某款主流的贝叶斯类检测器，跑出来的候选位点数量很漂亮，但真阳性率一直上不去。我建议他换一个思路，用Varscan2的mpileup2somatic模式重新跑一遍，然后把两个工具的结果取交集。结果交集里的位点，经过一轮人工核查，可信度明显高于任何单一工具的输出。这件事让我重新翻出了Varscan系列最原始的那篇文献，也就是Varscan1的发表文章，从头读了一遍。

Varscan1这篇文献，讲的是怎么从高通量测序数据里检测变异。它的核心贡献不是发明了某个复杂的统计模型，而是提出了一套非常务实的、基于"读段计数加阈值过滤"的检测框架。这套框架在今天的角度看可能显得朴素，但恰恰是这种朴素，让它在低频突变、肿瘤异质性样本、以及深度不均匀的靶向测序数据里，表现得比很多"高级"模型更稳。如果你正在做体细胞突变检测、肿瘤panel分析，或者只是想知道为什么Varscan至今还在很多流程里占一个位置，那这篇文献值得花时间读一读。

我下面会从文献的核心方法讲起，拆解它的统计逻辑，然后重点讲它在实际项目里怎么用、参数怎么调、和别的工具怎么配合，最后分享几个我踩过的坑。内容会涉及一些命令行操作和参数解释，但不需要你有很深的统计学背景，只要跑过比对和变异检测的基本流程就能看懂。

2. Varscan1文献到底提出了什么：读段计数框架的底层逻辑

2.1 文献的核心命题：用计数代替概率模型

Varscan1那篇文献的出发点很直接：当时主流的变异检测工具，比如基于贝叶斯推断的那一批，都在试图用一个概率模型去描述"这个位点是不是变异"。模型很优雅，但有个实际问题——它对测序深度、碱基质量、比对质量的分布假设很敏感。一旦实际数据的分布偏离假设，比如靶向测序里某些区域深度极高、某些区域深度极低，模型给出的置信度就会失真。

Varscan1换了个思路：我不去建模整个概率分布，我直接数。数这个位点上支持参考碱基的读段有多少，支持变异碱基的读段有多少，然后根据一个预设的阈值来判断。这个阈值包括最低覆盖深度、最低变异读段数、最低变异频率、以及一个基于Fisher精确检验或卡方检验的p值门槛。文献里把这个框架称为"read count-based"方法，翻译过来就是基于读段计数的检测。

这个思路的好处是透明。你看到一个位点被判定为变异，你可以直接去看它的深度是多少、变异读段占比是多少、p值是多少，每一个数字都能追溯到原始比对文件。相比之下，概率模型给出的后验概率，你很难直观地解释它为什么是这个值。对于临床相关的检测，这种可解释性非常重要。

2.2 统计检验的选择：Fisher精确检验在配对样本里的作用

Varscan1文献里一个关键的技术细节，是它在处理配对样本（比如肿瘤和正常配对）时，用了Fisher精确检验来比较两个样本在同一个位点上的等位基因计数。具体来说，它构建一个2x2的列联表：行是肿瘤和正常，列是变异读段和参考读段。然后检验这两个样本的变异频率是否有显著差异。

为什么用Fisher精确检验而不是卡方检验？文献里给的理由是，在低覆盖深度或者低变异频率的情况下，列联表里的期望频数会很小，卡方检验的近似条件不满足，而Fisher精确检验是基于超几何分布的精确检验，不依赖大样本近似。这一点在实际数据里很关键，因为肿瘤样本里很多位点的深度并不高，尤其是ctDNA或者FFPE样本，深度可能只有几十乘，变异读段可能只有几个，这时候卡方检验的p值会不准，Fisher检验更可靠。

我自己的经验是，在深度低于50x的位点上，Fisher检验和卡方检验给出的p值差异可以大到几个数量级。如果你用卡方检验的阈值去卡，会漏掉一批真实的低频突变。所以Varscan在配对模式下的这个选择，是有实际数据支撑的，不是随便选的。

2.3 对测序错误的处理：为什么它不依赖碱基质量重校准

Varscan1文献里还有一个容易被忽略的点：它对测序错误的处理方式。很多工具会先用碱基质量分数（BASE_QUAL）和比对质量分数（MAP_QUAL）对读段进行过滤或加权，然后再做检测。Varscan1的做法更简单：它设定一个最低碱基质量阈值和一个最低比对质量阈值，低于阈值的读段直接不计入计数，高于阈值的读段一视同仁。

这个做法看起来粗糙，但文献里解释了原因：在当时的测序平台上，碱基质量分数的校准并不完美，尤其是在读段末端和特定序列上下文里，质量分数会系统性偏高或偏低。如果过度依赖质量分数加权，反而会引入偏差。直接用一个硬阈值过滤，虽然损失了一些信息，但避免了偏差传播。

这个思路在今天依然有参考价值。我见过不少项目，花大量时间做碱基质量重校准（BQSR），结果变异检测的灵敏度并没有明显提升，反而因为重校准引入了新的不确定性。Varscan的这种"够用就行"的过滤策略，在某些场景下反而是更稳的选择。

3. 从文献到命令行：Varscan1核心参数的实操解读

3.1 最低覆盖深度与最低变异读段数：一对需要联合调整的参数

Varscan1文献里给出的默认参数是：最低覆盖深度8x，最低变异读段数2，最低变异频率0.2，p值阈值0.05。这四个参数里，前两个是最需要根据实际数据调整的。

最低覆盖深度8x这个默认值，在文献发表的年代是合理的，因为当时的测序成本高，很多样本的深度就在10x到30x之间。但现在做肿瘤panel或者全外显子，目标区域的深度普遍在100x以上，8x这个阈值就太宽松了。我一般会把--min-coverage设到20x甚至30x，具体取决于样本类型。如果是ctDNA样本，深度可能上千乘，但变异频率极低，这时候深度阈值可以设高一些，比如50x，同时把--min-var-freq降到0.01甚至0.005。

最低变异读段数2这个默认值，在深度高的时候问题不大，但在深度低的时候会出问题。比如一个位点深度只有10x，变异读段有2个，频率是0.2，刚好卡在阈值上。这种位点很可能是测序错误或者比对错误。我通常会把--min-var-count设到3或者4，尤其是在深度低于30x的位点上。但注意，这个参数不能设得太高，否则会漏掉真实的低频突变。在ctDNA场景里，变异读段数可能就只有2到3个，这时候你只能把阈值降到2，然后靠其他过滤条件来控假阳性。

这两个参数的联合调整，我的经验是：先根据样本类型确定一个深度下限，然后在这个深度下限的基础上，反推变异读段数的合理范围。比如深度下限是20x，变异频率下限是0.05，那变异读段数至少是1，但为了控假阳性，我会设到2。如果深度下限是100x，变异频率下限是0.01，那变异读段数至少是1，但我会设到3，因为深度高了，测序错误累积的绝对数量也多了。

3.2 p值阈值与变异频率阈值：谁先卡谁后卡

Varscan1文献里把p值阈值和变异频率阈值并列作为过滤条件，但实际使用的时候，这两个条件的先后顺序会影响结果。我的做法是先用变异频率阈值做粗筛，再用p值阈值做精筛。

原因很简单：变异频率是一个直观的、与深度无关的指标，它直接告诉你这个位点上变异读段的比例。如果一个位点的变异频率只有0.01，那不管p值多显著，它都很可能是噪声。先用频率卡一道，可以把大量低频率的噪声位点去掉，减少后续p值检验的计算量和多重检验校正的负担。

p值阈值在Varscan里是通过--p-value参数控制的，默认是0.05。但在配对样本模式下，Varscan会对每个位点做Fisher检验，然后输出p值。如果你做的是全基因组或者全外显子，位点数以百万计，0.05的阈值会导致大量的假阳性。这时候需要做多重检验校正，比如用Benjamini-Hochberg方法控制FDR。Varscan本身不直接做校正，但你可以把输出的p值拿出来，用R或者Python做校正，然后再卡一个校正后的阈值。

我通常的做法是：在Varscan里把--p-value设得宽松一些，比如0.1，先把候选位点都输出出来，然后在下游用FDR校正，把校正后p值小于0.05的位点留下。这样比直接在Varscan里卡0.05更灵活，也更容易控制最终的假阳性率。

3.3 参数调整的实操表格：不同场景下的推荐配置

下面这个表格是我在不同项目里总结出来的参数配置，供参考。注意这些不是绝对标准，你需要根据自己的数据特点做微调。

这个表格里的参数，我在多个项目里验证过，基本能覆盖大部分场景。但有一个例外：如果你的样本有严重的污染或者异质性很高，可能需要把min-var-freq再调低，同时把min-coverage调高，用深度来换灵敏度。

4. Varscan1在当代流程里的位置：和别的工具怎么配合

4.1 为什么我不建议单独用Varscan出最终结果

Varscan1的读段计数框架很稳，但它有一个明显的短板：它对局部比对错误的敏感性。如果一个位点附近有插入缺失或者重复序列，比对工具可能会把读段错误地比对到这个位点上，导致变异读段计数虚高。Varscan本身不做局部重比对，所以这种错误它会照单全收。

我的做法是：Varscan只作为流程里的一个环节，而不是终点。通常我会用Varscan先跑一遍，拿到候选位点，然后用一个基于单倍型的工具（比如GATK的HaplotypeCaller或者FreeBayes）再跑一遍，取两个工具的交集。交集里的位点，假阳性率会显著降低。如果两个工具的结果有差异，我会把差异位点单独拿出来，用IGV人工看一眼比对情况，再决定保留还是丢弃。

这个策略在肿瘤样本里尤其有用，因为肿瘤样本里经常有体细胞插入缺失和结构变异，这些区域的比对质量往往不高，单一工具很容易出错。两个工具的交集，相当于用两种不同的错误模型互相验证，可靠性更高。

4.2 和GATK的Mutect2对比：什么时候选Varscan

Mutect2是现在肿瘤体细胞检测的主流工具，它的模型比Varscan复杂得多，考虑了肿瘤纯度、倍性、测序错误的方向性等。在大多数情况下，Mutect2的灵敏度和特异性都优于Varscan。但有几个场景，我会优先选Varscan。

第一个场景是深度极高的靶向测序。Mutect2在深度超过500x的时候，计算量会急剧增加，而且它的模型对超深度的处理并不总是最优。Varscan的计数框架在超深度下反而更直接，因为它不依赖复杂的概率计算，只是数数。我做过一个对比，在1000x深度的panel数据上，Varscan跑完只要几分钟，Mutect2跑了几个小时，而且两者的真阳性结果高度重叠。

第二个场景是样本质量很差的情况。比如FFPE样本，DNA降解严重，测序错误模式很复杂。Mutect2的模型假设在降解样本里可能不成立，而Varscan的硬阈值过滤反而更鲁棒。当然，这并不意味着Varscan在FFPE上就没有假阳性，只是它的假阳性模式更可预测，更容易通过下游过滤去掉。

第三个场景是你需要快速筛查。Varscan的运行速度比Mutect2快一个数量级，如果你只是想做初步筛查，或者样本量很大需要快速过一遍，Varscan是更好的选择。

4.3 和Varscan2的关系：文献里的Varscan1和现在的版本有什么区别

Varscan1那篇文献发表的时候，工具还只支持单样本的变异检测。后来Varscan2增加了配对样本的体细胞检测、拷贝数变异检测、以及一些过滤功能。现在大家用的varscan somatic和varscan mpileup2somatic，都是Varscan2的功能。

但Varscan2的核心算法，还是Varscan1那套读段计数加Fisher检验的框架。Varscan2增加的功能，更多是在输入输出格式、过滤选项、以及多线程支持上。所以读Varscan1的文献，对于理解Varscan2的行为依然是有帮助的。比如Varscan2里的--min-coverage-normal和--min-coverage-tumor这两个参数，就是Varscan1里最低覆盖深度概念在配对样本里的延伸。

我建议的阅读顺序是：先读Varscan1的文献，理解计数框架和Fisher检验的逻辑，然后再去看Varscan2的文档，看它增加了哪些参数和功能。这样你在调参的时候，心里是有底的，知道每个参数背后的统计含义，而不是盲目试。

5. 那些文献里没写、但我踩过的坑

5.1 碱基质量阈值设太高会漏掉真实突变

Varscan有一个--min-avg-qual参数，默认是15。这个参数的意思是，一个碱基的质量分数如果低于15，就不计入计数。我一开始觉得这个阈值太宽松，就把它设到了30。结果发现，一批已知的阳性位点被漏掉了。

后来我去看原始比对文件，发现这些位点的变异碱基质量分数普遍在20到25之间。为什么？因为这些位点附近的序列上下文导致测序质量偏低，但变异是真实的。如果我把阈值设到30，这些位点就被过滤掉了。Varscan1文献里其实提到了这个问题，它建议的阈值是15，不是没有道理的。质量分数低的碱基里，确实有噪声，但也有真实的信号。阈值设太高，灵敏度会下降。

我现在的做法是：--min-avg-qual保持默认的15，但在下游过滤的时候，把平均质量分数作为一个参考指标。如果一个位点的变异读段平均质量分数低于20，我会标记出来，人工核查，而不是直接丢弃。

5.2 正常样本的深度不够会导致配对检验失效

Varscan的配对模式需要正常样本的深度足够，才能做Fisher检验。如果正常样本的深度太低，比如低于10x，Fisher检验的p值会不可靠，导致大量的假阳性或者假阴性。

我遇到过一个项目，正常样本用的是口腔拭子，深度只有8x左右。Varscan跑出来的体细胞突变数量异常多，而且很多位点在正常样本里也有变异读段，只是频率低。这就是因为正常样本深度不够，Fisher检验无法有效区分种系变异和体细胞变异。

解决办法有两个：一是提高正常样本的测序深度，至少到30x；二是如果深度实在提不上去，就用单样本模式跑肿瘤样本，然后和一个种系变异数据库（比如dbSNP或者gnomAD）做比对，把已知的种系变异过滤掉。第二种方法不如配对检验准确，但在深度受限的情况下是可行的替代方案。

5.3 重复标记对Varscan结果的影响被低估了

Varscan本身不做重复标记（duplicate marking），它假设你输入的BAM文件已经做过重复标记了。如果你忘了做这一步，PCR重复会被当成独立的读段计入计数，导致变异频率虚高。

我见过一个案例，一个位点的真实变异频率是0.05，但因为PCR重复，Varscan算出来的频率是0.15，直接超过了阈值，被误判为阳性。后来做了重复标记，这个位点就消失了。

所以流程里一定要在Varscan之前做重复标记，用Picard的MarkDuplicates或者samtools的markdup都可以。而且要注意，重复标记之后，最好把重复读段直接去掉，而不是只标记。Varscan默认会忽略标记为重复的读段，但如果你用的工具只标记不去除，需要确认Varscan是否正确处理了这些标记。

5.4 输出格式的坑：VCF里的p值不是校正后的

Varscan输出的VCF文件里，p值列是原始的p值，没有做多重检验校正。如果你直接拿这个p值去卡0.05，在全外显子或者全基因组数据里，假阳性会非常多。

我一开始没注意这一点，直接用Varscan输出的p值做过滤，结果候选位点数量是预期的好几倍。后来把p值拿出来，用R的p.adjust函数做BH校正，再卡0.05，位点数量就正常了。

Varscan的VCF里还有一个SOMATIC标志，表示这个位点被判定为体细胞突变。但这个标志是基于原始p值的，所以也不能完全信任。我的做法是：把Varscan的VCF导入到R或者Python里，提取p值，做校正，然后重新标记体细胞突变。这样得到的最终结果才可靠。

6. 把Varscan1文献读薄：我总结的三条核心原则

6.1 计数比建模更透明，透明比优雅更实用

Varscan1文献最核心的贡献，是它证明了在变异检测这个任务上，简单的读段计数加阈值过滤，可以做到和复杂概率模型相当甚至更好的效果。这个结论在今天依然成立。尤其是在临床相关的检测里，可解释性往往比统计优雅更重要。你不需要向别人解释一个复杂的后验概率是怎么算出来的，你只需要告诉他：这个位点深度是100x，变异读段有15个，频率是0.15，Fisher检验p值是0.001。这些数字任何人都能看懂，也都能验证。

6.2 参数没有绝对的最优，只有和场景匹配的合理

Varscan1文献里给出的默认参数，是在当时的数据条件下总结出来的。现在的测序深度、样本类型、建库方式都变了，默认参数不一定适用。你需要根据自己数据的深度分布、变异频率预期、假阳性容忍度来调整。我上面给的参数表格，也只是参考，不是标准答案。关键是要理解每个参数背后的统计含义，然后根据实际情况做取舍。

6.3 工具是流程的一部分，不是全部

Varscan1文献讲的是检测算法，但实际项目里，检测只是流程的一个环节。上游的比对质量、重复标记、碱基质量校正，下游的过滤、注释、人工核查，都会影响最终结果。我见过太多人把精力全花在调Varscan的参数上，却忽略了上游BAM文件的质量。一个比对质量差的BAM文件，用再好的检测工具也跑不出可靠的结果。所以读文献的时候，要把它放在整个流程的上下文里理解，而不是孤立地看算法本身。

最后分享一个我最近在用的技巧：如果你手头有已知的阳性对照样本，比如一个包含已知突变位点的细胞系或者参考样本，可以用它来快速验证Varscan的参数配置。把对照样本跑一遍，看已知位点是否都被检测到，同时看假阳性位点有多少。这个办法比盲目调参高效得多，通常跑两三轮就能找到适合你数据的参数组合。