一文搞懂———碱基质量值校正 BQSR

2023-03-17 本文已影响0人生信小书童

1 BQSR定义

Base Quality Score Recalibration(BQSR) 代表碱基质量评分再校准。简而言之，它是一个数据预处理步骤，检测测序仪在估计每个碱基调用的准确性时所产生的系统错误。

2 为什么要进行BQSR校正

测序仪给出的质量值不可靠。而后续的检测变异在很大程度上要依赖于变异质量值。测序仪给出的质量值受各种系统的影响，导致数据中的碱基质量分数过高或过低，其中一些错误是由测序反应的物理或者化学原理，一些可能是由于设备制造的缺陷

3 BQSR校准的原理

碱基质量值再校准(BQSR)是一个应用机器学习对这些误差进行经验建模并相应地调整质量评分的过程。例如，我们可以确定，对于给定的run，每当我们连续检出两个A核苷酸时，我们检出的下一个碱基的错误率会高1%。因此，在读操作中，任何在AA之后的基本调用的质量分数都应该降低1%。我们对几个不同的协变量(主要是碱基序列上下文、碱基在reads中的位置或cycle的个数)以一种可加的方式进行处理。因此，同一碱基的质量分数可能会因某种原因而增加，因另一原因而降低。碱基质量校正不会改变碱基的类型，只会降低或者升高碱基的质量值。测序仪厂商通常会倾向于将质量值提高。

此过程可应用于包含来自任何排序平台的数据的BAM文件，该测序平台在预期的尺度上输出基本质量分数。我们自己使用了来自Illumina、SOLiD、454、Complete Genomics和Pacific Biosciences测序仪的几代数据。

4 BQSR校准步骤

碱基再校准过程包括两个关键步骤:首先，BaseRecalibrator 命令 基于输入数据和一组已知突变（已知的突变位点就是一些dbsnp库啦）构建共变模型，生成重新校准文件; 然后 ApplyBQSR命令 根据模型调整数据中的基本质量分数，生成一个新的BAM文件。输出已知突变是被用来掩盖真实(预期)变化位点的碱基，以避免将真实变量计算为误差。在屏蔽站点之外，每个不匹配都被算作错误。剩下的主要是统计。

5 BQSR校准细节

（1）BaseRecalibrator 构建模型步骤

# baserecal
/bioapp/java8/bin/java -Xmx16g -Djava.io.tmpdir=${outdir}/${sample}_tmp -jar /bioapp/gatk-package-4.0.11.0-local.jar \
                BaseRecalibrator \
                -I=${outdir}/${sample}_rmdup.bam \
                --known-sites=$DBSNP_B37  \
                --known-sites=$PSINDEL \
                --known-sites=$GSINDEL \
                -O=${outdir}/${sample}_recal.grp  \
                -R ${ref}

为了构建重新校准模型，第一个工具将遍历输入BAM文件中的所有读取，并将关于碱基的以下特征的数据制成表格:
1、查看read属于哪个分组
2、机器报出的质量值
3、碱基所在的reads的位置
4、当前碱基+后一个碱基（二核苷酸）

对于每个bin，计算bin中的碱基数量以及这些碱基与参考库不匹配的频率，根据已知的dbSNP数据库，排除已知在总体中变化的位点。该信息以GATKReport格式输出到重新校准文件。

请注意，较准步骤应用“yates”校正低占用bins。而不是从# mismatches / # bases中推断真实的Q值，我们实际上是从(# mismatches + 1) / (# bases + 2)中推断出来的。这很好地解决了过拟合问题，这对拥有数十亿个碱基的数据集影响很小，但对于避免稀疏数据中稀有箱子的过度自信至关重要。

（2）ApplyBQSR 步骤

/bioapp/java8/bin/java -Xmx16g -Djava.io.tmpdir=$outdir/${sample}_tmp -jar /bioapp/gatk-package-4.0.11.0-local.jar \
                ApplyBQSR -bqsr=${outdir}/${sample}_recal.grp \
                -I=${outdir}/${sample}_rmdup.bam \
                -O=${outdir}/${sample}_rmdup_realigned_recal.bam

第二个工具再次遍历所有读取，使用重新校准文件根据它所在的箱子调整每个基数的分数。所以实际上新的质量分数是:
1、汇总报告质量值和经验质量值之前的差异
2、加上质量值的偏移
3、加上x循环的质量值和双核苷酸质量值的影响

经过重新校准，reads质量分数比以前更接近他们的经验分数。这意味着它们可以统计上可靠的方式用于下游数据的处理，例如变异检测。此外，通过考虑循环和序列上下文的质量变化，我们可以在读取中识别出真正高质量的碱基，经常找到一个碱基的子集，即使最初没有碱基被标记为Q30。

6 影响校正的因素

6.1 Read groups

重新校准系统是按照read-group分类的，这意味着它使用 @RG 标记按读组划分数据。这允许它对每个 read-group 执行重新校准，这反映了一个read属于哪个库，以及它在流单元上的哪个通道进行了排序。我们知道系统偏差可能发生在一个通道而不是另一个，或者一个库而不是另一个，所以能够在每个序列数据单元内重新校准使建模过程更加准确。因此，这意味着可以在具有多个读组的BAM文件上运行BQSR。但是，请注意，内存需求与文件中read-group的数量成线性关系，因此具有许多read-group的文件可能需要大量RAM来存储所有协变量数据。

6.2 数据量

质量值校准的一个关键决定因素是每个 bin 中观察到的碱基和不匹配的数量。当比对reads较少时，此过程将不能很好地工作。我们通常期望每个读组有超过100M个碱基;根据经验，数字越大效果越好。

7 校准理由

你几乎总是应该对你的测序数据进行重新校准。在人类数据中，考虑到我们所拥有的详尽的变异数据库，几乎所有剩余的不匹配——甚至在癌症中——都将是错误，所以为数据确定准确的错误模型非常容易，这对下游分析至关重要。对于非人类数据，如果你的生物没有可用的资源，你可能需要引导你自己的变体集，这可能会多一点工作，但这是值得的。

建立已知变异的方法：
1、首先，对原始的、未重新校准的数据进行第一轮变量检测。
2、然后，通过将其作为VCF文件提供给baserecaliator，将您最信任的变量集用作已知变量的数据库。
3、最后，使用重新校准的数据进行一轮真正的变量调用。这些步骤可以重复几次直到收敛。