MSBS7101 Lecture 6 · 变异检测与注释（一）：LoFreq 与 GATK

本讲介绍变异检测（variant calling）：遗传变异三大类型、区分真实变异与测序错误的四项标准、VCF 格式，并实操 LoFreq 与 GATK Mutect2 两条流程，从全外显子组 BAM 检出变异并在 IGV 中人工核查。

一、本讲信息

  • 课程 ：MSBS7101 Essential Skills in High-throughput Sequencing Data Analysis（高通量测序数据分析核心技能，课件首页标注课程码为 SBMS7201/MSBS7101）
  • 讲次 ：Variant Calling and Annotation 1（变异检测与注释·第一讲）
  • 日期 ：2026-10-06
  • 材料 ：讲义 Variant_calling_2026.pdf（39 页）+ Colab 实操 notebook Variant_Calling.ipynb
  • 前置知识 ：上一讲完成测序数据质控（FastQC）、比对（BWA）、比对文件处理（samtools）与可视化（IGV）

二、学习目标

课件第 2 页列出三条目标（课堂录音补充：Mutect2 流程步骤繁多，不要求逐步背诵，但要理解每一步为什么存在 ）：
  1. 描述 VCF（Variant Call Format）格式的特点；
  1. 理解变异检测流程（variant calling pipeline）的目的与各步骤的作用；
  1. 使用生物信息学工具实际生成变异检测结果。

三、核心内容

3.1 为什么要研究遗传变异

遗传变异（genetic variation） 指个体或细胞之间 DNA 序列的差异。它决定了生物体绝大多数可观察的性状（phenotype），也构成大多数人类疾病与健康问题的基础；环境因素与基因的共同作用同样是表型的重要来源。课堂录音补充：在下一代测序（Next-Generation Sequencing, NGS）的众多用途中（后续课程还会讲 RNA 测序等），检测遗传变异是最主要的用途之一——临床上对患者 DNA 进行测序，目的往往是获得疾病诊断的线索。

3.2 遗传变异的三大类型

遗传变异的三大类型
遗传变异的三大类型
图 1（课件 p5）：遗传变异按尺度分为拷贝数变异、结构变异与序列水平变异三类，右图为染色体示意。
  1. 拷贝数变异（Copy Number Variation, CNV） ：整条染色体或染色体臂的增益/丢失。典型例子是唐氏综合征（Down syndrome），患者有三条 21 号染色体（trisomy 21）。NGS 也能检测这类变异，但相关算法不在本课程范围内。
  1. 结构变异（Structural Variation, SV） ：涉及较大片段的重排，包括缺失（deletion）、重复（duplication）、倒位（inversion）与易位（translocation）。课件例子为面肩肱型肌营养不良症（Facioscapulohumeral Muscular Dystrophy, FSHD）。课堂录音详细解释：该病涉及 4 号染色体上一段串联重复（正常约 11–100 个拷贝），拷贝数减少时导致 DUX4 基因异常表达，其产物对肌肉细胞有害，患者在青少年期后肌肉逐渐受损；这类变异规模约 30–3000 kb，由于短读长（约 150 bp）无法数清重复拷贝数，几乎无法用常规 NGS 检测。
  1. 序列水平变异（Sequence level variation） ：本课程的重点，分为两类——
      • 单核苷酸多态性（Single Nucleotide Polymorphism, SNP） ：单个碱基的替换。人类是二倍体（每个位点有两条 DNA 拷贝），若仅一条拷贝发生改变称为杂合（heterozygous） ，两条都改变称为纯合（homozygous） 。
      • 插入/缺失（insertions/deletions, indels） ：少数碱基的插入或缺失。录音说明：indel 与结构变异之间没有严格界线，通常小于 100 bp（多数 50 bp 以内）称为 indel，达到数百碱基至千碱基则归为结构变异。
SNP 与 indel 图示
SNP 与 indel 图示
图 2（课件 p8）：以参考基因组（reference genome）对照二倍体基因组序列，展示杂合/纯合 SNP、缺失（CTT）与插入（CG）的序列形式，变异碱基以红色标出。
课件以镰刀型细胞贫血（sickle cell anaemia） 作为 SNP 致病的经典例子：β-珠蛋白基因上一个碱基替换（T→A）改变密码子，使 mRNA 上编码谷氨酸的密码子变为编码缬氨酸，血红蛋白结构异常，红细胞呈镰刀状、携氧能力下降。课堂录音补充：该等位基因在疟疾流行地区具有选择优势——寄生虫难以在镰状红细胞中生存，这解释了其在特定人群中的分布。

3.3 变异检测要解决的两个核心问题

课件第 10 页给出一个具体场景：一条读段（read）比对到人类参考基因组 hg38 的 2 号染色体 MSH6 基因区域后，与参考序列存在不匹配的碱基。由此引出两个问题：
  1. 如何区分测序错误与真实突变——在整个 BAM 文件中逐个人工判断不可行，需要高效的方法；
  1. 变异有什么功能影响——这部分留到下一讲的变异注释（variant annotation）。
本讲解决第一个问题。

3.4 真实变异的四个判断特征

变异检测算法区分测序错误（sequencing error）与真实变异，依据四类信息（课件 p11）：
真实变异的判断特征
真实变异的判断特征
图 3（课件 p11）：比对截图中同一区域的多个候选变异，有的只出现在单条读段、有的出现在多条读段、有的读段上密集出现多个变异，用以说明四项判断标准。
  1. 碱基质量分数（Phred/Q score） ：测序仪对每个碱基给出的质量值，Q 越高该碱基被读错的概率越低（上讲内容）；
  1. 支持变异的读段数量与正反链（strand） ：变异出现在越多读段、且同时出现在正链和负链读段上，越可信；若只出现在单条读段或全部同一方向，则可疑；
  1. 单条读段中变异的数量 ：一条读段上同时出现多个变异会增加可疑程度；
  1. 变异在读段中的位置 ：读段两端碱基质量通常偏低，变异若总落在读段末端，可信度下降。
课堂录音强调：没有任何单一标准能单独定论 ，算法会综合全部信息打分；即便如此，变异检测工具的 PASS 判定也不能完全相信——老师展示了一个变异正负链读段支持为 19:1 却被判 PASS 的例子。实际研究甚至临床生物信息学中，发现重要变异后必须在 IGV 等浏览器中人工核查，这是导师和临床平台都会要求的基本步骤。

3.5 基本变异检测流程

基本变异检测流程
基本变异检测流程
图 4（课件 p12）：从测序读段到注释结果的完整流程——质控（FastQC）、比对（BWA）、比对清理（samtools）、变异检测（LoFreq）与功能注释（VEP），各步之间以 FASTQ/SAM(BAM)/BAM/VCF 数据格式衔接。
本讲在上讲流程（FastQC → BWA → samtools → IGV）之后新增两步：LoFreq 推断 SNV 与 indel（本讲），VEP 判断功能意义（下一讲）。SNV（single nucleotide variant）指单核苷酸变异。

3.6 实操数据：全外显子组测序

本讲在 Google Colab 中运行 notebook：挂载 Google Drive、用 condacolab 安装 lofreq、在独立环境安装 GATK 4.6.2.0 与 openjdk 17，并安装 samtools/bcftools/htslib/bedtools/bamtools/igv-notebook。每次重新连接 Colab 都要重做这些设置 。
数据文件（Materials/Variant_Calling.ipynb）：
  • WXS_example_sorted.bam：来自（模拟的）同一结直肠癌患者血液 的全外显子组测序（Whole Exome Sequencing, WES）——外显子组测序只测定外显子（编码区），因此 IGV 中只有外显子区域有覆盖；
  • WXS_example_cancer_sorted.bam：来自同一患者结肠癌组织 的 WES 数据，供课后作业中与血液结果比较使用。
生殖系与体细胞突变（germline versus somatic） 是理解这组数据的关键（课件 p17）：
Germline 与 somatic 突变
Germline 与 somatic 突变
图 5（课件 p17）：生殖系突变遗传自亲代或发生在受精前，全身细胞普遍携带；体细胞突变则在个体发育过程中后天获得，仅存在于部分细胞（尤其是癌细胞）中。
  • 生殖系突变（germline mutation） ：遗传获得或发生在受精之前，全身所有细胞都有；
  • 体细胞突变（somatic mutation） ：特定细胞后天获得，只存在于部分细胞，在肿瘤中常见——例如结肠细胞获得突变后发展为结肠癌。
研究肿瘤时，标准做法是比较肿瘤 DNA 与正常（血液）DNA 的差异。课堂录音补充：体细胞突变随年龄不断累积；流行病学研究提示约三分之二的突变来源于内源因素或 DNA 复制本身。本讲先只处理血液（germline）样本。

3.7 LoFreq：快速变异检测与 VCF 格式

LoFreq 是一个简单直接的变异检测工具（录音说明：它并非行业主流，发表文献少用，但适合教学；它也支持并行计算和 tumor/normal 的 somatic 模式）。核心命令（notebook）：
  • --no-default-filter：只检出（call）变异、不做假阳性过滤——教学时先看全部候选；lofreq filter 再补上质量标签；--print-all 在打标签的同时保留所有变异，便于对比。
  • LoFreq 在检测的同时还会对读段做局部重新比对处理，因此其报告的深度可能与 IGV 直接显示的读段数略有出入（录音）。
VCF 文件结构
VCF 文件结构
图 6（课件 p22）：LoFreq 输出的 VCF——前 16 行为以 # 开头的文件头（记录命令、参数与字段含义），其后每行一个变异，含 CHROM/POS/ID/REF/ALT/QUAL/FILTER/INFO 列。
VCF（Variant Call Format） 是存储变异检测结果的标准文本格式：
  • 文件头 ：以 # 开头，记录所用命令、参考基因组与各字段缩写的定义；
  • 数据行 （每行一个变异）的主要列：
    • CHROM/POS：染色体与位置；ID：变异编号（尚未注释时为 .）；
    • REF/ALT：参考等位基因与替代等位基因（ALT 应与参考基因组一致核对）；
    • QUAL：综合碱基质量、链偏倚、读段内变异数、变异位置等信息得出的质量分，分高者标 PASS；
    • FILTER：PASS 或未通过的原因标签——课堂演示中，覆盖深度低于 10 的位点、链偏倚达显著水平的位点被分别标记；
    • INFO：DP（该位置原始深度）、AF（等位基因频率，即支持 ALT 的读段比例，AF=1 表示 100% 读段含变异）、SB（链偏倚）、DP4（参考/替代等位基因各在正/负链的读段数）等。
IGV 核查 ：利用 VCF 中的位置信息在 IGV 中直接跳转查看（如 notebook 中的 chr2:5876477-5876517）。课堂演示：第一个变异 T→C 共 24 条读段支持（15 正链 + 9 负链），与 VCF 的 DP=24 一致，基本确定为真；另一处 G→A 只有 5 条读段且全部来自负链，同时触发链偏倚与低深度两项标记，判为不可信。

3.8 GATK：社区标准流程

GATK（Genome Analysis Toolkit） 由 Broad Institute（哈佛与 MIT 附属机构）开发，是生物信息学界公认的变异检测标准工具，广泛用于研究与临床。GATK 的流程在「原始比对 BAM」与「变异检测」之间增加了预处理步骤，目的是生成「分析就绪（analysis-ready）」的 BAM——录音形象地总结为：想得到好的变异检测结果，需要先「退回一步」优化比对文件。
主要步骤（notebook 命令与课件 p28–p36）：
1. 添加读组（AddOrReplaceReadGroups） ：GATK 要求 BAM 文件头中有读组（read group, RG）元数据（样本、文库、平台、运行批次等）。这是历史设计——GATK 最初面向多样本混合的 BAM；即使只分析单样本，也必须补上 RG 才能运行（录音：这一步不改变分析内容，属于必要的格式步骤）。
2. 标记重复序列（MarkDuplicates） ：NGS 数据中有两类重复读段（duplicates）——
  • PCR 重复 ：建库时 DNA 量不足需 PCR 扩增，扩增偏好使某些片段被过度复制；
  • 光学重复（optical duplicates） ：测序芯片（flow cell）上同一 DNA 分子簇在信号扩增后，被相机误读为两个点，产生两条「来源不同」实则相同的读段（录音解释）。
MarkDuplicates 给这些重复读段打上标记，下游分析时将其忽略/降权，等价于把多条相同读段折叠为一条，防止同一分子被重复计数而虚增变异「证据」。课堂数据中约 10% 为重复读段，对外显子数据属正常；若达 80–90% 则说明实验有问题。
3. 碱基质量重校准（Base Quality Score Recalibration, BQSR） ：测序错误率并非在所有序列上下文（context）中均匀——录音举例：GG 语境错误率可能接近零，而 TG 语境明显更高。若直接沿用测序仪给出的原始质量分，会高估或低估真实错误率。分两步：
  • BaseRecalibrator：扫描比对读段，结合已知变异位点（dbSNP138、Mills 数据集作为「真值」参考），按循环周期、序列上下文、机器特性训练一个系统误差统计模型，输出重校准表；
  • ApplyBQSR：按该表逐碱基调整质量分，写入新 BAM（MD_BR.bam）。
课堂演示其效果：某碱基质量从 33 校准为 22，标记重复后峰值覆盖从 32 条降至 26 条——单点差异看似不大，但在全基因组百万级变异的规模上可显著降低假阳性。外显子数据需用 -L intervals.bed -ip 10 把分析限制在捕获区域并向两侧各扩 10 bp（padding，因 DNA 片段末端可能略超出探针捕获区）。
4. 比对质量指标（alignment metrics） ：notebook 运行 bamtools stats、CollectReadCounts、CollectMultipleMetrics、CollectHsMetrics 生成约 8 类 QC 文件（课件 p33 的表格：比对摘要、按循环的碱基分布与质量分布、插入片段长度分布、读长直方图、靶向区域覆盖等）。录音说明：这些指标在临床实验室质控中尤为重要，例如靶向数据中 off-target 读段比例过高提示实验异常。
5. 变异检测（Mutect2） ：GATK 提供两个主要检测器——HaplotypeCaller 只用于 germline；Mutect2（MuTect2） 可同时检测 germline 与 somatic 变异（本讲为简化只用 Mutect2 处理 germline 血液样本；录音：纯 germline 场景 HaplotypeCaller 噪声更低）。关键参数：-tumor 指定样本名、--germline-resource 提供人群等位基因频率参考（gnomAD，Broad 维护的大规模人群聚合数据库，含分人群频率，用于标注哪些变异在人群中已知）、--f1r2-tar-gz 输出读段方向信息、--min-base-quality-score 20 过滤低质量碱基。
6. 概率过滤（probabilistic filtering） （课件 p35）：
  • LearnReadOrientationModel：从 F1R2 输出学习读段方向偏差模型；
  • GetPileupSummaries + CalculateContamination：利用一组常见 SNP 估计样本间的交叉污染水平；
  • FilterMutectCalls：结合上述模型与污染估计做概率化 过滤（而非固定硬阈值），输出带 FILTER 标签的 VCF。
7. 规范化与进一步过滤 （课件 p36）：
  • bcftools norm：把多等位基因记录（如同一位置 A>C 与 A>G）拆分为多条双等位基因记录，并对 indel 做左对齐与规范化——不同工具对变异坐标的约定不同，规范化后记录才可比较、可合并；
  • bcftools filter 保留满足以下条件的变异：等位基因频率（Variant Allele Frequency, VAF）AF ≥ 2%、样本深度 DP ≥ 50、人群频率 POPAF ≥ 2 或在 gnomAD 中缺失。录音特别解释：POPAF 是 −log10 尺度 ，−log10(0.01)=2，故 POPAF≥2 等价于人群频率 ≤1%——即只保留罕见或人群中未见的变异；POPAF 缺失也接受，表示该变异在数据库中不存在。
  • 最后用 bcftools annotate 精简 INFO 字段（只保留等位基因深度 AD），得到最终文件，本例检出 3133 个变异 （录音）。
8. IGV 复核结果 （课件 p38，录音详解四个案例）：
  • slippage 标签 ：一处 G→GT 插入，其下游是连续长串 T——NGS 在同聚物（homopolymer）区域难以数清碱基个数，这类插入既可能真实也可能是测序滑动错误，需人工判断；
  • germline 标签 ：一处 G→C 变异，读段支持恰好约 50:50、正负链均衡，在 gnomAD 中可查到（罕见但存在）——germline 样本每个位点只有两份拷贝，VAF 期望约 50%；实际中 70:30 以上的比例即属可接受；
  • PASS 却实为数据库中已知变异 ：该变异是致病性（pathogenic）的无义突变（stop-gain），之所以未被标为 germline，是因为 af-only-gnomad 参考库刻意剔除了可能致病的变异——若不剔除，用户盲目只保留「不在人群数据库中」的变异时就会把致病变异一并滤掉。这是参考库设计上的重要细节；
  • clustered events ：多个变异密集出现在相邻位置，Mutect2 能将其识别为聚类事件（LoFreq 未明确标注），可疑度增加。
录音的总结论：Mutect2 远比 LoFreq 强大，今后自己的研究应使用前者；但无论用什么工具，都不能只筛 PASS 就直接使用 ，必须人工检查数据。

3.9 下一讲预告

课件 p39：下一讲将讨论变异过滤决策的更多细节（请先用 IGV 浏览 PASS 与被过滤的变异），并进行变异注释（variant annotation，使用 VEP）。

四、关键图解

见上文图 1–图 6：
图
页码
内容
图 1
p5
遗传变异三大类型总览
图 2
p8
SNP（杂合/纯合）与 indel 的序列图示
图 3
p11
真实变异的四项判断标准与比对截图
图 4
p12
基本变异检测流程图（FASTQ→BAM→VCF）
图 5
p17
germline 与 somatic 突变的区别
图 6
p22
LoFreq 输出的 VCF 文件结构

五、课前思考 / 自测

（材料中未附练习题，以下依据课件内容自拟）
  1. 某变异的全部支持读段都来自负链，为什么这会显著降低其可信度？这与测序过程中的哪类系统误差有关？
  1. PCR 重复与光学重复（optical duplicate）的成因有何不同？如果不做 MarkDuplicates，变异检测的假阳性为什么会升高？
  1. BQSR 为什么需要提供 dbSNP 等已知变异位点（known sites）？外显子数据为什么要用 -L intervals.bed -ip 10 限定区域并加 10 bp padding？
  1. 过滤条件 POPAF>=2 在 −log10 尺度下对应的人群频率是多少？为什么 POPAF 缺失（"."）的变异也被保留？
  1. germline 样本中 VAF≈50% 与 VAF≈5% 的变异分别提示什么？若一个致病变异出现在 gnomAD 中，为什么 af-only-gnomad 参考库要刻意排除它？

六、参考材料

  1. 主讲课件：Materials/Variant_calling_2026.pdf（39 页，Prof. Jason Wong）
  1. 实操 notebook：Materials/Variant_Calling.ipynb（Google Colab，HKU Practical Bioinformatics 2026 - Variant Calling）
  1. 课堂录音转写：Recordings/MSBS7101-1006.txt（英文语音识别转写，噪声较多，仅提取可靠信息）

七、课堂要点补充（据课堂录音整理）

  1. 背景动机 ：有同学反映上一讲不清楚「为什么要比对」，因此本讲先补遗传变异背景——NGS 数据最广泛的用途就是检测遗传变异，临床测序旨在为疾病提供诊断依据。
  1. indel 与结构变异的界线 ：没有严格定义，惯例是 <100 bp（多为 ≤50 bp）称 indel，数百 bp 至 kb 级称结构变异；凡无需专门 SV 算法、常规 NGS 流程即可检测的可视为 indel。
  1. 不要盲信工具的 PASS ：LoFreq 判 PASS 的一个变异正负链支持为 19:1，人工核查明显可疑；即使最好的变异检测工具也会出错，发现重要结果必须在 IGV 中人工检查——这是导师与临床平台的普遍要求。
  1. LoFreq 的定位 ：简单、适合教学，但非行业主流；其 somatic 模式可在作业中用于肿瘤/正常成对数据。
  1. 重校准的现实意义 ：单碱基质量分的调整（如 33→22）看似微小，但在全基因组数百万变异的规模上能显著降低假阳性率——这是 GATK 增加预处理步骤的意义。
  1. QC 指标的用途 ：插入片段长度分布、靶向覆盖、off-target 比例等指标在临床实验室是必查项；重复读段比例约 10% 对外显子数据正常，80–90% 则提示实验失败。
  1. gnomAD ：Broad Institute 维护的人群聚合等位基因频率数据库（汇总大规模人群数据、按人群分层），既用于 Mutect2 的 germline 资源，也可人工查询变异是否已知。查询时注意使用与参考基因组匹配的版本（如 hg38）。
  1. 致病变异与参考库设计 ：af-only-gnomad 剔除可能致病的变异，正是为了防止用户「只保留人群未见过变异」时把致病变异错误过滤掉——本讲数据中即出现这样一例 stop-gain 的 PASS 变异。
  1. WES 数据的覆盖形态 ：外显子覆盖可超过 200×，内含子区极低，且总存在少量 off-target 读段，属正常现象。
  1. 作业提示 ：肿瘤样本 BAM（WXS_example_cancer_sorted.bam）将在作业中与血液样本结果比较使用。
上一篇
MSBS7102 Lecture 6 · ggplot2 数据可视化
下一篇
MSBS7102 Lecture 5 · R 语言入门与 Tidyverse
Loading...
2026-11-3