01 基因组学实践教程

基因组学研究的是 DNA 序列本身：谁的基因组里有什么变异，这些变异是否落在功能区域，以及在群体里如何分布。即便今天大多数湿实验都在转录组或单细胞这一层，基因组重测序仍然是找到疾病相关变异、做 GWAS、做群体遗传分析的起点。

本专栏打算把"拿到一批测序数据之后怎么走到变异列表"这条主线讲清楚。重点放在重测序和短读长小变异，其他方向（长读长、结构变异、甲基化等）会作为延伸而不是主干。

基因组学项目的三条典型主线

新手容易把"基因组学"当成一回事，其实它是三种很不同的项目类型：

这三条主线分析流程差很多：胚系关注致病性、肿瘤关注 VAF 和驱动基因、群体关注频率和 LD。先想清楚自己在做哪一种，再选工具，比闷头跑流水线效率高。

一个项目大致是什么样子

一个典型的 WES 或 WGS 项目会拿到若干个体的双端 FASTQ 数据。从这里到"一份可以注释、可以过滤、可以做关联分析的 VCF"之间，要走的步骤大致是：

"BWA + GATK" 是短读长小变异的事实标准，文献、社区和官方 best practices 都围绕它展开。如果只做基因型分型而不做科研级变异检测，bcftools mpileup + bcftools call 也能跑出可用结果，代码更短。

常见工具栈

这套组合在教学和小到中等规模真实项目里基本够用：

长读长（PacBio、ONT）和结构变异会用到另一套工具链（minimap2、Sniffles、CuteSV 等），这里暂不展开。

推荐公开数据集

学这个专栏可以从三类数据入手：

NA12878 有 GIAB 发布的"黄金标准" VCF，可以直接用来评估你自己跑出的结果是否合理，是练习 variant calling 流程最方便的参照物。

最小可跑的例子

下面用 Bioconductor 的 VariantAnnotation 包做一次最简单的 VCF 读取和浏览。它自带一份 chr22.vcf.gz 示例，不需要联网下载：

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("VariantAnnotation", "TxDb.Hsapiens.UCSC.hg19.knownGene"))

library(VariantAnnotation)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

# 读入包里自带的 chr22 示例 VCF
fl <- system.file("extdata", "chr22.vcf.gz", package = "VariantAnnotation")
vcf <- readVcf(fl, "hg19")
vcf

# 看一下 VCF 结构
header(vcf)
rowRanges(vcf)[1:5]

# 提取基因型信息
geno(vcf)$GT[1:5, 1:5]

# 把变异定位到基因区域
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
loc <- locateVariants(vcf, txdb, CodingVariants())
head(loc)

# 统计每类变异的数量（内含子、外显子、UTR 等）
all_loc <- locateVariants(vcf, txdb, AllVariants())
table(all_loc$LOCATION)

locateVariants 会把每个 VCF 记录映射到最近的基因和区域类型。table(all_loc$LOCATION) 的输出是典型的"变异区域分布"结果——UTR、外显子、内含子、基因间各多少——这也是 variant calling 报告里必写的一栏。

真实重测序项目比这段要多两件事：前半段的 BWA 比对 + GATK 变异检测流水线（要跑几个小时、要有足够磁盘）、后半段的过滤策略和多样本联合分型。但 VCF 拿到手之后的操作思路，和上面这段是一致的。

专栏模块规划

所有 8 个模块已上线。03 和 04 带可跑脚本和真实数据图。

推荐前置知识

• 编程基础：R、Python、Bash 学习路径
• Jupyter 与交互式分析环境
• 公共数据库与数据检索

常见误区

误区 1：以为变异越多越好

WGS 跑出来 5M 变异，99% 是已知 SNP，对你研究问题毫无意义。真正有价值的是罕见变异 + 致病变异，需要靠 ClinVar / gnomAD 过滤。"我的样本有 5M 个变异" 不是亮点。

误区 2：用 SNV-only 思路做肿瘤分析

肿瘤里 CNV、SV、染色体重排同样关键 — 漏掉这些等于半残。WGS 项目要补 Manta / GRIDSS（SV）和 ASCAT / FACETS（CNV）。

误区 3：把胚系变异工具用在肿瘤上

HaplotypeCaller 假设变异在 50% 或 100% 的 reads 里（杂合或纯合）。肿瘤是异质性群体，VAF 可以是 5%、20%、80% 等任意值，必须用 Mutect2 或类似肿瘤专用 caller。

误区 4：不做 BQSR 直接 call 变异

BQSR 校正测序仪系统偏差，不做的话假阳性率提高 10-30%。GATK best practices 里这是必须步骤，不要图省事跳过。

误区 5：群体遗传项目不做 LD pruning 直接 PCA

PCA 会被高 LD 区域（HLA、着丝粒等）主导，主成分变成"局部 LD 模式"而不是真实的群体结构。PCA 前必须 LD prune（PLINK --indep-pairwise 50 5 0.2）。

参考资源

• GATK Best Practices
• Genome in a Bottle
• 1000 Genomes Project
• gnomAD
• Ensembl VEP
• samtools / bcftools 手册

02 数据类型：WGS vs WES vs Panel

基因组测序有三种主要策略，选择取决于研究目的和预算：

WGS vs WES 的分析差异

肿瘤 WES 的特殊性

肿瘤样本通常配对测序（tumor + matched normal），用 Mutect2 做体细胞突变检测。输出的 MAF 文件是 maftools 的标准输入。

# Mutect2 典型调用
gatk Mutect2 \
  -R reference.fa \
  -I tumor.bam \
  -I normal.bam \
  -normal normal_sample_name \
  -O somatic.vcf.gz

# VCF -> MAF 转换
vcf2maf.pl --input-vcf somatic.vcf.gz --output-maf somatic.maf \
  --tumor-id tumor --normal-id normal --ref-fasta reference.fa

怎么选：WGS / WES / Panel 的判断逻辑

实务上的选择标准比 textbook 里复杂。一份决策表：

新人最常见错误：用 WES 找结构变异。WES 的捕获 enrichment 让 break point 上的 reads 严重不均匀，CNV / SV 在 WES 上几乎不可靠。

常见坑

坑 1：tumor-only 没有 matched normal

没有配对正常样本时，无法区分体细胞突变和胚系变异。tumor-only 模式可以跑（Mutect2 --tumor-only），但要严格用 gnomAD 频率过滤掉所有 AF > 0.001 的位点，否则结果是肿瘤 + 胚系混合。

坑 2：WES 的 BED 文件版本不一致

不同 WES 试剂盒（Agilent SureSelect / Twist / IDT xGen）有自己的 target BED，用错 BED 会让 70% 的覆盖度统计是假的。BED 必须从厂商网站下载和试剂盒版本对应的那一份。

坑 3：Panel 数据用 WGS 的过滤参数

Panel 测序深度 500-2000x，WGS 30-50x。用 WGS 的 GQ ≥ 20 阈值在 Panel 上会几乎不过滤（Panel 的 GQ 普遍很高），需要用更严格的标准比如 VAF + 链平衡。

坑 4：用错 reference fasta

GATK 要求 reference 必须有 .dict 和 .fai 索引，且染色体顺序匹配。hg19 和 hg38 不能混用，混用会让所有变异坐标错位。

坑 5：Mutect2 跑 tumor-only 不加 PoN

Panel of Normals (PoN) 是从一批正常样本里收集的"看上去像变异但其实是测序伪影"的位点。没 PoN 的 Mutect2 假阳性率会很高，至少要做一份 50-100 个正常样本的 PoN。

下一步

接着深入：

• 02 比对与变异检测流程 — 从 FASTQ 到 VCF 的完整命令行流水线
• 03 VCF 注释与可视化 — VCF 拿到后第一步：注释功能影响

横向延伸：

• 04 maftools 肿瘤突变分析 — 如果是肿瘤项目，看完 02 直接跳这里
• 05 变异过滤与质量评估 — variant filtering 决定结果质量

参考资源

• GATK Best Practices
• Mutect2 文档
• vcf2maf

03 比对与变异检测流程

从 FASTQ 到 VCF 的标准流水线。这一章给出完整的 bash 命令序列，适合在服务器上跑。

BWA-MEM 比对

# 建索引（一次性）
bwa index reference.fa

# 比对（每个样本）
bwa mem -t 8 -R "@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA" \
  reference.fa sample1_R1.fq.gz sample1_R2.fq.gz \
  | samtools sort -@ 4 -o sample1.sorted.bam

samtools index sample1.sorted.bam

-R 参数加 read group 信息，GATK 后续步骤必须有。

去重复 + BQSR

# 标记 PCR 重复
gatk MarkDuplicates \
  -I sample1.sorted.bam \
  -O sample1.dedup.bam \
  -M sample1.dedup_metrics.txt

# Base Quality Score Recalibration
gatk BaseRecalibrator \
  -R reference.fa \
  -I sample1.dedup.bam \
  --known-sites dbsnp.vcf.gz \
  --known-sites mills_indels.vcf.gz \
  -O sample1.recal_table

gatk ApplyBQSR \
  -R reference.fa \
  -I sample1.dedup.bam \
  --bqsr-recal-file sample1.recal_table \
  -O sample1.recal.bam

HaplotypeCaller（胚系变异）

# 单样本模式（输出 gVCF）
gatk HaplotypeCaller \
  -R reference.fa \
  -I sample1.recal.bam \
  -O sample1.g.vcf.gz \
  -ERC GVCF

# 多样本联合分型
gatk CombineGVCFs -R reference.fa \
  -V sample1.g.vcf.gz -V sample2.g.vcf.gz \
  -O cohort.g.vcf.gz

gatk GenotypeGVCFs -R reference.fa \
  -V cohort.g.vcf.gz \
  -O cohort.vcf.gz

Mutect2（体细胞变异）

gatk Mutect2 \
  -R reference.fa \
  -I tumor.recal.bam \
  -I normal.recal.bam \
  -normal normal_sample \
  --germline-resource gnomad.vcf.gz \
  -O somatic_unfiltered.vcf.gz

gatk FilterMutectCalls \
  -R reference.fa \
  -V somatic_unfiltered.vcf.gz \
  -O somatic_filtered.vcf.gz

流水线工具

手动跑上面这些命令容易出错。推荐用 nf-core/sarek：

nextflow run nf-core/sarek \
  --input samplesheet.csv \
  --genome GRCh38 \
  --tools mutect2,haplotypecaller \
  --outdir results/

为什么 GATK Best Practices 这么麻烦

新手第一次跑这套常常想问"为什么要这么多步？" 每一步都是为了解决一个具体的假阳性来源：

每一步都不是装饰。理解每一步在防什么，比死记参数更有价值。

常见坑

坑 1：忘了 read group 信息

GATK 后续步骤强制要求 BAM 里有 @RG 标签。bwa mem -R "@RG\tID:..." 一开始就要加。否则 MarkDuplicates 直接报错。后期补 RG 用 gatk AddOrReplaceReadGroups。

坑 2：BWA 用了 -M 标志

老教程用 -M 让多比对 reads 标记为 secondary，但新版 GATK + Picard 期望 supplementary。用 -M 在新版 pipeline 里会让 MarkDuplicates 行为异常。不加 -M 是默认正确。

坑 3：GVCF 文件多样本时手动 cat

GVCF 是结构化文件，不能直接 cat。必须用 gatk CombineGVCFs 或 GenomicsDBImport，否则下游报错。

坑 4：Mutect2 不指定 germline-resource

--germline-resource gnomad.vcf.gz 让 Mutect2 把已知胚系变异从体细胞结果里去掉。不加这个参数 tumor-normal 配对的结果还能凑合，tumor-only 模式直接错。

坑 5：BAM 文件没 sort 就 dedup

MarkDuplicates 要求 BAM 已按坐标排序。bwa mem 输出的 SAM 是 read-name 顺序，要先 samtools sort -@ 4。直接 dedup 会报错或给出错误结果。

下一步

接着深入：

• 03 VCF 注释与可视化 — VCF 拿到手后的第一步
• 05 变异过滤与质量评估 — 决定结果可信度的关键

横向延伸：

• 04 maftools 肿瘤突变分析 — Mutect2 输出的 MAF 用 maftools 出图
• nf-core/sarek 教程 — 一键跑完整 pipeline

参考资源

• BWA 手册
• GATK Best Practices
• nf-core/sarek
• DeepVariant

04 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 不够用时换更强大的注释工具

参考资源

• VariantAnnotation
• Ensembl VEP
• ANNOVAR

05 maftools 肿瘤突变分析

maftools 是肿瘤外显子组（WES）分析里最常用的 R 包。它读入 MAF 格式的体细胞突变文件，一套函数出完整的突变景观图、驱动基因分析、突变签名等。

本章用 maftools 自带的 TCGA LAML（急性髓系白血病，193 个样本）数据演示。

MAF 格式

MAF（Mutation Annotation Format）是 TCGA 定义的标准格式，每行一个突变，关键列：

从 VCF 转 MAF 用 vcf2maf.pl（需要 VEP 注释）。

真实示例

配套脚本 genome04_maftools_sci.R 输出 6 张图：

Rscript scripts/genomics/genome04_maftools_sci.R

每张图看什么

图 1：MAF 总览 dashboard。左上：每个样本的突变数；右上：变异分类分布；左下：变异类型；右下：SNV 碱基替换类型。一张图看完整个队列的突变概况。

图 2：Oncoplot（突变景观图）。每列一个样本，每行一个基因（按突变频率排序）。颜色代表突变类型。这是肿瘤基因组文章里最核心的一张图 —— 一眼看出哪些基因在队列里反复突变。

LAML 里 FLT3、DNMT3A、NPM1 是 top 3 驱动基因，和文献完全一致。

图 3：DNMT3A 的 lollipop 图。横轴是蛋白结构域，每个棒棒糖代表一个突变位点，高度代表该位点在队列里出现的次数。R882 是 DNMT3A 的热点突变，在 AML 里反复出现。

图 4：基因间的共突变 / 互斥关系。绿色 = 共现（co-occurrence），粉色 = 互斥（mutual exclusivity）。互斥的基因对可能在同一条通路上（突变一个就够了）。

图 5：肿瘤突变负荷（TMB）分布。AML 是低 TMB 肿瘤（中位 ~10 个突变/样本），和黑色素瘤、肺癌（几百个）形成对比。TMB 是免疫治疗响应的预测指标之一。

图 6：转换/颠换比（按样本）。Ti/Tv 偏离正常范围可能提示特定的突变过程（如 APOBEC 活性会增加 C>T 转换）。

核心代码

library(maftools)

laml <- read.maf(
  maf = system.file("extdata", "tcga_laml.maf.gz", package = "maftools"),
  clinicalData = system.file("extdata", "tcga_laml_annot.tsv", package = "maftools")
)

plotmafSummary(laml, dashboard = TRUE)
oncoplot(laml, top = 10)
lollipopPlot(laml, gene = "DNMT3A", AACol = "Protein_Change")
somaticInteractions(laml, top = 15)

套到自己数据上

把 read.maf() 的路径换成自己的 MAF 文件即可。如果只有 VCF：

# VCF -> MAF（需要 VEP + vcf2maf）
vcf2maf.pl --input-vcf somatic.vcf.gz --output-maf somatic.maf \
  --tumor-id TUMOR --normal-id NORMAL --ref-fasta hg38.fa \
  --vep-path /path/to/vep --vep-data /path/to/vep_cache

下载资源

不想本地装环境？在 BioF3 上跑

06 变异过滤与质量评估

原始 VCF 里有大量假阳性。过滤策略决定了最终结果的可靠性。

VQSR vs 硬过滤

硬过滤典型参数

# SNP 过滤
gatk VariantFiltration -R ref.fa -V raw.vcf.gz \
  --filter-expression "QD < 2.0" --filter-name "LowQD" \
  --filter-expression "FS > 60.0" --filter-name "HighFS" \
  --filter-expression "MQ < 40.0" --filter-name "LowMQ" \
  --filter-expression "MQRankSum < -12.5" --filter-name "LowMQRS" \
  --filter-expression "ReadPosRankSum < -8.0" --filter-name "LowRPRS" \
  -O filtered_snps.vcf.gz

# Indel 过滤
gatk VariantFiltration -R ref.fa -V raw.vcf.gz \
  --filter-expression "QD < 2.0" --filter-name "LowQD" \
  --filter-expression "FS > 200.0" --filter-name "HighFS" \
  --filter-expression "ReadPosRankSum < -20.0" --filter-name "LowRPRS" \
  -O filtered_indels.vcf.gz

质量评估指标

• Ti/Tv ratio：WGS ~2.0-2.1，WES ~2.8-3.0。偏低 = 假阳性多
• Het/Hom ratio：人类 WGS ~1.5-2.0
• 已知变异比例：和 dbSNP 的重叠率应该 > 95%（WGS）
• Mendelian error rate：有 trio 数据时，子代不符合孟德尔遗传的比例应 < 1%

一份 VCF 的诊断 checklist

拿到一份新 VCF 时，按这个顺序快速判断质量：

1. 总变异数：WGS ~4-5M，WES ~50-100k，Panel <1000。数量异常意味着流水线有问题
2. Ti/Tv：上面已说，偏离正常意味着假阳性
3. dbSNP 重叠率：> 95% 健康，< 90% 有问题
4. Het/Hom：偏离 1.5-2.0 可能样本污染或近交
5. Singleton 比例：单次出现的变异比例。WGS 应 < 5%，过高说明假阳性
6. 覆盖度均匀性：变异 DP 分布的 CV 应该 < 0.5

这 6 个 5 分钟能查完，比单纯看 padj 信息量大得多。

常见坑

坑 1：VQSR 对小样本量直接套

VQSR 需要 30+ 样本才能学出可靠模型。< 30 样本必须用硬过滤，硬上 VQSR 结果不可信。

坑 2：硬过滤参数照搬 GATK 推荐值

GATK 推荐值是 WGS 的经验值，对 WES、Panel、低深度数据不一定合适。先看自己数据 QD / FS / MQ 的分布，再选阈值。

坑 3：过滤了 PASS 之外的所有

GATK Best Practices 把过滤标记写进 FILTER 列（PASS 或具体过滤名），不是删除。bcftools view -f PASS 才能拿到只有 PASS 的子集。直接读完整 VCF 后忘了筛 FILTER 列是常见错误。

坑 4：用 GQ 阈值过滤 INDEL

GQ（genotype quality）对 SNP 比较合理，INDEL 的 GQ 经常很低（call indel 本身难），用 GQ ≥ 20 会丢一半真 INDEL。INDEL 用 QD + FS 过滤更合适。

坑 5：trio 项目不查 Mendelian error

有父母和子代的 trio 数据，子代不符合孟德尔遗传的位点几乎都是假阳性。vcftools --mendel 可以查这个错误率，超过 1% 说明数据质量有问题。

下一步

接着深入：

• 06 群体遗传：PCA / 祖源 / LD — 过滤后的 VCF 进入群体分析
• 07 GWAS 入门 — GWAS 项目的过滤更严格（HWE、call rate 等）

横向延伸：

• 03 VCF 注释与可视化 — 过滤前后做注释对比
• GATK 硬过滤指南 — 官方阈值的来源和细节

参考资源

• GATK 硬过滤指南
• GATK VQSR 教程

07 群体遗传：PCA / 祖源 / LD

群体遗传分析用大量个体的基因型数据回答"这些人从哪来、彼此什么关系、哪些位点受到选择"。

常用工具

PCA 分析

# VCF -> PLINK 格式
plink2 --vcf cohort.vcf.gz --make-bed --out cohort

# LD pruning（去掉高 LD 的 SNP，避免 PCA 被局部 LD 主导）
plink2 --bfile cohort --indep-pairwise 50 5 0.2 --out pruned
plink2 --bfile cohort --extract pruned.prune.in --make-bed --out cohort_pruned

# PCA
plink2 --bfile cohort_pruned --pca 10 --out cohort_pca

输出 cohort_pca.eigenvec 就是每个样本的 PC1~PC10 坐标，用 R 画散点图。

ADMIXTURE 祖源分析

# 跑 K=2 到 K=6
for K in 2 3 4 5 6; do
  admixture --cv cohort_pruned.bed $K | tee log_K${K}.out
done

# 选最优 K：看 CV error 最低的那个
grep "CV error" log_K*.out

LD 衰减

# 计算 LD（r²）随距离的衰减
plink2 --bfile cohort --ld-window-r2 0 --ld-window 1000 --ld-window-kb 500 \
  --out ld_decay

LD 衰减速度反映群体的有效群体大小和重组率。

群体遗传分析的核心问题

在跑工具之前先想清楚要回答什么。这类项目大致是三类问题：

新人最常见错误：跑了 PCA 看到分群就完事，没思考"为什么这样分群"。PCA 只给坐标，不给解释。

常见坑

坑 1：PCA 不做 LD pruning

不 prune 时 PCA 会被 HLA、着丝粒等高 LD 区域主导，PC1 变成 "HLA 单体型聚类" 而不是真实群体结构。--indep-pairwise 50 5 0.2 是必须步骤。

坑 2：ADMIXTURE 跑一个 K 就报告

K 选择不是凭感觉。至少跑 K=2 到 K=10，看 CV error 最低的那个。CV error 平的话说明数据本身就没明显结构，硬选 K 没意义。

坑 3：用相关样本做群体遗传

第一代亲属的样本会扭曲所有群体结构估计。先用 PLINK --king-cutoff 0.0884 删掉 2 度以内亲属，再做下游分析。

坑 4：全样本一起 PCA 后才发现混杂

不同祖源的样本应该分别先做 QC，再一起跑分析。混在一起做 QC 时 MAF 阈值会被群体结构干扰。

坑 5：Fst 计算用了 SNV 全集

Fst 应该排除单态位点（MAF = 0）和高缺失率位点。vcftools --maf 0.05 --max-missing 0.95 是标准过滤。不过滤直接算 Fst 噪声很大。

下一步

接着深入：

• 07 GWAS 入门 — 群体结构 PC 是 GWAS 的标准协变量

横向延伸：

• 05 变异过滤与质量评估 — 群体遗传项目对过滤要求更严
• PLINK 2 GWAS 教程 — 群体遗传 + 关联分析的官方文档

参考资源

• PLINK 2 文档
• ADMIXTURE
• 1000 Genomes 群体结构教程

08 GWAS 入门

全基因组关联分析（GWAS）找的是"哪些基因组位点和某个表型（疾病、身高、药物响应）有统计学关联"。

基本流程

基因型数据（PLINK 格式）+ 表型文件
  → QC（MAF、HWE、缺失率、亲缘关系）
  → 关联检验（线性/逻辑回归，校正 PC）
  → Manhattan plot + QQ plot
  → 显著位点注释

PLINK 2 做关联

# 二分类表型（case/control）
plink2 --bfile cohort \
  --pheno phenotype.txt \
  --covar covariates.txt \
  --glm \
  --out gwas_results

# 输出 gwas_results.PHENO1.glm.logistic.hybrid

Manhattan plot

library(qqman)

results <- read.table("gwas_results.PHENO1.glm.logistic.hybrid",
                      header = TRUE)
manhattan(results, chr = "CHROM", bp = "POS", p = "P", snp = "ID",
          suggestiveline = -log10(1e-5), genomewideline = -log10(5e-8))

QQ plot

QQ plot 检查 p 值的整体膨胀（genomic inflation factor λ）。λ > 1.1 说明有群体分层或其他混杂没控制好。

qq(results$P)

常见坑

• 群体分层：不同祖源的人混在一起会产生假关联。用 PCA 的前几个 PC 做协变量
• 多重检验：全基因组显著性阈值是 5×10⁻⁸（Bonferroni 校正 ~1M 独立检验）
• LD 结构：一个显著信号可能对应一整个 LD block 里的几十个 SNP，真正的因果变异需要 fine-mapping

补充几条常见误判：

坑：样本量不足直接跑 GWAS

GWAS 需要大样本量（至少几千，理想几万）才能稳定检测到 OR ~1.1 的弱效应。< 1000 样本跑 GWAS 几乎不可能找到 5e-8 显著位点，先评估 power 再决定要不要做。

坑：不做 HWE 过滤

显著偏离 Hardy-Weinberg 平衡的位点通常是基因型分型错误。case/control 设计中只在 control 样本里检查 HWE（case 因为关联本身可能偏离 HWE）。

坑：QQ plot 严重 inflation 直接发表

λ > 1.05 通常说明群体分层没控制好。多加几个 PC 做协变量，重新跑直到 λ 回归 1.0 附近。

坑：找到一个显著 SNP 就当成 causal variant

一个显著信号可能对应一个 LD block 里的几十个相关 SNP，**真正的 causal 变异需要 fine-mapping（FINEMAP / SUSIE）**或功能验证（eQTL、CRISPR）。

坑：忽略 imputation 质量

INFO score < 0.3 的 imputed 变异不可信。imputation 后必须 --info 0.3 过滤，不然假阳性堆积。

下一步

接着深入：

• 08 临床变异解读与报告 — 找到的显著位点如何翻译成临床意义

横向延伸：

• 06 群体遗传：PCA / 祖源 / LD — GWAS 的群体结构 QC 是 06 章的延伸
• GWAS Catalog — 你找到的显著位点是不是已经被报道过

参考资源

• PLINK 2 GWAS 教程
• qqman R 包
• GWAS Catalog
• LocusZoom

09 临床变异解读与报告

从"一份 VCF"到"能给临床医生看的报告"，中间需要变异分级、数据库查询和标准化报告格式。

ACMG 变异分级

美国医学遗传学学会（ACMG）把胚系变异分为 5 级：

常用注释数据库

VEP 注释 + ClinVar

vep --input_file variants.vcf \
    --output_file annotated.vcf \
    --cache --dir_cache /path/to/vep_cache \
    --assembly GRCh38 \
    --everything \
    --plugin ClinVar,/path/to/clinvar.vcf.gz

报告模板

临床报告通常包含：

1. 患者信息 + 检测方法
2. 阳性发现（Pathogenic / Likely pathogenic 变异）
3. VUS 列表
4. 质量指标（覆盖度、Ti/Tv）
5. 方法学描述
6. 局限性声明

肿瘤报告的特殊性

肿瘤报告额外需要：

• 突变等位基因频率（VAF）
• 肿瘤突变负荷（TMB）
• 微卫星不稳定性（MSI）状态
• 可靶向突变（OncoKB Level 1-4）

临床报告 vs 科研分析的关键差异

科研分析容易把"显著"或"有趣"当成结论，临床报告则要回答"这个发现能不能影响治疗 / 诊断决策"。差异：

报告里最重要的句子：明确局限性。"本检测使用 ABC 试剂盒，覆盖 X 个基因 Y% 的 coding region，对 SV / CNV / 重复扩增不敏感"。这种话说清楚，临床医生才能正确解读"阴性结果"的含义。

常见坑

坑 1：把 dbSNP 收录当成"良性"

dbSNP 包含所有报道过的变异，不代表良性。判断良性要看 ClinVar 的 Benign / Likely benign 等明确分级 + gnomAD 群体高频。

坑 2：用错版本的 ClinVar

ClinVar 月更，临床报告必须冻结使用日期。不要每次跑都用 latest，那样无法回溯历史结论。

坑 3：VUS 当成"基本没事"沟通

VUS 字面意思是 "意义未明"，但实务上很多临床医生当成"良性偏倾向"。报告里要明确写"VUS 不应影响临床决策"，避免误解。

坑 4：不区分胚系和体细胞

有些肿瘤变异其实是胚系（如 BRCA1 致病变异），漏掉胚系判断会让患者亲属错过预防性筛查。Tumor + Normal 配对是必须的。

坑 5：报告里没有 raw QC 信息

只给变异列表不给覆盖度 / Ti/Tv / 比对率，复核者没法判断阴性结果是否可靠。报告附录至少包含每个基因的 mean coverage 和 callable region 比例。

下一步

基因组学专栏 9 篇到这里完成。

接着深入：

• 多组学整合 overview — 基因组数据和转录 / 表观整合
• 多组学整合 module05 MOFA2 — 临床多组学因子分析

横向延伸：

• 04 maftools 肿瘤突变分析 — 肿瘤报告的工具基础
• ACMG 2015 指南 — 临床变异分级的标准

参考资源

• ClinVar
• ACMG 2015 指南
• OncoKB
• InterVar
• Franklin（Genoox 自动分级）