变异检测对比：GATK/DeepVariant/Strelka2¶

# === GATK4 ===
conda install -c bioconda gatk4

# === DeepVariant ===
# 推荐Docker方式
docker pull google/deepvariant:1.6.1

# 或Singularity（HPC常用）
singularity pull deepvariant.sif docker://google/deepvariant:1.6.1

# === Strelka2 ===
conda install -c bioconda strelka

# === GATK最佳实践流程 ===
REF=ref.fa
BAM=aligned.bam

# 第1步：标记重复
gatk MarkDuplicates \
  -I $BAM \                                # 输入BAM
  -O marked.bam \                          # 输出BAM
  -M dup_metrics.txt                       # 重复统计

# 第2步：碱基质量重校准(BQSR)
gatk BaseRecalibrator \
  -R $REF \                                # 参考基因组
  -I marked.bam \                          # 输入BAM
  --known-sites dbsnp.vcf.gz \            # 已知变异位点
  --known-sites Mills_indels.vcf.gz \     # 已知InDel
  -O recal_table                           # 校准表

gatk ApplyBQSR \
  -R $REF \
  -I marked.bam \
  --bqsr-recal-file recal_table \          # 校准表
  -O recal.bam                             # 校准后的BAM

# 第3步：变异检测
gatk HaplotypeCaller \
  -R $REF \                                # 参考基因组
  -I recal.bam \                           # 校准后BAM
  -O sample.g.vcf.gz \                    # 输出GVCF
  -ERC GVCF                                # GVCF模式（推荐）

# 第4步：联合分型（多样本）
gatk CombineGVCFs \
  -R $REF \
  -V sample1.g.vcf.gz \                   # 样本1的GVCF
  -V sample2.g.vcf.gz \                   # 样本2的GVCF
  -O combined.g.vcf.gz

gatk GenotypeGVCFs \
  -R $REF \
  -V combined.g.vcf.gz \
  -O genotyped.vcf.gz                      # 联合分型结果

# 第5步：变异过滤(VQSR或硬过滤)
# 硬过滤（样本少时）
gatk SelectVariants -V genotyped.vcf.gz --select-type-to-include SNP -O snps.vcf.gz
gatk VariantFiltration \
  -V snps.vcf.gz \
  --filter-expression "QD<2.0 || FS>60.0 || MQ<40.0 || MQRankSum<-12.5 || ReadPosRankSum<-8.0" \
  --filter-name "FILTER" \
  -O filtered_snps.vcf.gz

# === DeepVariant（Docker方式，简单高效） ===
docker run \
  -v "$(pwd):/input" \                     # 挂载当前目录
  -v "$(pwd)/output:/output" \             # 输出目录
  google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \                       # 模型类型: WGS/WES/PACBIO/ONT
  --ref=/input/ref.fa \                    # 参考基因组
  --reads=/input/aligned.bam \             # 输入BAM
  --output_vcf=/output/deepvariant.vcf.gz \  # 输出VCF
  --output_gvcf=/output/deepvariant.g.vcf.gz \  # 输出GVCF
  --num_shards=8                            # 并行数（≈线程数）

# 就这一步！不需要BQSR、不需要硬过滤
# DeepVariant内部CNN已经学会了这些校正

# === Strelka2 胚系变异检测 ===
# 配置
configureStrelkaGermlineWorkflow.py \
  --bam aligned.bam \                      # 输入BAM
  --referenceFasta ref.fa \               # 参考基因组
  --runDir strelka_germline                # 运行目录

# 运行
strelka_germline/runWorkflow.py \
  -m local \                               # 本地运行
  -j 8                                     # 线程数

# 输出:
# strelka_germline/results/variants/variants.vcf.gz（所有变异）
# strelka_germline/results/variants/genome.S1.vcf.gz（基因组级VCF）

# === GATK Mutect2（体细胞变异，肿瘤分析必用） ===
gatk Mutect2 \
  -R ref.fa \
  -I tumor.bam \                           # 肿瘤样本
  -I normal.bam \                          # 配对正常样本
  -normal normal_sample_name \             # 正常样本名
  --germline-resource gnomad.vcf.gz \     # 人群变异数据库
  --panel-of-normals pon.vcf.gz \         # 正常样本面板
  -O mutect2.vcf.gz

# 过滤
gatk FilterMutectCalls \
  -V mutect2.vcf.gz \
  -R ref.fa \
  -O filtered_mutect2.vcf.gz

# === Strelka2 体细胞变异 ===
configureStrelkaSomaticWorkflow.py \
  --normalBam normal.bam \                 # 正常样本
  --tumorBam tumor.bam \                   # 肿瘤样本
  --referenceFasta ref.fa \
  --runDir strelka_somatic

strelka_somatic/runWorkflow.py -m local -j 8
# 输出: somatic.snvs.vcf.gz + somatic.indels.vcf.gz

场景是什么？
├── 胚系变异 + 少量样本 → DeepVariant（最简单最准）
├── 胚系变异 + 大队列 → GATK（联合分型优势）
├── 体细胞变异 → GATK Mutect2（金标准）
├── 快速初筛 → Strelka2（最快）
└── 临床级分析 → GATK + DeepVariant 双验证

# === 变异检测速查 ===
# GATK（标准流程）
gatk MarkDuplicates -I in.bam -O dedup.bam -M metrics
gatk HaplotypeCaller -R ref.fa -I dedup.bam -O out.g.vcf -ERC GVCF

# DeepVariant（一步到位）
run_deepvariant --model_type=WGS --ref=ref.fa --reads=in.bam --output_vcf=out.vcf

# Strelka2
configureStrelkaGermlineWorkflow.py --bam in.bam --referenceFasta ref.fa --runDir out/
out/runWorkflow.py -m local -j 8

# 体细胞: gatk Mutect2 -I tumor.bam -I normal.bam -normal NAME -O somatic.vcf
# 选择: 简单准确→DeepVariant | 大队列→GATK | 速度→Strelka2 | 肿瘤→Mutect2