MSBS7101 Lecture 7 · 变异检测与注释（二）：bcftools 与 VEP

本讲承接变异检测，用 bcftools 过滤与比较 VCF、bedtools 与 UCSC Table Browser 筛选编码区变异，再用 VEP 与 AlphaMissense 注释并优先化致病变异。

1. 本讲信息

  • 课程 ：MSBS7101 Essential Skills in High-throughput Sequencing Data Analysis（课件页脚亦标注 SBMS7201/MSBS7101）
  • 讲次 ：Variant Calling and Annotation 2（承接上一讲的变异检测）
  • 日期 ：2026-10-20
  • 教师 ：Prof. Jason Wong（School of Biomedical Sciences, HKU, jwhwong@hku.hk）
  • 材料 ：讲课 PDF《Variant_annotation_2026》（41 页）+ Google Colab 实操 notebook《Variant_Annotation.ipynb》
  • 定位 ：本讲是「找变异」到「理解变异」的过渡——上一讲用 LoFreq 和 GATK Mutect2 两个工具完成了变异检测（variant calling），本讲学习如何整理、比较、注释这些变异，并判断哪些才是真正值得关注的。

2. 学习目标

依据课件第 2 页：
  1. 使用 BCFtools 处理 VCF 文件：查看、过滤、取交集、添加自定义 FILTER 标签。
  1. 从 UCSC Genome Browser 获取基因组注释（annotation track、保守性、Table Browser、Session 管理）。
  1. 使用 Variant Effect Predictor（VEP）及其过滤器识别感兴趣的变异，并结合 AlphaMissense 评分评估错义变异的致病性。

3. 核心内容

3.1 流程回顾：本讲在全流程中的位置

高通量测序数据分析的标准流程是：原始 reads 质控（fastQC）→ 比对到参考基因组（bwa mem，得到 SAM）→ 转换排序为 BAM（samtools sort）→ 标记重复（gatk MarkDuplicates）→ 碱基质量校正（gatk BaseRecalibrator + ApplyBQSR）→ 收集比对质控指标（bamtools stats、gatk CollectMultipleMetrics 等）→ 变异检测（gatk Mutect2 / LoFreq）得到 VCF。本讲从最后一格「VCF」出发：VCF（Variant Call Format，变异调用格式）是记录每个变异位置、碱基改变和质量信息的文本文件，拿到它只代表「找到了变异」，后续还要判断哪些变异可靠、落在什么功能区、有什么功能影响。
p3 流程回顾图
p3 流程回顾图
课件第 3 页：从 fastq 质控、比对、BAM 处理到 Mutect2 产出 VCF 的完整流程图，本讲内容紧接 VCF 之后。

3.2 VCF 的「二进制化」：bgzip + tabix

类比 SAM→BAM 的关系：VCF 是文本格式，处理大文件时需转为压缩且带索引的形式。与 BAM 不同，VCF 不转成单一二进制，而是用 bgzip 压缩（blocked gzip，分块压缩格式）、再用 tabix 建立坐标索引：
两点注意（课件第 5 页 + 课堂录音）：
  • 不能用普通 gzip：tabix 只能索引 bgzip 格式；bgzip 和 tabix 都来自 htslib ，与 samtools、bcftools 出自同一作者、互相兼容。
  • 若 VCF 未按坐标排序，tabix 会报错，需先排序；由 LoFreq 等较老工具生成的 VCF 头部可能不符合最新 BCF 规范，bcftools 会给 warning，可忽略。
bcftools （课件第 4 页）是操作 VCF/BCF 的命令行工具，功能包括：查看（view/head）、过滤（view/filter）、取子集、改表头（reheader）、合并（merge）、标准化（normalize）、注释（annotate）、统计及格式转换。

3.3 查看与过滤 VCF：view 与 filter

常用查看命令：
过滤（filtering）有三类依据（课件第 7 页）：坐标区间（genome region）、预计算的 FILTER 列、INFO 列中的变异信息。FILTER 列 记录变异检测工具给出的质控判定（如 PASS、MinDP10 等）；INFO 列 存放测序深度（DP, depth）、等位基因频率（AF, allele frequency）、链偏倚等原始统计量。
  • bcftools view -f PASS：最直接的过滤，只保留 PASS 变异。本例 518 个变异过滤后剩 216 个 PASS。
  • bcftools view -i 'INFO/DP>=20' / -e 'INFO/DP<20'：按 INFO 字段过滤。-i（include）写保留条件，-e（exclude）写排除条件，两者互补，选哪种取决于哪种表达更直观。
  • 条件可以叠加：课件练习「只保留 AF < 0.2 且 PASS 的变异」，即把 -f PASS 与 AF 条件链在同一条命令里，可得少数罕见变异。
  • bcftools filter -m + -s 'LowDP' -e 'DP<5'：更复杂的过滤，并可改写 FILTER 列——给 DP<5 的变异追加一个名为 LowDP 的标签（-m + 追加而不覆盖原有标签，-s 指定标签名）。课堂演示命名为「MinDP5」；这样自定义的过滤器可以告诉他人你的过滤标准，便于结果共享与复现。

3.4 用 bcftools isec 比较两个 VCF：LoFreq vs Mutect2

不同变异检测工具的结果不会完全一致。bcftools isec（intersection）比较多个 VCF，输出重叠与独有的变异：
-f "PASS,germline" 表示只比较 FILTER 为 PASS 或 germline 的变异（germline 指与公共数据库已知胚系变异重合、被标为 germline 的位点，本质上也是可信变异）。结果生成 4 个文件加一个 README：
  • 0000.vcf：仅第一个文件（LoFreq）独有的变异——196 个
  • 0001.vcf：仅第二个文件（Mutect2）独有的变异——62 个
  • 0002.vcf / 0003.vcf：两者共有的变异——20 个 （内容相同，但分别保留来自两个 caller 的记录信息）
p13 isec 结果
p13 isec 结果
课件第 13 页：LoFreq 与 Mutect2 PASS/germline 变异比较结果——LoFreq 独有 196 个，共有 20 个，Mutect2 独有 62 个。
在 IGV（Integrative Genomics Viewer，基因组浏览器）中逐一检查后可解释差异的来源（课件第 14 页 + 课堂录音）：
  • LoFreq 独有的变异多为成簇突变 （clustered mutations）：相邻位置同时出现多个突变，往往是比对错误（misalignment）所致，可信度低。
  • Mutect2 独有的变异多有 DP<10（深度低）、与 indel 重叠、或位于同源多聚物（如连续 T 碱基）区域——短读长测序在重复碱基串上容易数错拷贝数（链滑移，strand slippage）。
  • 两者共有的变异则很干净，例如某位点支持碱基恰好 50:50，是典型的杂合（heterozygous）变异。
bedtools 复习 ：bedtools intersect（旧名 intersectBed）可做同样的事——支持 BED/GFF/VCF/BAM 输入，-v 参数输出 A 中不与 B 重叠的条目（独有变异），交换 -a/-b 分别得到两边的独有集。课堂验证：bedtools 得到与 isec 相同的数字（196/62/20）。区别在于 isec 一条命令同时给出全部四种组合，bedtools 则需多次调用或先提取 PASS 子集。

3.5 UCSC Genome Browser：注释的来源

UCSC Genome Browser 是基因组注释的重要在线来源。课堂以 hg38（人类参考基因组第 38 版）中的 FGFR2 基因为例复习操作，并新增两个为本讲和作业服务的功能：
  1. 从零搭建视图 ：默认 track 太多，可用「hide all tracks」清空，再逐个添加需要的 track（如 GENCODE 基因注释），点 refresh 刷新——做论文配图时常用。
  1. Track search 快速找 track ：直接搜索关键词（如 conservation）比在长列表里翻找更快。
  1. 保守性 track（conservation） ：基于约 30 种哺乳动物的多序列比对，某位置在越多物种中相同，分数越高。可见外显子（编码区）的保守性普遍高于内含子和基因间区——进化上重要的序列受选择压力，不易变动。
p21 保守性
p21 保守性
课件第 21 页：UCSC Genome Browser 保守性 track——外显子区域的保守性得分明显高于周边序列。
  1. OMIM 疾病变异 track ：ClinVar/OMIM 类 track 显示与疾病关联的已知变异。FGFR2 上放大可见某变异与一种颅缝早闭类发育综合征关联；同一氨基酸位置可对应多种不同碱基改变、不同疾病。手工在浏览器里比对自己的变异与疾病变异是可行但低效的方式——这正是后文用工具自动注释的动机。
  1. Session（会话）保存与分享 ：My Data → My Sessions 可把当前视图（位置、track 设置）存成 .txt 文件；对方加载该文件即可在任意浏览器完全复现你的视图。Assignment 2 将要求提交 session 文件 。

3.6 Table Browser：下载编码外显子坐标

UCSC Table Browser（Tools → Table Browser）可批量导出注释区间。课堂演示从 GENCODE v48（基因注释数据库）导出全基因组 （注意选 genome 而非当前显示的 position 范围）的编码外显子（coding exons），格式选 BED（基因组坐标区间格式），命名为 hg38_gencode_exon_coding.bed。BED 文件每一行是一个外显子区间，附带基因名、转录本名和 CDS（coding sequence，编码序列）标记；由于一个基因有多个转录本、外显子之间大量互相重叠，后续取交集时需注意去重（见 3.8 的 -u 参数）。

3.7 变异的分类与功能影响

拿到变异后首先关心它落在基因组哪类区域（课件第 26–28 页）：
  • 基因间区（intergenic region） ：约占基因组 2/3，全基因组测序检出的变异大多数在此。
  • 基因内非编码区（non-coding） ：如内含子，可能影响基因调控。
  • 编码区（coding） ：落在外显子内，进一步分为：
    • 同义变异（synonymous） ：碱基改变但氨基酸不变；
    • 非同义变异（non-synonymous） ：氨基酸改变，包括错义（missense，改变蛋白序列/结构） 和无义（nonsense，提前终止、截断蛋白） 。
p26 变异分类
p26 变异分类
课件第 26 页：按基因结构（TSS=转录起始位点）对 SNPs/突变的分类示意——基因间区、非编码区、编码区（再分同义/非同义）。
功能影响要点：
  • 序列变异是进化（自然选择）的原料，也用于法医 DNA 指纹；只有小部分直接致病或影响疾病风险与治疗反应。
  • 非同义变异改变蛋白；同义或非编码变异也可能通过改变转录/翻译效率、mRNA 稳定性、转录因子（TF）结合或 miRNA 结合而发挥作用。
  • 绝大多数变异是中性（neutral）的：不足 1% 会带来可察觉的表型改变。 因此「找到变异」远远不够，必须进一步筛选出「有后果」的变异。课程示例中一个 A→C 变异使谷氨酰胺（Glutamine）变为天冬氨酸（Aspartic acid），即是非同义变异。

3.8 用 bedtools 找编码区变异

有了外显子 BED 文件，即可用区间交集筛选编码变异：
-u（unique）参数很关键：每个变异无论与多少个外显子重叠都只输出一次——否则同一变异会因不同转录本的外显子重复出现而被打印多遍。结果：216 个 PASS 变异中只有 24 个 落在编码区，直观体现了「绝大多数变异不在编码区」。
课堂练习将命令用管道（pipe）串联，筛选纯合（homozygous，AF=1）的编码变异 ：
得到 24 个编码变异中的 16 个 纯合变异。两个细节（课堂录音）：
  • 管道中给 -a 传 - 表示从标准输入读取上一条命令的输出，可省去中间文件。
  • 此时不能加 -H ：去掉表头后 bedtools 无法识别这是 VCF 格式；-H 只适合配合 wc -l 计数等场景。

3.9 VEP：在线注释变异的功能后果

以上都是「手工」区间层面的注释，真正的功能注释交给专用工具。VEP（Variant Effect Predictor） 是 Ensembl（欧洲生物信息学联盟）提供的在线工具（www.ensembl.org/Tools/VEP），可连接 ClinVar、OMIM、PharmGKB、UniProt 等数据库为变异提供上下文。使用流程（课件第 30–34 页）：
  1. 先清洗数据 ：用 bcftools 过滤出 PASS（或纯合/杂合）变异再上传——把未通过质控的变异交给 VEP 没有意义。
  1. 上传前必须选对基因组版本 ：本例选 GRCh38（即 UCSC 的 hg38，VEP 默认配 GENCODE 50 转录本集）。版本选错则所有坐标都错，注释完全无效。
  1. 参数选择：变异命名法选 HGVS （Human Genome Variation Society，人类基因组变异协会命名标准，临床通用，便于交流）；人群频率可加 gnomAD/ExAC（显示每个变异在数据库人群中的携带情况）；也可勾选表型（phenotype）注释。
  1. 结果是一个带后果（consequence）、基因、HGVS 名称、频率等列的变异列表，可用过滤器（filter）按分子后果筛选，如只看 stop gained（提前终止密码子获得） 或 missense 变异。
p31 VEP 界面
p31 VEP 界面
课件第 31 页：VEP 在线工具——选择 hg38 后上传 WXS_germline_filter_PASS.vcf.gz，并连接 ClinVar、OMIM 等数据库注释。
课堂示例（MSH6 基因）：VEP 在结果中共标出 4 个 stop gained 变异。回 IGV 检查（右键 Sort by base 按碱基排序，chr2:47800395–47800415）：其中 2 个频率很低且相邻、来自软剪切读段（soft-clipped reads），判断为假阳性 ；另有 1 个真实的 CGA→TGA 变异——TGA 是终止密码子，该变异会使蛋白提前截断。课堂上老师对同一基因出现多个 stop gained 表示怀疑并逐一核查，示范了「工具注释 → IGV 人工验证」的完整判断链。VEP 结果界面功能有限，正式分析建议下载完整结果文件。

3.10 AlphaMissense：错义变异致病性预测

AlphaMissense 是 DeepMind 开发的 AI 模型，在 AlphaFold2（蛋白质结构预测模型）基础上整合预测的结构信息与进化信息，对人类蛋白组中所有可能的错义变异预测致病性（课件第 38–39 页）：
  • 输出概率得分 0–1：> 0.8 判为有害（damaging），< 0.2 大概率良性（benign） ，介于两者之间为不确定。
  • 只适用于错义变异 （missense）；无义突变本身就截断蛋白、必然有害，无需该评分。
  • 课堂演示：在 VEP 结果中过滤 missense 变异并叠加 AlphaMissense > 0.8 条件，得到少数高分变异；其中一个位于 MSH2 基因、人群频率极低、预测影响高——稀有 + 高影响 + 重要基因，是典型的需要优先分析的候选。
p39 AlphaMissense 评分
p39 AlphaMissense 评分
课件第 39 页：AlphaMissense 概率得分解释——0 到 1，>0.8 视为 damaging，<0.2 视为 probably benign。

3.11 下讲预告

下一讲进入 RNA-seq 数据分析（差异表达）；Assignment 2 将于下周发布 ，内容基于过去三周的材料（不含 RNA-seq）。

4. 关键图解

插图
内容
assets/p3-pipeline-recap.png
全流程回顾：QC → 比对 → BAM 处理 → Mutect2 → VCF，本讲从 VCF 继续
assets/p13-isec-venn.png
LoFreq 与 Mutect2 变异交集：196 独有 / 20 共有 / 62 独有
assets/p21-conservation-exons.png
UCSC 保守性 track：外显子区保守性显著升高
assets/p26-snp-types.png
变异按基因组位置分类：基因间区 → 非编码 → 编码 → 同义/非同义
assets/p31-vep-upload.png
VEP 上传界面：选 hg38、上传 PASS VCF
assets/p39-alphamissense-score.png
AlphaMissense 得分阈值：>0.8 damaging，<0.2 benign

5. 课前思考 / 自测

  1. 为什么 VCF 压缩必须用 bgzip 而不是普通 gzip？tabix 索引在整个流程中起什么作用？
  1. bcftools view -i 与 -e 有何区别？写出「保留 AF<0.2 且 PASS 的变异」的命令。
  1. bcftools isec 输出的 0002.vcf 与 0003.vcf 都含共有变异，两者有何不同？
  1. intersectBed 加 -u 与不加 -u 的输出有何差异？在什么场景下必须加？
  1. 同义（synonymous）变异不改变氨基酸，是否就一定没有功能影响？举例说明其他可能的作用机制。
  1. VEP 提交时选错基因组版本（如把 hg38 数据提交到 hg19）会发生什么？AlphaMissense 评分为什么不适用于 stop gained 变异？

6. 参考材料

  • 讲课 PDF：Materials/Variant_annotation_2026.pdf（41 页，主课件）
  • 实操 Notebook：Materials/Variant_Annotation.ipynb（Google Colab，含全部命令与练习）
  • 课堂录音：Recordings/MSBS7101-1020-summary.md、Recordings/MSBS7101-1020-transcript.md

7. 课堂要点补充（基于课堂录音）

  • 工具族谱 ：bgzip、tabix 来自 htslib，与 samtools、bcftools 同一作者出品，格式互操作；这是 VCF 压缩必须用 bgzip 的原因。
  • LoFreq 旧版表头 warning 可忽略 ：bcftools 提示表头不完全符合最新 BCF 规范，属警告而非错误。
  • 假阳性判读案例 （IGV 实操）：成簇出现的相邻突变多为比对错误；位于连续相同碱基（同源多聚物）区的 indel 附近变异可能源于测序链滑移；与软剪切读段相关的低频 stop gained 不可靠；真实变异表现为支持碱基比例清晰（如杂合约 50:50）。
  • 命令组织习惯 ：过滤条件用 -i/-e 哪种写法取决于哪个更直观；条件可链式叠加；-H 只在计数/查看时用，管道传给 bedtools 时必须保留表头；-a - 可让 bedtools 直接读管道输入。
  • UCSC 实操提示 ：做图前先 hide all 再按需加 track；用 track search 快速定位 conservation 等 track；Table Browser 导出时务必选 genome 全基因组范围；session 文件本质是记录位置与 track 设置的文本，Assignment 2 需要提交。
  • 数据规模直觉 ：人类个体与参考基因组相比约有数百万个变异，其中编码变异只占很小部分（本例 216 个 PASS 中仅 24 个编码）；不足 1% 的变异有可察觉表型效应。
  • VEP 当天网页故障 ：课堂现场 VEP 网页一度无法运行（界面数月前刚更新、报错无提示），老师口头讲解了完整流程与预期结果；课后可自行重跑，也可对 Mutect2 的 VCF 重复练习。
  • 选变异的优先级思路 （老师总结）：人群频率极低 + 预测影响高 + 位于重要基因（如 MSH2/MSH6，DNA 错配修复基因）的变异应优先深入分析。
上一篇
MSBS7102 Lecture 7 · R 迷你项目：达尔文实验与配对 t 检验
下一篇
MSBS7102 Lecture 6 · ggplot2 数据可视化
Loading...
2026-11-3