全外显子组测序的数据质控是个特别容易被走过场的环节。跑个 FastQC看到一排绿勾截图贴进报告收工。我刚开始接手流程搭建的时候也是这么干的直到有一批样本在变异检出阶段冒出一堆说不清的假阳性回头翻质控记录才发现那批样本的 on-target rate 只有 40% 出头fold-80 接近 3.0等于把一大半测序数据白扔在了目标区域外面。重新捕获来不及整个项目延期两周。从那以后我就把全外显子组的数据质控从一个看一眼的环节改成了流程里真正会拦人的闸门。流程走到第四步前面实验设计、建库测序、原始数据交接这些都干完了接下来要决定的是这批数据到底能不能往下走、哪些样本要重做、哪些指标可以睁一只眼闭一只眼。说白了质控这一步不产生任何科学结论但它决定了后面所有结论的可信度。这篇讲的是我在实操里怎么设计这一环三层指标分别看什么、每个数字背后的含义、参数为什么这么设、什么情况下必须叫停、以及那些看着吓人其实没事的数字怎么区分。不管你是刚接触生信流程、只会照着脚本敲命令的新手还是手里已经有一套管线、想把质控做成自动放行的老手应该都能从里面找到能直接抄走的东西。1. 质控在整个WES流程里的位置与设计思路1.1 为什么质控不是跑个FastQC就完事很多人对质控的理解停留在原始数据干净不干净这一层实际上 WES 的质控要分三个完全不同的层次来看每一层对应的问题类型、工具、判断标准都不一样。第一层是read 层面也就是刚下机的 fastq 文件。这一层关心的是碱基质量、接头残留、N 碱基比例、GC 分布、是否有其他物种污染。工具是 FastQC、fastp、MultiQC。这一层出问题通常是测序或者建库环节的锅。第二层是比对层面也就是比对到参考基因组之后的 BAM 文件。这一层关心的是比对率、重复率、插入片段长度分布、捕获效率、目标区域覆盖度和均一性。工具是 BWA-MEM、samtools、Picard/GATK、mosdepth、Qualimap。这一层出问题可能是捕获探针批次、文库复杂度、参考基因组版本不匹配。第三层是样本层面也就是把所有样本放一起横向看。这一层关心的是样本间污染、性别是否和记录一致、有没有串样、批次效应。工具是 verifybamid2、SNP 指纹比对、MultiQC 汇总。这一层出问题最要命因为它不会在单样本报告里显形只有横向比对才看得出来。我见过不少流程只做了第一层FastQC 全绿就放行结果第二层、第三层的问题一直到变异过滤阶段才暴露。所以我在自己的流程里把三层都做成强制卡点任何一层亮红灯就中断不允许先跑下去看看。1.2 三层质控的分工与放行标准为了让流程能自动判断我给每个指标都定了三种状态绿自动放行、黄记录并继续但要在报告里标注、红中断等人工确认。这个分级很关键因为如果把所有指标都设成硬阈值稍微偏一点就中断流程根本跑不下去反过来全部放行质控就失去意义。下面这张表是我目前用的核心指标和阈值人类样本、双端 150bp、常规 WES 捕获芯片的场景其他物种或者特殊建库方式要相应调整层级指标绿放行黄记录红中断readQ30 比例 90%85%–90% 85%read接头残留占比 0.5%0.5%–2% 2%readN 碱基比例 0.5%0.5%–1% 1%read修剪后保留 read 比例 95%90%–95% 90%比对比对率 98%95%–98% 95%比对重复率 20%20%–35% 35%比对on-target rate 70%55%–70% 55%比对平均目标深度 100x80–100x 80x比对≥20x 覆盖比例 95%90%–95% 90%比对fold-80 penalty 1.51.5–2.0 2.0样本污染率 FREEMIX 1%1%–3% 3%样本性别一致性一致—不一致这张表不是拍脑袋定的。Q30 90% 是 Illumina 平台近几年的常规水平低于 85% 基本可以认为测序环节出了问题。比对率 98% 对人类样本来说不算苛刻因为参考基因组质量好正常样本轻松能到 99% 以上。on-target rate 和捕获芯片的设计、探针批次强相关不同批次的芯片表现差异不小所以我把黄色区间设得比较宽。注意阈值一定要跟自己的平台、芯片、建库方式绑定。别人流程里的 95% 搬到你这儿可能是灾难因为建库方式、测序仪型号、甚至 flowcell 批次都会影响基线。我的建议是先用二十个以上已知合格的样本跑一遍统计出自己实验室的分布再定阈值。1.3 工具链选型的取舍逻辑质控工具的选择上我踩过的坑比想象中多。这里把几个关键选择和背后的理由说清楚省得你重复试错。FastQC 这个没得选它就是行业事实标准报告模块固定、解读资料多、MultiQC 原生支持。它的问题是不做任何处理只出报告所以修剪还得靠别的工具。修剪工具我在 fastp、Trimmomatic、cutadapt 之间反复横跳过。最后长期用 fastp原因有三个一是速度快C 写的多线程同样数据量比 Trimmomatic 快三到五倍这对大批量样本很关键二是双端接头能自动识别不需要你手动去填 adapter 序列省掉了接头序列填错导致修剪不干净这类低级错误三是单次运行同时输出 json 和 htmlMultiQC 直接能读报告链路顺。Trimmomatic 的优势是参数是 step 式的、每一步做什么非常明确可控性强老流程里常见如果你想精细控制每一步的剪切逻辑它仍然值得用。cutadapt 更适合处理非 Illumina 平台或者特殊接头。比对工具用 BWA-MEM这个和 GATK 的官方推荐一致。核心原因是 WES 里 indel 附近会有大量 soft-clip 的 readBWA-MEM 对局部比对和 soft-clip 的处理更成熟后续 indel 检出更稳。bowtie2 在短 read 场景更快但对长 indel 的支持一般我一般不用于 WES。覆盖度统计这块早期用 samtools depth后来全换成 mosdepth。原因是 samtools depth 单线程一个 30x 的 WES BAM 跑下来要几十分钟mosdepth 多线程加算法优化同样的活儿一两分钟出结果而且原生支持 BED 区间、窗口统计、阈值量化输出格式也规整。这个替换对整条流程的耗时改善非常明显。汇总报告用 MultiQC把所有样本的 FastQC、fastp、Picard、mosdepth 结果一次性拉到一个页面上横向比对特别直观。批次效应、离群样本往往就是在这个页面上第一眼发现的。2. 原始下机数据的质控从FastQC报告里读出真问题2.1 FastQC七个关键模块到底怎么看FastQC 一页报告十几个模块新手容易逐条看绿勾其实真正需要盯的就那么几个。我按重要性排一下。Per base sequence quality是最核心的一张图。横轴是 read 上的碱基位置纵轴是 Phred 质量值。Phred 分数 Q 和错误率的换算关系是Q -10 × log10(P)也就是 Q20 对应 1% 的错误率Q30 对应 0.1%。正常的双端 150bp 数据图形应该是前几个碱基稍微低一点这是随机引物起始带来的偏好性属于正常中间平稳末端缓慢下滑。如果末端出现断崖式下跌说明测序后期的信号衰减明显需要在修剪时加尾部滑窗剪切。如果整条 read 都低那就是测序本身的问题修剪救不回来。Per sequence GC content值得单独说。理论上这条曲线应该是一个接近正态的单峰峰的位置和参考基因组的 GC 含量吻合。如果出现双峰或者整个分布明显右移八成是污染——可能是细菌、真菌也可能是另一个物种的样本串进来了。这时候不要急着往下跑先做物种溯源。Sequence Duplication Levels里蓝线是理论分布、红线是实际分布。WES 因为要经过 PCR 扩增重复率天然比 PCR-free 的 WGS 高15%–25% 很常见不用紧张。真正要警惕的是红线远远跑到蓝线右边重复率超过 40%那说明文库复杂度太低可能是起始 DNA 量不够、PCR 循环数过多这种样本即使测了有效信息量也有限。Overrepresented sequences列出占比超过 0.1% 的序列。这里出现接头序列、引物序列、rRNA 片段都很正常但如果出现一段你完全陌生的序列且占比很高就要警惕外源污染。Per base N content一般应该贴着底部超过 1% 说明测序仪在某个位置读不出来通常是硬件或者试剂问题。Adapter Content是累积曲线理想状态是全贴底。这条曲线一旦抬头就必须修剪否则接头序列会被当成 read 比对到参考基因组上制造假阳性变异。Per tile sequence quality是按 flowcell 的 tile 分块看质量。如果某一块区域整片泛红说明那个位置的测序出了问题通常是气泡、堵塞或者试剂分布不均。这类问题属于硬件层面你没法处理但必须在报告里记录下来作为后续判断依据。2.2 接头污染与低质量碱基的处理策略接头污染的来源主要有两个一是插入片段长度小于测序读长导致 read 末端读进了接头序列二是接头二聚体。前者在插入片段分布偏短的文库比如 FFPE 样本、降解样本里特别常见。处理思路要分情况。如果只是零星的接头残留fastp 的自动检测加默认修剪就够了。如果 Adapter Content 曲线抬头很明显或者插入片段整体偏短就需要打开滑窗质量剪切把末端低质量区一并切掉。这里有个很实际的取舍修剪得越狠数据越干净但有效 read 也越少。我在实际操作中会看修剪后保留 read 比例这个指标如果低于 90%就要回头想想是不是参数设得过严而不是无脑加码。有些同行为了追求 FastQC 全绿把参数调得极狠结果丢了 20% 的数据得不偿失。还有一个容易忽略的点双端 read 的重叠区纠错。fastp 默认会对双端 read 的重叠区域做一致性校正把低质量的碱基用高质量的互补碱基覆盖掉。这个功能对降低整体错误率有帮助但对 WES 来说重叠区的校正会轻微改变碱基质量值的分布如果你后续要做严格的 BQSR需要留意一下。我的做法是默认开启因为这个校正带来的收益通常大于那点分布扰动。2.3 fastp的完整参数配置与逐项解释下面这套是我目前稳定在用的配置人类双端 150bp WES 场景可以直接抄fastp \ -i raw_R1.fastq.gz -I raw_R2.fastq.gz \ -o clean_R1.fastq.gz -O clean_R2.fastq.gz \ --detect_adapter_for_pe \ --cut_front --cut_tail \ --cut_window_size 4 --cut_mean_quality 20 \ -q 20 -u 30 -n 5 -l 50 \ --correction \ --thread 16 \ --json qc/sample.json --html qc/sample.html \ --report_title sample_fastp_report逐个拆一下为什么这么设。--detect_adapter_for_pe是双端接头自动检测。虽然 fastp 内置了 Illumina 常见接头序列但不同建库试剂盒的接头序列可能有细微差别自动检测比死填序列稳。这个参数一定要开尤其是你不确定建库试剂盒型号的时候。--cut_front --cut_tail --cut_window_size 4 --cut_mean_quality 20这一组是滑窗质量剪切。逻辑是从 read 两端各取一个 4bp 的窗口如果窗口内平均质量低于 Q20就把这个窗口切掉然后往前推进。开启cut_front是为了处理 read 开头那几个碱基的随机引物偏好性带来的低质量cut_tail处理末端衰减。窗口设 4 而不是更大是因为窗口越大越容易误伤4bp 的颗粒度对 150bp 的 read 比较合适。-q 20定义合格碱基的门槛-u 30表示一条 read 里不合格碱基超过 30% 就整条丢弃。这两个配合使用。为什么不用更严的 Q30因为 Q30 作为丢弃阈值太狠了会丢掉大量只是末端稍微差一点的 read而 Q20 对应的 1% 错误率已经能满足 WES 变异检出的需求。-n 5表示 N 碱基超过 5 个就丢弃。-l 50表示修剪后长度短于 50bp 的 read 直接丢弃。这个 50 是经验值正常 WES 插入片段在 150–250bp修剪后还能留 100bp 以上所以 50 基本不会误伤。但如果你做的是 FFPE 或者降解严重的样本插入片段本身就短这时候要降到 30否则会丢掉大量有效 read。我吃过这个亏一批 FFPE 样本用了-l 50直接丢了快两成数据。--correction开启重叠区纠错。前面说了默认开。--thread 16线程数按机器核数调一般设到物理核数就行设太高反而因为内存带宽成为瓶颈。还有两个参数我特意没用--dedup和--trim_poly_g。--dedup是 fastp 自带的序列级去重它会直接丢掉重复的 read 对。这个做法我强烈不推荐原因在下一章展开。--trim_poly_g是修剪 polyG 尾巴这是 NextSeq/NovaSeq 双色荧光平台的特有现象如果你用的是这类平台需要加上如果是四色平台的旧机型就不用。提示修剪完一定要重新跑 FastQC和修剪前的报告并排放确认接头曲线下去了、末端质量抬上来了。我见过有人改完参数没复检等于白改。3. 比对层面的质控比对率、重复率、捕获效率3.1 BWA-MEM比对与read group的坑比对这一步重点不是跑通而是跑对。BWA-MEM 本身的参数不多但 read group 这个事必须说透。bwa mem -t 16 -K 100000000 \ -R RG\tID:sample01\tSM:sample01\tPL:ILLUMINA\tLB:lib01\tPU:unit01 \ ref.fa clean_R1.fastq.gz clean_R2.fastq.gz \ | samtools sort - 8 -m 2G -o sample01.sorted.bam - samtools index sample01.sorted.bam-R里的SM字段是样本名ID是 read group 标识LB是文库名PL是测序平台。这几个字段里SM最关键因为 GATK 后续会按SM把 read 分组判断哪些 read 属于同一个样本。如果SM写错了会出现两种典型灾难一是同一个样本的多个文库被当成不同样本导致变异检出时把它当混合样本处理二是不同样本的SM被写成一样导致所有样本被合并最后出的 VCF 只有一条样本记录。这个坑我在早期项目里踩过一次当时是批量脚本里的变量替换出了错几百个样本的SM全成了一个值。表面上看流程跑得顺顺利利直到变异检出的 VCF 打开只有一列样本才发现不对全部重跑。-K参数控制每次批处理加载到内存的碱基数默认是 1000 万。这个值会影响内存占用和比对的批次划分设大了内存吃得凶设小了批次多、速度稍慢。我一般设成 1 亿配合足够的线程数速度比较理想。-t 16是线程数。这里有个实操建议BWA-MEM 的线程效率不是线性的超过 16 线程后提升不明显而 samtools sort 是另一个瓶颈。我的做法是 BWA 用 16 线程管道后接 samtools sort 用 8 线程、2G 内存缓冲整体吞吐比较平衡。3.2 重复率的正确解读与去重时机重复率这个指标特别容易被误读。先说清楚它是怎么来的。WES 的重复 reads 有两个来源一是 PCR 扩增偏倚同一个原始 DNA 分子被扩增了多份测出来是同一条序列二是测序深度足够高时随机取样的碰撞。对于 100x 以上的深度第二种来源占比不低这部分重复是数学上不可避免的不代表文库质量差。所以看到一个 25% 的重复率不要慌。判断标准要看同批次比较如果同批次样本普遍在 15%–20%某个样本突然 40%那说明这个样本的文库复杂度确实低值得回头看建库记录。如果全批次都是 35%那更可能是芯片/测序深度/建库方式带来的系统性特征不是单样本问题。去重的工具我用 GATK 的 MarkDuplicates底层就是 Picardgatk MarkDuplicates \ -I sample01.sorted.bam \ -O sample01.md.bam \ -M qc/sample01.md_metrics.txt \ --CREATE_INDEX true \ --VALIDATION_STRINGENCY LENIENT注意这里是Mark标记而不是Remove移除。标记的意思是给重复 read 打上 flag但不从 BAM 里删掉。为什么不删因为后续的变异检出环节需要这些重复 read 来提供支持证据尤其是低等位频率的位点。GATK 的 HaplotypeCaller 默认会忽略标记为重复的 read但保留它们在文件里便于回溯和可视化检查。直接删掉就失去了这层信息。去重的时机也很关键。标准顺序是比对 → 排序 → 标记重复 → BQSR碱基质量值重校准→ 变异检出。去重必须在 BQSR 之前因为重复 read 会让同一个位置的碱基被重复计数干扰质量值的经验估计。这个顺序如果搞反了BQSR 的校正结果会有系统性偏差虽然不一定报错但会影响后续的变异质量。还有一个常见误区试图在 fastq 层面用 fastp 的--dedup去重。这个做法的问题是fastp 只能基于序列完全一致来判断重复而真正的重复判断要基于比对坐标——两条序列不完全一样但比对到同一个起始位置的 read才是真正的 PCR 重复。所以 fastq 层去重会漏掉大量真重复同时误删一些序列相同但实际来自不同基因组位置的 read比如重复区域。我的做法是 fastq 层完全不去重全部交给 BAM 层的 MarkDuplicates。3.3 捕获效率与覆盖度的计算实操WES 和 WGS 最大的区别就在这一步你花钱测的是全基因组但真正有用的只是那 1%–2% 的编码区。所以捕获效率on-target rate直接决定了这笔钱花得值不值。计算 on-target rate 最直接的方式是用 GATK 的 CalculateHsMetricsgatk CalculateHsMetrics \ -I sample01.md.bam \ -O qc/sample01.hs_metrics.txt \ -R ref.fa \ --BAIT_INTERVALS bait.interval_list \ --TARGET_INTERVALS target.interval_list \ --PER_TARGET_COVERAGE qc/sample01.per_target_coverage.txt这里要区分两个概念bait是探针覆盖的区域target是你真正关心的区域。探针设计时 bait 通常比 target 略大因为要留出捕获的余量。计算 on-target rate 用的是 bait计算覆盖度用的是 target。把这两个搞混是新手常见错误会导致指标口径不一致不同批次没法横向比。输出文件里有几个关键字段值得重点看字段含义参考范围PCT_SELECTED_BASESon-target rate 0.70MEAN_TARGET_COVERAGE目标区平均深度 100PCT_TARGET_BASES_20X≥20x 覆盖的靶碱基比例 0.95FOLD_80_BASE_PENALTY覆盖均一性 1.5PCT_PF_READS_ALIGNED比对率 0.98AT_DROPOUT/GC_DROPOUTAT/GC 偏好导致的覆盖缺失越低越好ZERO_CVG_TARGETS_PCT完全无覆盖的靶区比例 0.02FOLD_80_BASE_PENALTY这个指标值得单独解释因为它最容易被忽略但信息量最大。它的计算方式是把所有靶碱基按覆盖深度排序取第 20 百分位的深度值用平均深度除以这个值。理想情况是 1.0意思是 80% 的靶区域深度都达到了平均值。如果算出来是 2.5意味着 80% 的区域深度只达到平均值的 40%覆盖严重不均——少数区域深度极高大量区域深度很低。这种样本即使平均深度看着漂亮比如 150x实际有效覆盖可能很差。AT_DROPOUT和GC_DROPOUT反映的是捕获探针的序列偏好。GC 含量极高的区域比如启动子区、第一外显子捕获效率天然低这是技术限制。如果这两个值异常高说明该样本在高 GC 或高 AT 区域的覆盖缺失明显后续做 CNV 分析时要特别小心。除了 CalculateHsMetrics我还会用 mosdepth 做一次交叉验证因为它速度更快而且能输出窗口级的深度分布便于画图mosdepth \ -t 8 \ -b target.bed \ -n \ --by 200 \ --quantize 0:1:10:20:30:50:100: \ qc/sample01 \ sample01.md.bam--by 200表示按 200bp 窗口统计--quantize定义了深度分档输出文件里会直接给出每一档覆盖的碱基比例。-n表示不输出 per-base 文件那个文件几十 GB一般不需要。跑完会得到sample01.mosdepth.summary.txt和几个 regions 文件summary 里有总平均深度、各染色体的平均深度regions 文件里是每个窗口的深度。如果你想自己算一份定制化的报告用 mosdepth 的窗口输出配合一段 Python 就够了import gzip def coverage_summary(regions_gz, thresholds(1, 10, 20, 30, 50)): depths [] with gzip.open(regions_gz, rt) as fh: for line in fh: parts line.rstrip(\n).split(\t) depths.append(float(parts[-1])) n len(depths) mean sum(depths) / n res {windows: n, mean_depth: round(mean, 2)} for t in thresholds: hit sum(1 for d in depths if d t) res[fpct_ge_{t}x] round(hit / n * 100, 2) ordered sorted(depths) p20 ordered[int(n * 0.2)] res[fold80] round(mean / p20, 2) if p20 0 else None return res print(coverage_summary(qc/sample01.regions.bed.gz))这段脚本按窗口算不是按碱基算所以算出来的百分比和 CalculateHsMetrics 会有一点差异但趋势一致用来做快速筛查足够了。如果要严格按碱基算得读 per-base 文件那个开销大得多我一般只在对某个样本有疑问的时候才做。注意用samtools depth算平均深度有一个经典陷阱。不加-a参数时它只输出有 read 覆盖的位置没覆盖的位置压根不出现在结果里。这时候你拿总深度除以行数算出来的平均值会系统性偏高因为零覆盖的位置被排除了。必须加-a把零覆盖位置也输出或者干脆换 mosdepth。我第一次踩这个坑的时候一个实际只有 60x 的样本算出来是 95x白高兴一场。4. 样本层面的质控污染、性别与一致性4.1 污染率估算看不见的杀手样本间污染是质控里最隐蔽的问题因为它不会让任何单个指标明显异常但会系统性地污染下游所有变异检出结果。轻度污染的样本低等位频率的假阳性会明显增多重度污染的样本基因型判断直接失真纯合会被判成杂合。估算污染我用 verifybamid2它基于一组公共的常见 SNP 位点通过比对 read 上的等位基因频率分布来反推污染比例verifybamid2 \ --Reference ref.fa \ --BamFile sample01.md.bam \ --SVDPrefix 1000g.phase3.100k.b38.vcf.gz.dat \ --Output qc/sample01.selfSM \ --Verbose \ --NumThread 8输出文件里的FREEMIX就是污染率估计值。我的判断标准是低于 1% 视为干净1%–3% 要记录并在报告里标注超过 3% 就要考虑样本降级或者重抽。超过 5% 基本可以判定污染这个样本的变异检出结果不能直接用于临床或科研结论。CHIPMIX是另一个字段用来估计是否存在同一个体的另一种组织混入比如血液样本里混了口腔上皮。这个值异常的情况比较少见但如果出现值得回头查样本采集记录。这里要说清楚一点verifybamid2 需要一份参考的 SNP 位点集合和对应的 SVD 模型文件这些属于公开资源具体来源按你所在环境的规范获取不要随便网上下一个来源不明的文件。模型文件和参考基因组的版本必须严格对应版本不一致算出来的污染率会是错的。4.2 性别核查与样本指纹比对性别核查是个成本极低但价值极高的检查。原理很简单男性只有一条 X 染色体Y 染色体有覆盖女性有两条 X 染色体X 上的杂合位点比例接近正常常染色体Y 染色体基本无覆盖。实操上有两种做法。一种是看覆盖度统计 X 和 Y 染色体的平均深度除以常染色体的平均深度做归一化。正常男性 X 的相对深度约 0.5Y 约 0.5正常女性 X 约 1.0Y 接近 0。这个方法快但对覆盖度不均一的样本会有波动。另一种是看杂合度统计 X 染色体非拟常染色体区域上的杂合 SNP 比例。女性应该在 0.3%–0.5% 这个量级男性接近 0。这个方法更准但需要先跑一遍变异检出。我一般两个都做交叉验证。如果两个结果都指向同一个性别但和登记记录不符那就是样本错配必须整批回溯。样本指纹比对是另一个必做项尤其在多样本同批上机的情况下。做法是挑一组高多态性的 SNP 位点从每个样本的 BAM 里提取基因型两两计算一致性。同批次内除了已知的重复样本其他样本之间的一致性应该在 50%–60% 左右随机两个无关个体的一致性就在这个量级。如果两个样本的一致性超过 90%那基本可以确定是同一个样本被测了两次或者发生了串样。这一步的价值在于它能发现前两层质控完全看不出来的问题。两个样本各自的数据都很漂亮FastQC 全绿、比对率 99%、覆盖度 120x但它们其实是同一个人的 DNA——这种问题只有横向比对才能发现。4.3 质控汇总把散落的指标收成一张表单样本的报告再详细也没法替代横向汇总。我的做法是每个样本跑完质控后把关键指标写进一个统一的 TSV最后用 MultiQC 或者自建的脚本做可视化。字段大致是这样import pandas as pd fields [ sample_id, raw_reads, clean_reads, clean_ratio, q30_before, q30_after, adapter_pct, mapping_rate, dup_rate, insert_size_median, on_target_rate, mean_depth, pct_20x, pct_30x, fold80, gc_dropout, at_dropout, freemix, sex_check, fingerprint_cluster ] df pd.DataFrame(columnsfields) df.to_csv(qc/batch_summary.tsv, sep\t, indexFalse)这张表的作用是让离群样本无处遁形。我会按mean_depth和fold80画散点图正常样本会聚成一团离群的点一眼就能看出来。也会按on_target_rate排序看看有没有明显低于批次中位数的。这里有个经验批次效应比单样本异常更常见也更容易被忽略。如果某一批样本的 on-target rate 集体比上一批低 10 个百分点那大概率是捕获探针批次或者实验操作的问题而不是每个样本各自出了状况。这种情况下单看每个样本的报告都是黄灯但不严重横向一看才知道是系统性问题。所以汇总这一步不是可选项是必做项。提示把历史批次的汇总表也留着做成累积视图。新批次进来的时候和过去半年的数据比一比能提前发现试剂批次、操作人员、设备状态的变化。这个习惯帮我躲过好几次批量返工。5. 常见问题排查与踩坑实录5.1 问题速查表下面这张表是我这些年攒下来的基本覆盖了 WES 质控里八成以上的异常情况现象可能原因排查动作处理方式Adapter Content 曲线抬头插入片段短于读长、接头二聚体看插入片段分布、Overrepresented sequences开启自动接头检测 尾部剪切必要时缩短-lPer sequence GC 双峰外源物种污染、样本串入抽样比对统计物种分布、看比对率确认污染源样本降级或重测Q30 低于 85%测序试剂问题、flowcell 老化看 Per tile 图定位区域联系测序方核查考虑重测比对率低于 95%参考基因组版本错、接头未去净、污染核对参考版本、复检接头、抽 read 看比对情况修正参考版本重跑补修剪重复率超过 35%文库复杂度低、起始量不足、PCR 循环过多查建库记录、同批次横向比标记去重后评估有效深度必要时重测on-target rate 低于 55%捕获失败、探针批次问题查捕获实验记录、对照同批次重新捕获不要硬跑fold-80 超过 2.0GC 偏好、捕获不均看 GC_DROPOUT、画深度分布图记录并在下游用 callable loci 过滤目标区平均深度不足数据量不够、捕获效率低算理论数据量需求补测数据或降低检出灵敏度预期FREEMIX 超过 3%样本交叉污染查样本处理记录、做 SNP 指纹重抽样本整批回溯操作流程性别与记录不符样本错配、标签贴错SNP 指纹两两比对整批回溯重新采样两样本指纹一致性 90%同一份样本重复上机、串样核对样本编号、查实验记录剔除重复样本重新计算批次5.2 几个印象深刻的坑第一个坑是平均深度够了但覆盖不够。有一批样本平均深度都在 130x 以上看着很漂亮但做变异检出的时候某些样本的检出位点数明显偏少。回头算PCT_TARGET_BASES_20X发现有几个样本只有 88%比批次中位数低了 8 个百分点。原因就是 fold-80 偏高深度分布极度不均——少数区域被过度测序大量区域深度不足。这说明平均深度这个指标有欺骗性必须配合覆盖比例和均一性一起看。后来我在流程里加了一条规则PCT_TARGET_BASES_20X低于 92% 就自动标黄必须在报告里说明。第二个坑是插入片段的长度阈值。前面提过-l 50在 FFPE 样本上丢了近两成数据。那次之后我改成按样本类型动态配置新鲜组织/血液样本用 50FFPE 或降解样本用 30同时看 fastp 报告里的插入片段峰值峰值低于 120bp 就自动切到宽松模式。这个改动看起来不起眼但对降解样本的数据利用率提升非常明显。第三个坑是 downsample 的必要性。同一个项目里有的样本测了 12Gb有的测了 8Gb深度差了 50%。这种情况下深度高的样本在变异检出时会有系统性优势导致样本间比较失真。我的做法是在 BAM 层面做 downsample把所有样本统一到相同的有效深度gatk DownsampleSam \ -I sample01.md.bam \ -O sample01.ds.bam \ -P 0.75 \ --RANDOM_SEED 42 \ --CREATE_INDEX true-P是保留比例按目标深度除以当前深度算。--RANDOM_SEED一定要固定否则每次跑结果不一样没法复现。这个操作在多样本联合分析比如 CNV、群体频率比较里尤其重要。第四个坑是不要用 FastQC 的PASS/FAIL直接做流程判断。FastQC 的判定逻辑很粗某些模块对 WES 数据天然就是 FAIL比如 Sequence Duplication Levels 和 Per base sequence content 的前几位这些 FAIL 完全正常不影响后续分析。如果流程按 FastQC 的 fail 数量来判断放行会导致大量合格样本被拦下。我的做法是忽略 FastQC 的判定只读原始数值用自己定义的阈值判断。第五个坑是质控报告的可追溯性。早期我没存中间文件只留了最终报告。后来有一次要复查半年前的一个样本发现当时用的 fastp 版本和参数都没记录只能凭记忆推断。现在我的做法是每次运行都把完整命令行、软件版本、参考基因组 md5 写进一个 log 文件和报告一起归档。这个习惯在复现问题和应对审计时价值极高。最后分享一个我一直在用的判断原则质控的目的不是让所有指标都变成绿色而是知道每个偏离正常值的数字意味着什么风险。有些偏离可以接受有些必须叫停区别在于这个偏离会不会影响下游结论的可靠性。搞不清这一点质控就会变成一种仪式看着很规范实际上什么都没防住。 
长河君 · 2026/11/3 6:20:25 · 阅读 1.2k 全外显子组测序数据质控：三层指标、阈值与实战排查
全外显子组测序的数据质控是个特别容易被走过场的环节。跑个 FastQC看到一排绿勾截图贴进报告收工。我刚开始接手流程搭建的时候也是这么干的直到有一批样本在变异检出阶段冒出一堆说不清的假阳性回头翻质控记录才发现那批样本的 on-target rate 只有 40% 出头fold-80 接近 3.0等于把一大半测序数据白扔在了目标区域外面。重新捕获来不及整个项目延期两周。从那以后我就把全外显子组的数据质控从一个看一眼的环节改成了流程里真正会拦人的闸门。流程走到第四步前面实验设计、建库测序、原始数据交接这些都干完了接下来要决定的是这批数据到底能不能往下走、哪些样本要重做、哪些指标可以睁一只眼闭一只眼。说白了质控这一步不产生任何科学结论但它决定了后面所有结论的可信度。这篇讲的是我在实操里怎么设计这一环三层指标分别看什么、每个数字背后的含义、参数为什么这么设、什么情况下必须叫停、以及那些看着吓人其实没事的数字怎么区分。不管你是刚接触生信流程、只会照着脚本敲命令的新手还是手里已经有一套管线、想把质控做成自动放行的老手应该都能从里面找到能直接抄走的东西。1. 质控在整个WES流程里的位置与设计思路1.1 为什么质控不是跑个FastQC就完事很多人对质控的理解停留在原始数据干净不干净这一层实际上 WES 的质控要分三个完全不同的层次来看每一层对应的问题类型、工具、判断标准都不一样。第一层是read 层面也就是刚下机的 fastq 文件。这一层关心的是碱基质量、接头残留、N 碱基比例、GC 分布、是否有其他物种污染。工具是 FastQC、fastp、MultiQC。这一层出问题通常是测序或者建库环节的锅。第二层是比对层面也就是比对到参考基因组之后的 BAM 文件。这一层关心的是比对率、重复率、插入片段长度分布、捕获效率、目标区域覆盖度和均一性。工具是 BWA-MEM、samtools、Picard/GATK、mosdepth、Qualimap。这一层出问题可能是捕获探针批次、文库复杂度、参考基因组版本不匹配。第三层是样本层面也就是把所有样本放一起横向看。这一层关心的是样本间污染、性别是否和记录一致、有没有串样、批次效应。工具是 verifybamid2、SNP 指纹比对、MultiQC 汇总。这一层出问题最要命因为它不会在单样本报告里显形只有横向比对才看得出来。我见过不少流程只做了第一层FastQC 全绿就放行结果第二层、第三层的问题一直到变异过滤阶段才暴露。所以我在自己的流程里把三层都做成强制卡点任何一层亮红灯就中断不允许先跑下去看看。1.2 三层质控的分工与放行标准为了让流程能自动判断我给每个指标都定了三种状态绿自动放行、黄记录并继续但要在报告里标注、红中断等人工确认。这个分级很关键因为如果把所有指标都设成硬阈值稍微偏一点就中断流程根本跑不下去反过来全部放行质控就失去意义。下面这张表是我目前用的核心指标和阈值人类样本、双端 150bp、常规 WES 捕获芯片的场景其他物种或者特殊建库方式要相应调整层级指标绿放行黄记录红中断readQ30 比例 90%85%–90% 85%read接头残留占比 0.5%0.5%–2% 2%readN 碱基比例 0.5%0.5%–1% 1%read修剪后保留 read 比例 95%90%–95% 90%比对比对率 98%95%–98% 95%比对重复率 20%20%–35% 35%比对on-target rate 70%55%–70% 55%比对平均目标深度 100x80–100x 80x比对≥20x 覆盖比例 95%90%–95% 90%比对fold-80 penalty 1.51.5–2.0 2.0样本污染率 FREEMIX 1%1%–3% 3%样本性别一致性一致—不一致这张表不是拍脑袋定的。Q30 90% 是 Illumina 平台近几年的常规水平低于 85% 基本可以认为测序环节出了问题。比对率 98% 对人类样本来说不算苛刻因为参考基因组质量好正常样本轻松能到 99% 以上。on-target rate 和捕获芯片的设计、探针批次强相关不同批次的芯片表现差异不小所以我把黄色区间设得比较宽。注意阈值一定要跟自己的平台、芯片、建库方式绑定。别人流程里的 95% 搬到你这儿可能是灾难因为建库方式、测序仪型号、甚至 flowcell 批次都会影响基线。我的建议是先用二十个以上已知合格的样本跑一遍统计出自己实验室的分布再定阈值。1.3 工具链选型的取舍逻辑质控工具的选择上我踩过的坑比想象中多。这里把几个关键选择和背后的理由说清楚省得你重复试错。FastQC 这个没得选它就是行业事实标准报告模块固定、解读资料多、MultiQC 原生支持。它的问题是不做任何处理只出报告所以修剪还得靠别的工具。修剪工具我在 fastp、Trimmomatic、cutadapt 之间反复横跳过。最后长期用 fastp原因有三个一是速度快C 写的多线程同样数据量比 Trimmomatic 快三到五倍这对大批量样本很关键二是双端接头能自动识别不需要你手动去填 adapter 序列省掉了接头序列填错导致修剪不干净这类低级错误三是单次运行同时输出 json 和 htmlMultiQC 直接能读报告链路顺。Trimmomatic 的优势是参数是 step 式的、每一步做什么非常明确可控性强老流程里常见如果你想精细控制每一步的剪切逻辑它仍然值得用。cutadapt 更适合处理非 Illumina 平台或者特殊接头。比对工具用 BWA-MEM这个和 GATK 的官方推荐一致。核心原因是 WES 里 indel 附近会有大量 soft-clip 的 readBWA-MEM 对局部比对和 soft-clip 的处理更成熟后续 indel 检出更稳。bowtie2 在短 read 场景更快但对长 indel 的支持一般我一般不用于 WES。覆盖度统计这块早期用 samtools depth后来全换成 mosdepth。原因是 samtools depth 单线程一个 30x 的 WES BAM 跑下来要几十分钟mosdepth 多线程加算法优化同样的活儿一两分钟出结果而且原生支持 BED 区间、窗口统计、阈值量化输出格式也规整。这个替换对整条流程的耗时改善非常明显。汇总报告用 MultiQC把所有样本的 FastQC、fastp、Picard、mosdepth 结果一次性拉到一个页面上横向比对特别直观。批次效应、离群样本往往就是在这个页面上第一眼发现的。2. 原始下机数据的质控从FastQC报告里读出真问题2.1 FastQC七个关键模块到底怎么看FastQC 一页报告十几个模块新手容易逐条看绿勾其实真正需要盯的就那么几个。我按重要性排一下。Per base sequence quality是最核心的一张图。横轴是 read 上的碱基位置纵轴是 Phred 质量值。Phred 分数 Q 和错误率的换算关系是Q -10 × log10(P)也就是 Q20 对应 1% 的错误率Q30 对应 0.1%。正常的双端 150bp 数据图形应该是前几个碱基稍微低一点这是随机引物起始带来的偏好性属于正常中间平稳末端缓慢下滑。如果末端出现断崖式下跌说明测序后期的信号衰减明显需要在修剪时加尾部滑窗剪切。如果整条 read 都低那就是测序本身的问题修剪救不回来。Per sequence GC content值得单独说。理论上这条曲线应该是一个接近正态的单峰峰的位置和参考基因组的 GC 含量吻合。如果出现双峰或者整个分布明显右移八成是污染——可能是细菌、真菌也可能是另一个物种的样本串进来了。这时候不要急着往下跑先做物种溯源。Sequence Duplication Levels里蓝线是理论分布、红线是实际分布。WES 因为要经过 PCR 扩增重复率天然比 PCR-free 的 WGS 高15%–25% 很常见不用紧张。真正要警惕的是红线远远跑到蓝线右边重复率超过 40%那说明文库复杂度太低可能是起始 DNA 量不够、PCR 循环数过多这种样本即使测了有效信息量也有限。Overrepresented sequences列出占比超过 0.1% 的序列。这里出现接头序列、引物序列、rRNA 片段都很正常但如果出现一段你完全陌生的序列且占比很高就要警惕外源污染。Per base N content一般应该贴着底部超过 1% 说明测序仪在某个位置读不出来通常是硬件或者试剂问题。Adapter Content是累积曲线理想状态是全贴底。这条曲线一旦抬头就必须修剪否则接头序列会被当成 read 比对到参考基因组上制造假阳性变异。Per tile sequence quality是按 flowcell 的 tile 分块看质量。如果某一块区域整片泛红说明那个位置的测序出了问题通常是气泡、堵塞或者试剂分布不均。这类问题属于硬件层面你没法处理但必须在报告里记录下来作为后续判断依据。2.2 接头污染与低质量碱基的处理策略接头污染的来源主要有两个一是插入片段长度小于测序读长导致 read 末端读进了接头序列二是接头二聚体。前者在插入片段分布偏短的文库比如 FFPE 样本、降解样本里特别常见。处理思路要分情况。如果只是零星的接头残留fastp 的自动检测加默认修剪就够了。如果 Adapter Content 曲线抬头很明显或者插入片段整体偏短就需要打开滑窗质量剪切把末端低质量区一并切掉。这里有个很实际的取舍修剪得越狠数据越干净但有效 read 也越少。我在实际操作中会看修剪后保留 read 比例这个指标如果低于 90%就要回头想想是不是参数设得过严而不是无脑加码。有些同行为了追求 FastQC 全绿把参数调得极狠结果丢了 20% 的数据得不偿失。还有一个容易忽略的点双端 read 的重叠区纠错。fastp 默认会对双端 read 的重叠区域做一致性校正把低质量的碱基用高质量的互补碱基覆盖掉。这个功能对降低整体错误率有帮助但对 WES 来说重叠区的校正会轻微改变碱基质量值的分布如果你后续要做严格的 BQSR需要留意一下。我的做法是默认开启因为这个校正带来的收益通常大于那点分布扰动。2.3 fastp的完整参数配置与逐项解释下面这套是我目前稳定在用的配置人类双端 150bp WES 场景可以直接抄fastp \ -i raw_R1.fastq.gz -I raw_R2.fastq.gz \ -o clean_R1.fastq.gz -O clean_R2.fastq.gz \ --detect_adapter_for_pe \ --cut_front --cut_tail \ --cut_window_size 4 --cut_mean_quality 20 \ -q 20 -u 30 -n 5 -l 50 \ --correction \ --thread 16 \ --json qc/sample.json --html qc/sample.html \ --report_title sample_fastp_report逐个拆一下为什么这么设。--detect_adapter_for_pe是双端接头自动检测。虽然 fastp 内置了 Illumina 常见接头序列但不同建库试剂盒的接头序列可能有细微差别自动检测比死填序列稳。这个参数一定要开尤其是你不确定建库试剂盒型号的时候。--cut_front --cut_tail --cut_window_size 4 --cut_mean_quality 20这一组是滑窗质量剪切。逻辑是从 read 两端各取一个 4bp 的窗口如果窗口内平均质量低于 Q20就把这个窗口切掉然后往前推进。开启cut_front是为了处理 read 开头那几个碱基的随机引物偏好性带来的低质量cut_tail处理末端衰减。窗口设 4 而不是更大是因为窗口越大越容易误伤4bp 的颗粒度对 150bp 的 read 比较合适。-q 20定义合格碱基的门槛-u 30表示一条 read 里不合格碱基超过 30% 就整条丢弃。这两个配合使用。为什么不用更严的 Q30因为 Q30 作为丢弃阈值太狠了会丢掉大量只是末端稍微差一点的 read而 Q20 对应的 1% 错误率已经能满足 WES 变异检出的需求。-n 5表示 N 碱基超过 5 个就丢弃。-l 50表示修剪后长度短于 50bp 的 read 直接丢弃。这个 50 是经验值正常 WES 插入片段在 150–250bp修剪后还能留 100bp 以上所以 50 基本不会误伤。但如果你做的是 FFPE 或者降解严重的样本插入片段本身就短这时候要降到 30否则会丢掉大量有效 read。我吃过这个亏一批 FFPE 样本用了-l 50直接丢了快两成数据。--correction开启重叠区纠错。前面说了默认开。--thread 16线程数按机器核数调一般设到物理核数就行设太高反而因为内存带宽成为瓶颈。还有两个参数我特意没用--dedup和--trim_poly_g。--dedup是 fastp 自带的序列级去重它会直接丢掉重复的 read 对。这个做法我强烈不推荐原因在下一章展开。--trim_poly_g是修剪 polyG 尾巴这是 NextSeq/NovaSeq 双色荧光平台的特有现象如果你用的是这类平台需要加上如果是四色平台的旧机型就不用。提示修剪完一定要重新跑 FastQC和修剪前的报告并排放确认接头曲线下去了、末端质量抬上来了。我见过有人改完参数没复检等于白改。3. 比对层面的质控比对率、重复率、捕获效率3.1 BWA-MEM比对与read group的坑比对这一步重点不是跑通而是跑对。BWA-MEM 本身的参数不多但 read group 这个事必须说透。bwa mem -t 16 -K 100000000 \ -R RG\tID:sample01\tSM:sample01\tPL:ILLUMINA\tLB:lib01\tPU:unit01 \ ref.fa clean_R1.fastq.gz clean_R2.fastq.gz \ | samtools sort - 8 -m 2G -o sample01.sorted.bam - samtools index sample01.sorted.bam-R里的SM字段是样本名ID是 read group 标识LB是文库名PL是测序平台。这几个字段里SM最关键因为 GATK 后续会按SM把 read 分组判断哪些 read 属于同一个样本。如果SM写错了会出现两种典型灾难一是同一个样本的多个文库被当成不同样本导致变异检出时把它当混合样本处理二是不同样本的SM被写成一样导致所有样本被合并最后出的 VCF 只有一条样本记录。这个坑我在早期项目里踩过一次当时是批量脚本里的变量替换出了错几百个样本的SM全成了一个值。表面上看流程跑得顺顺利利直到变异检出的 VCF 打开只有一列样本才发现不对全部重跑。-K参数控制每次批处理加载到内存的碱基数默认是 1000 万。这个值会影响内存占用和比对的批次划分设大了内存吃得凶设小了批次多、速度稍慢。我一般设成 1 亿配合足够的线程数速度比较理想。-t 16是线程数。这里有个实操建议BWA-MEM 的线程效率不是线性的超过 16 线程后提升不明显而 samtools sort 是另一个瓶颈。我的做法是 BWA 用 16 线程管道后接 samtools sort 用 8 线程、2G 内存缓冲整体吞吐比较平衡。3.2 重复率的正确解读与去重时机重复率这个指标特别容易被误读。先说清楚它是怎么来的。WES 的重复 reads 有两个来源一是 PCR 扩增偏倚同一个原始 DNA 分子被扩增了多份测出来是同一条序列二是测序深度足够高时随机取样的碰撞。对于 100x 以上的深度第二种来源占比不低这部分重复是数学上不可避免的不代表文库质量差。所以看到一个 25% 的重复率不要慌。判断标准要看同批次比较如果同批次样本普遍在 15%–20%某个样本突然 40%那说明这个样本的文库复杂度确实低值得回头看建库记录。如果全批次都是 35%那更可能是芯片/测序深度/建库方式带来的系统性特征不是单样本问题。去重的工具我用 GATK 的 MarkDuplicates底层就是 Picardgatk MarkDuplicates \ -I sample01.sorted.bam \ -O sample01.md.bam \ -M qc/sample01.md_metrics.txt \ --CREATE_INDEX true \ --VALIDATION_STRINGENCY LENIENT注意这里是Mark标记而不是Remove移除。标记的意思是给重复 read 打上 flag但不从 BAM 里删掉。为什么不删因为后续的变异检出环节需要这些重复 read 来提供支持证据尤其是低等位频率的位点。GATK 的 HaplotypeCaller 默认会忽略标记为重复的 read但保留它们在文件里便于回溯和可视化检查。直接删掉就失去了这层信息。去重的时机也很关键。标准顺序是比对 → 排序 → 标记重复 → BQSR碱基质量值重校准→ 变异检出。去重必须在 BQSR 之前因为重复 read 会让同一个位置的碱基被重复计数干扰质量值的经验估计。这个顺序如果搞反了BQSR 的校正结果会有系统性偏差虽然不一定报错但会影响后续的变异质量。还有一个常见误区试图在 fastq 层面用 fastp 的--dedup去重。这个做法的问题是fastp 只能基于序列完全一致来判断重复而真正的重复判断要基于比对坐标——两条序列不完全一样但比对到同一个起始位置的 read才是真正的 PCR 重复。所以 fastq 层去重会漏掉大量真重复同时误删一些序列相同但实际来自不同基因组位置的 read比如重复区域。我的做法是 fastq 层完全不去重全部交给 BAM 层的 MarkDuplicates。3.3 捕获效率与覆盖度的计算实操WES 和 WGS 最大的区别就在这一步你花钱测的是全基因组但真正有用的只是那 1%–2% 的编码区。所以捕获效率on-target rate直接决定了这笔钱花得值不值。计算 on-target rate 最直接的方式是用 GATK 的 CalculateHsMetricsgatk CalculateHsMetrics \ -I sample01.md.bam \ -O qc/sample01.hs_metrics.txt \ -R ref.fa \ --BAIT_INTERVALS bait.interval_list \ --TARGET_INTERVALS target.interval_list \ --PER_TARGET_COVERAGE qc/sample01.per_target_coverage.txt这里要区分两个概念bait是探针覆盖的区域target是你真正关心的区域。探针设计时 bait 通常比 target 略大因为要留出捕获的余量。计算 on-target rate 用的是 bait计算覆盖度用的是 target。把这两个搞混是新手常见错误会导致指标口径不一致不同批次没法横向比。输出文件里有几个关键字段值得重点看字段含义参考范围PCT_SELECTED_BASESon-target rate 0.70MEAN_TARGET_COVERAGE目标区平均深度 100PCT_TARGET_BASES_20X≥20x 覆盖的靶碱基比例 0.95FOLD_80_BASE_PENALTY覆盖均一性 1.5PCT_PF_READS_ALIGNED比对率 0.98AT_DROPOUT/GC_DROPOUTAT/GC 偏好导致的覆盖缺失越低越好ZERO_CVG_TARGETS_PCT完全无覆盖的靶区比例 0.02FOLD_80_BASE_PENALTY这个指标值得单独解释因为它最容易被忽略但信息量最大。它的计算方式是把所有靶碱基按覆盖深度排序取第 20 百分位的深度值用平均深度除以这个值。理想情况是 1.0意思是 80% 的靶区域深度都达到了平均值。如果算出来是 2.5意味着 80% 的区域深度只达到平均值的 40%覆盖严重不均——少数区域深度极高大量区域深度很低。这种样本即使平均深度看着漂亮比如 150x实际有效覆盖可能很差。AT_DROPOUT和GC_DROPOUT反映的是捕获探针的序列偏好。GC 含量极高的区域比如启动子区、第一外显子捕获效率天然低这是技术限制。如果这两个值异常高说明该样本在高 GC 或高 AT 区域的覆盖缺失明显后续做 CNV 分析时要特别小心。除了 CalculateHsMetrics我还会用 mosdepth 做一次交叉验证因为它速度更快而且能输出窗口级的深度分布便于画图mosdepth \ -t 8 \ -b target.bed \ -n \ --by 200 \ --quantize 0:1:10:20:30:50:100: \ qc/sample01 \ sample01.md.bam--by 200表示按 200bp 窗口统计--quantize定义了深度分档输出文件里会直接给出每一档覆盖的碱基比例。-n表示不输出 per-base 文件那个文件几十 GB一般不需要。跑完会得到sample01.mosdepth.summary.txt和几个 regions 文件summary 里有总平均深度、各染色体的平均深度regions 文件里是每个窗口的深度。如果你想自己算一份定制化的报告用 mosdepth 的窗口输出配合一段 Python 就够了import gzip def coverage_summary(regions_gz, thresholds(1, 10, 20, 30, 50)): depths [] with gzip.open(regions_gz, rt) as fh: for line in fh: parts line.rstrip(\n).split(\t) depths.append(float(parts[-1])) n len(depths) mean sum(depths) / n res {windows: n, mean_depth: round(mean, 2)} for t in thresholds: hit sum(1 for d in depths if d t) res[fpct_ge_{t}x] round(hit / n * 100, 2) ordered sorted(depths) p20 ordered[int(n * 0.2)] res[fold80] round(mean / p20, 2) if p20 0 else None return res print(coverage_summary(qc/sample01.regions.bed.gz))这段脚本按窗口算不是按碱基算所以算出来的百分比和 CalculateHsMetrics 会有一点差异但趋势一致用来做快速筛查足够了。如果要严格按碱基算得读 per-base 文件那个开销大得多我一般只在对某个样本有疑问的时候才做。注意用samtools depth算平均深度有一个经典陷阱。不加-a参数时它只输出有 read 覆盖的位置没覆盖的位置压根不出现在结果里。这时候你拿总深度除以行数算出来的平均值会系统性偏高因为零覆盖的位置被排除了。必须加-a把零覆盖位置也输出或者干脆换 mosdepth。我第一次踩这个坑的时候一个实际只有 60x 的样本算出来是 95x白高兴一场。4. 样本层面的质控污染、性别与一致性4.1 污染率估算看不见的杀手样本间污染是质控里最隐蔽的问题因为它不会让任何单个指标明显异常但会系统性地污染下游所有变异检出结果。轻度污染的样本低等位频率的假阳性会明显增多重度污染的样本基因型判断直接失真纯合会被判成杂合。估算污染我用 verifybamid2它基于一组公共的常见 SNP 位点通过比对 read 上的等位基因频率分布来反推污染比例verifybamid2 \ --Reference ref.fa \ --BamFile sample01.md.bam \ --SVDPrefix 1000g.phase3.100k.b38.vcf.gz.dat \ --Output qc/sample01.selfSM \ --Verbose \ --NumThread 8输出文件里的FREEMIX就是污染率估计值。我的判断标准是低于 1% 视为干净1%–3% 要记录并在报告里标注超过 3% 就要考虑样本降级或者重抽。超过 5% 基本可以判定污染这个样本的变异检出结果不能直接用于临床或科研结论。CHIPMIX是另一个字段用来估计是否存在同一个体的另一种组织混入比如血液样本里混了口腔上皮。这个值异常的情况比较少见但如果出现值得回头查样本采集记录。这里要说清楚一点verifybamid2 需要一份参考的 SNP 位点集合和对应的 SVD 模型文件这些属于公开资源具体来源按你所在环境的规范获取不要随便网上下一个来源不明的文件。模型文件和参考基因组的版本必须严格对应版本不一致算出来的污染率会是错的。4.2 性别核查与样本指纹比对性别核查是个成本极低但价值极高的检查。原理很简单男性只有一条 X 染色体Y 染色体有覆盖女性有两条 X 染色体X 上的杂合位点比例接近正常常染色体Y 染色体基本无覆盖。实操上有两种做法。一种是看覆盖度统计 X 和 Y 染色体的平均深度除以常染色体的平均深度做归一化。正常男性 X 的相对深度约 0.5Y 约 0.5正常女性 X 约 1.0Y 接近 0。这个方法快但对覆盖度不均一的样本会有波动。另一种是看杂合度统计 X 染色体非拟常染色体区域上的杂合 SNP 比例。女性应该在 0.3%–0.5% 这个量级男性接近 0。这个方法更准但需要先跑一遍变异检出。我一般两个都做交叉验证。如果两个结果都指向同一个性别但和登记记录不符那就是样本错配必须整批回溯。样本指纹比对是另一个必做项尤其在多样本同批上机的情况下。做法是挑一组高多态性的 SNP 位点从每个样本的 BAM 里提取基因型两两计算一致性。同批次内除了已知的重复样本其他样本之间的一致性应该在 50%–60% 左右随机两个无关个体的一致性就在这个量级。如果两个样本的一致性超过 90%那基本可以确定是同一个样本被测了两次或者发生了串样。这一步的价值在于它能发现前两层质控完全看不出来的问题。两个样本各自的数据都很漂亮FastQC 全绿、比对率 99%、覆盖度 120x但它们其实是同一个人的 DNA——这种问题只有横向比对才能发现。4.3 质控汇总把散落的指标收成一张表单样本的报告再详细也没法替代横向汇总。我的做法是每个样本跑完质控后把关键指标写进一个统一的 TSV最后用 MultiQC 或者自建的脚本做可视化。字段大致是这样import pandas as pd fields [ sample_id, raw_reads, clean_reads, clean_ratio, q30_before, q30_after, adapter_pct, mapping_rate, dup_rate, insert_size_median, on_target_rate, mean_depth, pct_20x, pct_30x, fold80, gc_dropout, at_dropout, freemix, sex_check, fingerprint_cluster ] df pd.DataFrame(columnsfields) df.to_csv(qc/batch_summary.tsv, sep\t, indexFalse)这张表的作用是让离群样本无处遁形。我会按mean_depth和fold80画散点图正常样本会聚成一团离群的点一眼就能看出来。也会按on_target_rate排序看看有没有明显低于批次中位数的。这里有个经验批次效应比单样本异常更常见也更容易被忽略。如果某一批样本的 on-target rate 集体比上一批低 10 个百分点那大概率是捕获探针批次或者实验操作的问题而不是每个样本各自出了状况。这种情况下单看每个样本的报告都是黄灯但不严重横向一看才知道是系统性问题。所以汇总这一步不是可选项是必做项。提示把历史批次的汇总表也留着做成累积视图。新批次进来的时候和过去半年的数据比一比能提前发现试剂批次、操作人员、设备状态的变化。这个习惯帮我躲过好几次批量返工。5. 常见问题排查与踩坑实录5.1 问题速查表下面这张表是我这些年攒下来的基本覆盖了 WES 质控里八成以上的异常情况现象可能原因排查动作处理方式Adapter Content 曲线抬头插入片段短于读长、接头二聚体看插入片段分布、Overrepresented sequences开启自动接头检测 尾部剪切必要时缩短-lPer sequence GC 双峰外源物种污染、样本串入抽样比对统计物种分布、看比对率确认污染源样本降级或重测Q30 低于 85%测序试剂问题、flowcell 老化看 Per tile 图定位区域联系测序方核查考虑重测比对率低于 95%参考基因组版本错、接头未去净、污染核对参考版本、复检接头、抽 read 看比对情况修正参考版本重跑补修剪重复率超过 35%文库复杂度低、起始量不足、PCR 循环过多查建库记录、同批次横向比标记去重后评估有效深度必要时重测on-target rate 低于 55%捕获失败、探针批次问题查捕获实验记录、对照同批次重新捕获不要硬跑fold-80 超过 2.0GC 偏好、捕获不均看 GC_DROPOUT、画深度分布图记录并在下游用 callable loci 过滤目标区平均深度不足数据量不够、捕获效率低算理论数据量需求补测数据或降低检出灵敏度预期FREEMIX 超过 3%样本交叉污染查样本处理记录、做 SNP 指纹重抽样本整批回溯操作流程性别与记录不符样本错配、标签贴错SNP 指纹两两比对整批回溯重新采样两样本指纹一致性 90%同一份样本重复上机、串样核对样本编号、查实验记录剔除重复样本重新计算批次5.2 几个印象深刻的坑第一个坑是平均深度够了但覆盖不够。有一批样本平均深度都在 130x 以上看着很漂亮但做变异检出的时候某些样本的检出位点数明显偏少。回头算PCT_TARGET_BASES_20X发现有几个样本只有 88%比批次中位数低了 8 个百分点。原因就是 fold-80 偏高深度分布极度不均——少数区域被过度测序大量区域深度不足。这说明平均深度这个指标有欺骗性必须配合覆盖比例和均一性一起看。后来我在流程里加了一条规则PCT_TARGET_BASES_20X低于 92% 就自动标黄必须在报告里说明。第二个坑是插入片段的长度阈值。前面提过-l 50在 FFPE 样本上丢了近两成数据。那次之后我改成按样本类型动态配置新鲜组织/血液样本用 50FFPE 或降解样本用 30同时看 fastp 报告里的插入片段峰值峰值低于 120bp 就自动切到宽松模式。这个改动看起来不起眼但对降解样本的数据利用率提升非常明显。第三个坑是 downsample 的必要性。同一个项目里有的样本测了 12Gb有的测了 8Gb深度差了 50%。这种情况下深度高的样本在变异检出时会有系统性优势导致样本间比较失真。我的做法是在 BAM 层面做 downsample把所有样本统一到相同的有效深度gatk DownsampleSam \ -I sample01.md.bam \ -O sample01.ds.bam \ -P 0.75 \ --RANDOM_SEED 42 \ --CREATE_INDEX true-P是保留比例按目标深度除以当前深度算。--RANDOM_SEED一定要固定否则每次跑结果不一样没法复现。这个操作在多样本联合分析比如 CNV、群体频率比较里尤其重要。第四个坑是不要用 FastQC 的PASS/FAIL直接做流程判断。FastQC 的判定逻辑很粗某些模块对 WES 数据天然就是 FAIL比如 Sequence Duplication Levels 和 Per base sequence content 的前几位这些 FAIL 完全正常不影响后续分析。如果流程按 FastQC 的 fail 数量来判断放行会导致大量合格样本被拦下。我的做法是忽略 FastQC 的判定只读原始数值用自己定义的阈值判断。第五个坑是质控报告的可追溯性。早期我没存中间文件只留了最终报告。后来有一次要复查半年前的一个样本发现当时用的 fastp 版本和参数都没记录只能凭记忆推断。现在我的做法是每次运行都把完整命令行、软件版本、参考基因组 md5 写进一个 log 文件和报告一起归档。这个习惯在复现问题和应对审计时价值极高。最后分享一个我一直在用的判断原则质控的目的不是让所有指标都变成绿色而是知道每个偏离正常值的数字意味着什么风险。有些偏离可以接受有些必须叫停区别在于这个偏离会不会影响下游结论的可靠性。搞不清这一点质控就会变成一种仪式看着很规范实际上什么都没防住。