跳转至

Python

  • 分类于 Python
  • 需要 3 分钟阅读时间

1. 原文参考

● 利用贝叶斯定理判断SNV 一个位点所有可能的genotype如下图： ["AA", "CC", "TT", "GG", "AC", "AT", "AG", "CT", "CG", "TG"]

利用贝叶斯定理判断SNV的方法就是根据测序数据(D)计算每一种genotype(G)的概率，即Pr(G|D)。

也就是计算所有十种可能性的概率， 那么，概率最高的就是该位点的genotype。

根据贝叶斯定理计算Pr(G|D)公式如下：

\[Pr(G | D)=Pr(G)*Pr(D | G)/Pr(D)\]

我们的目的是比较不同genotype的概率，因此计算Pr(G)*Pr(D|G)/Pr(D)即可。

\[Pr(G|D)\propto Pr(G)*Pr(D|G)\]

现在有两个问题：

  1. 如何计算Pr(D|G)?
  2. 如何计算P(G)?

1.如何计算Pr(D|G)? 下面通过一个例子来说明计算Pr(D|G)的方法。data如下：

假设G为A1A2，如下：

Pr(D|G)的计算公式如下：

\[Pr(\mathcal{D}|G)=\prod_{i=1..n}Pr(b_{i}|\{A_{1},A_{2}\})=\prod_{i=1..n}\frac{Pr(b_{i}|A_{1})+Pr(b_{i}|A_{2})}{2}\]
\[b_{j}\in\{\mathrm{A},\mathrm{G},\mathrm{A}\}\]

P(b|A)的概率如下， 其中e为对应碱基的质量值：

\[\text{当 }b_i=A_j\quad Pr(b_i|A_j)=1-e_i\]
\[\text{当 }b_i\neq A_j\quad Pr(b_i|A_j)=e_i/3\]

当Genotype=AG时，套用上面的公式，首先分别计算不同碱基对应的结果：

\[Pr(b_1=\mathtt A|\mathtt A\mathtt G)=\quad\frac{1}{2}\left((1-10^{-2})+\frac{10^{-2}}{3}\right)=\quad0.49667\]
\[Pr(b_2=\mathtt G|\mathtt A\mathtt G)=\quad\frac{1}{2}\left(\frac{10^{-1}}{3}+(1-10^{-1})\right)=\quad0.466667\]
\[Pr(b_{3}=\mathtt A|\mathtt A\mathtt G)=\quad\frac{1}{2}\left((1-10^{-5})+\frac{10^{-5}}{3}\right)=\quad0.499997\]

最终的结果如下：

\[Pr({D}|{AG})=Pr(b_1={A}|{AG}) * Pr(b_2={G}|{AG}) * Pr(b_3={A}|{AG}) = 0.49667 * 0.466667 * 0.499997=0.115888\]
  1. 如何计算P(G)？

计算P(G)， 就是计算10种可能的genotype的概率， 即Pr(AA)，Pr(TT)，Pr(CC)等等。

假定参考碱基为G，如果杂合突变率为0.001， 纯合突变率为0.0005。genotype的概率如下： 有了Pr(D|AG)和Pr(G)的概率， Pr(AG|D)就可以计算出来了， 如下：

\[Pr({AG}|{D})= Pr({D}|{AG})Pr({AG})=0.115888*0.001=0.000116\]

同理可以计算出其他genotype的概率， 如下：

因此Pr(AG)是10种genotype的概率最高的， 该位点的genotype为AG。

实际SNV calling的方法，可能还有很多细节， 或者其他方法，我这里分享的东西也可能存在不少错误，如果你发现了， 欢迎留言。

2. 代码实现

Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
import math
import itertools as it

# 计算基因型的先验概率
def calc_prior_prob(ref_base, g1, g2):
    # 假设已知的参数
    # REF_BASE = 'G'  # 参考碱基
    HET_RATE = 0.001  # 杂合突变率
    HOM_RATE = 0.0005  # 纯合突变率
    if g1 == g2 == ref_base:
        return 1.0 - HET_RATE - HOM_RATE
    elif g1 != g2:
        return HET_RATE
    else:
        return HOM_RATE

# 根据PHRED质量值计算错误率
def get_error_rate(qual):
    # 将PHRED分数转换为概率
    error_rate = 10 ** (-qual / 10.0)
    return error_rate

# 计算读取数据的条件概率
def calc_likelihood2(data, g1, g2, base_qualities):
    likelihood = 1.0
    for base, qual in zip(data, base_qualities):
        p_err = get_error_rate(qual)
        p_A1 = get_genotype_rate(base, p_err, g1)
        p_A2 = get_genotype_rate(base, p_err, g2)
        likelihood *= (p_A1 + p_A2) / 2
    return likelihood


def get_genotype_rate(base, p_err, genotype):
    if base == genotype:
        return 1.0 - p_err
    else:
        return p_err / 3

# 计算后验概率
def calc_posterior_prob(ref_base, data, base_qualities):
    probs = {}
    genotype_list = get_genotype(ref_base, data)
    for g1,g2 in genotype_list:
        prior = calc_prior_prob(ref_base, g1, g2)
        likelihood = calc_likelihood2(data, g1, g2, base_qualities)
        probs[(g1, g2)] = prior * likelihood
    return probs


def get_genotype(ref_base, data):
    res = []
    tmp = list(set(data))
    if ref_base not in tmp:
        tmp.append(ref_base)
    for e in it.combinations_with_replacement(tmp, 2):
        res.append(e)
    return res


# 使用示例
REF_BASE = 'G'  # 参考碱基
data = ['A', 'G', 'A']
base_qualities = [20, 10, 50]
posterior_probs = calc_posterior_prob(REF_BASE, data, base_qualities)

# 输出每种基因型的后验概率
for genotype, prob in sorted(posterior_probs.items(), key=lambda x: x[1], reverse=True):
    print(f"Genotype {genotype}: {prob:.6f}")