windows ubuntu 子系统，肿瘤全外篇，3. gatk中的BaseRecalibrator，HaplotypeCaller，ApplyVQSR

1.碱基质量校正。

2.单样本变异调用

3.变异质控

2中，我们对测序数据进行了比对，bam排序，标记重复和建立索引。这次我们就直接可以进入gatk流程了。

1.碱基质量校正。

BaseRecalibrator 是 GATK（Genome Analysis Toolkit）工具集中的一个工具，用于校正测序数据中碱基质量的偏差，从而提高变异检测的准确性。

在测序数据分析中，由于测序平台、化学试剂、PCR等因素的影响，碱基质量可能存在偏差，例如一些碱基可能比其他碱基更容易出错。BaseRecalibrator 通过分析测序数据中的已知变异位点，统计每个碱基的观察频率和预期频率，然后计算出修正因子，以校正碱基质量的偏差。

BaseRecalibrator 的工作流程大致如下：

收集变异位点信息： 从已知的变异位点（如 dbSNP 数据库中的位点）中提取出来，并根据这些位点的信息统计每个碱基的观察频率和质量分数。
建立模型： 基于收集到的信息，建立一个模型，用于描述碱基质量与观察频率之间的关系。
计算校正因子： 使用建立的模型，计算每个碱基的校正因子，用于修正碱基质量的偏差。
生成校正后的 BAM 文件： 应用计算得到的校正因子，对测序数据中的碱基质量进行校正，生成校正后的 BAM 文件。

1.1生成校正表

gatk BaseRecalibrator \
        -R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
        -I ./L1.markdup.bam \
        -known-sites /mnt/h/db/gatk/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
        -known-sites /mnt/h/db/gatk/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
        -O L1.IndelRealigner.intervals

gatk BaseRecalibrator: 这是运行 GATK 工具集中 BaseRecalibrator 工具的命令。

-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta: -R 参数指定参考基因组的 FASTA 文件路径。这里使用了 Homo sapiens 的 hg38 版本作为参考基因组。

-I ./L1.markdup.bam: -I 参数指定输入的 BAM 文件，其中 L1.markdup.bam 是标记去重后的测序数据文件。

-known-sites /mnt/h/db/gatk/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz: -known-sites 参数指定已知的变异位点文件，这里使用了 1000 Genomes 项目中高置信度的单核苷酸多态性（SNPs）的 VCF 文件。

-known-sites /mnt/h/db/gatk/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz: -known-sites 参数再次指定已知的变异位点文件，这次使用了 Mills & 1000 Genomes 项目中的金标准插入/缺失（Indels）的 VCF 文件。

-O L1.IndelRealigner.intervals: -O 参数指定输出文件的路径和文件名，这里生成的文件名为 L1.IndelRealigner.intervals，用于存储间插片段重对齐所需的间隔信息。

1.2 应用校正表

gatk ApplyBQSR \
        -R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
        -I ./L1.markdup.bam \
        -bqsr L1.IndelRealigner.intervals \
        -O L1_recalibrated_reads.bam

gatk ApplyBQSR: 这是运行 GATK 工具集中 ApplyBQSR 工具的命令。

-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta: -R 参数指定参考基因组的 FASTA 文件路径，这里使用了 Homo sapiens 的 hg38 版本作为参考基因组。

-I ./L1.markdup.bam: -I 参数指定输入的 BAM 文件，其中 949743-T_L2_1.markdup.bam 是标记去重后的测序数据文件。

-bqsr L1.IndelRealigner.intervals: -bqsr 参数指定之前生成的校正因子文件，这里使用了L1.IndelRealigner.intervals 文件。

-O L1_recalibrated_reads.bam: -O 参数指定输出文件的路径和文件名，这里生成的文件名为 L1_recalibrated_reads.bam，用于存储校正后的测序数据。

2.单样本变异调用

HaplotypeCaller是GATK（Genome Analysis Toolkit）中的一个重要工具，用于对测序数据进行变异（SNP和Indel）的检测和分析。它采用了全局的分析方法，利用多个样本的信息来识别可能的变异位点，并生成高质量的变异呼叫。用这个方法不用提前局部重比对。
主要特点：

基于组装的变异检测：HaplotypeCaller通过对碱基片段进行组装来识别变异位点，而不是简单地对每个碱基位置进行独立的分析。这种方法可以更准确地区分变异和测序错误。
局部重比对：在变异检测之前，HaplotypeCaller会进行局部重比对，以更好地处理插入缺失（Indels）。
多样本模式：HaplotypeCaller能够同时处理多个样本的测序数据，以提高变异检测的准确性。
基于贝叶斯方法的变异呼叫：HaplotypeCaller使用贝叶斯方法来计算每个可能的变异的概率，并在确定最终呼叫时考虑这些概率。
局部重比对：首先，HaplotypeCaller会对比对过的测序数据进行局部重比对，以纠正插入缺失导致的比对错误，并提高变异检测的准确性。
组装候选片段：接下来，HaplotypeCaller会从测序数据中组装出可能的等位基因片段（haplotypes），这些片段是由候选变异引起的，或者是由于测序错误而引入的。
变异检测：然后，对每个组装出的候选片段进行变异检测，通过比较片段与参考基因组之间的差异来识别SNP和Indel。
变异呼叫：最后，使用贝叶斯方法计算每个候选变异的概率，并基于这些概率进行变异呼叫，以确定最终的变异集合。

gatk HaplotypeCaller \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-I L1_recalibrated_reads.bam \
--dbsnp /mnt/h/db/gatk/hg38/dbsnp_146.hg38.vcf.gz \
-O L1_raw.vcf \
1>L1_log.HC 2>&1

-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta: -R 参数指定参考基因组的 FASTA 文件路径，这里使用了 Homo sapiens 的 hg38 版本作为参考基因组。

-I L1_recalibrated_reads.bam: -I 参数指定输入的 BAM 文件，其中 L1_recalibrated_reads.bam 是经过校正后的测序数据文件。

--dbsnp /mnt/h/db/gatk/hg38/dbsnp_146.hg38.vcf.gz: --dbsnp 参数指定已知的变异位点文件，这里使用了 dbSNP 数据库中的 hg38 版本的 VCF 文件。

-O L1_raw.vcf: -O 参数指定输出文件的路径和文件名，这里生成的文件名为 L1_raw.vcf，用于存储未经过过滤的原始变异结果。

1>L1_log.HC 2>&1: 这部分是将标准输出和标准错误输出重定向到文件 L1_log.HC 中，以便记录日志信息。

GATK HaplotypeCaller之后，需要进行质控，可以用GATK Variant Quality Score Recalibration（VQSR）。VQSR的原理是利用自身数据和已知变异位点集之间的重叠，通过高斯混合模型（GMM）构建一个分类器来对变异数据进行打分，从而评估每个位点的可靠性。具体而言，VQSR使用已知变异位点集作为训练集，通过统计特征（如变异质量、测序深度等）和已知位点的质量得分，构建一个GMM模型。然后，该模型会根据自身数据的特征和已知位点的质量得分，为每个变异位点分配一个新的质量得分。最终，VQSR根据这个新的质量得分，对变异位点进行过滤和分类，以确定哪些位点是高可靠度的变异。通过这一方法，VQSR可以帮助识别那些在自身数据中可能存在假阳性或假阴性的变异位点，并提供更可靠的变异集合，从而提高变异检测的准确性和可信度。

3.变异质控

接下来的命令需要 R: conda insatll R

3.1生成校正模型

gatk VariantRecalibrator \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-V L1_raw.vcf \
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 /mnt/h/db/gatk/hg38/hapmap_3.3.hg38.vcf.gz \
-resource:omini,known=false,training=true,truth=false,prior=12.0 /mnt/h/db/gatk/hg38/1000G_omni2.5.hg38.vcf.gz \
-resource:1000G,known=false,training=true,truth=false,prior=10.0 /mnt/h/db/gatk/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 /mnt/h/db/gatk/hg38/dbsnp_146.hg38.vcf.gz \
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum -mode SNP -tranche 100.0 \
-tranche 99.9 -tranche 99.0 -tranche 95.0 -tranche 90.0 \
-O ./snprecal/L1_WES.snp.recal.vcf \
--tranches-file ./snprecal/L1_WES.snp.tranches \
--rscript-file ./snprecal/L_WES.snp.plots.R

-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta：指定参考基因组的FASTA文件路径。

-V L1_raw.vcf：指定原始的VCF文件路径，这是待校正的文件。

-resource:hapmap,known=false,training=true,truth=true,prior=15.0 /mnt/h/db/gatk/hg38/hapmap_3.3.hg38.vcf.gz：指定用于训练和校正的参考资源。这里包括hapmap、omini、1000G和dbsnp四个资源。每个资源都有不同的权重（prior），hapmap和omini的truth值为true，而1000G和dbsnp的known值为true。这些资源通常用于帮助确定变异的真实性和可靠性。

-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum：指定用于校正的注释标签（annotations）。这些标签包括深度（DP）、质量/深度比值（QD）、分段比例（FS）、strand偏差（SOR）、读取位置秩和（ReadPosRankSum）以及映射质量秩和（MQRankSum）等。

-mode SNP：指定模式为SNP，表示这个命令用于对SNP进行校正。

-tranche：指定用于切分数据的阈值。这里有四个不同的阈值：100.0、99.9、99.0和95.0。这些阈值用于确定变异的可靠性等级，例如，阈值100.0表示最高可靠性。

-O ./snprecal/L1_WES.snp.recal.vcf：指定输出的校正后的VCF文件路径。

--tranches-file ./snprecal/L1_WES.snp.tranches：指定输出的阈值文件路径，其中包含了根据给定阈值进行的变异分级。

--rscript-file ./snprecal/L_WES.snp.plots.R：指定输出的R脚本文件路径，用于绘制校正过程中生成的各种图表。

这个命令通过对原始的VCF文件进行训练和校正，生成了一个校正后的VCF文件，其中包含了经过过滤和分级的SNP数据，以及相应的阈值文件和R脚本用于可视化分析。

3.2 运用校正模型

gatk ApplyVQSR \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-V ../L_1_raw.vcf \
--ts-filter-level 99.0 \
--tranches-file L1_WES.snp.tranches \
--recal-file L1_WES.snp.recal.vcf \
-mode SNP \
-O 19P0126636WES.snps.VQSR.vcf.gz

-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta：指定参考基因组的FASTA文件路径。

-V ../L1_raw.vcf：指定原始的VCF文件路径，这是待应用VQSR的文件。

--ts-filter-level 99.0：指定过滤阈值的水平，这里设定为99.0，表示将满足99.0%可靠性阈值的变异保留下来。

--tranches-file L1_WES.snp.tranches：指定之前生成的阈值文件路径，用于对变异进行过滤。

--recal-file L1_WES.snp.recal.vcf：指定之前生成的校正文件路径，用于对变异进行校正。

-mode SNP：指定模式为SNP，表示这个命令用于对SNP数据进行VQSR。

-O 19P0126636WES.snps.VQSR.vcf.gz：指定输出的VCF文件路径，这是应用VQSR后的结果文件，经过过滤和校正的SNP数据。

这个命令将会根据之前建立的校正模型，对原始的VCF文件中的SNP进行过滤和校正，并生成一个经过VQSR处理的结果文件。

这就生成了一个VCF文件，流程跑完了，但是感觉还是啥也不会，先不说了，还是有空补充下基础知识吧。跑流程的过程中遇到了无数大坑，一坑更比一坑坑，大家想学，还是要自己跑一下。