想象你是一座超级邮局的质检主任。每天,测序仪这台“基因邮筒”会吐出上亿封基因密信——每一封信都有四个部分:以 @ 开头的信封编号,由 A、T、C、G 写成的正文,一行 + 分隔线,以及一串像乱码一样的“邮戳质量符”。这些密信必须经过质检员抽查(FASTQ 质控)、分拣员按参考地图归档(比对生成 BAM)、侦探逐字核对找出拼写差异(变异检测),最后整理成一份标准问题清单(VCF)。本文会带你以这条“信使流水线”为线索,读懂 FASTQ、BAM、VCF 三种数据格式,并掌握从原始测序数据到可靠变异检测的完整质控流程。

岗位关联

掌握这套标准格式与质控流程,是进入生信分析与精准医疗领域的入场券。学会它,你就从“只会调用工具”升级为“能诊断数据质量、设计可靠流水线”的数据诊断专家。

目录
  1. 可视化演示
  2. 1. 测序与FASTQ格式
  3. 2. FASTQ质控
  4. 3. 比对与BAM
  5. 4. BAM质控
  6. 5. 变异与VCF
  7. 6. 变异过滤流程
  8. AI 应用场景
  9. Python 演示
  10. 常见误区
  11. 面试题

可视化演示

展示从FASTQ原始测序数据到BAM比对再到VCF变异检测的完整质控流程,交互切换查看各阶段关键指标。

本节配有交互动效(进入视口后自动加载,约 3KB)。

'+modeData.fastq.caption+'

1. 测序与FASTQ格式

邮局每天都会收到海量“基因密信”,每一封信都严格遵循四行格式。第一行 @ 后面写着信封编号;第二行是正文,由 A、T、C、G 四种字母排列而成;第三行是一个 + 分隔线;第四行是一串与正文长度完全相等的字符,像是信封上的邮戳——但这里的邮戳不是日期,而是每个字母的可信度。

在生信世界里,这种“密信”格式叫做 FASTQ。第二行的字母串是测得的 DNA 序列,第四行的每个字符则对应一个 Phred 质量分数。字符的 ASCII 码减去 33 就是质量值 Q。Q 值越高,说明这个位置被读错的可能性越小。测序仪在读出每个碱基时都会附加这样一条质量信息,让下游分析者知道哪些位置“字迹模糊”。

设想开信时,某个位置的字母 A 可能被误读成 C。质量分数 Q 与错误概率 P 之间有清晰的数学关系:Q 越大,错误概率越小。业内常说的 Q20、Q30,正是测序仪出厂质检的通行证。

读懂这封信,才能判断它是否值得被拆开、被归档。

Q=−10×log⁡10P\displaystyle Q = -10 \times \log_{10} P
通俗解释:质量分数 Q 就像是邮戳上的星级评级。Q30 代表这封信的这个字母有 99.9% 的概率没被读错,只有 0.1% 的可能性出错。
  1. QPhred 质量分数,衡量单个碱基可信度的数值。
    Q\displaystyle Q
  2. -10负号将概率取负对数后的值转为正数,10 是为了让分数落在常用量级。
    −10\displaystyle -10
  3. log10 P以 10 为底对错误概率 P 取对数;错误概率小,Q 值就大。
    log⁡10P\displaystyle \log_{10} P

2. FASTQ质控

质检员收到密信后不会直接归档,而是先抽样检查:信封有没有破损、地址是否模糊、邮票是否贴歪。测序数据的“身份证检查”同样由 FastQC 这类工具完成,它会把上百万封信的邮戳质量汇总成一份九宫格体检报告。

报告中最常看的两张图是 Per base sequence quality(每个碱基位置的质量分布)和 Per base sequence content(每个位置的碱基比例)。如果前端质量明显下滑,说明测序仪前几个循环不稳定;如果序列末端出现大量同一个 adapter 片段,说明“信封贴了多余的邮票”——发生了接头污染。

对于质量差的读段,我们需要用 Trimmomatic 或 fastp 进行修剪:去掉 3' 端低质量尾巴、去除接头序列、丢弃整条过短或过差的 read。这个动作就像把模糊不清的信件退回或裁掉破损的边角。质量阈值常用 Q30,因为 Q30 对应的错误概率只有千分之一。

只有通过初筛的信件,才允许被送去分拣大厅。

P=10−Q10=10−3010=0.001\displaystyle P = 10^{-\frac{Q}{10}} = 10^{-\frac{30}{10}} = 0.001
通俗解释:Q30 的碱基平均错误率是 0.1%,也就是 1000 个字母里预计只有 1 个被读错。
  1. P单碱基错误概率,越小说明碱基越可信。
    P\displaystyle P
  2. 10^{-Q/10}这是 Phred 公式的反解,把 Q 值换算成错误概率。
    10−Q10\displaystyle 10^{-\frac{Q}{10}}
  3. = 0.001代入 Q=30 后的计算结果,等于千分之一。

3. 比对与BAM

合格的密信进入分拣大厅。分拣员人手一份“参考地图”(参考基因组),任务是把每一封信的正文定位到地图的具体街道和门牌号——这个动作称为比对(Alignment)。归档时,管理员会把同一街道的信件按门牌号排好,并压缩保存,这就是 BAM 文件。

BAM 的文本前辈叫 SAM。SAM 每一行记录一条读段的比对结果,包括标志位(FLAG)、参考坐标、比对质量(MAPQ)、CIGAR 字符串等。CIGAR 描述了这封信在到达过程中哪些字母发生了插入、删除或错配。BAM 是 SAM 的二进制压缩版,占用空间更小、程序读取更快,几乎成了实际分析中的标准存储格式。

常用的比对工具是 BWA,它先为参考地图建立索引,再通过 seeding 和 extension 策略快速找到每个 read 的最佳落点。比对完成后,还需要按坐标排序(samtools sort)、标记并去除 PCR 重复(markdup)、建立索引(index),就像档案室给每个抽屉编号、剔除复印副本。

归档柜不是装完就完事,还要检查是否有信件摆错了格子。

MAPQ=−10×log⁡10P(read 被比对到错误位置)\displaystyle MAPQ = -10 \times \log_{10} P(\text{read 被比对到错误位置})
通俗解释:MAPQ 是“分拣员的自信程度”。分数低说明这封信可能被放到了错误的街道,下游变异检测时不应过度信任它。
  1. MAPQ比对质量值,描述一个 read 落点正确的概率。
    MAPQ\displaystyle MAPQ
  2. -10 × log10与 Phred 相同的负对数转换,让概率变成易读的分数。
    −10×log⁡10\displaystyle -10 \times \log_{10}
  3. P(被比对到错误位置)该 read 被放在参考基因组错误位置的概率。
    P(read 被比对到错误位置)\displaystyle P(\text{read 被比对到错误位置})

4. BAM质控

档案管理员不会立刻宣布归档完成,而是会抽查几个数据:某些街道是不是几乎没信(低覆盖度)?某个抽屉里是不是堆满了复印副本(PCR 重复)?信封长度是否落在合理范围(插入片段分布)?这些检查统称为 BAM 质控。

第一项看的是比对率。对全基因组测序,通常要求 95% 以上的 read 能比对到参考基因组;比对率过低说明样本可能被污染或参考基因组不匹配。第二项看测序深度与覆盖度,即基因组上每个位置被多少封信覆盖。第三项看插入片段分布:双端测序时两个 read 之间的 DNA 片段长度应呈正态分布,若出现异常尖峰,提示文库构建异常。第四项看重复率:由 PCR 扩增引入的重复 read 会被标记,比例过高会降低有效深度。

覆盖度能用一条简洁公式来估算:已知读段数、平均读长和基因组大小,就能算出期望深度。samtools flagstat 可以统计比对率,bedtools 可以计算每个位点的覆盖深度,samtools markdup 可以标记重复。

档案室的整理质量,直接决定了之后侦探能否准确找到线索。

Coverage=N×LG\displaystyle Coverage = \frac{N \times L}{G}
通俗解释:覆盖度是“平均每封信被读了几遍”。30× 表示基因组每个位置大约被 30 封不同的信覆盖。
  1. Coverage测序深度/覆盖度,衡量目标区域被重复测量的次数。
    Coverage\displaystyle Coverage
  2. N × LN 是可比对的 read 总数,L 是平均读长,两者相乘得到总测序数据量。
    N×L\displaystyle N \times L
  3. ÷ G除以参考基因组长度 G,得到单位长度的读段数。
    G\displaystyle G

5. 变异与VCF

现在轮到侦探出场。侦探会拿着每个样本的“归档信件”与参考地图逐字比对,把不同的字母圈出来,记录到一张标准问题清单上,这就是 VCF(Variant Call Format)。VCF 里的每一行,代表基因组上的一个变异位点。

VCF 是制表符分隔的文本文件,关键列有:CHROM(染色体)、POS(位置)、ID(rs 号)、REF(参考碱基)、ALT(替代碱基)、QUAL(质量值)、FILTER(过滤标记)和 INFO(附加信息)。变异类型主要包括 SNP(单个碱基替换)和 Indel(插入/缺失)。有些 VCF 还带 FORMAT 和样本列,记录每个样本的基因型(GT)、等位基因深度(AD)、测序深度(DP)、基因型质量(GQ)等。

调用变异的神器是 GATK HaplotypeCaller,它会在每个区域重新组装候选单倍型,再输出候选位点。此时每个位点的 QUAL 值代表“该位点不是纯合参考(存在变异)的证据强度”。

问题清单记录得越全,后续分析就越省力。

QUAL=−10×log⁡10P(该位点无变异)\displaystyle QUAL = -10 \times \log_{10} P(\text{该位点无变异})
通俗解释:QUAL 是侦探在“问题清单”上给每条记录打的星级。星级越高,代表这里真的存在差异的证据越强。
  1. QUAL变异位点的质量分数。
    QUAL\displaystyle QUAL
  2. -10 × log10与 Phred 一致的负对数转换。
    −10×log⁡10\displaystyle -10 \times \log_{10}
  3. P(该位点无变异)如果这个概率越小,QUAL 就越高,越支持存在变异。
    P(该位点无变异)\displaystyle P(\text{该位点无变异})

6. 变异过滤流程

侦探圈出的候选变异不能直接写入终版报告——因为其中混着不少“笔误”和“噪声”。经验丰富的复核员会逐条把关,这就是 GATK 最佳实践 中的变异过滤流程。

首先用 HaplotypeCaller 得到原始 VCF,然后根据多个统计指标做硬过滤:QD(QUAL 除以深度)代表每个证据的独立支持度;FS(Fisher 链偏置)衡量正负链上的变异是否均衡;MQ 是覆盖该位点的 read 的比对质量。常见的硬过滤条件如 QD < 2.0 或 FS > 60.0 会直接判为可疑。对于人群样本,还可以使用 VQSR(Variant Quality Score Recalibration)——用高斯混合模型学习已知位点的特征分布,给每个变异打一个更精准的置信分数。

过滤后还要用已知的 dbSNP、ClinVar 位点做交叉验证,看是否找回已知致病位点;或者用 Sanger 测序对关键位点做实验复核。整个流程像给问题清单增加“复核人签名”,只有通过层层把关的变异,才最终写入交付报告。

只有通过层层把关的变异,才配得上“可靠”二字。

QD=QUALDP\displaystyle QD = \frac{QUAL}{DP}
通俗解释:QD 是“每个证明的平均含金量”。同样 100 分证据,被 10 条 read 支持比被 100 条 read 支持更经不起推敲。
  1. QD质量值与深度的比值,衡量每个单位深度的变异证据强度。
    QD\displaystyle QD
  2. QUAL位点的原始变异质量分数。
    QUAL\displaystyle QUAL
  3. DP该位点覆盖的总 read 深度。
    DP\displaystyle DP

AI 应用场景

肿瘤靶向测序变异检测

阶段:数据预处理 → 变异检测 → 变异过滤

解决的问题:肿瘤组织测序会产生大量低质量数据和 PCR 重复,传统手动处理耗时长且阈值不统一。AI 辅助流水线能自动完成 FASTQ 质控、BAM 比对质控和 VCF 过滤,将人工调整阈值的时间从数天压缩到数小时。

具体案例:某肿瘤靶向测序项目中,分析团队构建了从 FASTQ 到 VCF 的全自动流水线:fastp 去接头 → BWA 比对 → GATK HaplotypeCaller → 硬过滤 → 注释。流水线每天处理 50 个样本,变异检出的一致率从 90% 提升到 98%。

遗传病家系外显子组分析

阶段:数据质控 → 比对 → 变异注释与过滤

解决的问题:家系外显子组测序需要从大量候选变异中筛选致病突变。通过严格的 BAM 质控和 VQSR 校准,可显著降低假阳性,再用 ANNOVAR 或 VEP 注释变异,最后结合家系共分离模式锁定可疑位点。

具体案例:一例罕见遗传病家系中,分析流程先排除低质量变异(FS > 60、QD < 2),再保留罕见错义突变,最终在候选基因中发现了与表型共分离的致病位点,辅助临床确诊。

Python 演示

python37 行
import io  def parse_fastq(text):    lines = text.strip().split('\n')    for i in range(0, len(lines), 4):        yield lines[i], lines[i + 1], lines[i + 2], lines[i + 3]  def phred_to_prob(qchar):    q = ord(qchar) - 33    return 10 ** (-q / 10.0)  def mean_quality(quals):    qs = [ord(c) - 33 for c in quals]    return sum(qs) / float(len(qs))  def filter_fastq(text, min_mean_q=20):    kept = []    for name, seq, plus, quals in parse_fastq(text):        if mean_quality(quals) >= min_mean_q:            kept.append((name, seq, plus, quals))    return kept  if __name__ == '__main__':    sample = (        "@read1\nACGTACGT\n+\nIIIIIIII\n"        "@read2\nTTTTCCCC\n+\n####!!!!\n"        "@read3\nAACCGGTT\n+\nHHHHHHHH\n"    )    kept = filter_fastq(sample, min_mean_q=20)    print("保留读段数:", len(kept))    for name, seq, plus, quals in kept:        print(name, seq)
目的:演示一个最简 FASTQ 质控过滤流程——解析四行格式、计算平均 Phred 质量、丢弃低质量读段。
  1. def parse_fastq(text)按四行一组解析 FASTQ 文本,生成器逐个产出读段。
  2. ord(qchar) - 33将 ASCII 质量字符转成 Phred 质量分 Q。
  3. mean_quality(quals)计算该读段所有碱基的平均质量分。
  4. if mean_quality(quals) >= min_mean_q保留平均质量达到阈值的读段,滤掉低质量“密信”。
  5. __name__ == '__main__'让脚本既能被导入也能直接运行。

常见误区

误区Q30 是不是代表这个碱基有 30% 的概率是正确的?所以 Q30 的数据很差,Q100 才叫好。

正解Q30 代表该碱基的错误概率只有 0.1%,也就是准确率 99.9%。Phred 分数与错误概率是对数关系,Q 值每升高 10,错误率就下降 10 倍。Q20、Q30 是常见的合格线,不是“30% 正确”。

为何常见:人们容易把“数值小 = 质量差”的直觉迁移过来,忽略了 Phred 分数是负对数转换后的结果。

误区BAM 和 SAM 是两种不同的文件格式,SAM 是文本,BAM 是二进制,内容完全不同,不能互相转换。

正解BAM 本质上是 SAM 的二进制压缩版,记录的信息完全一致,只是压缩后体积更小、读取更快。两者可以通过 samtools view 自由转换,且不会丢失信息。分析中为了节省磁盘空间,几乎统一使用 BAM。

为何常见:用户看到文件名后缀不同、打开方式不同,容易误以为是两种独立格式。

误区QUAL 分数高的变异,临床上一定更致命或多态效应更强。

正解QUAL 只代表“这个位点真实存在变异的统计证据强度”,与变异对蛋白质功能或疾病的影响大小毫无关系。一个 QUAL=99 的同义突变可能无害,而 QUAL=50 的无义突变可能致病。

为何常见:把“质量值”和“严重性/致病性”两个概念混淆了,临床解读需要参考注释和 ACMG 指南。

误区只要把 `QD < 2.0`、`FS > 60.0` 这些阈值套到任何 VCF 上,就能得到可靠的变异集。

正解GATK 官方阈值来自人类 WGS 外显子组的最佳实践;靶向测序、RNA-seq、细菌基因组、不同测序深度等场景的统计分布差异很大,硬编码阈值会导致过多假阳或假阴。应根据数据特征调整,或用 VQSR 学习经验分布。

为何常见:官方文档给出了“推荐值”,很多人当成“标准值”直接使用,忽略了统计阈值的条件性。

面试题

Q1. FASTQ 文件中的第二行记录的是什么?

单选题 · easy

  • A.碱基序列← 正确答案
  • B.质量分数
  • C.序列标识
  • D.染色体位置

答案:A。FASTQ 四行结构依次是标识符、碱基序列、+分隔符、质量字符串。第二行是测得的 DNA 序列;质量分数在第四行,标识符在第一行,位置信息不属于 FASTQ。

Q2. Phred 质量分数 Q=30 对应的单碱基错误概率约为?

单选题 · easy

  • A.0.1%← 正确答案
  • B.1%
  • C.10%
  • D.0.01%

答案:A。根据 P=10^(-Q/10),Q=30 时 P=10^(-3)=0.001,即 0.1%。这相当于每 1000 个碱基中预计有 1 个可能读错;1% 对应 Q20,0.01% 对应 Q40。

Q3. 关于 SAM/BAM 文件,下列说法正确的有?

多选题 · easy

  • A.压缩二进制的SAM← 正确答案
  • B.与SAM文本内容相同← 正确答案
  • C.可互相无损转换← 正确答案
  • D.只能存比对标志

答案:ABC。BAM 是 SAM 的二进制压缩版本,二者记录的信息一致,可通过 samtools view 互相转换。D 错误,BAM 存储完整比对信息,包括序列、质量、CIGAR、FLAG 等。

Q4. FastQC 报告中的 "Per base sequence quality" 图,横轴通常代表什么?

单选题 · mid

  • A.读段的位置← 正确答案
  • B.染色体坐标
  • C.碱基质量
  • D.测序平台

答案:A。该图展示每条 read 从 5' 到 3' 每个位置的质量分布,横轴是读段上的碱基位置。染色体坐标属于比对后概念,FastQC 处理的是原始 FASTQ 数据。

Q5. samtools sort 命令的主要作用是什么?

单选题 · mid

  • A.按坐标排序← 正确答案
  • B.按质量排序
  • C.按名称排序
  • D.去除重复

答案:A。比对后的 SAM/BAM 默认按 read 名称顺序排列,下游变异检测需要按参考坐标顺序读取,所以必须先按坐标排序。质量排序无此用途,去重由 markdup 完成。

Q6. 以下属于 BAM 质控关键指标的有?

多选题 · mid

  • A.整体比对百分比← 正确答案
  • B.有效覆盖深度← 正确答案
  • C.插入片段长度分布← 正确答案
  • D.参考碱基构成比例

答案:ABC。BAM 质控关注比对率、测序深度、覆盖度、插入片段分布和重复率。参考碱基构成比例由参考基因组决定,与 BAM 数据质量无关,不属于常规质控指标。

Q7. VCF 文件中 ALT 列记录的是?

单选题 · mid

  • A.替代型等位基因← 正确答案
  • B.参考型等位基因
  • C.基因功能注释信息
  • D.样本基因型信息

答案:A。VCF 的核心列中 REF 是参考基因组碱基,ALT 是检出的替代等位基因。基因注释信息由 ANNOVAR/VEP 等工具补充,样本基因型在 FORMAT/样本列中。

Q8. 关于 VCF 中 QUAL 值的正确理解是?

单选题 · hard

  • A.通常QUAL越高变异证据越强← 正确答案
  • B.通常QUAL越高变异耐性越差
  • C.QUAL与测序深度完全无关
  • D.QUAL等于所有碱基质量数值之和

答案:A。QUAL 是 -10 log10 P(该位点无变异),越高代表真实变异的证据越强。它并不代表临床严重性;虽然深度会影响 QUAL 的计算,但并不能说完全无关,也不是所有碱基质量之和。

Q9. GATK 硬过滤时常用的统计指标包括?

多选题 · hard

  • A.质量比QD← 正确答案
  • B.链偏置FS← 正确答案
  • C.比对质量MQ← 正确答案
  • D.碱基GC含量

答案:ABC。GATK 硬过滤经典条件围绕 QD、FS、MQ、MQRankSum、ReadPosRankSum 等统计量;GC 含量是 FASTA/FASTQ 质控或文库偏好分析指标,不属于 GATK 过滤因子。

Q10. VQSR 变异质量重校准主要依赖的统计方法是?

单选题 · extreme

  • A.高斯混合模型← 正确答案
  • B.线性回归模型
  • C.决策树模型
  • D.支持向量机模型

答案:A。VQSR 利用已知的 1000 Genomes、HapMap、dbSNP 等位点作为训练真集,用高斯混合模型刻画高质量与低质量变异的特征分布,从而输出每个位点的 VQSLOD 分数。线性回归、决策树和 SVM 均不是 GATK 官方采用的方法。