焦磷酸测序（2）Quality分值和Error Model

terence

测序仪得到信号后，第一步即是将信号转换为碱基序列。即是用同样的样本，测序时试剂的剂量和批次不同，出来的结果也会不同。不同的实验室出来的测序结果也不全是一样的。哇，那么怎么评价一次测序的质量呢？有没有一个统一的标准？这还得感谢Phil Green老爷爷了（他的实验室全是一帮老头，好像还在不断的coding，老当益壮呀。Jun Liu有次说过，D. Rubin向他抱怨，现在的学生都不喜欢写代码:) ）。Phil Green提出了Phred score，也是现在新一代测序计算Quality的一般思路。

什么叫质量好，什么叫质量差呢？直观上来说，如果一个base出错的概率很高，那么这个测序结果的可靠性就差，自然质量就低了。如果一个base的出错概率很低，那么质量就高。但是，说到这也还没有解决问题呀。对于测序，我们又不知道一个base背后的真正碱基是什么，又怎么计算出错的概率？ok，先写下Phil Green计算quality的公式，再来看看有什么方法可以解决这个问题。

Quality Score = -10 log( P( true base | observed base ) ), true base != observed base

P( true base | observed base ) 是出错的概率，举个例子，observed base 是 A， 出错的概率就是 P( C | A ) + P( G | A ) + P( T | A )，意思是产生A的真实碱基{C,G,T}的概率。Green很巧妙的用了Bayes公式，

P( true base | observed base ) ~ p( observed base | true base ) p( true base )

哇，变成这个形式，就可以轻而易举的计算quality了。这个公式指出了一条明路，可以人为造出一堆序列（test fragment），然后送去测序，再把测出来的读段map回这堆序列，这样我们不但知道真实的序列，而且知道读段和真实序列的关系，就可以统计p( observed base | true base )了。这个称作测序仪的Error Model。

现在，新一代测序仪都是按照这个思路计算quality分值。

但是，但是，但是，这种做法有一个致命的缺陷，假设Mapping是靠谱的，总能正确的找到读段是从测试序列的某个地方产生出来的。序列比对（alignment）对打分函数（scoring function）的依赖性很高，如果要mapping靠谱，打分函数就要反映测序仪的错误特点。但是，要知道测序仪的错误的特点，又必须先mapping才能统计。都想口出闽南三字经了，把人都绕晕了。实际上，这是一个典型的数据缺失问题，可以用Rubin发明的EM算法来解决。缺失的数据是alignment，而Durbin提出的HMM可以穷举所有可能的alignment，并且给每个alignment一个概率，这样就可以迭代求解测序仪的Error Model了。

Error Model最早是水侠提出的，后来USC的Lei Li把这个model给推广了。

有一次我问老板，是不是计算机的人都喜欢自己写程序，不喜欢用别人的代码。老板说，写程序是基本功。我只好默默的写自己的程序了。呵呵。

zouquan

赞原创！

我觉得计算机人都爱自己写程序，就像每个闲的人都想自己设计装修房子一样，纯粹是压力小时间多，呵呵