想象你是一座超级邮局的质检主任。每天,测序仪这台“基因邮筒”会吐出上亿封基因密信——每一封信都有四个部分:以 @ 开头的信封编号,由 A、T、C、G 写成的正文,一行 + 分隔线,以及一串像乱码一样的“邮戳质量符”。这些密信必须经过质检员抽查(FASTQ 质控)、分拣员按参考地图归档(比对生成 BAM)、侦探逐字核对找出拼写差异(变异检测),最后整理成一份标准问题清单(VCF)。本文会带你以这条“信使流水线”为线索,读懂 FASTQ、BAM、VCF 三种数据格式,并掌握从原始测序数据到可靠变异检测的完整质控流程。
掌握这套标准格式与质控流程,是进入生信分析与精准医疗领域的入场券。学会它,你就从“只会调用工具”升级为“能诊断数据质量、设计可靠流水线”的数据诊断专家。
可视化演示
展示从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,正是测序仪出厂质检的通行证。
读懂这封信,才能判断它是否值得被拆开、被归档。
- QPhred 质量分数,衡量单个碱基可信度的数值。
- -10负号将概率取负对数后的值转为正数,10 是为了让分数落在常用量级。
- log10 P以 10 为底对错误概率 P 取对数;错误概率小,Q 值就大。
2. FASTQ质控
质检员收到密信后不会直接归档,而是先抽样检查:信封有没有破损、地址是否模糊、邮票是否贴歪。测序数据的“身份证检查”同样由 FastQC 这类工具完成,它会把上百万封信的邮戳质量汇总成一份九宫格体检报告。
报告中最常看的两张图是 Per base sequence quality(每个碱基位置的质量分布)和 Per base sequence content(每个位置的碱基比例)。如果前端质量明显下滑,说明测序仪前几个循环不稳定;如果序列末端出现大量同一个 adapter 片段,说明“信封贴了多余的邮票”——发生了接头污染。
对于质量差的读段,我们需要用 Trimmomatic 或 fastp 进行修剪:去掉 3' 端低质量尾巴、去除接头序列、丢弃整条过短或过差的 read。这个动作就像把模糊不清的信件退回或裁掉破损的边角。质量阈值常用 Q30,因为 Q30 对应的错误概率只有千分之一。
只有通过初筛的信件,才允许被送去分拣大厅。
- P单碱基错误概率,越小说明碱基越可信。
- 10^{-Q/10}这是 Phred 公式的反解,把 Q 值换算成错误概率。
- = 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比对质量值,描述一个 read 落点正确的概率。
- -10 × log10与 Phred 相同的负对数转换,让概率变成易读的分数。
- P(被比对到错误位置)该 read 被放在参考基因组错误位置的概率。
4. BAM质控
档案管理员不会立刻宣布归档完成,而是会抽查几个数据:某些街道是不是几乎没信(低覆盖度)?某个抽屉里是不是堆满了复印副本(PCR 重复)?信封长度是否落在合理范围(插入片段分布)?这些检查统称为 BAM 质控。
第一项看的是比对率。对全基因组测序,通常要求 95% 以上的 read 能比对到参考基因组;比对率过低说明样本可能被污染或参考基因组不匹配。第二项看测序深度与覆盖度,即基因组上每个位置被多少封信覆盖。第三项看插入片段分布:双端测序时两个 read 之间的 DNA 片段长度应呈正态分布,若出现异常尖峰,提示文库构建异常。第四项看重复率:由 PCR 扩增引入的重复 read 会被标记,比例过高会降低有效深度。
覆盖度能用一条简洁公式来估算:已知读段数、平均读长和基因组大小,就能算出期望深度。samtools flagstat 可以统计比对率,bedtools 可以计算每个位点的覆盖深度,samtools markdup 可以标记重复。
档案室的整理质量,直接决定了之后侦探能否准确找到线索。
- Coverage测序深度/覆盖度,衡量目标区域被重复测量的次数。
- N × LN 是可比对的 read 总数,L 是平均读长,两者相乘得到总测序数据量。
- ÷ G除以参考基因组长度 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 × log10与 Phred 一致的负对数转换。
- P(该位点无变异)如果这个概率越小,QUAL 就越高,越支持存在变异。
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质量值与深度的比值,衡量每个单位深度的变异证据强度。
- QUAL位点的原始变异质量分数。
- DP该位点覆盖的总 read 深度。
AI 应用场景
肿瘤靶向测序变异检测
阶段:数据预处理 → 变异检测 → 变异过滤
解决的问题:肿瘤组织测序会产生大量低质量数据和 PCR 重复,传统手动处理耗时长且阈值不统一。AI 辅助流水线能自动完成 FASTQ 质控、BAM 比对质控和 VCF 过滤,将人工调整阈值的时间从数天压缩到数小时。
具体案例:某肿瘤靶向测序项目中,分析团队构建了从 FASTQ 到 VCF 的全自动流水线:fastp 去接头 → BWA 比对 → GATK HaplotypeCaller → 硬过滤 → 注释。流水线每天处理 50 个样本,变异检出的一致率从 90% 提升到 98%。
遗传病家系外显子组分析
阶段:数据质控 → 比对 → 变异注释与过滤
解决的问题:家系外显子组测序需要从大量候选变异中筛选致病突变。通过严格的 BAM 质控和 VQSR 校准,可显著降低假阳性,再用 ANNOVAR 或 VEP 注释变异,最后结合家系共分离模式锁定可疑位点。
具体案例:一例罕见遗传病家系中,分析流程先排除低质量变异(FS > 60、QD < 2),再保留罕见错义突变,最终在候选基因中发现了与表型共分离的致病位点,辅助临床确诊。
Python 演示
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)
def parse_fastq(text)按四行一组解析 FASTQ 文本,生成器逐个产出读段。ord(qchar) - 33将 ASCII 质量字符转成 Phred 质量分Q。mean_quality(quals)计算该读段所有碱基的平均质量分。if mean_quality(quals) >= min_mean_q保留平均质量达到阈值的读段,滤掉低质量“密信”。__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 文件中的第二行记录的是什么?
- A.碱基序列← 正确答案
- B.质量分数
- C.序列标识
- D.染色体位置
答案:A。FASTQ 四行结构依次是标识符、碱基序列、+分隔符、质量字符串。第二行是测得的 DNA 序列;质量分数在第四行,标识符在第一行,位置信息不属于 FASTQ。
Q2. Phred 质量分数 Q=30 对应的单碱基错误概率约为?
- 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 文件,下列说法正确的有?
- A.压缩二进制的SAM← 正确答案
- B.与SAM文本内容相同← 正确答案
- C.可互相无损转换← 正确答案
- D.只能存比对标志
答案:ABC。BAM 是 SAM 的二进制压缩版本,二者记录的信息一致,可通过 samtools view 互相转换。D 错误,BAM 存储完整比对信息,包括序列、质量、CIGAR、FLAG 等。
Q4. FastQC 报告中的 "Per base sequence quality" 图,横轴通常代表什么?
- A.读段的位置← 正确答案
- B.染色体坐标
- C.碱基质量
- D.测序平台
答案:A。该图展示每条 read 从 5' 到 3' 每个位置的质量分布,横轴是读段上的碱基位置。染色体坐标属于比对后概念,FastQC 处理的是原始 FASTQ 数据。
Q5. samtools sort 命令的主要作用是什么?
- A.按坐标排序← 正确答案
- B.按质量排序
- C.按名称排序
- D.去除重复
答案:A。比对后的 SAM/BAM 默认按 read 名称顺序排列,下游变异检测需要按参考坐标顺序读取,所以必须先按坐标排序。质量排序无此用途,去重由 markdup 完成。
Q6. 以下属于 BAM 质控关键指标的有?
- A.整体比对百分比← 正确答案
- B.有效覆盖深度← 正确答案
- C.插入片段长度分布← 正确答案
- D.参考碱基构成比例
答案:ABC。BAM 质控关注比对率、测序深度、覆盖度、插入片段分布和重复率。参考碱基构成比例由参考基因组决定,与 BAM 数据质量无关,不属于常规质控指标。
Q7. VCF 文件中 ALT 列记录的是?
- A.替代型等位基因← 正确答案
- B.参考型等位基因
- C.基因功能注释信息
- D.样本基因型信息
答案:A。VCF 的核心列中 REF 是参考基因组碱基,ALT 是检出的替代等位基因。基因注释信息由 ANNOVAR/VEP 等工具补充,样本基因型在 FORMAT/样本列中。
Q8. 关于 VCF 中 QUAL 值的正确理解是?
- A.通常QUAL越高变异证据越强← 正确答案
- B.通常QUAL越高变异耐性越差
- C.QUAL与测序深度完全无关
- D.QUAL等于所有碱基质量数值之和
答案:A。QUAL 是 -10 log10 P(该位点无变异),越高代表真实变异的证据越强。它并不代表临床严重性;虽然深度会影响 QUAL 的计算,但并不能说完全无关,也不是所有碱基质量之和。
Q9. GATK 硬过滤时常用的统计指标包括?
- A.质量比QD← 正确答案
- B.链偏置FS← 正确答案
- C.比对质量MQ← 正确答案
- D.碱基GC含量
答案:ABC。GATK 硬过滤经典条件围绕 QD、FS、MQ、MQRankSum、ReadPosRankSum 等统计量;GC 含量是 FASTA/FASTQ 质控或文库偏好分析指标,不属于 GATK 过滤因子。
Q10. VQSR 变异质量重校准主要依赖的统计方法是?
- A.高斯混合模型← 正确答案
- B.线性回归模型
- C.决策树模型
- D.支持向量机模型
答案:A。VQSR 利用已知的 1000 Genomes、HapMap、dbSNP 等位点作为训练真集,用高斯混合模型刻画高质量与低质量变异的特征分布,从而输出每个位点的 VQSLOD 分数。线性回归、决策树和 SVM 均不是 GATK 官方采用的方法。