DNA 测序分析全流程教程

对应大纲第三阶段 3.1 DNA 分析方向目标：从零理解 WGS/WES/Panel 测序数据是怎么被分析出"结果"的

第一步：质控

FASTQ 格式解析，Phred 质量分的含义（Q30 = 错误率 0.1%）
FastQC 报告怎么看，Trimmomatic 修剪接头

第二步：比对

参考基因组 hg38 是什么，为什么需要"地图"
BWA-MEM2 比对命令及读组（Read Group）的作用
SAMtools flagstat 看比对率，MarkDuplicates 去重复

第三步：变异检测

SNP/Indel/CNV/SV 的区别
GATK HaplotypeCaller → VQSR 过滤假阳性的完整逻辑
肿瘤体细胞变异（Mutect2）和胚系变异的关键区别，TMB/MSI 临床意义

第四步：注释

ANNOVAR 把坐标翻译成"这个基因、这个变异、健康人频率、临床意义"

先建立直觉：测序在做什么？

DNA 测序仪不能一口气读完 30 亿个碱基。它的做法是：

把 DNA 打碎成几百个碱基的短片段
对每个短片段独立测序，得到大量"读段"（Read）
把这些读段拼回参考基因组上，推断原始序列

这就是为什么分析流程的核心是"比对"——把碎片拼回已知地图上。

原始 DNA（30 亿碱基）
    ↓ 打碎
几亿个短读段（每条 150bp）
    ↓ 测序仪输出
FASTQ 文件（原始数据）
    ↓ 比对到参考基因组
BAM 文件（每条读段定位到基因组坐标）
    ↓ 找差异
VCF 文件（和参考基因组不同的位置 = 变异）

第一步：原始数据质量控制

FASTQ 文件长什么样？

拿到数据的第一步，永远是看数据质量。

@SRR001666.1 071112_SLXA-EAS1_s_7:5:1:817:345 length=72
GGGTGATGGCCGCTGCCGATGGCGTCAAATCCCACCAAGTTACCCTTAACAACTTAAGGGTTTTCAAATAGA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIII9IG9ICIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

四行是一条读段：

第 1 行：@ 开头，读段 ID
第 2 行：碱基序列（A/T/C/G）
第 3 行：+ 分隔符
第 4 行：每个碱基的质量分（ASCII 字符，与碱基一一对应）

Phred 质量分的含义：

质量分	错误概率	ASCII 字符
Q10	10%（1/10）	+
Q20	1%（1/100）	5
Q30	0.1%（1/1000）	?
Q40	0.01%（1/10000）	I

Q30 是行业常用的质量门槛，低于 Q20 的碱基通常要过滤掉。

用 FastQC 做质量报告

fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o ./qc_output/

FastQC 会生成 HTML 报告，重点看这几个指标：

检查项	正常状态	异常时意味着
Per base sequence quality	读段末端 Q 值下降是正常的	整体 Q<20 说明测序质量差
Per sequence GC content	接近正态分布	双峰说明有污染
Adapter Content	无或极少	需要 Trim 接头序列
Sequence Duplication Levels	PCR 扩增产物会高	Panel 测序高重复率正常，WGS 异常

质量修剪（Trimming）

如果原始数据质量不好，需要先修剪低质量碱基和接头：

# Trimmomatic：最常用的修剪工具
trimmomatic PE \
    sample_R1.fastq.gz sample_R2.fastq.gz \
    sample_R1_clean.fastq.gz sample_R1_unpaired.fastq.gz \
    sample_R2_clean.fastq.gz sample_R2_unpaired.fastq.gz \
    ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \  # 去除 Illumina 接头
    LEADING:3 \                             # 去除开头质量<3 的碱基
    TRAILING:3 \                            # 去除末尾质量<3 的碱基
    SLIDINGWINDOW:4:15 \                    # 滑动窗口过滤
    MINLEN:36                               # 过滤后长度至少 36bp

为什么是双端（PE = Paired-End）？ 现代测序通常对同一个 DNA 片段从两端分别测序，得到 R1 和 R2 两个文件。双端数据能提供片段长度信息，显著提高比对准确率——尤其是在结构变异检测中。

第二步：比对到参考基因组

参考基因组是什么？

参考基因组是科学界"认定的标准人类 DNA 序列"，由千人基因组等大型项目整理而来。

常用版本：

GRCh38/hg38：当前主流版本（2013 年发布）
GRCh37/hg19：旧版本，很多临床数据库仍用此版本，切换时要注意坐标不通用

下载参考基因组：

# 从 NCBI 下载 hg38
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/\
GCA_000001405.15_GRCh38_assembly_structure/Primary_Assembly/assembled_chromosomes/FASTA/

# 建立 BWA 索引（只需做一次，需要约 1 小时）
bwa index hg38.fa

用 BWA-MEM2 做比对

BWA-MEM2 是短读段比对的行业标准，MEM（Maximal Exact Match）算法适合 70-1000bp 的读段：

bwa mem \
    -t 16 \                        # 使用 16 个线程
    -R "@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA\tLB:lib1" \  # 读组信息（GATK 必须）
    hg38.fa \                      # 参考基因组
    sample_R1_clean.fastq.gz \
    sample_R2_clean.fastq.gz \
    | samtools sort -@ 16 -o sample.sorted.bam  # 直接排序输出 BAM

读组（Read Group）为什么重要？

GATK 在变异检测时用读组区分来自不同测序批次的数据，避免批次效应。@RG 行至少要有：

ID：读组 ID
SM：样本名
PL：测序平台（ILLUMINA/ONT/PACBIO）
LB：文库名

BAM 文件处理

# 建立 BAM 索引（后续工具定位坐标需要）
samtools index sample.sorted.bam

# 查看比对统计
samtools flagstat sample.sorted.bam

flagstat 的关键输出：

500000000 + 0 in total          # 总读段数
498500000 + 0 mapped (99.70%)   # 比对上的比例，WGS 正常应 >95%
250000000 + 0 properly paired   # 正确配对的双端读段

去重（Mark Duplicates）

PCR 扩增会产生多份完全相同的读段副本，它们不携带额外信息，保留会导致变异频率计算偏差：

# GATK MarkDuplicates
gatk MarkDuplicates \
    -I sample.sorted.bam \
    -O sample.markdup.bam \
    -M sample.markdup.metrics.txt

samtools index sample.markdup.bam

理解重复读段：如果两条读段起始位置完全相同，认为它们来自同一个 DNA 分子的 PCR 拷贝，只保留质量最高的一条。WES/Panel 因为捕获步骤的 PCR 放大，重复率通常 20-40%；WGS 通常 <10%。

第三步：变异检测

这是整个流程的核心——找出你的样本和参考基因组的差异在哪里。

什么是变异？

参考基因组：...ATCG[G]TAAGC...
你的样本：  ...ATCG[A]TAAGC...
                  ↑
               这个位置 G→A 就是一个 SNP

主要变异类型：

类型	全称	定义	示例
SNP	Single Nucleotide Polymorphism	单碱基替换	G→A
Indel	Insertion/Deletion	插入或缺失几个碱基	插入 ATG
CNV	Copy Number Variation	某段基因组拷贝数异常	某区域扩增为 4 倍
SV	Structural Variation	大片段结构变异	染色体易位、倒位

GATK HaplotypeCaller：胚系变异检测

胚系变异（Germline Variant）= 生来就有的变异，存在于体内所有细胞，与遗传病相关。

# 第一步：生成每个样本的 GVCF（中间格式，可多样本合并）
gatk HaplotypeCaller \
    -R hg38.fa \
    -I sample.markdup.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF \                    # 输出 GVCF 格式
    --dbsnp dbsnp_138.hg38.vcf.gz  # 已知 SNP 数据库，用于标注

# 第二步（多样本时）：合并 GVCF
gatk CombineGVCFs \
    -R hg38.fa \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    -O cohort.g.vcf.gz

# 第三步：联合基因分型
gatk GenotypeGVCFs \
    -R hg38.fa \
    -V cohort.g.vcf.gz \
    -O cohort.vcf.gz

为什么要用 GVCF 中间格式？

直接从 BAM 出 VCF 只记录"有变异的位置"，GVCF 还记录"没有变异的位置的证据"。这样当你后续加入新样本时，只需重新做第三步，不用重新跑每个样本的全流程。

GATK VQSR：变异质量校正

原始变异中有很多假阳性（噪音），VQSR 用机器学习模型过滤：

# SNP 的 VQSR
gatk VariantRecalibrator \
    -R hg38.fa \
    -V cohort.vcf.gz \
    --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.hg38.vcf.gz \
    --resource:omni,known=false,training=true,truth=true,prior=12.0 1000G_omni2.5.hg38.vcf.gz \
    --resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G_phase1.snps.high_confidence.hg38.vcf.gz \
    --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_138.hg38.vcf.gz \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
    -mode SNP \
    -O cohort.snp.recal \
    --tranches-file cohort.snp.tranches

gatk ApplyVQSR \
    -R hg38.fa \
    -V cohort.vcf.gz \
    -O cohort.snp.filtered.vcf.gz \
    --truth-sensitivity-filter-level 99.0 \
    --tranches-file cohort.snp.tranches \
    --recal-file cohort.snp.recal \
    -mode SNP

VQSR 使用真实已知变异数据库（HapMap、1000 Genomes）作为训练集，给每个变异打一个"像真变异"的概率分。

样本数不够时（< 30）用 Hard Filter 代替 VQSR：

gatk VariantFiltration \
    -V cohort.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "QD2" \
    --filter-expression "FS > 60.0" --filter-name "FS60" \
    --filter-expression "MQ < 40.0" --filter-name "MQ40" \
    -O cohort.hard_filtered.vcf.gz

第四步：变异注释

过滤后的 VCF 只知道"chr1:925952 G→A"，但这个变异有什么生物学意义？需要注释。

ANNOVAR 注释

# 下载注释数据库
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar refGene humandb/
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar gnomad211_genome humandb/
perl annotate_variation.pl -buildver hg38 -downdb -webfrom annovar clinvar_20220320 humandb/

# 运行注释
perl table_annovar.pl \
    cohort.hard_filtered.vcf.gz \
    humandb/ \
    -buildver hg38 \
    -out cohort.annotated \
    -remove \
    -protocol refGene,gnomad211_genome,clinvar_20220320,dbnsfp35a \
    -operation g,f,f,f \
    -nastring . \
    -vcfinput

注释后每个变异会增加：

Gene.refGene：变异在哪个基因
ExonicFunc.refGene：变异类型（同义突变/错义突变/终止密码子...）
gnomAD_genome_ALL：在健康人群中的频率（>1% 通常是良性多态）
CLNSIG：ClinVar 临床意义（Pathogenic/Benign/Uncertain significance）
SIFT/PolyPhen2：预测氨基酸替换对蛋白功能的影响

注释结果解读示例

Chr  Start   Ref  Alt  Gene    ExonicFunc    gnomAD   CLNSIG
chr1 925952  G    A    SAMD11  missense_SNV  0.0001   Pathogenic

这条变异：

位于 SAMD11 基因
是错义突变（氨基酸改变）
在健康人群频率仅 0.01%（罕见）
ClinVar 标注为致病性

WGS vs WES vs Panel 的关键区别

维度	WGS	WES	Panel
覆盖范围	全基因组 100%	外显子 ~1%	特定基因 <0.1%
测序深度	30×	100×	500-2000×
每样本成本	~1500 元	~300 元	~500 元
原始数据量	~100 GB	~10 GB	<1 GB
能发现的变异	SNP/Indel/SV/CNV	SNP/Indel	SNP/Indel（已知热点）
主要用途	科研、罕见病	遗传病诊断	肿瘤精准用药
合规要求	中等	中等	最高（临床报告）

深度（Depth/Coverage）的含义

参考基因组某位置：chr1:100000

读段1:  ----ATCG----
读段2:      ATCG----
读段3:    --ATCG--
读段4:  ----ATCG----
读段5:      ATCG----

该位置被 5 条读段覆盖 → 深度 = 5×

深度越高，越能可靠地识别低频变异。肿瘤 Panel 需要 500-2000× 的深度，才能检测到肿瘤中仅 1% 细胞携带的突变。

肿瘤样本的特殊之处

肿瘤分析和遗传病分析的核心差别：肿瘤变异是体细胞变异（后天获得），不在正常细胞中。

配对分析（Tumor-Normal Pair）

# GATK Mutect2：肿瘤体细胞变异检测
gatk Mutect2 \
    -R hg38.fa \
    -I tumor.markdup.bam -tumor-sample tumor_name \
    -I normal.markdup.bam -normal-sample normal_name \  # 同一患者的正常组织
    --germline-resource af-only-gnomad.hg38.vcf.gz \
    --panel-of-normals pon.vcf.gz \                    # 正常人群背景噪音
    -O tumor_vs_normal.vcf.gz

# 过滤假阳性
gatk FilterMutectCalls \
    -R hg38.fa \
    -V tumor_vs_normal.vcf.gz \
    -O tumor_vs_normal.filtered.vcf.gz

为什么要配对正常样本？

肿瘤样本中的"变异"可能是：

真正的体细胞突变（肿瘤特有）✓
这个人本来就有的胚系变异（正常细胞也有）✗
测序噪音 ✗

正常组织样本（通常取血液）用来区分 1 和 2。

肿瘤特有指标

TMB（Tumor Mutational Burden，肿瘤突变负荷）

定义：每兆碱基的体细胞突变数量
临床意义：TMB 高的肿瘤对免疫检查点抑制剂（PD-1/PD-L1）响应更好
计算：TMB = 过滤后体细胞 SNV 数量 / 外显子区域大小（Mb）

MSI（Microsatellite Instability，微卫星不稳定性）

定义：基因组中短串联重复序列（微卫星）长度的异常变化
临床意义：MSI-H 肿瘤对 PD-1 抑制剂有效（FDA 批准的泛癌种适应症）
检测工具：MSIsensor2、MANTIS

完整流程一览

# ============ 质控 ============
fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o ./qc/
trimmomatic PE sample_R1.fastq.gz sample_R2.fastq.gz \
    sample_R1_clean.fastq.gz /dev/null sample_R2_clean.fastq.gz /dev/null \
    ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

# ============ 比对 ============
bwa mem -t 16 -R "@RG\tID:s1\tSM:sample1\tPL:ILLUMINA\tLB:lib1" \
    hg38.fa sample_R1_clean.fastq.gz sample_R2_clean.fastq.gz \
    | samtools sort -@ 16 -o sample.sorted.bam
samtools index sample.sorted.bam
samtools flagstat sample.sorted.bam  # 检查比对率

# ============ 去重 ============
gatk MarkDuplicates \
    -I sample.sorted.bam -O sample.markdup.bam -M sample.markdup.metrics.txt
samtools index sample.markdup.bam

# ============ 碱基质量校正（BQSR）============
gatk BaseRecalibrator \
    -R hg38.fa -I sample.markdup.bam \
    --known-sites dbsnp_138.hg38.vcf.gz \
    --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
    -O sample.recal.table
gatk ApplyBQSR \
    -R hg38.fa -I sample.markdup.bam \
    --bqsr-recal-file sample.recal.table \
    -O sample.bqsr.bam

# ============ 变异检测 ============
gatk HaplotypeCaller \
    -R hg38.fa -I sample.bqsr.bam \
    -O sample.g.vcf.gz -ERC GVCF \
    --dbsnp dbsnp_138.hg38.vcf.gz
gatk GenotypeGVCFs -R hg38.fa -V sample.g.vcf.gz -O sample.vcf.gz

# ============ 注释 ============
perl table_annovar.pl sample.vcf.gz humandb/ -buildver hg38 \
    -out sample.annotated -remove \
    -protocol refGene,gnomad211_genome,clinvar_20220320,dbnsfp35a \
    -operation g,f,f,f -nastring . -vcfinput

常见问题排查

现象	可能原因	排查方法
比对率 < 90%	物种污染、接头未去除	FastQC 检查接头；BLAST 随机读段确认物种
重复率 > 50%（WGS）	测序量不足/DNA 降解	查看文库插入片段大小分布
变异数量异常少	深度不足、参考基因组版本错误	检查平均覆盖深度；确认 hg19/hg38 版本
VQSR 失败	样本数太少（<30）	改用 Hard Filter
某基因区域覆盖极低	高 GC 区域、重复序列区域	bedtools coverage 检查区域深度

下一步：选择你的细分方向

跑完这条流程之后，你的方向选择决定了接下来的深度：

遗传病方向：学变异解读框架（ACMG 指南五级分类），理解 ClinVar/OMIM 数据库逻辑
肿瘤方向：学 Mutect2 体细胞变异检测，TMB/MSI 计算，CNV 分析（CNVkit）
群体遗传方向：学 GWAS（全基因组关联分析），PLINK 工具集，连锁不平衡
平台工程方向：把上面这条流程用 Nextflow 封装成可重复、可扩展的 pipeline

DNA 测序分析全流程教程

先建立直觉：测序在做什么？​

第一步：原始数据质量控制​

FASTQ 文件长什么样？​

用 FastQC 做质量报告​

质量修剪（Trimming）​

第二步：比对到参考基因组​

参考基因组是什么？​

用 BWA-MEM2 做比对​

BAM 文件处理​

去重（Mark Duplicates）​

第三步：变异检测​

什么是变异？​

GATK HaplotypeCaller：胚系变异检测​

GATK VQSR：变异质量校正​

第四步：变异注释​

ANNOVAR 注释​

注释结果解读示例​

WGS vs WES vs Panel 的关键区别​

深度（Depth/Coverage）的含义​

肿瘤样本的特殊之处​

配对分析（Tumor-Normal Pair）​

肿瘤特有指标​

完整流程一览​

常见问题排查​

下一步：选择你的细分方向​

先建立直觉：测序在做什么？

第一步：原始数据质量控制

FASTQ 文件长什么样？

用 FastQC 做质量报告

质量修剪（Trimming）

第二步：比对到参考基因组

参考基因组是什么？

用 BWA-MEM2 做比对

BAM 文件处理

去重（Mark Duplicates）

第三步：变异检测

什么是变异？

GATK HaplotypeCaller：胚系变异检测

GATK VQSR：变异质量校正

第四步：变异注释

ANNOVAR 注释

注释结果解读示例

WGS vs WES vs Panel 的关键区别

深度（Depth/Coverage）的含义

肿瘤样本的特殊之处

配对分析（Tumor-Normal Pair）

肿瘤特有指标

完整流程一览

常见问题排查

下一步：选择你的细分方向