基于de novo似然方法（2）

方法：

概率评估的理论基础：

在本节中，我们对拼接质量的概率模型和测序过程的模型进行形式化，测序过程模型允许我们计算一个reads集合的任意一个特定拼接的似然度。我们证实我们的概率打分方法是正确的，因为真实的基因序列的得分最高。

拼接的似然度：

让A 表示一个给定的拼接是真实的基因拼接，让R表示观察到一系列reads的事件，我们在下文中将使用相同的符号表示拼接的序列和观察到这个拼接的事件。我们也用 相同的符号表示reads的集合和观察到这个reads集合的事件。

根据贝叶斯准则，给定一个观察到的reads的集合，产生相应拼接的概率可以被写成：

Pr[A|R]=Pr[R|A]Pr[A]/Pr[R]

Pr[A]是观察到基因序列A的先验概率，任意关于被拼接基因的先验知识都可以被包括进Pr[A]；但是，为了文章的目的，我们假设对于相同的read集合，任何合理的拼接的先验概率是常数。给出基因的一般的可以得到的信息，为Pr[A]构建一个精确的数学框架是超出本文范围的额外的扩展。

相似的，Pr[R]是观察到reads集合R的先验概率，因为我们的主要目标是比较相同read集的不同拼接，而不是获得一个对拼接质量的普遍的精确度量，我们也可以把Pr[R]认为是常数。因此，为了比较拼接器，Pr[A|R]和Pr[R|A]是等价的。后面，一个reads集合的后验概率，给定一个该数据的特定拼接，可以根据基于测序过程的大致模型很快的算出，并在本文中被当作拼接质量的代表。

我们假设各个reads是相互独立的，（在mate-pair实验中这个假设不成立，我们会在下文中讨论），那么Pr[P|A]=.......。如果reads的集合未排序，我们需要考虑产生相同reads集合的各种可能排列。因为这个值对于任何给定的reads集合是常数，我们在下文中都忽略它。

Pr[r|A]，在之后称为pr，可以用对测序过程的合适模型来计算。之后，我们将讨论复杂度递增的模型以及它们对似然度分值的影响。

假设我们有一个来自真实基因的reads的集合R，此集合是从基因中每个位置单端read，且长度固定，没有错误。给定一个reads的集合R，某一特定的read从真实基因中产生的概率正好等于read出现在R中的次数除以R的长度。（请注意当从repeats中产生时，多种read可以有相同的序列）。让Ns表示序列s在R中出现的次数，qs=Ns/|R|表示序列s从真实基因产生的概率。为了表示真实基因使得似然得分最高，我们假设有拼接A，ps表示序列s从拼接A产生的概率。

给定拼接A，我们的似然打分就是S中所有序列s的。。。。的乘积，可以被重写成。现在，请注意，因为|R|是个固定的常数，最大化似然得分等同于最大化。。。。

似然值可以被重写为。。。。。。。

其中Dkl（Q||P）是分布Q和R的kl-差异值，H（Q）是Q的香农熵。因为kl-差异值总是非负的，且当且仅当Q=P时，kl-差异值=0，因此，当拼接等同于真实的基因时，平均概率最大。

虽然真实基因使得模型中的似然值最大，仍有可能存在其他的拼接能达到一样的最优值，只要这些拼接对与每一个序列s，概率ps都等同于概率qs。比如，虽然发生了拼接错误，但是拼接错误之后却与产生的reads一致了，这时，就会发生上述的情况。这种情况凸现出现代测序实验固有的信息缺失，没有了额外的远程信息，reads本身提供的信息不足以区分一个基因多种可能的重构。

片段测序的无错模型：

测序过程的最基本模型是无错模型。在此模型中，我们假设reads的长度固定（一个更一般的read长度分布可被包含在这个模型里而不会影响从相同read集合获得的各个拼接的评估）。我们进一步假设reads是从基因上均匀采样的，也就是，基因上的每一个位点都可能是一个read的起始点。几乎所有的拼接概率模型都作了这个简化假设，尽管所有的现代测序技术都存在误差。一个更精确的、独立于技术的的模型可以通过增加额外的因素来考虑比如DNA构成误差。

为了一般性的目的，我们的讨论局限在在统一的采样模型上。而且，为了简单性的目的，我们假设（1）真实基因包含单独的成环的连续序列（2）我们的拼接也是一个单独的contig (3)每个read都可以被map到拼接上。我们将讨论对模型的扩展来放宽这些假设条件。

在这些假设条件下，给定一个拼接序列，我们可以计算一个read的概率为：pr=nr/2L

其中，nr表示read在长度为L的拼接中出现的位置的数目，除以2是因为reads是以相同的概率从一个DNA的正向和反向采样的，这个公式先前被Medvedev用来定义基因拼接的目标函数。

一个测序过程的实际模型：

上述的无错模型做出了许多与真实数据集特征不符的假设。这里我们显示如何扩展这个模型来考虑诸如测序错误、mate-pair信息、包含多个contig的拼接以及无法map的read的情况。

测序错误：

当前DNA的测序技术都存在小但是不能忽略的错误率。这里我们关注三种类型的错误：插入、缺失和核苷酸的替换。

在无错模型中，read从序列中位置j处产生的概率是1，如果read与此位置的reference精确匹配，否则是零。现在我们扩展这个模型，即read从reference上任意位置产生的概率是一个0和1之间的实数，代表着一种测序仪器从reference的特定位置产生特定read的概率。这个值明显依赖于read序列和reference序列在位置j处的差异。给定一个拼接序列，一个read的概率是此read在基因中所有位置出现的累计概率。

特别的，让我们用pr,j表示观察到read r在位置j处终止的概率。那么，read r的总的概率就是pr=pr，j（forward）+pr，j（reverse）/2L

单独的概率pr，j可以被计算，如果我们不考虑插入和缺失，而且只允许替代以e的概率发生。可以根据测序仪器生成read的质量来设定每个碱基发生替换错误的概率。那么，pr，j=e^s*(1-e)^(l-s)。其中s是将read r匹配到reference序列j位置所需要的替换次数。在更一般的情况下，pr，j可以使用动态规划算法来计算。

通过动态规划来精确计算概率：