跳转至

752. 全基因组测序WGS分析流程

一句话概述:从原始测序数据到变异检测的完整WGS分析管线——覆盖质控、比对、去重、变异检测、过滤和注释全流程,是基因组学分析的"黄金标准"流程。


核心知识点速查表

步骤工具输入 → 输出
质控FastQC+fastpFASTQ → 清洁FASTQ
比对BWA-MEM2FASTQ+参考基因组 → SAM
排序去重samtools+MarkDuplicatesSAM → 去重BAM
BQSRGATK BaseRecalibratorBAM → 校准BAM
变异检测GATK HaplotypeCallerBAM → GVCF
联合分型GATK GenotypeGVCFs多个GVCF → VCF
过滤GATK VQSR/硬过滤原始VCF → 过滤VCF
注释VEP/ANNOVARVCF → 注释VCF

一、流程总览(白话版)

WGS分析像一条"生产线": - 原材料检验(质控) → 加工(比对) → 去瑕疵(去重+BQSR) → 找差异(变异检测) → 质检(过滤) → 贴标签(注释)


二、完整流程代码

# ===== WGS标准分析流程 =====
SAMPLE="sample1"  # 样本名
REF="/data/reference/hg38.fa"  # 参考基因组
DBSNP="/data/reference/dbsnp_146.hg38.vcf.gz"  # dbSNP数据库
KNOWN_INDELS="/data/reference/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz"  # 已知indel
THREADS=16  # 线程数

# ===== 1. 质控 =====
fastp \
  -i ${SAMPLE}_R1.fastq.gz \  # 正向reads
  -I ${SAMPLE}_R2.fastq.gz \  # 反向reads
  -o ${SAMPLE}_clean_R1.fastq.gz \  # 输出清洁reads
  -O ${SAMPLE}_clean_R2.fastq.gz \  # 输出清洁reads
  -j ${SAMPLE}_fastp.json \  # JSON报告
  -h ${SAMPLE}_fastp.html \  # HTML报告
  -w ${THREADS} \  # 线程数
  --detect_adapter_for_pe  # 自动检测接头

# ===== 2. 比对到参考基因组 =====
# 建索引(只需一次)
bwa-mem2 index ${REF}  # 为参考基因组建BWA-MEM2索引

# 比对
bwa-mem2 mem \
  -t ${THREADS} \  # 线程数
  -R "@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tPL:ILLUMINA\tLB:lib1" \  # Read Group信息
  ${REF} \  # 参考基因组
  ${SAMPLE}_clean_R1.fastq.gz \  # 正向reads
  ${SAMPLE}_clean_R2.fastq.gz |  \  # 反向reads
  samtools sort -@ ${THREADS} -o ${SAMPLE}.sorted.bam  # 排序并输出BAM

samtools index ${SAMPLE}.sorted.bam  # 建立BAM索引

# ===== 3. 标记PCR重复 =====
gatk MarkDuplicates \
  -I ${SAMPLE}.sorted.bam \  # 输入BAM
  -O ${SAMPLE}.dedup.bam \  # 输出去重BAM
  -M ${SAMPLE}.dedup_metrics.txt \  # 重复率统计
  --REMOVE_DUPLICATES false  # 标记但不删除(推荐)

samtools index ${SAMPLE}.dedup.bam  # 重新建索引

# ===== 4. 碱基质量重校准(BQSR) =====
gatk BaseRecalibrator \
  -R ${REF} \  # 参考基因组
  -I ${SAMPLE}.dedup.bam \  # 输入BAM
  --known-sites ${DBSNP} \  # 已知SNP位点
  --known-sites ${KNOWN_INDELS} \  # 已知indel位点
  -O ${SAMPLE}.recal_data.table  # 输出重校准表

gatk ApplyBQSR \
  -R ${REF} \  # 参考基因组
  -I ${SAMPLE}.dedup.bam \  # 输入BAM
  --bqsr-recal-file ${SAMPLE}.recal_data.table \  # 重校准表
  -O ${SAMPLE}.recal.bam  # 输出校准后BAM

# ===== 5. 变异检测(HaplotypeCaller) =====
gatk HaplotypeCaller \
  -R ${REF} \  # 参考基因组
  -I ${SAMPLE}.recal.bam \  # 输入校准BAM
  -O ${SAMPLE}.g.vcf.gz \  # 输出GVCF(中间格式)
  -ERC GVCF \  # 输出GVCF模式(方便多样本联合分析)
  --native-pair-hmm-threads ${THREADS}  # HMM线程数

# ===== 6. 变异过滤(硬过滤示例) =====
# SNP过滤
gatk SelectVariants -V ${SAMPLE}.vcf.gz --select-type-to-include SNP -O snps.vcf.gz
gatk VariantFiltration -V snps.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 "SOR > 3.0" --filter-name "HighSOR" \
  -O snps_filtered.vcf.gz

# Indel过滤
gatk SelectVariants -V ${SAMPLE}.vcf.gz --select-type-to-include INDEL -O indels.vcf.gz
gatk VariantFiltration -V indels.vcf.gz \
  --filter-expression "QD < 2.0" --filter-name "LowQD" \
  --filter-expression "FS > 200.0" --filter-name "HighFS" \
  -O indels_filtered.vcf.gz

三、常见报错与解决

报错信息原因解决方案
BWA: missing index参考基因组未建索引bwa-mem2 index ref.fa
GATK: read group errorBAM缺少RG信息比对时加-R参数
BQSR: no known sites已知位点文件路径错误检查dbSNP文件路径
HaplotypeCaller: OOM内存不足增加Java堆内存--java-options "-Xmx32g"
VCF: no PASS variants过滤太严格放宽过滤阈值

四、面试高频问题

Q1: BQSR的原理和目的?

A: 测序仪报告的碱基质量值有系统性偏差(特别是在reads末端和某些序列context中)。BQSR用已知的SNP位点作为"标准答案",对质量值进行经验校正,提高变异检测准确性。

Q2: 为什么要用GVCF模式?

A: GVCF不仅记录变异位点,还记录"没有变异"的位点的置信度。这使得后续多样本联合分型时不需要重新处理BAM文件,效率高且可以增量添加新样本。


五、速查表

# WGS分析核心命令
fastp -i R1.fq.gz -I R2.fq.gz -o clean_R1.fq.gz -O clean_R2.fq.gz
bwa-mem2 mem -t 16 -R "@RG..." ref.fa R1.fq.gz R2.fq.gz | samtools sort -o sorted.bam
gatk MarkDuplicates -I sorted.bam -O dedup.bam -M metrics.txt
gatk BaseRecalibrator -R ref.fa -I dedup.bam --known-sites dbsnp.vcf -O recal.table
gatk ApplyBQSR -R ref.fa -I dedup.bam --bqsr-recal-file recal.table -O recal.bam
gatk HaplotypeCaller -R ref.fa -I recal.bam -O sample.g.vcf.gz -ERC GVCF

# 流程管理推荐: nf-core/sarek