03 VCF 注释与可视化
03 VCF 注释与可视化
VCF 文件拿到手之后，第一步是"每个变异落在什么基因、什么区域、有什么功能影响"。本章用 R 的 VariantAnnotation 包在内置 chr22 VCF 上演示。
真实示例
配套脚本 genome03_variant_anno_sci.R 输出 6 张图：
Rscript scripts/genomics/genome03_variant_anno_sci.R
每张图看什么
图 1：SNV vs Indel 的数量分布。WGS/WES 里 SNV 通常占 90%+。
图 2：变异落在哪些基因组区域。大部分在内含子和基因间区（非编码区占基因组 98%+）。
图 3：转换/颠换比。WGS 期望 ~2.0-2.1，WES 编码区 ~3.0。偏低可能说明假阳性多。
图 4：等位基因频率谱。经典 L 形：大部分变异是稀有的（低频）。
图 5：chr22 上的变异密度分布。某些区域密集可能对应基因密集区或重复序列。
图 6：编码区变异的功能后果（同义/错义/无义）。
下载资源
常见坑
坑 1：用错了 TxDb 物种 / 版本
TxDb.Hsapiens.UCSC.hg19.knownGene 和 TxDb.Hsapiens.UCSC.hg38.knownGene 不能混用。VCF header 里 ##reference 字段告诉你版本，先看再选 TxDb。
坑 2：变异类型只看 SNP / INDEL，忽略 MNP / SV
VCF 里还有 MNP（连续多碱基替换）、symbolic alleles（<DEL> / <INS> / <DUP> 这种 SV）等。isSNV(vcf) 只筛 SNP，简单脚本会漏掉一半。要按需判断每种类型。
坑 3：Ti/Tv 全队列算一个值
不同基因组区域 Ti/Tv 不一样：编码区高（~3.0），非编码区低（~2.0）。应该分区域看，混在一起的 Ti/Tv 没诊断价值。
坑 4：等位基因频率谱用 INFO/AF 直接画
VCF 的 INFO/AF 是 caller 估计值，不是真实群体频率。要看群体 AF 必须从 gnomAD / 1KG 数据库查，不能用 VCF 自带的 AF 当群体频率用。
坑 5：locateVariants 把 NA 当成"基因间"
部分变异 LOCATION 是 NA（坐标超出 TxDb 覆盖范围、未知 transcript），不是"基因间"。统计时要单独处理 NA，不然分布图被 NA 拉偏。
下一步
接着深入：
- 05 变异过滤与质量评估 — 注释完后过滤掉假阳性
- 08 临床变异解读与报告 — 用 ClinVar 注释做临床分级
横向延伸：
- 04 maftools 肿瘤突变分析 — 肿瘤项目的 VCF 走另一条专用注释路径
- VEP 在线注释 — VariantAnnotation 不够用时换更强大的注释工具