Hybrid error correction and de novo assembly of single-molecule sequencing reads

本文链接：https://blog.csdn.net/weixin_42472706/article/details/88951572

今天介绍的文献是关于三代测序的拼接算法的研究，或者准确来讲是一个利用二代测序和三代测序结果来拼接基因组序列的方法，当然这篇文章其实只是做了一个部分的创新，大部分拼接的原理是基于已有的拼接算法。当然，这篇文献其实有点旧了，是2012年发表在nature子刊biotechnology上的。
标题里的de novo是一个拉丁语貌似，文献中比较喜欢用这个词表示重新、从头开始的意思。single-molecular sequencing就是我们常说的三代测序。read在生物测序中指的是一个测序片段，称为一个read。因为，现在一次测序能检测的序列长度有限，所以真实测序中我们对一个完整基因组随机打断成若干碎片，碎片的长度小于我们一次测序能检测的长度，最后测出来的这些小片段就是read。然后我们基于这些碎片进行拼接，由于打断是随机的，所以我们可以把一些有重叠区域的两个read拼接起来。

背景

先看一下文章的背景知识，总的来讲就是二代测序每次得到的reads比较短，而三代测序得到的reads比较长。但是，三代测序准确度比二代测序差太远。read片段长，对于转录组或者基因组的拼接就越容易，有模型表明1000bp的read长度比100bp的序列长度拼接结果的连续性要至少强6倍（这里连续性是指，对于转录组或者基因组的拼接在现在的算法水平基本不可能拼出来完整的结果）。但是，三代测序的测序精度实在是太差了，就是拼出来了也和真实的基因组差得远。所以这个这篇文章提出了针对这个问题的算法。

摘要

概括的将一下这个文章的工作，这篇文献的主要想法就是先用short read（二代测序结果）去矫正三代测序的结果，也就是先把三代测序的测得不对碱基给改过来，再用矫正过的三代测序的long read去进行de novo assembly。没有明白的可以结合上面的图理解一下。（上面的图不是原文献的图片，是我从其他文献找到的图，然后自己P过以后的一个文献流程示意图）原文献的NCBI链接
总结一下就是，文献提出了一个算法叫做PBcR（PacBio是太平洋生物公司的名字，c表示corrected矫正的Reads）这个算法的名字真是太low了，还把测序公司的名字带上了 。这个算法只是Celera Assember11这个拼接算法软件的一个环节，最终把准确度从80%提高到了99.9%。

具体方法

首先看第一步：long reads来源于PacBio公司的测序结果，short reads来自于Illumina公司的454二代测序和PacBio CCS的二代测序（PacBio公司现在已经被Illimina公司兼并了，回头看这个文献貌似埋下了伏笔哈哈哈）
然后关于能匹配的上的要求有如下几条：

必须共有一段长度大于14bp的种子序列，我觉得意思就是如上图，并且short read必须能全长aligne（比对上），其他的不予考虑。（这里更新一下关于seed sequence的一个补充的内容，之前看到这里不是很能理解为什么称为seed序列。最近学习了BLAST多序列比对的算法，也就是NCBI我们常用的那个BLAST系列功能，感觉对seed的认识刷新了。因为现在来讲multiple sequence alignment本身还是Pairwise Sequence Alignment，只不过把所有的相似度高的结果放到一起了，但是不可能把你要查询的序列和数据库里的所有序列都pair alignment一下，所以BLAST就会把你的查询序列分成小段，叫做seed sequence。一般来讲蛋白质是3个氨基酸长度（eg.下图），核酸是11个碱基长度。然乎这个seed会去寻找数据库的匹配序列，如果完全贴合就会向两侧延伸，如果延伸的长度没有达到阈值就会脱落，放弃这条序列。只有一定延伸到一定阈值长度的靶向序列才会被留下，进行下一步的Pairwise Sequence Alignment进行打分，把查询序列切成小段洒向数据库的这个过程感觉很像播种，所以把这个序列叫做seed，我猜测作者这里使用的这个软件，让short reed匹配到long read上也可能使用了类似的方法来缩短时间）

如上图，只有部分重叠区域在末尾的不予考虑。
重叠区域不考虑两个来自相同测序技术的片段，也就是说只考虑short-short这种不考虑，当然long-long也不考虑。

第二步是这个算法最核心的一步，也就是对匹配上的序列进行一个筛选，也就是对一些匹配程度不是很高的序列筛选掉。我自己根据文献，绘制了一个原理图可以帮助理解为什么要选择Top C hits的reads。

绿色的bar表示whole genome，理论上来讲，一个short read能匹配上的long read如左侧所示，基本就和这个位置的测序深度（coverage / depth）差不多。（测序深度简单来说就是这一段序列最后拼接的时候，这个区域被覆盖了多少层，详情见维基百科链接。
但是，如果像右侧这个long read，正确的应该是蓝色区域，但是由于错误比较多，所以变成了上面黄绿色的彩色区域，正好就和左侧的short read匹配上了，本来如果匹配一个正确的short read（右图蓝色short bar图形），就会被正确矫正，但是现在匹配了这个蓝紫色的short bar，就会被错误矫正。所以选择匹配度高的前C个long read对于提高准确度是有必要的。

我们再来看一下，这个文献提供的算法流程示意图，对于pink bar表示long read，上面的black vertical(垂直) bars表示错误的位点。blue bar表示short bar。在b图，blue bar上的black vertical bars表示匹配不上的位点，并不是错误的位点。可以看出来，很多时候pink bar上的错误位点，基本都是blue bar上匹配不上的位点。然后，将高匹配度的保留，剩下的抛弃，最后进行矫正。在c图，blue bar上的black vertical bars表示匹配测序错误的位点，毕竟short也不是完全正确的，有些很罕见的情况，这个错误最后被代入到了校正后的purple bar里面，当然两端有些gap区域，因为前面第一步匹配的时候，末端匹配上的根本那就不被考虑，所以两端空白也是正常的，不过不用担心，再下一步的long reads的拼接，可能就又拼上去了。最后e图的灰色区域，本身是two-copy inexact repeat，只使用short read这样的片段是拼不出来的。这里补充一下生物基因组中的知识，基因组中有高度重复序列（2-10核苷酸），出现频率106-108之间；中度重复序列（100~300核苷酸），出现频率102；单一序列（1000bp），即结构基因，大多单拷贝。而二代测序的片段中值是100bp最大150bp，而三代测序的长度是700bp中值，最大1000bp。二代测序对于中度重复序列刚好长度不足，拼接困难。而三代测序就可以跨越这个瓶颈。

然后我们具体的看一下，top C hits是怎么筛选的。首先绘制直方图，怎么绘制呢，我们看我标记图中的那个点，横坐标9表示能匹配上9条long read，也就是ni=9，而这样的short read有~20000左右。红线表示我们切割的Top C的位置，也就是short read最多匹配到C条long reads。根据作者的说法，这个图的一般形状就类似于泊松分布，coverage的分布也是类似的。

具体红线被确定出来的公式在上面，首先选定一个比例底线T（后文提到这个T的选择一般是大于最高峰）；slope(K)就是K和K-1间的斜率。total(K)指的是K前所有short read占的比例。最后C的求解是，遵照最小的K，满足total(K)大于T比例，斜率大于上一点的斜率，并且对应的值小于上一点。其实这个条件很明显可以看出来，是指一个峰的右侧（下降一侧）。我觉得之所以设定T，是防止下面这种情况（蓝线所示），而要达到红线所示。

第三步，用AMOS consensus module来矫正long reads。
第四步，用Celera Assembler来拼接得到结果。
结果
下面是文章的一些结果，当然这个结果的展示也是一步步铺开的，有时候也是为了选择一些更优的方法而去比较得到结果。

首先看Result的第一部分，采用的数据集是酿酒酵母的基因组，然后生成的没有错误的模拟数据集，通过每次变化reads的长度来看拼接结果的好坏。如何衡量拼接结果呢？这里引入一个衡量标准，N50 contig size。什么意思呢，contig是重叠群的意思，也就是测序结束的reads片段拼接之后形成的一段段的序列（因为whole genome现在的技术还不能完全拼接，所以结果就是一段段的contig）；而N50的意思，可以根据上图理解，这是我在网络上找到的别人的图片链接在这里N50 contig解释，什么意思呢，我们得到contig以后先把他们按照从大到小的顺序排序，也就是1a图，然后我们从最大的开始往后累加，直到长度到达所有contig总长度的一半，我们最后添加上的这个contig长度就呗称为N50 contig size途中可以看出这个值等于60。我们也就是用N50来评判contigs的一个平均长度，然后用来衡量拼接结果的优劣。

上面这个图表示：作者比较了两种现成的软件他们基于不同的数学原理，the OLC assembler Celera Assembler11 和 the de Brujin assembler SOAPdenovo. 图中紫色的小星星表示基准，是用来自酵母菌的50× of real 76-bp Illumina paired-end (300 bp) reads数据使用SOAPdenovo assembly 拼接的结果。（简单解释一下这个数据的意思，50×是指的测序深度前文已经写过了，Illumina paired-end是这个公司的一种测序方法，对每一段序列两端测序，这里300bp就是两端各测150bp，一共300bp。76bp是指的最终测序的reads长度的中值。）
实际的效果是图中的几条曲线，模拟了一次10×的测序但是没有错误，随机打乱的长度分布跟别根据了PacBio的几种仪器(Q1, 2011)(Q2, 2011)和(Q1, 2012)。所以最后作者选择了Celera Assembler11进行下面的工作，因为他随着reads长度增加测序的效果会上升，而SOAPdenovo变化不明显。

我们再看一下结果的第二部分，数据集用了三个lambda噬菌体，大肠杆菌和酿酒酵母。

首先，我们对不同测序深度（short reads）的大肠杆菌基因组进行拼接，然后作者认为50×是一个最优的选择，测序深度再增加，矫正结果准确度的变化不太大了。
然后这三个数据集的结果如表格所示，矫正后的序列和参考序列进行了比较。使用了一个现成的软件Nucmer 3.23 (ref. 50). 只有reads大于500 bp才被考虑。
%TP (通量),校正后的非嵌合、正确修剪序列中未校正的原始碱基的百分比;
%Idy (同一性), 对比参考序列平均的准确度;
%Cov (coverage),以参考序列为标准，得到平均的覆盖度;
%Chimer, 矫正的碱基分开匹配到参考序列上，也就是嵌合体;
%Trim,匹配到参考序列上，匹配小于99.5% 长度的比例。

左侧图这个结果展示了不同测序结技术对大肠杆菌基因组的拼接结果。
右侧图展示了PBcR（矫正后的long reads）的深度增加，N50所体现的拼接效果性能的显著提升。但是10×已经可以到达一个较好的结果。

这一部分结果，作者拼接了鹦鹉的基因组。但是由于在当时鹦鹉没有可供参考的基因组数据，所以作者提供了一个网址。这个网址打开是一个拼接算法的比赛，我猜测作者也参加了这个比赛，所以才会在结果这里说自己做了鹦鹉的基因组拼接。然后作者怎么衡量拼接的结果呢，就是把这个比赛里面除了自己别人拼接的鹦鹉数据库作为模板，然后把自己矫正过的long read去匹配，然后99.6%的reads至少一次匹配上去了，93.7%的全长匹配上并且相似度在99.5%。

选用了玉米的幼苗样品，拼接结果错误率很低。

这个图示化的结果如上图，可以看到这个物种还有两个亚型。矫正过的结果有很大的改善。