ngs-quality-control

NGS 质量控制

概述

此技能将一套可复现的 QC 工作流应用于 Illumina 短读长测序数据的完整生命周期：原始 reads → 修剪后的 reads → 比对后的 BAM。它会生成确定性指标，适用于 WGS、WES 和靶向 panel 项目的 QC 门控。

流水线阶段：

1. 原始 reads QC — 对每个 FASTQ 运行 FastQC，并使用 MultiQC 汇总。
2. 接头/质量修剪 — 使用 fastp 检测双端接头。
3. 修剪后 QC — 对修剪后的 reads 运行 FastQC + MultiQC。
4. 比对 QC — samtools stats/flagstat、Picard 指标套件。
5. 覆盖度 QC — mosdepth（基因组窗口或目标 BED）。
6. 汇总 — 对所有工具输出运行 MultiQC，生成单页 HTML 摘要。

用法

假设 THREADS=8、OUT=qc_out，双端 FASTQ 文件为 sample_R1.fastq.gz / sample_R2.fastq.gz，坐标排序的 sample.bam 及其索引、参考基因组 ref.fasta，以及（对于 WES/panel）目标区间。

1. 原始 FASTQ QC — FastQC

mkdir -p $OUT/fastqc_raw
fastqc -t $THREADS -o $OUT/fastqc_raw *.fastq.gz

需要检查的 FastQC 模块：

• Per base sequence quality — 3′ 端的 Q30 下滑表示需要修剪。
• Per sequence quality scores — 双峰分布可能表示某个 tile/lane 质量不佳。
• Per base sequence content — 前约 12 bp 的偏差属于正常情况（随机六聚体引物）；其他位置的大幅偏差可能表示接头/污染。
• Per sequence GC content — GC 含量偏离物种的预期值 → 污染或偏差。
• Sequence Duplication Levels — 非扩增子文库中的高重复率 → 文库复杂度问题。
• Overrepresented sequences — 接头残留、rRNA、PhiX 残留。
• Adapter Content — 确认需要修剪的接头类型。

2. 接头与质量修剪 — fastp

mkdir -p $OUT/fastp
fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o $OUT/fastp/sample_R1.trim.fastq.gz -O $OUT/fastp/sample_R2.trim.fastq.gz \
  --detect_adapter_for_pe \
  --thread $THREADS \
  --html $OUT/fastp/sample.fastp.html \
  --json $OUT/fastp/sample.fastp.json

有用的附加参数：

• --qualified_quality_phred 20 — 碱基质量阈值（默认值为 15）。
• --length_required 50 — 丢弃修剪后短于 50 bp 的 reads。
• --cut_tail --cut_tail_window_size 4 --cut_tail_mean_quality 20 — 使用滑动窗口修剪 3′ 端低质量序列。
• --dedup — 在内存中去除 PCR 重复序列（按 read pair 身份判断；大多数工作流更适合在比对后使用 Picard/MarkDuplicates）。

3. 修剪后 QC + MultiQC

mkdir -p $OUT/fastqc_trim
fastqc -t $THREADS -o $OUT/fastqc_trim $OUT/fastp/*.trim.fastq.gz

multiqc $OUT -o $OUT/multiqc

MultiQC 会自动检测扫描目录中的 FastQC、fastp、Picard、samtools 和 mosdepth 输出。

4. 比对 QC — samtools

samtools flagstat -@ $THREADS sample.bam > $OUT/samtools/sample.flagstat.txt
samtools stats    -@ $THREADS sample.bam > $OUT/samtools/sample.stats.txt
samtools idxstats sample.bam               > $OUT/samtools/sample.idxstats.txt

从 samtools stats 提取：raw total sequences、reads mapped、error rate、insert size average、insert size standard deviation。

5. Picard 指标套件

Picard 需要在 ref.fasta 旁有序列字典（picard CreateSequenceDictionary R=ref.fasta）；对于杂交捕获指标，还需要 Picard 格式的区间列表（.interval_list，其头部 = BAM @SQ 行；不能是普通 BED 文件）。

# Convert BED → interval_list (required for CollectHsMetrics)
picard BedToIntervalList \
  I=targets.bed O=targets.interval_list SD=ref.dict

# General alignment metrics (all library types)
picard CollectAlignmentSummaryMetrics \
  R=ref.fasta I=sample.bam O=$OUT/picard/sample.alignment_metrics.txt

# Insert size distribution (paired-end)
picard CollectInsertSizeMetrics \
  I=sample.bam O=$OUT/picard/sample.insert_size.txt \
  H=$OUT/picard/sample.insert_size.pdf

# WGS coverage metrics
picard CollectWgsMetrics \
  R=ref.fasta I=sample.bam O=$OUT/picard/sample.wgs_metrics.txt

# WES / targeted panel hybrid-selection metrics
picard CollectHsMetrics \
  R=ref.fasta I=sample.bam O=$OUT/picard/sample.hs_metrics.txt \
  BAIT_INTERVALS=baits.interval_list \
  TARGET_INTERVALS=targets.interval_list

6. 覆盖度 — mosdepth

# WGS: 500 bp fixed windows
mosdepth --by 500 --threads $THREADS --fast-mode --no-per-base \
  $OUT/mosdepth/sample_wgs sample.bam

# WES / panel: per-target coverage
mosdepth --by targets.bed --threads $THREADS --no-per-base \
  $OUT/mosdepth/sample_wes sample.bam

需要关注的输出：

• *.mosdepth.summary.txt — 每条染色体和总体的平均深度。
• *.mosdepth.region.dist.txt — 覆盖度达到 ≥X 的区域累计比例（用于计算 % targets ≥20x）。
• *.regions.bed.gz — 每个窗口/目标的平均深度。

添加 --quantize 0:1:10:30: 可按覆盖度分箱，或添加 --thresholds 10,20,30,50，输出这些深度下每个区域的通过计数。

响应格式

• 先给出用户需要的命令或代码，再进行解释
• 结构为：确认输入 → 可运行代码 → 解释关键参数 → 注意事项
• 每项任务提供一个完整的可运行示例；不要展示所有替代方案
• 代码注释应简洁且实用（说明做什么，而非为什么存在）
• 目标：50–100 行代码，并附简短说明

核心概念

按测序类型选择 QC 工具

测序类型	覆盖度工具	覆盖度指标	Picard 套件	关键阈值
WGS	mosdepth --by 500	CollectWgsMetrics	PCT_20X ≥ 95%	平均深度 ≥ 30×
WES	mosdepth --by targets.bed	CollectHsMetrics	PCT_TARGET_BASES_20X ≥ 80%	靶向区域比例 ≥ 60%
靶向面板	mosdepth --by targets.bed	CollectHsMetrics	PCT_TARGET_BASES_20X ≥ 95%	平均深度 ≥ 500×（体细胞）
RNA-seq	—（使用 featureCounts/STAR 日志）	CollectRnaSeqMetrics	PCT_CODING_BASES	rRNA < 10%

所有类型： FastQC（原始数据 + 修剪后数据）→ fastp（修剪）→ samtools flagstat → MultiQC（汇总）。

质量编码检测

现代 Illumina（CASAVA ≥1.8，2011 年起）输出 Phred+33（Sanger）。旧版数据（Illumina 1.3–1.7、Solexa）使用 Phred+64。比对工具需要 Phred+33。

• FastQC 会在 Basic Statistics → Encoding 中报告编码。查看 Sanger / Illumina 1.9。

• 快速 CLI 检查 — 查看质量行的 ASCII 范围：

  zcat sample_R1.fastq.gz | head -4000 \
    | awk 'NR%4==0' \
    | od -An -c | tr -s ' \n' '\n' \
    | awk 'length==1' | sort -u
  

  字符位于 !–I 范围内 → Phred+33；字符为 @–h → Phred+64。

• 转换旧版数据：seqtk seq -Q64 -V in.fq > out.fq（将 +64 转换为 +33）。

修剪参数的用途

• 接头修剪 用于读段延伸到接头序列的情况（插入片段长度 < 读长）。--detect_adapter_for_pe 利用读段重叠信息，对于双端测序，其可靠性高于基于序列的检测。
• 质量修剪 用于双通道化学体系下 3′ 端质量衰减的情况（NextSeq/NovaSeq 在暗周期会输出 poly-G；这些平台请使用 --trim_poly_g，fastp 会根据平台检测自动启用该选项）。
• 除非了解文库设计（UMI、内嵌条形码），否则不要按固定长度硬截断。

重复：光学重复与 PCR 重复

• samtools flagstat 中的 “duplicates” 反映了 BAM 中设置的标记（通常由 Picard MarkDuplicates 或 samtools markdup 设置）。应在比对后运行 MarkDuplicates，而不是对修剪后的 FASTQ 运行。
• CollectWgsMetrics / CollectHsMetrics 会报告 PCT_EXC_DUPE — 因重复而排除的碱基比例。这是基于覆盖度的 QC 所采用的标准重复率。

覆盖度：平均深度与均一性

仅报告平均深度并不足够。应同时报告：

• 靶向区域平均覆盖度（HsMetrics 中的 MEAN_TARGET_COVERAGE；WgsMetrics 中的 MEAN_COVERAGE）。
• 覆盖广度 — % targets ≥20x / % genome ≥20x（Picard 中的 PCT_TARGET_BASES_20X、PCT_20X；mosdepth 的 region.dist.txt）。

平均深度为 200x，但只有 70% 的靶向区域达到 20x 的面板存在覆盖均一性问题（捕获效率、GC 偏倚）。

快速参考

质量阈值（常见门控标准）

指标	WGS	WES	靶向面板	来源
平均覆盖度	≥ 30×	靶向区域 ≥ 50×	≥ 500×（体细胞），≥ 100×（胚系）	WgsMetrics / HsMetrics
20× 覆盖广度	基因组 ≥ 95%	靶向区域 ≥ 80%	靶向区域 ≥ 95%	Picard PCT_20X / mosdepth
10× 覆盖广度（胚系）	≥ 90%	靶向区域 ≥ 90%	靶向区域 ≥ 98%	Picard PCT_10X
重复率	< 30%	< 40%	< 50%（建议使用 UMI）	MarkDuplicates / Picard PCT_EXC_DUPE
比对率	> 95%	> 95%	> 95%	samtools flagstat / AlignmentSummaryMetrics
平均碱基质量	≥ 80% 的碱基达到 Q30	同左	同左	FastQC / fastp
插入片段长度（PE150）	300–500 bp，SD < 平均值的 20%	150–300 bp	依文库而定	CollectInsertSizeMetrics
靶向区域比例（捕获）	不适用	≥ 60%（PCT_SELECTED_BASES）	≥ 70%	CollectHsMetrics
污染率	< 1%	< 1%	< 1%	VerifyBamID（不在此技能范围内）

单行命令

# Full raw → trimmed → FastQC → MultiQC
fastqc -t 8 -o qc/raw *.fastq.gz && \
fastp -i R1.fq.gz -I R2.fq.gz -o R1.trim.fq.gz -O R2.trim.fq.gz \
      --detect_adapter_for_pe --thread 8 --html fastp.html --json fastp.json && \
fastqc -t 8 -o qc/trim *.trim.fq.gz && \
multiqc qc -o qc/multiqc

# BAM QC bundle
samtools flagstat sample.bam > qc/sample.flagstat.txt
picard CollectAlignmentSummaryMetrics R=ref.fasta I=sample.bam O=qc/aln.txt
picard CollectInsertSizeMetrics I=sample.bam O=qc/ins.txt H=qc/ins.pdf
mosdepth --by 500 --threads 8 --fast-mode --no-per-base qc/sample sample.bam

工具角色一览

工具	输入	主要输出	用途
FastQC	FASTQ	每个文件对应的 HTML/zip	原始和剪切后 reads 的质量
fastp	FASTQ 双端数据	剪切后的 FASTQ + JSON/HTML	接头和质量剪切
samtools stats/flagstat	BAM	文本	比对计数、错误率
Picard AlignmentSummary	BAM	文本	链平衡、PF reads、Q20/Q30 比对碱基
Picard InsertSize	BAM	文本 + PDF	PE 插入片段分布
Picard WgsMetrics	BAM	文本	WGS 覆盖度 + 覆盖广度
Picard HsMetrics	BAM	文本	WES/Panel 靶区捕获率和覆盖度
mosdepth	BAM	BED + dist	快速窗口化／逐靶区深度
MultiQC	目录	单个 HTML	队列汇总

常见错误

• 错误： 只在剪切后运行 FastQC 正确： 对原始和剪切后的 reads 都运行 FastQC 原因： 剪切前 QC 可诊断文库质量；剪切后 QC 可确认剪切有效且没有过度剪切

• 错误： 将原始 BED 文件传给 Picard CollectHsMetrics 正确： 通过 picard BedToIntervalList I=x.bed O=x.interval_list SD=ref.dict 转换为 Picard 格式的 .interval_list 原因： BAIT_INTERVALS/TARGET_INTERVALS 需要包含完整 @SQ 头部的 .interval_list；原始 BED 会导致失败，或悄无声息地产生错误指标

• 错误： 使用 contig 名称或排序顺序不匹配的 interval-list 文件 正确： 确保 contig 名称（chr1 与 1）和 @SQ 顺序与 BAM 头部完全一致 原因： 不匹配会导致指标错误而不报错，或造成工具运行失败

• 错误： 对剪切后的 FASTQ 或按名称排序的 BAM 运行 MarkDuplicates 正确： 仅对按坐标排序的已比对 BAM 运行 MarkDuplicates 原因： MarkDuplicates 需要按坐标排序的已比对 reads，才能识别光学重复和 PCR 重复

• 错误： 仅将平均覆盖度作为 QC 指标报告 正确： 始终将平均深度与覆盖广度（% ≥20x）及重复率一起报告 原因： 均一性较差时，较高的平均值会掩盖未达标的靶区

• 错误： 对 2012 年之前的 FASTQ 数据信任自动编码检测 正确： 对古老数据进行比对前，明确验证 Phred 编码 原因： 将 Phred+64 数据传给现代比对器会导致碱基质量值悄无声息地出错

• 错误： 在运行每个 lane 的 QC 前合并 lane／样本 正确： 先对每个 lane／FASTQ 运行 FastQC；每个 lane 的 QC 通过后再合并 原因： 合并后，单个质量不佳的 lane 将无法识别

• 错误： 将 fastp 中的 --dedup 作为主要去重步骤 正确： 比对后使用 Picard MarkDuplicates 或 samtools markdup 去重 原因： fastp 的 --dedup 仅能在比对前捕获序列完全相同的重复，无法识别基于比对结果的重复

• 错误： 在大型基因组上运行 mosdepth 时不使用 --no-per-base 正确： 对大型基因组进行 QC 时始终使用 --no-per-base 原因： 逐碱基输出量极大，且很少需要；窗口／区域输出足以用于 QC

参考资料

• FastQC：https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
• fastp：Chen 等，Bioinformatics 2018，https://doi.org/10.1093/bioinformatics/bty560
• Picard：https://broadinstitute.github.io/picard/
• mosdepth：Pedersen 和 Quinlan，Bioinformatics 2018，https://doi.org/10.1093/bioinformatics/btx699
• 覆盖度建议：Sims 等，Nat Rev Genet 2014，https://doi.org/10.1038/nrg3642