GATK4 CNV — 基于读段深度的拷贝数变异检测流程¶

# GATK CNV是GATK4内置工具，安装GATK4即可
conda create -n gatk_cnv -c bioconda -c conda-forge gatk4=4.6.0.0
conda activate gatk_cnv

# 验证相关工具可用
gatk CollectReadCounts --help     # 收集读段深度
gatk CreateReadCountPanelOfNormals --help  # 创建PoN
gatk DenoiseReadCounts --help     # 去噪
gatk ModelSegments --help         # 分割建模

# 下载必要资源文件（hg38）
# common germline SNP位点（计算等位基因频率用）
# 来源：GATK资源包 af-only-gnomad.hg38.vcf.gz

# 步骤1：对每个正常样本收集读段计数
# (WES需要区间文件，WGS可用自动生成的全基因组区间)

# 生成区间列表（WES/靶向测序）
gatk PreprocessIntervals \
    -R reference.fasta \
    -L targets.interval_list \    # 捕获区域
    --bin-length 0 \              # WES: bin长度设0（每个外显子为一个bin）
    --padding 250 \               # 两侧各延伸250bp
    -O preprocessed_intervals.interval_list

# 对每个正常样本收集读段深度
for normal_bam in normal1.bam normal2.bam normal3.bam; do
    sample="${normal_bam%.bam}"
    gatk CollectReadCounts \
        -I $normal_bam \
        -L preprocessed_intervals.interval_list \
        --interval-merging-rule OVERLAPPING_ONLY \
        -O ${sample}.counts.hdf5   # 输出HDF5格式（比VCF快）
done

# 步骤2：构建PoN（至少需要10-20个正常样本，越多越好）
gatk CreateReadCountPanelOfNormals \
    --input normal1.counts.hdf5 \
    --input normal2.counts.hdf5 \
    --input normal3.counts.hdf5 \
    --minimum-interval-median-percentile 5.0 \  # 过滤低覆盖bin
    --output cnv_pon.hdf5           # 输出PoN文件

# 步骤1：收集肿瘤样本读段深度
gatk CollectReadCounts \
    -I tumor.bam \
    -L preprocessed_intervals.interval_list \
    --interval-merging-rule OVERLAPPING_ONLY \
    -O tumor.counts.hdf5

# 步骤2：收集等位基因计数（用于分析杂合性丢失LOH）
gatk CollectAllelicCounts \
    -I tumor.bam \
    -R reference.fasta \
    -L gnomad_common_sites.vcf.gz \   # 常见SNP位点
    -O tumor.allelicCounts.tsv

# 步骤3：去噪（对比PoN）
gatk DenoiseReadCounts \
    -I tumor.counts.hdf5 \
    --count-panel-of-normals cnv_pon.hdf5 \
    --standardized-copy-ratios tumor.standardizedCR.tsv \   # 标准化后的拷贝比
    --denoised-copy-ratios tumor.denoisedCR.tsv             # 去噪后的拷贝比

# 步骤4：分割建模
gatk ModelSegments \
    --denoised-copy-ratios tumor.denoisedCR.tsv \
    --allelic-counts tumor.allelicCounts.tsv \    # 等位基因计数（LOH分析）
    --output cnv_results \
    --output-prefix tumor                          # 输出前缀

# 步骤5：调用CNV片段
gatk CallCopyRatioSegments \
    --input cnv_results/tumor.cr.seg \            # ModelSegments输出
    --output tumor.called.seg                      # 最终CNV调用结果

# 完整GATK CNV体细胞流程（WGS，hg38）

REF="hg38.fasta"
# WGS区间：全基因组常染色体（1kb bin）
INTERVALS="hg38.autosomal.interval_list"
PON="cnv_panel_of_normals.hdf5"      # 预建的PoN
GNOMAD_SITES="small_exac_common_3.hg38.vcf.gz"

# 1. 收集肿瘤读段深度
gatk CollectReadCounts \
    -I tumor.bam -L $INTERVALS \
    --interval-merging-rule OVERLAPPING_ONLY \
    -O tumor.counts.hdf5

# 2. 收集等位基因深度（LOH分析）
gatk CollectAllelicCounts \
    -I tumor.bam -R $REF \
    -L $GNOMAD_SITES \
    -O tumor.allelicCounts.tsv

# 3. 去噪
gatk DenoiseReadCounts \
    -I tumor.counts.hdf5 \
    --count-panel-of-normals $PON \
    --standardized-copy-ratios tumor.stdCR.tsv \
    --denoised-copy-ratios tumor.denoisedCR.tsv

# 4. 建模分割
gatk ModelSegments \
    --denoised-copy-ratios tumor.denoisedCR.tsv \
    --allelic-counts tumor.allelicCounts.tsv \
    --output cnv_out --output-prefix tumor

# 5. 调用CNV
gatk CallCopyRatioSegments \
    -I cnv_out/tumor.cr.seg \
    -O tumor.called.seg

# 6. 可视化
gatk PlotDenoisedCopyRatios \
    --standardized-copy-ratios tumor.stdCR.tsv \
    --denoised-copy-ratios tumor.denoisedCR.tsv \
    --sequence-dictionary $REF.dict \
    --output-prefix tumor_cnv \
    --output cnv_plots/

gatk PlotModeledSegments \
    --denoised-copy-ratios tumor.denoisedCR.tsv \
    --allelic-counts cnv_out/tumor.hets.tsv \
    --segments cnv_out/tumor.modelFinal.seg \
    --sequence-dictionary $REF.dict \
    --output-prefix tumor_modeled \
    --output cnv_plots/