Software-污染识别-VerifyBamID

发表于 | 分类于 ， ，

VerifyBamID 是一款软件，用于验证特定文件中的读取是否与个体（或一组个体）先前已知的基因型匹配，并检查读取是否因两个样本混合而受到污染。当外部基因型可用时， verifyBamID 可以检测样本污染和交换。当外部基因型不可用时， verifyBamID 仍然可以稳健地检测样本交换。
github
verifyBamID2 文章：Ancestry-agnostic estimation of DNA sample contamination from sequence reads

软件简介

使用方法

1
2
3
4
5
VerifyBamID --DisableSanityCheck \
    --SVDPrefix 1000g.phase3.10k.b37.exome.vcf.gz.dat \
    --Reference hg19.genome.fa \
    --BamFile $final_bam \
    --Output $qc_dir/$sampleID

以下内容由AI解析完善

检测原理

一、整体执行逻辑

  1. 载入参考/基因型信息（GenMatrixBinary）

    根据 --SVDPrefix / 指定的 VCF/二进制前缀载入变异位点和/或已处理的矩阵/频率信息。实现里会：
    读取 VCF（或预处理数据），按位点提取等位基频（AF）、位点基本信息及（当可用）参考样本的基因型/个体 ID。
    过滤非常染色体位点（只保留 autosomes）。
    应用阈值过滤（minAF、minCallRate 等）。
    为后续步骤建立每个位点的等位基因频率或 genotype matrix。

  2. 在 BAM 上做 pileup（BamPileBases）

    在上一步的位点集合上，对输入 BAM 做 pileup，按位点收集观测到的碱基、对应的碱基质量（Phred）、mapping quality、read group、以及各类过滤（SAM flags、minMapQ、minQ、maxDepth、是否忽略 overlapping pair 等）。
将 Phred 质量转换成错误率（e = 10^(-Q/10)），并对有上限的 Q（maxQ）做截断；实现中有 fPhred2Err[] 数组用于快速查表。

  3. 每个位点计算观测数据的似然（基于基因型+测序误差+参考偏差+污染模型）

    对每个位点，基于该位点的等位基因频率/参考基因型信息和 pileup 的碱基/质量，计算在不同假设下（例如样本基因型为 REF/REF、REF/ALT、ALT/ALT；或在含污染 f 的混合模型下）观测到这些碱基的概率（位点似然）。
关键构件：
    基于 Phred 得到单个碱基的“发射概率”：P(observed base | sampled allele) = 1 - e 当匹配，否则 ≈ e/3（或实现中对“其他碱基/插入”等有特殊处理）。
    基于基因型的“抽样等位基因概率”：当基因型为 het 时，读到 Ref 或 Alt 的概率通常为 0.5/0.5，但实现允许“reference bias”参数来调整这些采样概率（源码中用 pSN 矩阵保存 Pr(sampled allele | true genotype)）。
    当考虑污染（freemix）时，位点的观测被建模为两个个体的混合：观测 = (1 - f) * 源样本 + f * 污染者。实现里对污染者可能使用群体等位基频（按 AF 假设 Hardy-Weinberg）或在有参考个体时用具体个体基因型来对污染进行建模/积分。
位点似然是对该位点所有读的发射概率（通常相乘，取对数求和以避免下溢）。

  4. 全基因组（全位点集合）上对污染率 f、参考偏差等参数的估计与假设检验

    通过在 f 的取值网格（VerifyBamIDArgs.grid，命令行默认 0.05）上计算全局 log-likelihood，先做粗略的网格搜索找出较好区域，再做局部优化（或 EM）来得到最大似然估计（MLE），这就是 freemix（估计的污染比例）。
若程序允许估计 reference-bias 或其他参数（由 bFreeMixOnly、bFreeRefBiasOnly、bFreeFull 等控制），会联合或交替估计这些参数以最大化总体似然。
计算对比统计量，例如：
    LLR（log-likelihood ratio）或 delta log-likelihood：比较 f = 0 (无污染假设) 与 f = mle 的对数似然差，用于衡量污染信号强度。
    p-value / z-score（若实现提供）或其他基于似然比的显著性度量。

    “自检 / 匹配参考个体”（self / best）

    当 VCF 中包含个体基因型时，VerifyBamID 能评估该 BAM 是否来自 VCF 中给定或任一参考个体：
    逐位点将测序似然与参考个体基因型进行比对，计算样本与参考个体相同/不同的似然或 posterior，从而输出“是否与某个 VCF 个体匹配”以及最佳匹配（best match）的 ID 和相应分数。
这一步与污染估计可以同时给出：“这个 BAM 看起来像哪个参考样本？” 以及“是否有显著污染”。

二、主要输出（摘要）与每个输出的计算/含义 下面列出 VerifyBamID 的常见主要结果项及其计算逻辑（不同版本在命名/格式上可能有差异，但基本统计是如下几类）：

  1. FREEMIX（或 mlefMix / estimated contamination f）

    意义：估计的外源污染比例 f（0-1），表示测序读中有多少比例来自非目标样本的 DNA。
计算逻辑：最大化全位点总体似然得到的污染参数 MLE。总体似然为各位点似然的乘积（或对数和），每个位点的似然在污染模型下为： P(obs | f) = sum_{g_source} P(g_source | prior) * [ (1 - f) * P(obs | g_source) + f * P(obs | contaminant) ]
— 其中 P(obs | g) 又是各读的发射概率乘积；contaminant 的成分可用群体 AF 折算成基因型分布并对其积分求和（即把污染者视为从群体中随机抽取）。
实现上常用网格搜索（指定 grid）+ 局部优化；输出通常是 MLE 的数值（如 0.023 表示 ~2.3% 污染）。

  2. LLK / deltaLLK / LLR（对数似然/似然比）

    意义：度量模型拟合改善程度，例如将 f = 0 与 f = mle 的对数似然差，或将“样本与某参考个体为同一人”的假设对比其它假设。
计算逻辑：对数似然 = sum_over_sites log P(obs_site | model_params)。LLR = 2*(logL(mle) - logL(null)) 或直接给出 logL 差值，用于判断是否显著偏离无污染假设。

  3. BEST MATCH / SELF 判定（bestOut / selfOut）

    意义：如果你传入了含基因型的 VCF，可得出样本最可能对应的 VCF 个体（best match），以及与给定样本 ID 的 self-check（样本是否与给定 ID 匹配）。
计算逻辑：对于每个参考个体 i，计算 P(obs | genotype_i)（或含污染模型下的观测似然），选取使似然最大的个体并输出相应统计（例如匹配概率、位点不一致数、覆盖位点数等）。

  4. 位点/总体覆盖深度与位点数（nMarkers、nBases、DP、mean depth）

    意义：用于说明用于估计的位点数量、总读数或平均深度，及过滤后实际参与计算的位点数。
计算逻辑：在 pileup 步骤统计，考虑过滤条件（minMapQ、minQ、maxDepth、exclude flags、非 autosome 的跳过等），并累积有效读数与有效位点。

  5. Call Rate / site call statistics

    意义：在 VCF 与 BAM 的交叉位点中有多少位点可用（达到最小覆盖/质量），以及用于匹配/估计的有效位点比例。
计算逻辑：在读取 VCF 时按 minCallRate 过滤位点；在 BAM pileup 时按 minQ/minMapQ/maxDepth 等过滤读并决定该位点是否被接受为“可调用”。

  6. Reference-bias 参数（若估计：pRefRef、pRefHet、pRefAlt）

    意义：描述同一真实基因型下，测序/对齐导致读到参考/替代等位基因的偏向；若存在系统性偏向（例如对参考等位基因更容易比对/呼出），会影响基因型似然计算与污染估计。
计算逻辑：在模型中把“读到哪个等位基因”的概率从均等（het 情况下通常 0.5/0.5）改为可调参数 pSN；程序可以在自由/固定的设置下用最大似然估计这些参数，使得总体似然最大。源码中 pSN[i*3+j] 存储 Pr(sampled=i | true=j)。

  7. 质量控制/诊断项（例如不一致位点数、各等位基因读数分布、各 ReadGroup 统计等）

    意义：帮助判断污染是否受某一文库/运行影响、是否需要按 read-group 分析等。
计算逻辑：由 pileup 阶段按 read group 分组统计碱基计数、深度、质量分布等；并在模型里可输出按 RG 的 contamination estimate（若实现支持）。

三、关键数学细节（更具体的概率表达式）

  1. 单碱基发射概率（简化表达）： e = 10^{-Q/10} （实现中 Q 超过 maxQ 时有截断） 若 sampled allele = a，observed base = b： P(b | sampled allele = a) = (1 - e) if b equals a ≈ e/3 if b different（或按 implementation 细节分配到“其他”类别）
    基因型到 sampled allele 的概率（含 reference-bias）： pSN 表示 Pr(sampled allele = i | true genotype = j)；例如
    若 true genotype = homozygous ref，则 sampled allele 非 ref 的概率很低（由 pRefRef/pRefAlt 参数控制）
若 het，则默认 0.5/0.5，但可被 pRefHet 调整为偏向 ref 或 alt
    含污染模型的位点似然（概念式）： P(obs_site | f) = sum_{g_source} P(g_source) prod_{reads r} [ (1 - f) P(read_r | g_source) + f * P(read_r | contaminant_model) ]
    contaminant_model 可以是基于群体 AF 的基因型分布积分（若不知道污染者是谁），或基于某参考个体的基因型（若做 pairwise 匹配）。
    全局对数似然 = sum_over_sites log P(obs_site | params)。最大化这个值得到参数估计。

四、命令中特别选项的影响

--DisableSanityCheck：跳过一些输入/文件一致性检查，程序会更快直接运行，但可能忽略格式/匹配错误。
--SVDPrefix 1000g.phase3.10k.b37.exome.vcf.gz.dat：表明使用预处理的 1000G（或对应）位点/矩阵信息作为参考（用于 AF、样本匹配或 SVD/投影以找 best match）。
--Reference hg19.genome.fa：参考基因组用于 alignment 相关的位点定位或某些比对相关过滤。
其它参数（默认为 VerifyBamIDArgs 中的值）：grid、minQ、maxDepth、minMapQ、genoError、minAF、minCallRate 等，都会影响位点筛选与似然计算。

五、给用户的实用说明（如何解读结果）

FREEMIX（估计污染率）是最重要的输出之一：一个小值（如 < 0.01）通常表明污染可以忽略；较大值（> ~0.02–0.05）需要关注和可能的样本重做/去污处理。
LLR（或 delta log-likelihood）说明污染估计相较于无污染模型的支持度。数值越大，越支持存在污染。
BEST MATCH / SELF check：如果你有 VCF 的个体 ID，查看最佳匹配和匹配概率可以检验样本 ID 对应性（样本错配会在这里显著表现）。
覆盖深度/位点数：小的有效位点数或低平均深度会降低污染估计的精度（不可靠）。
若程序同时估计 reference-bias，应检查该参数是否显著偏离中性（0.5），大偏差说明技术偏向可能需要考虑在下游分析中纠正。

六、典型输出文件/行（示例性的说明）

程序通常会在输出文件（你用 --Output sampleID 指定的前缀）或控制台打印一行或几行摘要，包括：sample ID / FREEMIX / LLR / #sites used / mean depth / bestMatchID / bestMatchScore 等。不同发布版具体字段名会有差别，但上面提到的统计项是核心。
-------------本文结束感谢您的阅读-------------