全外显子组测序数据质控：从FastQ到捕获区间的关键指标-平芜编程栈

1. 把质控放在流程第4步，到底在质控什么

很多人第一次搭全外显子组分析流程的时候，前三步跑得飞快，到了第4步数据质控就开始糊弄——fastp用默认参数过一遍，samtools flagstat看一眼 mapping rate 是 99.x%，截图往报告里一贴，就直奔变异检测去了。我自己最开始也是这么干的，直到有一次一批样本 call 出来的变异位点里，同一家族的两个成员在几十个已知位点上互相矛盾，回头一层层扒，最后发现根源在第4步被跳过的那张覆盖度表上：样本的平均深度有 120x，漂亮得很，但关键基因的几个外显子实际只有 3x 覆盖。这种"平均数骗人"的情况在全外显子数据里一点都不罕见。

这篇讲的是全外显子组生信分析流程的第四个环节：数据质控。它承接上游的原始数据下机、格式整理、参考基因组比对，往下游喂给变异检测、注释、过滤。质控做得好不好，直接决定了后面所有结果的可用性，因为变异检测这一步是没有"纠错"能力的——它只会忠实地把覆盖不均、接头残留、样本污染这些东西转化成假阳性或假阴性，然后你在下游花三天时间去解释一堆根本不该存在的位点。

适合谁看？如果你是刚接手 WES 分析流程的分析人员、做遗传病方向的研究生、或者需要给实验室搭一套可复现流程的工程师，这套内容可以直接抄作业。如果你已经跑过几十个样本但从来没认真读过HsMetrics那张表，那更值得往下看——里面藏着的问题，可能已经影响了你之前发的文章。

1.1 前面三步留下的隐患，全在这一步结算

按照典型的 WES 流程顺序，第1步是实验端的样本与文库设计，第2步是下机数据整理（把 BCL 转成 fastq、按样本拆分、生成样本表），第3步是参考基因组比对（BWA-MEM 或 BWA-MEM2 比对到 GRCh38，输出 BAM）。

这三步里，每一步都会埋雷，而雷都在第4步爆出来：

文库端的雷：建库时 PCR 循环数给多了，重复率飙升；片段化没控制好，插入片段偏短，测序时读穿到接头上去；input DNA 量不够，文库复杂度低。
测序端的雷：Q30 比例偏低、index 拆分错误、同一个 lane 的样本之间发生 index hopping。
比对端的雷：参考基因组版本和捕获芯片的区间文件版本对不上，导致大量 reads 被标成 off-target；BAM 排序、去重、碱基质量重校正的顺序搞错。

这些问题在原始 fastq 上看不出来，在 flagstat 上也看不出来，只有把 FastQ 层质控和 BAM 层质控两条线都跑通，才能定位到具体是哪一环出的问题。

我一直坚持一个原则：质控不是为了出一份合格报告，而是为了在花大代价做变异检测之前，用最小的成本判断"这个样本还能不能救"。所以我在流程里给质控这一步设了两个出口：一个是 pass 的样本自动流向下一环节，另一个是 fail 的样本带上诊断信息退回重跑或者直接废弃。

1.2 WES 和 WGS 在质控指标上的分歧点

如果直接把 WGS 的质控阈值套到 WES 上，会踩一堆坑。原因在于两者的数据分布根本不同：

维度	全基因组（WGS）	全外显子（WES）
测序深度要求	30x 左右即可	目标区域通常 100x 起，临床场景 150x
覆盖均一性	主要受 GC 影响	受捕获探针效率、GC、重复序列三重影响
重复率典型值	一般低于 5%	常见 5%~20%，捕获本身会富集
off-target 概念	不存在（全基因组都是 target）	核心指标，通常 20%~40% 属于正常
插入片段	通常 350~500bp	通常 150~250bp，更容易发生读穿
关键质控指标	平均深度、覆盖均一性	on-target 率、fold-80、target 区间覆盖度分布

举个具体例子：WGS 的样本重复率超过 10% 基本要怀疑建库，但 WES 的重复率到 15% 我一般不动它——因为捕获过程本身会富集特定片段，加上测序深度高，重复率天然就高一些。拿 WGS 的阈值去卡 WES，会导致大量正常样本被误判。

1.3 两层质控的边界：哪些问题必须在 FastQ 层解决

质控分两层，分层逻辑很清晰：FastQ 层解决"数据本身脏不脏"，BAM 层解决"数据落没落到该落的地方"。

FastQ 层能解决的问题：碱基质量差、接头残留、N 碱基比例高、overrepresented 序列、双端 reads 质量不一致。这些问题如果留到 BAM 层才发现，代价是比对效率低、假阳性比对增多，而且没法回头修——因为 BAM 是比对结果，你不可能在 BAM 上把接头切掉再重新比对（理论上可以，但实践中不如在 fastq 阶段就处理干净）。

BAM 层才能解决的问题：捕获效率、目标区域覆盖度、覆盖均一性、插入片段分布、重复率、样本污染程度。这些指标必须依赖比对位置信息，FastQ 层完全无能为力。

所以流程设计上，我一般把这两个层次做成两个独立的模块，中间用样本 ID 串联，任意一层 fail 都不会污染下游就是设计目标。下面第2节讲 FastQ 层，第3、4节讲 BAM 层。

2. FastQ 层的过滤：FastQC 报表怎么看才不算白看

FastQC几乎人人都在跑，但真正会读它报表的人不多。多数人的操作是——打开 HTML，看到一堆绿色 PASS 就关掉，看到黄色 WARN 也不管，看到红色 FAIL 才紧张一下，然后照样往下跑。这套使用方法基本等于没做质控。

我的习惯是把 FastQC 的 12 个模块按"能不能救命"重新排个优先级，重点看那几个真正会连锁影响下游的，其余的快速扫过即可。

2.1 十二个模块按"能救命"的顺序重排

按我的经验，WES 数据里真正需要盯住的模块只有五个，排个序：

第一优先：Per Base Sequence Quality。这个模块看的是每个测序循环位置上的质量值分布。正常情况下，双端 150bp 测序，前 100bp 质量应该稳定在 Q30 以上，末端会自然下降。如果从头到尾质量都很低，或者中间某一段突然塌陷，那基本是测序仪出了问题，这种数据无论怎么过滤都救不回来，直接要求重测。

第二优先：Adapter Content。WES 的重灾区。插入片段短的文库，read 2 会读穿进入接头区域，接头序列就会在 fastq 里成规模出现。这个指标是后面fastp要重点处理的。

第三优先：Per Sequence GC Content。正常情况下应该是一条以基因组平均 GC 为中心的正态分布曲线。如果出现明显的双峰，或者分布严重偏移，要么是样本被其他物种污染，要么是文库里混进了非目标序列。WES 数据因为捕获富集，GC 分布本就会略微偏离正态，但出现两个明显的峰，我基本会直接判定为污染。

第四优先：Sequence Duplication Levels。注意 FastQC 算的重复是"完全相同的序列"，在测序深度高的 WES 里这个数字天然偏高，看到 30% 以上的 duplication 不要立刻惊慌，它和后面 Picard 算的重复率口径不同。

第五优先：Overrepresented Sequences。这一项经常被忽略，但它是最直接的问题指示器。如果 top 序列全是接头，说明接头污染严重；如果 top 序列是某段 rRNA 或者线粒体序列，说明捕获特异性差或者样本降解。点开看清楚到底是什么序列，比看那个百分比有用得多。

剩下的模块里，Per tile sequence quality 能帮你判断是不是某个 flowcell 区域出了问题，Per base N content 在 N 比例超过 5% 时要警惕，其余几个我基本只在排查疑难样本时才会细看。

2.2 read-through 接头是 WES 最常见的假阳性来源

接头污染这件事值得单独说，因为它在 WES 里的危害被严重低估了。

WES 的插入片段中位数通常在 180bp 左右，而测序策略是双端 150bp。这意味着当插入片段小于 150bp 时，read 1 会一直读到接头区域，read 2 也会读到接头区域。这部分 reads 如果不清洗掉，比对软件会尝试把含接头的序列比对到参考基因组上，结果就是产生大量低质量比对，甚至在某些同源区域形成假阳性堆积。

更隐蔽的问题是：接头序列本身是人工合成的，在某些样本里会导致 reads 之间高度相似，从而提高重复率、降低有效深度。你在HsMetrics里看到的重复率异常高、覆盖度不均，源头可能就在接头没洗干净。

我的处理方式是分两步：先用fastp的--detect_adapter_for_pe自动检测双端接头（它内部会做 overlap 分析，对短插入片段的检测比传统工具更敏感），然后再加一次基于插入片段长度的显式过滤。

2.3 fastp 参数逐条解释：从默认值改起

下面这套参数是我在多个项目里反复调过的版本，适合双端 150bp 的 WES 数据：

fastp \ -i sample_R1.fastq.gz -I sample_R2.fastq.gz \ -o sample_clean_R1.fastq.gz -O sample_clean_R2.fastq.gz \ --detect_adapter_for_pe \ --qualified_quality_phred 15 \ --unqualified_percent_limit 40 \ --length_required 36 \ --cut_front --cut_tail --cut_mean_quality 20 \ --correction \ --overrepresentation_analysis \ --thread 8 \ --json sample_fastp.json \ --html sample_fastp.html

逐条说为什么这么设：

--qualified_quality_phred 15而不是默认的 15 以上某个值。Q15 对应错误率约 3%，这是很多流程里被验证过的平衡点。设成 Q20 会过度削减数据量，尤其对末端质量本来就在下降的 read 2 不友好；设成 Q10 又太宽松，保留的错配会影响后续比对精度。

--unqualified_percent_limit 40的意思是：一条 read 里如果有超过 40% 的碱基低于 Q15，整条丢弃。这个阈值比默认的 40% 稍微严格一点的版本效果更好，实测丢弃率在 1%~3% 之间，属于可接受范围。如果你的数据质量本来就差，可以放宽到 50%，但别更高。

--length_required 36是过滤后保留的最短长度。低于 36bp 的 reads 在比对时会产生大量的多重比对，对 WES 来说是有害无益。有些流程会设成 50，我不推荐——尤其是在插入片段偏短的批次里，会把大量有效数据也一起扔掉。

--cut_front --cut_tail --cut_mean_quality 20是滑窗式的质量剪切，从两端切掉平均质量低于 Q20 的碱基。这一步主要针对 read 末端质量衰减，对提升比对效率有明显帮助。

--correction开启双端校正：当一对 reads 在某位置上互相矛盾时，用高质量的那个去纠正低质量的那个。这个功能在低深度数据上有争议（可能过度纠正），但在 WES 这种高深度场景下是很划算的。

--overrepresentation_analysis把 overrepresented 序列分析也做一遍，输出的 json 里会带上序列内容，方便你判断到底富集的是什么。

跑完之后，我的习惯是把过滤前后的数据量、Q30 比例、接头含量做成一张对比表，挂在报告里。这张表的用处是：当某个样本在后续 BAM 层质控出现异常时，你能快速判断问题是出在过滤前还是过滤后。如果过滤后 Q30 依然只有 70%，那说明原始数据质量就差，责任不在过滤环节。

提示：不要对fastp的过滤后数据量下降感到恐慌。WES 数据正常过滤后保留 95% 以上是常态，如果降到 90% 以下，说明原始数据本身有结构性问题（比如接头污染极重或者质量普遍偏低），这时候应该回头看原始数据，而不是调参数硬撑。

3. 比对结果里被误读最严重的三个指标

bwa mem跑完，用samtools flagstat一看，99.5% mapped，很多人这时候心就放下来了。我要泼一盆冷水：mapping rate 是全流程里信息量最低的质控指标之一。

原因很简单：人类基因组里充斥着重复序列和低复杂度区域，比对软件哪怕把一条错误来源的 read 硬塞到某个重复区域上去，它也会算作 mapped。99% 的 mapping rate 只能说明"数据能被贴上参考基因组"，完全不能说明"贴得对不对、贴的位置是不是你要的目标区域"。

3.1 mapping rate 高不等于数据干净

真正有诊断价值的比对指标有这么几个：

samtools flagstat sample.markdup.bam samtools stats sample.markdup.bam > sample.stats.txt

看flagstat输出时，我关注的是properly paired 的比例。双端数据如果 proper pair 比例低于 90%，说明要么插入片段分布异常（比如片段太短，两端比对到不同位置），要么数据里有大量嵌合 reads。这个指标比 mapping rate 敏感得多。

flagstat里还有两个数字容易被忽略：secondary 和 supplementary。如果 secondary 比例异常高（超过 5%），说明数据里有大量多重比对，多半是重复序列含量高或者 reads 太短；supplementary 高则通常出现在结构变异密集的区域，也跟插入片段偏短有关。

samtools stats的SN段落信息更丰富，几个关键字段：

字段	含义	经验参考值
`insert size average`	平均插入片段长度	与文库设计值偏差不超过 15%
`insert size standard deviation`	插入片段标准差	小于平均值的 20%
`error rate`	整体测序错误率	低于 0.5%，超过 1% 要警惕
`average quality`	平均碱基质量	高于 30
`bases mapped (cigar)`	实际比对上的碱基数	与 mapped reads 数量互相印证
`reads properly paired`	正常配对数	占比应高于 90%

这里有个实操细节：samtools stats的error rate是基于 CIGAR 里的错配数估算的，如果你的样本和参考基因组有较大的个体差异，这个值会偏高，属于正常现象。所以看这个指标要结合样本背景，不要死扣阈值。

3.2 插入片段分布是建库质量的 X 光片

插入片段分布藏着的信息量极大，我基本靠这一张图就能判断一批文库的建库质量。

用Picard CollectInsertSizeMetrics能直接输出分布图和统计表：

java -jar picard.jar CollectInsertSizeMetrics \ I=sample.markdup.bam \ O=sample.insert_size_metrics.txt \ H=sample.insert_size_histogram.pdf \ M=0.5

几种典型的异常形态，我都在实际项目里遇到过：

第一种：分布整体左移，中位数掉到 120bp 以下。说明片段化过头了。这类数据的后果是 read 2 大量读穿接头，有效比对长度下降，虽然 mapping rate 还是很高，但实际可用的覆盖度分布很差。

第二种：分布出现双峰。两个峰意味着文库里有两种不同长度的片段群体，通常是因为片段化步骤不均匀，或者混入了不同批次的文库。这种情况在多重样本混批的时候要特别小心。

第三种：分布右侧拖了一条长尾巴。一直到 1000bp 都有信号，说明文库里有基因组 DNA 残留（捕获前没有彻底片段化干净）。这类数据的 off-target 率会明显升高，因为长片段更容易被非特异性地捕获下来。

第四种：分布异常窄，标准差极小。一个尖锐的单峰，往往是 PCR 过度扩增导致文库复杂度降低的信号，这种样本的重复率一定很高。

我的做法是把这批样本的插入片段中位数和标准差做成一张散点图，如果某几个样本明显偏离群体，直接标记出来，作为后续判断的依据。这比逐个看 PDF 高效得多。

3.3 重复率：高不一定是 PCR 问题

Picard MarkDuplicates输出的重复率（PERCENT_DUPLICATION）是另一个被误读的重灾区。看到 25% 就认为建库有问题，这个判断经常是错的。

重复率的解读必须结合两个前提：测序深度和文库复杂度。

测序深度越高，重复率越高——这是数学必然。假设文库真实复杂度是 1 亿个独立片段，你测 100x 深度可能只抽到其中一小部分，重复自然就上来了。所以在比较不同样本的重复率之前，先确认它们的测序深度是否可比。

文库复杂度低导致的重复，才是真正的问题。判断方法很简单：看重复率是否随测序深度饱和。如果测序量翻了一倍，重复率依然线性增长，说明文库复杂度不足；如果重复率开始饱和，说明只是深度带来的假性重复。

WES 的重复率参考区间我一般这么划：

低于 10%：优秀，文库多样性好。
10%~20%：正常，多数 WES 样本在这个区间。
20%~30%：需要留意，检查一下插入片段分布和文库起始量记录。
高于 30%：大概率有问题，优先怀疑文库起始 DNA 量不足或者 PCR 循环数过多。

还有一个细节：去重前后的一致性。我习惯在去重前算一次估算重复率（samtools stats能提供粗略值），去重后再用 Picard 算一次精确值，两个数字差太多说明 BAM 的 flag 有问题（比如之前的流程把 duplicate flag 标错了）。这种问题在拼接不同来源的 BAM 时经常出现。

4. 捕获区间指标才是 WES 质控的真正命门

前面三节讲的所有指标，在 WGS 里同样要看。而真正把 WES 和 WGS 区分开的，是捕获区间相关的这一整套指标——这也是我认为"数据质控"这一步在 WES 流程里价值最高的地方。

因为 WES 的本质是"用一个探针集把基因组的 1%~2% 捞出来"，捞得准不准、捞得匀不匀，决定了后面所有分析的上限。捞偏了，你测 200x 也没用——目标区域只覆盖了 50x。

4.1 Picard HsMetrics 里必须盯住的字段

CollectHsMetrics是这一层的核心工具，但它输出的字段有四十多个，第一次看会懵。我按重要性把必须看的字段挑出来：

java -jar picard.jar CollectHsMetrics \ I=sample.markdup.bam \ O=sample.hs_metrics.txt \ R=GRCh38.fa \ BAIT_INTERVALS=capture_bait.interval_list \ TARGET_INTERVALS=capture_target.interval_list \ TARGET_TERRITORY=60000000 \ COVERAGE_CAP=200

注意BAIT_INTERVALS和TARGET_INTERVALS是两个不同的文件：bait 是探针实际覆盖的区间（通常比 target 大一圈），target 是你真正关心的分析区间（通常是编码外显子加上少量侧翼）。很多人只提供 target，导致 on-bait 率算不准。

关键字段的解读：

字段	含义	我的判断标准
`PCT_SELECTED_BASES`	on+near bait 碱基占比	高于 80% 良好，低于 70% 需排查
`PCT_OFF_BAIT`	完全落在 bait 之外的碱基占比	低于 25% 可接受
`FOLD_ENRICHMENT`	相对全基因组的富集倍数	一般 15~30 倍，与探针集设计相关
`MEAN_TARGET_COVERAGE`	目标区域平均深度	按实验设计定，常规 100x 起
`PCT_TARGET_BASES_20X`	目标区域达到 20x 的碱基占比	遗传病场景要求 ≥95%
`PCT_TARGET_BASES_100X`	达到 100x 的占比	高置信度变异检测看这个
`FOLD_80_BASE_PENALTY`	覆盖均一性指标（fold-80）	≤2 良好，>2.5 说明均一性差
`AT_DROPOUT`	AT 富集区域的覆盖丢失程度	<10 可接受，越小越好
`GC_DROPOUT`	GC 富集区域的覆盖丢失程度	<10 可接受，越小越好
`PCT_CHIMERAS`	嵌合 reads 比例	<5%，偏高说明文库构建有问题

FOLD_80_BASE_PENALTY这个指标特别值得讲，因为它的含义经常被误解。它的定义是：要达到目标区域 80% 的碱基都达到某个覆盖度，需要的最小覆盖深度是多少倍的均值。理想情况是 1，意味着覆盖完全均匀，平均 100x 就处处 100x。实际 WES 数据通常在 1.5~2.0 之间。

这个值的意义在于，它直接告诉你"平均深度有多不可信"。如果 fold-80 是 2.5，那你平均 100x 的样本，实际上有 20% 的目标碱基深度低于 40x，这些区域的变异检测可靠性就很成问题。所以在看MEAN_TARGET_COVERAGE的时候，必须同时看 fold-80，两个数字一起才有意义。

4.2 平均深度、均一性、低覆盖比例三者的关系

这三个指标经常被混为一谈，我用一个具体的例子来说明它们的差别。

假设两个样本：

样本 A：平均深度 120x，PCT_TARGET_BASES_20X= 96%，fold-80 = 1.8
样本 B：平均深度 130x，PCT_TARGET_BASES_20X= 88%，fold-80 = 2.6

从平均值看 B 更好，但从可用性看 A 明显更优。B 有 12% 的目标碱基低于 20x，在遗传病分析里这个比例意味着平均每个基因都有若干外显子无法可靠判断变异，实际交付价值大打折扣。

这就是为什么我在质控表里把PCT_TARGET_BASES_20X设成一票否决项，而不是只看平均深度。经验上，这个比例低于 90% 的样本，我会标记为"可分析但有缺口"，需要在报告里明确列出低覆盖区域；低于 85% 的样本，我会直接建议补测或者重做。

还有一层更细的判断：低覆盖区域是不是集中在特定类型的区域。如果缺口集中在 GC 极高（>70%）或极低（<25%）的区域，那是捕获探针的固有限制，属于行业共性问题，可以通过补充特殊探针或者改用其他方法解决，不算样本失败。如果缺口随机散布在各个正常 GC 区域，那更可能是捕获过程本身出了问题。

4.3 用 mosdepth 把覆盖度落到区间级别

HsMetrics给的是汇总数字，但真正定位问题需要区间粒度的覆盖度。这一层我用mosdepth，它比bedtools genomecov快得多，内存占用也友好：

mosdepth \ -t 8 \ -n \ -b capture_target.bed \ -Q 20 \ -T 1,10,20,30,50,100 \ --by capture_target.bed \ sample \ sample.markdup.bam

参数解释：-b指定区间文件做 on-target 统计；-Q 20只统计质量值 ≥20 的碱基（这一步很重要，否则低质量比对会虚高覆盖度）；-T指定多个深度阈值，输出每个区间达到这些阈值的碱基比例；--by让 mosdepth 按区间分块输出统计。

跑完后得到sample.regions.bed.gz（每个目标区间的平均深度）、sample.thresholds.bed.gz（每个区间达到各阈值的比例）、sample.mosdepth.summary.txt（汇总）。有了区间级数据，你可以做几件很有用的事：

第一件：找出深度低于 20x 的区间清单。直接从 thresholds 文件里筛，输出成一张表挂在报告里。这张表的价值在于，当后续变异检测在某区域没有结果时，你能立刻回查"是本来就没覆盖，还是覆盖了但没检出"。

第二件：做深度分布直方图。把regions.bed.gz的深度列抽出来画直方图，正常样本应该是一个以均值附近为中心的分布。如果出现明显的长尾（大量区间深度接近 0），说明捕获失败；如果分布异常窄，说明探针效率过于均匀（这种情况很罕见，出现了要考虑是不是数据经过了某种平滑处理）。

第三件：按基因聚合。用bedtools intersect把目标区间和基因注释做交集，然后按基因汇总覆盖度。这个粒度对临床报告特别有用，因为它直接对应"哪些基因可以出报告"。

zcat sample.regions.bed.gz | \ awk '{OFS="\t"; print $1,$2,$3,$4}' | \ bedtools intersect -a - -b genes.bed -wa -wb | \ awk '{sum[$8]+=$4; cnt[$8]++} END {for (g in sum) print g, sum[g]/cnt[g]}' \ > gene_level_coverage.txt

这段脚本比较糙，实际用的时候要处理区间重叠和基因多转录本的问题，但思路是这样：把碱基级覆盖度聚合成基因级，才有临床解释价值。

5. 质控结论怎么变成流程里的一个开关

前面讲的都是"怎么测"，这一节讲"怎么判"。质控如果不能自动化地给出 pass/fail 结论，那它就只是一个装饰品——因为人工判断不可复现，也无法规模化。

5.1 一套能落地的阈值表

下面这张表是我在多个项目里迭代出来的，不是教科书标准，是根据实际交付经验调的。使用时要注意，不同应用场景（科研 vs 临床、生殖系 vs 体细胞）的严格程度不一样，体细胞突变检测对深度的要求要高得多。

指标	层面	合格线	警戒线	一票否决
过滤后 Q30 比例	FastQ	≥85%	80%~85%	<80%
过滤后保留率	FastQ	≥95%	90%~95%	<88%
接头残留比例	FastQ	<0.5%	0.5%~2%	>2%
properly paired 比例	BAM	≥95%	90%~95%	<88%
平均插入片段	BAM	150~250bp	130~150bp	<120bp
重复率	BAM	<20%	20%~30%	>35%
目标区域平均深度	Capture	≥100x	80~100x	<60x
`PCT_TARGET_BASES_20X`	Capture	≥95%	90%~95%	<85%
fold-80	Capture	≤2.0	2.0~2.5	>2.8
`PCT_OFF_BAIT`	Capture	<20%	20%~30%	>40%
GC_DROPOUT	Capture	<5	5~10	>15
污染估计值	混合	<0.01	0.01~0.03	>0.05

这张表的关键在于分档而不是一刀切。警戒线区间的样本不直接废弃，而是标记出来，在下游分析时对该样本的结果保持审慎态度，并且在报告里注明。只有一票否决的项目才触发自动退回。

污染估计这一项我单独提一下。常用做法是用群体频率信息做估计（比如基于常见位点的等位基因频率偏离），输出一个混合比例。这个值超过 0.03 基本可以确认样本间存在交叉污染，此时无论其他指标多好，这个样本的变异检测结果都不可信——因为污染会产生大量低比例的假阳性位点，而且这些位点在过滤时很难和真正的嵌合体区分开。

5.2 不合格样本的排查顺序

样本 fail 之后不要急着下结论，按下面的顺序排查，效率最高：

第一步：先确认不是流程 bug 再怀疑样本。检查区间文件版本是否与参考基因组一致、BAM 是否真的做了去重、samtools stats里 CIGAR 是否正常。我遇到过好几次"覆盖度异常"，最后发现是 interval_list 用的是 hg19 的坐标而 BAM 是 hg38 的，这种情况下目标区间整体错位，覆盖度当然是零。

第二步：看 FastQ 层是否已经报警。如果 FastQ 层质量就差，那 BAM 层的问题基本都是连带反应，不用往下查了，直接退回重测序。

第三步：看插入片段分布和重复率。这两个指标同时异常，指向文库构建环节。翻出建库记录看起始 DNA 量、PCR 循环数、片段化条件。

第四步：看 off-target 比例和 GC_DROPOUT。这两个异常但前几步正常，指向捕获环节。可能是捕获试剂批次问题、杂交温度偏差、或者样本 DNA 质量差导致杂交效率低。

第五步：看污染估计。前面都正常但变异检测结果混乱，多半是污染。这时候要回头检查同一批次的样本，看是不是配对样本或者相邻孔位发生了交叉。

这个顺序的逻辑是：从流程到数据、从上游到下游、从整体到局部，每一步排除一类可能性，避免在错误的层面积累排查成本。

5.3 报告聚合与留痕

质控做好了，最后一步是把结果聚合成一份可归档的报告。这里我用MultiQC：

multiqc . \ --config multiqc_config.yaml \ -n wes_qc_report \ -f

MultiQC会自动识别fastp、samtools stats、Picard、mosdepth的输出格式并聚合。但默认配置不够用，我一般会写一个配置文件，做三件事：

第一，自定义阈值线的颜色，让 pass/warn/fail 三档在图上直接可见，不用去对表。

第二，把关键指标提取成一张总表，放在报告第一页，包含样本 ID、数据量、Q30、重复率、平均深度、PCT_TARGET_BASES_20X、fold-80、最终判定。这张表是整个质控输出的核心，其他图表都是它的支撑材料。

第三，固定报告版本号和数据指纹。每次质控输出都记录使用的参考基因组版本、区间文件版本、各工具版本号。这一步看起来繁琐，但在追溯历史样本时价值巨大——半年后你发现某批样本结果可疑，能不能复现当时的质控条件，全看这份记录。

留痕这件事我踩过坑。早期项目没记录工具版本，后来fastp升级后默认参数变了，同一批数据重跑出来的过滤结果和当初不一样，导致无法复现。从那之后我在流程里强制写版本信息到报告的 metadata 段落。

6. 返工过几次之后总结的几个坑

下面这几个问题都是我在实际项目里真实踩过的，每个都值得单独拎出来讲。

6.1 bait 与 target 区间文件版本错配的连锁反应

这个问题最隐蔽，因为它不会报错——流程照样跑完，指标照样输出，只是数字难看。

具体表现是：PCT_SELECTED_BASES只有 40% 左右，PCT_OFF_BAIT高到 50% 以上，但 mapping rate 是 99.8%，FastQ 层质控全绿。第一次遇到的时候我怀疑是捕获失败，差点把整批样本判废。

排查方法很朴素：把 bait 区间文件的前几行坐标打印出来，和 BAM 里某个已知基因的坐标做对比。

head -3 capture_bait.interval_list samtools view sample.bam chr1:1000000-1001000 | head -2

如果区间文件里的坐标和 BAM 头部的参考序列版本不一致（比如一个用chr1，另一个用1；或者一个是 hg19 的坐标量级，另一个是 hg38 的），立刻就能看出来。

这个问题的高发场景是：区间文件从供应商那里拿到的版本，参考基因组是用公共资源下载的版本，两者中间隔了一次版本升级。我的做法是在流程入口强制校验染色体命名规则：读 BAM header 的@SQ行，检查染色体名格式，再和区间文件的第一列做集合比对，不一致就中断流程。

samtools view -H sample.bam | awk '/^@SQ/ {print $2}' | sed 's/SN://' | sort > bam_chroms.txt cut -f1 capture_bait.interval_list | grep -v '^@' | sort -u > interval_chroms.txt diff bam_chroms.txt interval_chroms.txt

这十几行脚本能省掉半天排查时间，非常值得加进流程。

6.2 平均深度及格但外显子边界漏覆盖

第二个坑更刁钻：样本的MEAN_TARGET_COVERAGE有 110x，PCT_TARGET_BASES_20X也有 94%，看起来很健康，但下游变异检测在某些基因的外显子边界处反复出现低质量位点。

原因在于外显子边界区域的覆盖度天然偏低。捕获探针的设计是围绕外显子展开的，但探针在边界处的效率会下降，加上可变剪接的存在，边界附近本来就有覆盖度梯度。如果只看汇总指标，这些局部问题会被均值掩盖。

解决办法是把覆盖度分析下沉到外显子级别的首尾各 20bp：

# 提取每个外显子的首尾20bp作为重点关注区域 bedtools flank -i exons.bed -g genome.txt -l 20 -r 0 > exon_start_flank.bed bedtools flank -i exons.bed -g genome.txt -l 0 -r 20 > exon_end_flank.bed cat exon_start_flank.bed exon_end_flank.bed | sort -k1,1 -k2,2n > exon_boundary.bed mosdepth -b exon_boundary.bed -Q 20 exon_bd sample.markdup.bam zcat exon_bd.regions.bed.gz | awk '$4 < 20 {print}' > low_cov_boundary.txt

得到的low_cov_boundary.txt就是需要警惕的区域清单。在实际报告中，我会把这些区域和最终检出的变异位点做交集，如果某个变异正好落在低覆盖边界上，它的可信度要打折扣，需要 Sanger 验证。

这个做法看起来麻烦，但它解决了一个根本问题：WES 的"合格"是一个统计意义上的概念，而临床决策需要的是区域级别的确定性。二者之间必须靠这种细粒度分析来桥接。

6.3 样本交叉污染的早期信号

污染这件事，指标层面能看到的信号其实很晚。等到污染估计值超过 0.05，样本基本已经废了。但有一些早期信号可以提前发现：

信号一：杂合位点比例异常。在常染色体上，正常样本的杂合/纯合比例大致在一个范围内。如果某样本的杂合位点比例显著偏高，可能是污染引入了额外的等位基因。这个要在变异检测之后才能算，属于事后指标，但它能反过来提示你回查质控数据。

信号二：PCT_CHIMERAS异常。嵌合 reads 比例升高有时是污染的间接表现，因为来自不同样本的片段在文库构建时可能形成非自然配对。

信号三：重复率在样本间出现规律性差异。如果同一批次的样本重复率呈现某种梯度（比如按孔位顺序递增），要警惕建库过程中的交叉。

信号四：FastQ 层的 overrepresented 序列在不同样本间高度相似。如果两个本应无关的样本，其 top overrepresented 序列完全一致，那很可能是污染。这个信号在 FastQ 层就能看到，是最早的预警。

我在流程里加了一步自动检查：把同一批次所有样本的 top 20 overrepresented 序列取出来做两两比对，如果两个样本的序列集合重合度超过 80%，就标记为"疑似污染关联"，人工复核。

这个检查的成本极低，但抓出过两次真实的样本混淆问题——其中一次是两个样本在文库构建时标签贴反了，数据层面完全看不出来，只有通过这条线索才发现。

还有一个我个人的习惯：每批样本跑完之后，把质控指标按样本在板上的位置画成热图。如果某个区域出现指标聚集性异常（比如某一列的样本普遍重复率高），那基本可以确定是操作环节的问题，而不是样本本身的问题。这种空间维度的分析，比逐个看数字有效得多。

最后分享一个观念上的转变。我早期做质控，追求的是"所有指标都过线"，后来发现这个目标本身有问题——因为指标之间是有取舍的，过度清洗会损失数据，过度保留会引入噪音。现在我的判断标准变成了：这批数据能不能支撑我要做的那个具体分析。做生殖系变异检测，我需要的是目标区域 20x 以上的覆盖比例；做体细胞突变检测，我需要的是更严格的均一性和更高的深度；做拷贝数分析，我需要的是 reads 分布的均匀性而不是绝对深度。目标不同，质控的侧重点就不同，阈值表也应该跟着调整。这套思路没法照抄，但值得每个做流程的人自己走一遍。