Pacbio SMRT君的逆袭—算法篇（三代测序那些事 第四期）

可靠消息，国内某测序巨头准备引进Pacbio SMRT测序平台PacBio RSII系统开始进军三代测序市场，按照惯例，其他小巨头也会紧随其后，所以，小编深深觉得，今年国内肯定要刮起一阵不小的三代测序风，三代测序市场的竞争也将变得激烈起来。

从09年PacBio RS首批测试数据公布暴露出15%的原始错误率（主要为InDel错误）后开始遭人诟病到13年PacBio SMRT首次进入中国市场（Nextomics首家推出）而又普遍不被看好再到如今被国内主要NGS测序公司争相追捧，在小编看来，这算得上是PacBio SMRT君的一次华丽逆袭。

这期，小编将为大家盘点那些让PacBio SMRT君华丽逆袭的那些算法。

显然，直接使用原始错误率为15%且大多数为InDel错误的Reads进行基因组拼装是不可行的，因为大多数基因组组装软件所能忍受的Reads上限错误率<~5%-10%。因此对PacBio RS平台所产生的Reads（错误随机分布）进行比对 & 校正成为了此类数据应用于基因组组装的第一步。

然而大多数比对软件主要是针对高准确率短读长的NGS测序数据设计，比如 SOAP、Bowite、BWA、Maq、SHRIMP、ELAND等，无法对这种读长为数kb、原始错误率15%的三代数据进行比对，虽然在BWA基础上修改得到的可容忍较高错误率且能够进行长读长数据进行比对的BWA-SW却有着比对率不高的缺点。

以下两款专门针对PacBio数据的新型比对软件的出现改变了PacbioSMRT君的这一窘境，也开启了PacBio君的逆袭之路。

NO.1 BLASR （主要针对微生物基因组）

BLASR是基于经典动态规划思想设计的局部比对软件，使用了BWT（Burrows-wheeler Transform）格式的索引结构进行匹配区域快速定位，对于候选匹配区域与最终匹配区域的确定使用了运算速度更快的稀疏动态规划算法（Sparse Dynamic Programing）。

加利福尼亚大学的软件开发人员首先使用产生于大肠杆菌E.coli OH104：H4基因组的48X的PacBio数据（10.7% insertion、4.3% deletion、0.9% substitution）对BLASR的mapping 率、mapping速度进行了评估，结果显示，BLASR mapping率为90%，高于BWA-SW的50%，运行时间为20min 54 S ，小于BWA-SW的434 min 5S。

为了评估该软件的mapping准确率，软件开发人员使用了大肠杆菌基因组的Pacbio 模拟数据对基因组进行mapping，结果显示90%以上的reads mapping 准确度在99.99%以上。见图8^[1]。

        BLASR凭借其优秀的mapping品质与多款Consensus算法软件（AMOS、Quiver、PBDAG-Con等）一起构成了多款后来出现的针对Pacbio数据的基因组组装软件的核心校正算法，这些软件包括二三代混合组装软件PBcR-BLASR、AHA，纯三代组装软件HGAP。

其中纯三代组装软件HGAP（Hierarchical Genome-Assembly Process）是一款基于分级组装思想的基因组组装软件。其大致流程为：

1）挑选较长的reads作为seed reads（>6kb）

2）使用BLASR将较短的的reads mapping 的seed reads上，使用PBDAG -Con一致性算法对reads进行校正并进行预组装

3）使用CA算法对预组装得到的准确率较高的长Reads进行组装

4）将原始reads mapping回组装好的基因组，使用新型一致性算法Quiv er对所得基因组进行进一步校正，最终得到准确率大于99.9999%（QV60） 的高质量的微生物基因组图谱。大致流程见图1。

        Stephen Turner等人分别使用了100X、90X、100X的大肠杆菌、栖热菌、肝素黄杆菌PacBio数据对HGAP的组装效果进行了评估，组装结果中Contig数分别为2、3、1^[1]。

随着Pacbio SMRT的读取长度、运行通量的增加，之前出现的二三代混合组装软件相较于纯三代组装软件HGAP，无论是在测序成本还是组装指标均没有优势，因此，这种二三代混合组装的策略在微生物基因组组装中慢慢被淘汰，对于之前提到的那些二三代混合组装小编就不做过多赘述^[3][4][5]。

NO.2 MinHash（针对大型动植物基因组）

MinHash也是一款基于经典的动态规划思想设计局部序列比对软件，与BLASR不同的是，它采用了最小哈希算法（MinHash）实现了匹配区域的快速定位，该过程如图1所示。

        大致流程为：1）reads Kmer化， 2）将Kmer通过 Hash方程转化为整数格式的fingerprints，3）挑选各自fingerprint最小的Kmer组成用于比对的Kmer集合Sketch 4）使用Jaccard相似度计算Kmer相似度 5）若相似度超过阈值，则返回基因组区域使用动态规划算法详细比对 6）找出匹配区域。

包括Pacbio 的Chen-Shan Chin在内的软件开发人员使用拟南芥、果蝇、人类的PacBio测序数据评估了MinHash算法在大型基因组测序数据比对过程中的性能。结果显示Mapping率在80%左右，准确率均在90%以上，而运行时间仅为15-21 CPU h ，而另一款主要应用于微生物基因组数据比对的BLASR的在大型基因组测序数据比对时运行时间高达上百 CPU h [6]。

之后软件开发人员将MinHash结合新型Consensus软件FalconSense，再整合到OLC组装算法软件Celera Assembler（CA）中得到了大型基因组纯三代组装算法PBcR-MinHash。软件开发人员分别使用了121X、144X、54X的果蝇、拟南芥、人类葡萄胎的Pacbio数据对PBcR-MinHash的组装性能进行了评估，三个物种的ContigN50分别达到了11Mb、20Mb、4Mb。

至此，PacBio SMRT君通过内在的修为弥补了表面缺陷，完成了自己的逆袭之路。

未来组生物（Nextomics Biosciences）基于HGAP已经推出了多款微生物基因组完成图产品，在动植物基因组方面也成功召开了基于三代测序技术的杜仲基因组新闻发布会，感兴趣的小伙伴可以电话或邮箱联系我们。

Paper：

[1] Mark J Chaisson et al. Mapping single molecule sequencing readsusing basic local alignment with successive refinement (BLASR): application andtheory. BMC Bioinformation . 2012

[2] Stephen W Turner et al. Nonhybrid, finished microbial genomeassemblies from long-read SMRT sequencing data. Nature Mehods. 2013

[3] Sergey Koren et al. Hybrid error correction an de novo assembly ofsingle molecule sequencing reads. NatBiotechnol . 2012

[4] Ali Bashir et al. A Hybrid Approach for the Automated Finishing ofBacterial Genomes. Nat Biotechnol .2012

[5] Filipe J Ribeiro et al. Finished bacterial genomes from shotgunsequence data. Genome Research. 2012

[6] Konstantin Berlin et al. Assembling large genomes with singlemolecule sequencing and locality sensitive hasing. Bio Rxiv. 2014

声明：本文原创，若转载，请注明来源