● 利用贝叶斯定理判断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)\]
现在有两个问题：
- 如何计算Pr(D|G)?
- 如何计算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\]
- 如何计算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的方法，可能还有很多细节， 或者其他方法，我这里分享的东西也可能存在不少错误，如果你发现了， 欢迎留言。
| 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}")
|
