GATK 的 Germline mutation 流程--肿瘤基因组测序数据分析专栏

生信菜鸟团

发布于 2021-11-15 17:32:52

3.8K0

发布于 2021-11-15 17:32:52

文章被收录于专栏：生信菜鸟团

简介

基因组测序，最重要的就是检测变异位点，对于家系数据、遗传病研究，更多的是关心 Germline mutation 生殖突变。当然，部分肿瘤研究也会关注 Germline mutation。GATK 对这类变异的检测有一整套流程，主要用到的工具是：HaplotypeCaller 、GenomicsDBImport、GenotypeGVCFs、VariantRecalibrator、 ApplyVQSR 等工具

流程图

Germline mutation 分析，对样本没有太多的要求，肿瘤非配对样本也可以分析。不过方法上有两种，单个样本和多个样本（队列）略有不同。流程图是：对于多样本或队列样本：

对于单个样本：

在这里，仅介绍多样本或者队列样本的 GATK Germline mutation 流程，基于 BQSR 得到的 BAM 进行分析，主要有以下几步：

HaplotypeCaller

第一步先对每个样本 call 突变，用到了 HaplotypeCaller ，而且是在 GVCF 模式下，代码是：

${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" HaplotypeCaller  \
      -R ${GENOME}             \
      -ERC GVCF                \
      -I 5.gatk/${id}_bqsr.bam \
      -O 6.gvcf/${id}.g.vcf.gz \
      --intervals ${bed}       \
      1>./6.gvcf/${id}_HC.log 2>&1

GenomicsDBImport

第二步是对第一步结果的整合，如果有50个样本，就有 50 个 *g.vcf.gz 文件，最后得到一个类似数据库的文件夹。代码是：

${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" GenomicsDBImport \
        -R ${GENOME}                                      \
        $(ls 6.gvcf/*g.vcf.gz | awk '{print "-V "$0" "}') \
        -L ${bed}                                         \
        --merge-input-intervals TRUE                      \
        --genomicsdb-workspace-path ./6.gvcf/gvcfs_db     \
        1>./6.gvcf/GenomicsDBImport.log 2>&1

需要注意的是，对于外显子数据，除了指定 bed 文件，还要加上一个参数 --merge-input-intervals TRUE 。如果不加，对于每一个 bed 文件上的坐标（即bed文件的每一行），程序就会循环一次，并在 ./6.gvcf/gvcfs_db 文件夹中生成一个子文件夹，如果 bed 文件有 20W 行，就会有 20W 个文件夹，每个文件夹大小~100M，这个数据量是非常恐怖的。不仅如此，运行时间也大大增加。而加上参数 --merge-input-intervals TRUE 后，程序会对 bed 文件中的坐标进行整合，同一条染色体会整合到一起运行，并将结果保存到同一个文件夹中。

GenotypeGVCFs

第三步成为联合基因分型，按照官方文档所述：在这一步，收集所有每个样本的 GVCF 结果（已经在上一步整合为数据库）并将它们一起传递给联合基因分型工具 GenotypeGVCF。这会产生一组联合调用的 SNP 和 indel ，准备进行过滤。这种队列范围的分析即使在困难的位点也能灵敏地检测变异，并产生一个平方的基因型矩阵，该矩阵提供有关所考虑的所有样本中所有感兴趣位点的信息，这对于许多下游分析很重要。此步骤运行速度非常快，可以在将样本添加到队列时随时重新运行，从而解决所谓的 N+1 问题。代码是：

${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" GenotypeGVCFs \
        -R ${GENOME} \
        -V gendb://6.gvcf/gvcfs_db \
        -O ./6.gvcf/germline_merge.vcf.gz \
        1>./6.gvcf/GenotypeGVCFs.log 2>&1

VariantRecalibrator

在数据预处理有一步是 BQSR 碱基质量重矫正，而这一步则是变异质量重矫正，两者是两回事，用到的工具也不同，需要区分开。这一步实际上是基于机器学习的方法，对原始的 vcf 文件进行变异质量重矫正并且进行过滤。不过存在一个的缺点：该算法需要高质量的已知变体集作为训练和真实资源，而对于许多生物来说，这些资源尚不可用。它还需要相当多的数据来了解好与坏变体的概况，因此在仅涉及一个或几个样本的小数据集、靶向测序数据、RNAseq 上使用可能很困难甚至不可能使用，以及非模式生物。对于上述提到的情况，需要改用硬过滤的方法，可以参考：Hard-filtering germline short variants 代码是：

    ## 首先是对 SNP 位点运行 VariantRecalibrator
    ${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" VariantRecalibrator  \
    -R ${GENOME}                                                                 \
    -V 6.gvcf/germline_merge.vcf.gz                                              \
    --max-gaussians 4                                                            \
    --resource:hapmap,known=false,training=true,truth=true,prior=15.0  ${hapmap} \
    --resource:omni,known=false,training=true,truth=false,prior=12.0   ${omni}   \
    --resource:1000G,known=false,training=true,truth=false,prior=10.0  ${phase1} \
    --resource:dbsnp,known=true,training=false,truth=false,prior=2.0   ${dbsnp}  \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR                \
    -mode SNP                           \
    -O 6.gvcf/snp.recal                 \
    --tranches-file 6.gvcf/snp.tranches \
    --rscript-file 6.gvcf/snp.plots.R   \
    1>./6.gvcf/SNP_VariantRecalibrator.log 2>&1

    ## 接着是对 INDEL 位点运行 VariantRecalibrator
    ${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" VariantRecalibrator  \
    -R ${GENOME}                                                                 \
    -V 6.gvcf/germline_merge.vcf.gz                                              \
    --resource:mills,known=false,training=true,truth=true,prior=12.0    ${indel} \
    --resource:dbsnp,known=true,training=false,truth=false,prior=2.0    ${dbsnp} \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR                \
    --resource:axiomPoly,known=false,training=true,truth=false,prior=10.0  ${Axiom_Exome} \
    -mode INDEL                           \
    -O 6.gvcf/indel.recal                 \
    --tranches-file 6.gvcf/indel.tranches \
    --rscript-file 6.gvcf/indel.plots.R   \
    1>./6.gvcf/INDEL_VariantRecalibrator.log 2>&1

ApplyVQSR

    ## 然后就是对 SNP 运行 ApplyVQSR
    ${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" ApplyVQSR \
    -V 6.gvcf/germline_merge.vcf.gz       \
    --recal-file 6.gvcf/snp.recal         \
    --tranches-file 6.gvcf/snp.tranches   \
    --truth-sensitivity-filter-level 99.0 \
    --create-output-variant-index true    \
    -mode SNP                             \
    -O 6.gvcf/snp.vqsr.vcf.gz             \
    1>./6.gvcf/SNP_ApplyVQSR.log 2>&1

    ## 对 INDEL 运行 ApplyVQSR，注意，这里输入的时候，采用上一步校正过 SNP 质量后得到的结果
    ${GATK} --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" ApplyVQSR \
    -V 6.gvcf/snp.vqsr.vcf.gz             \
    --recal-file 6.gvcf/indel.recal       \
    --tranches-file 6.gvcf/indel.tranches \
    --truth-sensitivity-filter-level 99.0 \
    --create-output-variant-index true    \
    -mode INDEL                           \
    -O 6.gvcf/snp.indel.vqsr.vcf.gz       \
    1>>./6.gvcf/INDEL_ApplyVQSR.log 2>&1

最后经过过滤的位点在 FILTER 会标记上 PASS 标签。