GATK4 — 业界标准变异检测流程工具包¶

# 方法1：conda安装（推荐，版本4.6.0.0）
conda create -n gatk4 -c bioconda -c conda-forge \
    gatk4=4.6.0.0 \           # 指定最新稳定版本
    samtools \
    picard \
    bwa
conda activate gatk4

# 验证安装
gatk --version   # 应显示 4.6.0.0

# 方法2：Docker运行（最干净，无依赖冲突）
docker pull broadinstitute/gatk:4.6.0.0
docker run -v $(pwd):/data broadinstitute/gatk:4.6.0.0 \
    gatk --version

# 方法3：下载预编译包
wget https://github.com/broadinstitute/gatk/releases/download/4.6.0.0/gatk-4.6.0.0.zip
unzip gatk-4.6.0.0.zip
./gatk-4.6.0.0/gatk --version

# GATK需要Java 17
java -version   # 确认Java版本

# 第一步：统计已知变异位点的碱基质量偏差
gatk BaseRecalibrator \
    -I aligned.bam \                         # 输入比对BAM
    -R reference.fasta \                     # 参考基因组
    --known-sites dbsnp.vcf.gz \             # 已知SNP数据库
    --known-sites mills_indels.vcf.gz \      # 已知Indel数据库
    -O recal_data.table                      # 输出校正表

# 第二步：应用校正，生成校正后BAM
gatk ApplyBQSR \
    -I aligned.bam \
    -R reference.fasta \
    --bqsr-recal-file recal_data.table \
    -O recal.bam                             # 校正后的BAM文件

# 模式一：直接输出VCF（少量样本）
gatk HaplotypeCaller \
    -R reference.fasta \
    -I recal.bam \
    -O output.vcf.gz \
    --native-pair-hmm-threads 4    # GPU加速的HMM线程数

# 模式二：输出GVCF（多样本联合检测推荐）
gatk HaplotypeCaller \
    -R reference.fasta \
    -I recal.bam \
    -O sample.g.vcf.gz \
    -ERC GVCF                      # 输出每个位置信息的GVCF

# 合并多个GVCF到GenomicsDB数据库
gatk GenomicsDBImport \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    -V sample3.g.vcf.gz \
    --genomicsdb-workspace-path genomicsdb \
    -L chr1 -L chr2                # 按染色体分批运行（推荐）

# 从数据库联合基因分型
gatk GenotypeGVCFs \
    -R reference.fasta \
    -V gendb://genomicsdb \        # 输入GenomicsDB路径
    -O cohort.vcf.gz

# 完整GATK4 GVCF流程（适用于人类基因组hg38）

REF="hg38.fasta"
DBSNP="dbsnp_138.hg38.vcf.gz"
MILLS="Mills_and_1000G_gold_standard.indels.hg38.vcf.gz"

# 1. BQSR预处理
gatk BaseRecalibrator \
    -I sample.bam -R $REF \
    --known-sites $DBSNP --known-sites $MILLS \
    -O recal.table

gatk ApplyBQSR -I sample.bam -R $REF \
    --bqsr-recal-file recal.table -O sample_recal.bam

# 2. HaplotypeCaller生成GVCF
gatk HaplotypeCaller \
    -R $REF -I sample_recal.bam \
    -O sample.g.vcf.gz -ERC GVCF \
    --native-pair-hmm-threads 8

# 3. 变异过滤（硬过滤示例，无足够样本数时替代VQSR）
gatk VariantFiltration \
    -R $REF -V cohort.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "QD2" \
    --filter-expression "MQ < 40.0" --filter-name "MQ40" \
    -O filtered.vcf.gz

# 4. 选取通过过滤的SNP
gatk SelectVariants \
    -V filtered.vcf.gz \
    --select-type-to-include SNP \
    --exclude-filtered \
    -O snps_pass.vcf.gz