758. 多样本联合基因分型¶

单独分型的问题：
  样本A在chr1:12345位置，测序深度只有3x
  → 单独看，不确定这里是纯合参考(0/0)还是杂合(0/1)

联合分型的优势：
  样本A: 3x → 不确定
  样本B: 30x → 明确是杂合(0/1)
  样本C: 20x → 明确是杂合(0/1)
  ...100个样本中有50个在这个位置有变异
  → 集体信息帮助确认：样本A大概率也是杂合(0/1)

关键优势：
  ① 低深度位点借助群体信息提高准确性
  ② 区分"真的没变异"和"没测到"
  ③ 获得统一的多样本VCF，方便下游GWAS等分析

传统联合检测（已淘汰）：
  100个样本一起跑HaplotypeCaller → 极慢
  新来1个样本 → 101个全部重跑 → 这就是N+1问题

GVCF方案（现在的标准）：
  Step 1: 每个样本独立生成GVCF（一次性，可并行）
  Step 2: 合并所有GVCF（GenomicsDBImport）
  Step 3: 联合分型（GenotypeGVCFs）

  新来1个样本？
  → 只需Step 1生成它的GVCF
  → Step 2增量更新GenomicsDB
  → Step 3重新联合分型
  → 不需要重跑老样本的Step 1！

# ===== 每个样本独立生成GVCF =====
# 前提：已完成BWA比对 + MarkDuplicates + BQSR

# 对每个样本运行HaplotypeCaller的GVCF模式
gatk HaplotypeCaller \
  -R reference.fa \                     # 参考基因组
  -I sample1_recal.bam \               # 校准后的BAM文件
  -O sample1.g.vcf.gz \                # 输出GVCF文件（.g.vcf.gz）
  -ERC GVCF \                          # 关键：GVCF模式
  --native-pair-hmm-threads 4 \        # HMM线程数
  -L intervals.bed                      # 目标区间（WES必需，WGS可选）

# GVCF vs VCF的区别：
# VCF只记录有变异的位点
# GVCF记录所有位点（包括参考纯合位点），用"块"压缩参考区域
# 这样联合分型时才能区分"参考纯合"和"没有数据"

# 批量处理（bash循环）
for sample in sample1 sample2 sample3 sample4 sample5; do
    gatk HaplotypeCaller \
      -R reference.fa \
      -I ${sample}_recal.bam \          # 每个样本的BAM
      -O ${sample}.g.vcf.gz \           # 每个样本的GVCF
      -ERC GVCF \
      --native-pair-hmm-threads 4 &     # 后台并行运行
done
wait  # 等待所有任务完成
echo "所有GVCF生成完毕"

# ===== GenomicsDBImport：合并多个GVCF =====

# 方法一：直接指定每个GVCF文件
gatk GenomicsDBImport \
  -V sample1.g.vcf.gz \                # 样本1的GVCF
  -V sample2.g.vcf.gz \                # 样本2的GVCF
  -V sample3.g.vcf.gz \                # 样本3的GVCF
  --genomicsdb-workspace-path my_database \  # 输出GenomicsDB目录
  -L intervals.list \                   # 区间列表（必需！）
  --reader-threads 4 \                  # 读取线程数
  --batch-size 50                       # 每批处理的样本数

# 方法二：用sample map文件（推荐，样本多时更方便）
# 先创建sample_map.txt（制表符分隔）
# 格式：样本名\tGVCF路径
cat > sample_map.txt << 'EOF'
sample1 /path/to/sample1.g.vcf.gz
sample2 /path/to/sample2.g.vcf.gz
sample3 /path/to/sample3.g.vcf.gz
sample4 /path/to/sample4.g.vcf.gz
sample5 /path/to/sample5.g.vcf.gz
EOF

gatk GenomicsDBImport \
  --sample-name-map sample_map.txt \    # 样本映射文件
  --genomicsdb-workspace-path my_database \  # 输出目录
  -L chr1 -L chr2 -L chr3 \            # 按染色体指定区间
  --reader-threads 4 \                  # 读取线程
  --merge-input-intervals \             # 合并小区间提高效率
  --bypass-feature-reader \             # GATK4.2.3+加速选项(10-15%提速)
  --tmp-dir /tmp/gatk_tmp              # 临时目录（需要大空间）

# ===== 增量更新：添加新样本 =====
# 已有GenomicsDB后，新样本直接追加
gatk GenomicsDBImport \
  --genomicsdb-update-workspace-path my_database \  # 更新已有数据库
  -V new_sample.g.vcf.gz \             # 新样本的GVCF
  -L intervals.list                     # 区间（需与原始一致）

# ===== GenotypeGVCFs：联合分型 =====
gatk GenotypeGVCFs \
  -R reference.fa \                     # 参考基因组
  -V gendb://my_database \             # 输入GenomicsDB（注意gendb://前缀）
  -O joint_called.vcf.gz \             # 输出联合分型VCF
  --tmp-dir /tmp/gatk_tmp \            # 临时目录
  -G StandardAnnotation \              # 标准注释
  -G AS_StandardAnnotation              # 等位基因特异性注释

# 按染色体并行加速
for chr in chr{1..22} chrX chrY; do
    gatk GenotypeGVCFs \
      -R reference.fa \
      -V gendb://my_database \
      -O joint_${chr}.vcf.gz \          # 每条染色体一个VCF
      -L ${chr} &                       # 后台并行
done
wait

# 合并所有染色体的结果
gatk GatherVcfs \
  $(for chr in chr{1..22} chrX chrY; do echo "-I joint_${chr}.vcf.gz"; done) \
  -O joint_called_all.vcf.gz            # 最终合并VCF

# ===== 过滤策略一：VQSR（推荐，需要≥30个样本）=====
# SNP的VQSR
gatk VariantRecalibrator \
  -R reference.fa \
  -V joint_called.vcf.gz \             # 联合分型VCF
  --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap.vcf.gz \
  --resource:omni,known=false,training=true,truth=true,prior=12.0 omni.vcf.gz \
  --resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G.vcf.gz \
  --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp.vcf.gz \
  -an QD -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an MQ \  # 注释特征
  -mode SNP \                           # SNP模式
  -O snp_recal.recal \                  # 输出重校准文件
  --tranches-file snp_recal.tranches    # 分层文件

gatk ApplyVQSR \
  -R reference.fa \
  -V joint_called.vcf.gz \
  -O joint_snp_filtered.vcf.gz \        # 输出过滤后VCF
  --recal-file snp_recal.recal \
  --tranches-file snp_recal.tranches \
  --truth-sensitivity-filter-level 99.5 \  # 敏感性阈值
  -mode SNP                             # SNP模式

# ===== 过滤策略二：硬过滤（样本少时用）=====
# 提取SNP
gatk SelectVariants -V joint_called.vcf.gz -select-type SNP -O raw_snps.vcf.gz

# SNP硬过滤
gatk VariantFiltration \
  -R reference.fa \
  -V raw_snps.vcf.gz \
  -O filtered_snps.vcf.gz \
  --filter-expression "QD < 2.0" --filter-name "LowQD" \          # 质量/深度
  --filter-expression "FS > 60.0" --filter-name "HighFS" \        # 链偏好
  --filter-expression "SOR > 3.0" --filter-name "HighSOR" \       # 链偏好比
  --filter-expression "MQ < 40.0" --filter-name "LowMQ" \         # 比对质量
  --filter-expression "MQRankSum < -12.5" --filter-name "LowMQRS" \
  --filter-expression "ReadPosRankSum < -8.0" --filter-name "LowRPRS"

# 提取INDEL
gatk SelectVariants -V joint_called.vcf.gz -select-type INDEL -O raw_indels.vcf.gz

# INDEL硬过滤
gatk VariantFiltration \
  -R reference.fa \
  -V raw_indels.vcf.gz \
  -O filtered_indels.vcf.gz \
  --filter-expression "QD < 2.0" --filter-name "LowQD" \
  --filter-expression "FS > 200.0" --filter-name "HighFS" \
  --filter-expression "SOR > 10.0" --filter-name "HighSOR" \
  --filter-expression "ReadPosRankSum < -20.0" --filter-name "LowRPRS"

# ===== CombineGVCFs：小队列的备选 =====
# 当GenomicsDBImport有问题时使用（不推荐大队列）
gatk CombineGVCFs \
  -R reference.fa \
  -V sample1.g.vcf.gz \                # GVCF文件1
  -V sample2.g.vcf.gz \                # GVCF文件2
  -V sample3.g.vcf.gz \                # GVCF文件3
  -O combined.g.vcf.gz                  # 输出合并GVCF

# 然后联合分型
gatk GenotypeGVCFs \
  -R reference.fa \
  -V combined.g.vcf.gz \               # 注意：这里不用gendb://前缀
  -O joint_called.vcf.gz

# GenomicsDBImport vs CombineGVCFs：
# GenomicsDBImport：高效，支持增量更新，推荐
# CombineGVCFs：简单，但大队列时慢且耗内存

# ===== 内存管理 =====
# GenomicsDBImport使用C/C++底层库
# Java内存不要设满，留20%给C库
gatk --java-options "-Xmx32g -Xms32g" GenomicsDBImport \
  --sample-name-map sample_map.txt \
  --genomicsdb-workspace-path my_database \
  -L intervals.list \
  --reader-threads 5 \                  # 读取线程（不要太多）
  --batch-size 50 \                     # 每批50个样本
  --tmp-dir /scratch/tmp                # 临时目录放在快速存储上

# ===== 分区并行策略 =====
# 将基因组分成多个区间，每个区间独立处理
# 最后合并结果
# scatter-gather模式
for i in $(seq 1 24); do
    gatk GenomicsDBImport \
      --sample-name-map sample_map.txt \
      --genomicsdb-workspace-path db_shard_${i} \
      -L shard_${i}.intervals \
      --reader-threads 2 &
done
wait

# ===== 联合分型速查 =====

# Step 1: 生成GVCF（每个样本独立）
gatk HaplotypeCaller -R ref.fa -I sample.bam \
    -O sample.g.vcf.gz -ERC GVCF

# Step 2: 合并GVCF
gatk GenomicsDBImport \
    --sample-name-map sample_map.txt \
    --genomicsdb-workspace-path my_db -L intervals.list

# Step 3: 联合分型
gatk GenotypeGVCFs -R ref.fa \
    -V gendb://my_db -O joint.vcf.gz

# Step 4: 过滤
# 大队列(≥30): VQSR
# 小队列(<30): 硬过滤

# 增量更新（新样本加入）
gatk GenomicsDBImport \
    --genomicsdb-update-workspace-path my_db \
    -V new_sample.g.vcf.gz -L intervals.list

# 关键参数：
# -ERC GVCF → 生成GVCF（必须！）
# gendb:// → GenotypeGVCFs的输入前缀
# --bypass-feature-reader → GATK4.2.3+加速选项
# --batch-size 50 → 控制内存使用