变异检测流程（GATK Best Practices）¶

BAM（比对后） → MarkDuplicates → BQSR → HaplotypeCaller → GVCF
                                                              ↓
                                          GenomicsDBImport（合并多样本）
                                                              ↓
                                          GenotypeGVCFs（联合基因分型）
                                                              ↓
                                          VQSR / 硬过滤 → 最终 VCF

# 如果 BAM 未排序，先排序（通常 BWA 输出经过 samtools sort 已排好）
# -@ 8 表示用 8 个线程加速
# -o 输出文件名
samtools sort -@ 8 -o sample_sorted.bam sample.bam

# 建索引（后续所有工具都需要 .bai 索引文件）
samtools index -@ 8 sample_sorted.bam

# GATK4 MarkDuplicates（已整合原 Picard 功能）
# --INPUT           输入的排序后 BAM
# --OUTPUT          输出标记重复后的 BAM
# --METRICS_FILE    重复率统计报告
# --REMOVE_DUPLICATES false  标记但不删除（保留信息）
# --CREATE_INDEX true  自动创建索引

gatk MarkDuplicates \
    --INPUT sample_sorted.bam \
    --OUTPUT sample_markdup.bam \
    --METRICS_FILE sample_dup_metrics.txt \
    --REMOVE_DUPLICATES false \
    --CREATE_INDEX true

# 正常重复率：WGS 10-30%，WES 10-40%
# 如果 >50%，说明文库复杂度低，数据质量有问题
grep "PERCENT_DUPLICATION" sample_dup_metrics.txt

# 第一步：建立重校准模型（BaseRecalibrator）
# -R                参考基因组 FASTA
# -I                输入 BAM
# --known-sites     已知变异位点（用来排除真实变异，只看测序误差）
#                   需要下载 Broad 提供的资源包（dbSNP + Mills + 1000G）
# -O                输出重校准表

gatk BaseRecalibrator \
    -R ref/hg38.fa \
    -I sample_markdup.bam \
    --known-sites ref/dbsnp_146.hg38.vcf.gz \
    --known-sites ref/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
    --known-sites ref/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
    -O sample_recal_data.table

# 第二步：应用重校准（ApplyBQSR）
# 把校准表应用到 BAM 文件，生成校准后的 BAM
gatk ApplyBQSR \
    -R ref/hg38.fa \
    -I sample_markdup.bam \
    --bqsr-recal-file sample_recal_data.table \
    -O sample_bqsr.bam

# 再跑一次 BaseRecalibrator（用校准后的 BAM）
gatk BaseRecalibrator \
    -R ref/hg38.fa \
    -I sample_bqsr.bam \
    --known-sites ref/dbsnp_146.hg38.vcf.gz \
    --known-sites ref/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
    -O sample_recal_after.table

# 生成对比报告
gatk AnalyzeCovariates \
    -before sample_recal_data.table \
    -after sample_recal_after.table \
    -plots recalibration_plots.pdf

# 单样本模式（输出 GVCF 格式，便于后续多样本联合分析）
# -R            参考基因组
# -I            校准后的 BAM
# -O            输出 GVCF（含变异位点 + 非变异位点的置信度信息）
# -ERC GVCF    启用 GVCF 模式（面试重点：为什么用 GVCF？）
# --native-pair-hmm-threads 4  Pair-HMM 计算线程数

gatk HaplotypeCaller \
    -R ref/hg38.fa \
    -I sample_bqsr.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF \
    --native-pair-hmm-threads 4

# 如果只分析特定区域（如外显子组），加 -L 参数
# gatk HaplotypeCaller ... -L exome_targets.bed

# GenomicsDBImport：把多个样本的 GVCF 导入一个 GenomicsDB 数据库
# --genomicsdb-workspace-path  数据库存储路径
# -V                          每个样本的 GVCF 文件
# -L                          分析区域（按染色体拆分可并行加速）

gatk GenomicsDBImport \
    --genomicsdb-workspace-path cohort_db_chr1 \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    -V sample3.g.vcf.gz \
    -L chr1 \
    --reader-threads 4

# 样本多时用 sample map 文件代替多个 -V 参数
# sample_map.txt 格式：样本名\t文件路径
# gatk GenomicsDBImport --sample-name-map sample_map.txt ...

# 对 GenomicsDB 进行联合基因分型
# -V gendb://  指向 GenomicsDBImport 创建的数据库
# -O           输出多样本联合 VCF

gatk GenotypeGVCFs \
    -R ref/hg38.fa \
    -V gendb://cohort_db_chr1 \
    -O cohort_chr1_raw.vcf.gz

# 第一步：对 SNP 建立 VQSR 模型
# -V                    原始 VCF
# --resource            训练资源（truth=true 表示高置信真阳性集）
# -an                   用于训练的注释特征
# --mode SNP            只处理 SNP
# -O                    输出重校准文件
# --tranches-file       输出灵敏度梯度文件

gatk VariantRecalibrator \
    -R ref/hg38.fa \
    -V cohort_raw.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=false,prior=12.0 1000G_omni2.5.hg38.vcf.gz \
    --resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G_phase1.snps.hg38.vcf.gz \
    --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_146.hg38.vcf.gz \
    -an QD -an MQRankSum -an ReadPosRankSum -an FS -an MQ -an SOR \
    --mode SNP \
    -O cohort_snp.recal \
    --tranches-file cohort_snp.tranches

# 第二步：应用 VQSR 过滤
# --truth-sensitivity-filter-level 99.5  保留 99.5% 的真阳性
gatk ApplyVQSR \
    -R ref/hg38.fa \
    -V cohort_raw.vcf.gz \
    -O cohort_snp_filtered.vcf.gz \
    --recal-file cohort_snp.recal \
    --tranches-file cohort_snp.tranches \
    --truth-sensitivity-filter-level 99.5 \
    --mode SNP

# 第三步：对 InDel 重复上述步骤（使用不同参数和资源）
gatk VariantRecalibrator \
    -R ref/hg38.fa \
    -V cohort_snp_filtered.vcf.gz \
    --resource:mills,known=false,training=true,truth=true,prior=12.0 Mills_and_1000G.hg38.vcf.gz \
    --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_146.hg38.vcf.gz \
    -an QD -an ReadPosRankSum -an FS -an SOR \
    --mode INDEL \
    --max-gaussians 4 \
    -O cohort_indel.recal \
    --tranches-file cohort_indel.tranches

gatk ApplyVQSR \
    -R ref/hg38.fa \
    -V cohort_snp_filtered.vcf.gz \
    -O cohort_final_filtered.vcf.gz \
    --recal-file cohort_indel.recal \
    --tranches-file cohort_indel.tranches \
    --truth-sensitivity-filter-level 99.0 \
    --mode INDEL

# 先把 SNP 和 InDel 分开
gatk SelectVariants \
    -V cohort_raw.vcf.gz \
    --select-type-to-include SNP \
    -O cohort_raw_snps.vcf.gz

gatk SelectVariants \
    -V cohort_raw.vcf.gz \
    --select-type-to-include INDEL \
    -O cohort_raw_indels.vcf.gz

# SNP 硬过滤（GATK 推荐阈值）
# QD < 2.0        质量除以深度，过低说明支持变异的证据不够
# FS > 60.0       Fisher链偏差，过高说明变异只出现在一条链上
# MQ < 40.0       比对质量均值，过低说明reads比对不可靠
# MQRankSum < -12.5  参考/变异allele比对质量差异
# ReadPosRankSum < -8.0  变异出现在read的位置偏向末端（测序误差高）
# SOR > 3.0       链偏差的另一个指标

gatk VariantFiltration \
    -V cohort_raw_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 "MQRankSum < -12.5" --filter-name "LowMQRankSum" \
    --filter-expression "ReadPosRankSum < -8.0" --filter-name "LowReadPosRankSum" \
    --filter-expression "SOR > 3.0" --filter-name "HighSOR" \
    -O cohort_filtered_snps.vcf.gz

# InDel 硬过滤（阈值和 SNP 不同）
gatk VariantFiltration \
    -V cohort_raw_indels.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 "LowReadPosRankSum" \
    --filter-expression "SOR > 10.0" --filter-name "HighSOR" \
    -O cohort_filtered_indels.vcf.gz

# 第一步：生成 pileup 并计算基因型似然值
# mpileup 逐位点统计每个碱基的支持reads数
# -f         参考基因组
# -q 20      最低比对质量（过滤比对差的reads）
# -Q 20      最低碱基质量
# --threads 8  多线程
# -a "FORMAT/AD,FORMAT/DP"  输出等位基因深度和总深度

bcftools mpileup \
    -f ref/hg38.fa \
    -q 20 -Q 20 \
    --threads 8 \
    -a "FORMAT/AD,FORMAT/DP,FORMAT/ADF,FORMAT/ADR" \
    sample_bqsr.bam | \
# 第二步：变异调用
# call 基于似然值判断基因型
# -m         使用多等位基因调用模型（推荐，替代旧的 -c 模型）
# -v         只输出变异位点（不输出非变异位点）
# --ploidy 2 二倍体
# -Oz        输出 gzip 压缩的 VCF
bcftools call \
    -m -v \
    --ploidy 2 \
    --threads 8 \
    -Oz -o sample_bcftools_raw.vcf.gz

# 第三步：基础过滤
bcftools filter \
    -i 'QUAL>=30 && DP>=10 && DP<=200 && MQ>=40' \
    -Oz -o sample_bcftools_filtered.vcf.gz \
    sample_bcftools_raw.vcf.gz

# 建索引
bcftools index sample_bcftools_filtered.vcf.gz

# 拉取 DeepVariant Docker 镜像
docker pull google/deepvariant:1.10.0

# 运行 DeepVariant（一条命令完成全部流程）
# --model_type  WGS/WES/PACBIO/ONT（根据数据类型选择）
# --ref         参考基因组
# --reads       输入 BAM（需经过排序和建索引）
# --output_vcf  输出 VCF
# --output_gvcf 输出 GVCF（可选）
# --num_shards  并行数（=CPU核心数）

docker run -v "$(pwd):/input" -v "$(pwd)/output:/output" \
    google/deepvariant:1.10.0 \
    /opt/deepvariant/bin/run_deepvariant \
    --model_type=WGS \
    --ref=/input/ref/hg38.fa \
    --reads=/input/sample_bqsr.bam \
    --output_vcf=/output/sample_deepvariant.vcf.gz \
    --output_gvcf=/output/sample_deepvariant.g.vcf.gz \
    --num_shards=16

# 综合过滤：质量 + 深度 + 链偏差 + 缺失率
bcftools filter \
    -i 'QUAL>=30 && INFO/DP>=10 && INFO/DP<=500 && MQ>=40' \
    -Oz -o step1.vcf.gz input.vcf.gz

# 按样本级别过滤（GQ 和 DP）
bcftools filter \
    -i 'FMT/GQ>=20 && FMT/DP>=8' \
    -S . \
    -Oz -o step2.vcf.gz step1.vcf.gz

# 过滤缺失率（>20% 样本缺失的位点去掉）
bcftools view \
    -i 'F_MISSING < 0.2' \
    -Oz -o final.vcf.gz step2.vcf.gz

# 统计过滤前后变异数量
bcftools stats input.vcf.gz | grep "^SN"
bcftools stats final.vcf.gz | grep "^SN"

# 第一步：VCF 转 ANNOVAR 输入格式
# convert2annovar.pl 是 ANNOVAR 的格式转换脚本
convert2annovar.pl -format vcf4 input.vcf > input.avinput

# 第二步：基因注释 + 数据库注释
# table_annovar.pl 一条命令完成多数据库注释
# -buildver hg38          参考基因组版本
# -protocol              注释数据库列表
#   refGene              基因区域注释（外显子/内含子/UTR）
#   gnomad41_genome      人群频率（gnomAD v4.1）
#   clinvar_20240917     临床意义（致病/良性/意义不明）
#   dbnsfp47a            功能预测分数（SIFT/PolyPhen2/CADD等）
# -operation            每个数据库对应的注释类型（g=基因，f=过滤）

table_annovar.pl input.avinput humandb/ \
    -buildver hg38 \
    -out sample_annotated \
    -protocol refGene,gnomad41_genome,clinvar_20240917,dbnsfp47a \
    -operation g,f,f,f \
    -nastring . \
    -csvout

# 一条命令完成注释（速度最快）
# GRCh38.p14 是预构建的人类基因组数据库
# -v         详细输出
# -stats     生成 HTML 统计报告

java -Xmx8g -jar snpEff.jar \
    -v GRCh38.p14 \
    -stats snpeff_report.html \
    input.vcf.gz > annotated.vcf

# 常见输出字段解读：
# missense_variant     错义突变（氨基酸改变）
# synonymous_variant   同义突变（氨基酸不变）
# frameshift_variant   移码突变（InDel导致阅读框改变，通常有害）
# stop_gained          无义突变（提前出现终止密码子）
# splice_region_variant 剪接区域变异

# Ensembl VEP（功能最全面）
# --cache           使用本地缓存（需提前下载）
# --assembly GRCh38 指定基因组版本
# --everything      输出所有可用注释
# --vcf             输出 VCF 格式
# --plugin CADD     加载 CADD 评分插件（预测有害性）

vep --input_file input.vcf.gz \
    --output_file annotated_vep.vcf \
    --cache \
    --assembly GRCh38 \
    --everything \
    --vcf \
    --plugin CADD,whole_genome_SNVs.tsv.gz \
    --fork 8