CN113300720A - 长dna序列存储的插入删节分段识别方法 - Google Patents

Abstract

本发明公开了长DNA序列存储的插入删节分段识别方法。针对由叠加水印的长DNA分子经过三代测序、组装得到的存在插入删节错误的长序列在采用硬判决纠错时性能存在损失，以及采用软判决算法时，状态转移网格图过大导致中间度量计算复杂度高问题，首先对长序列建立隐马尔科夫模型，利用硬判决前向‑后向算法计算每个碱基距离原始位置的偏移，给每个碱基建立偏移索引；然后对水印序列分段并按照分段后的水印序列边界的偏移索引将读出序列划分为若干长度不同的片段；最后依次对短序列片段采用软判决前向‑后向算法估计符号似然概率，完成纠错。本发明提出的分段式纠错方法可以有效地减少网格图的大小，降低纠错复杂度。

Description

长DNA序列存储的插入删节分段识别方法

技术领域

本发明属于利用DNA的数据存储领域，尤其涉及长DNA序列存储的插入删节分段识别方法。

背景技术

随着信息技术的飞速发展，全球产生的数据量快速增长。然而，当前存储技术的发展速度严重滞后于数据量的快速增长，并且光盘、硬盘和磁带等传统标准存储介质寿命有限，维护成本高，难以满足日益增长的数据存储需求。随着合成技术、测序技术以及组装技术的发展，合成DNA凭借密度高，能耗低，介质保存时间久等优点为海量数据的档案存储提供了另一种选择，引起了研究者们的广泛兴趣。近年来，研究者尝试将信息数据存储在DNA中，并证明其在“冷”数据存储方面可行性。然而，由于在DNA合成、样本扩增、测序以及最终碱基识别中的缺陷，多种类型的错误会对测序读段造成破坏，包括核苷酸插入、删除和替换错误。

DNA信道中引入的错误过程与通信系统中的信息传输过程遇到错误的过程类似。例如，高速通信系统中的时钟抖动、偏移造成插入、删节错误；无线光通信系统中的差分脉冲位置调制(Differential Pulse-Position Modulation，DPPM)系统中，码片因噪声发生跳表引起插入、删节；高密度磁存储系统中的比特图形化介质(Bit-Patterned Media，BPM)存储系统中，介质缺陷、读写电路不完善或时钟抖动、磁头震动问题都会造成插入和删节的发生。最常见的DNA合成错误是单碱基的删节，在大规模并行寡核苷酸合成过程中，替代和插入错误也很常见。Heckel等人通过分析之前的研究结果，发现错误主要来自合成和测序，也受到存储DNA降解和聚合酶链反应(polymerase chain reaction，PCR)的影响，在测序时，替代错误比删节和插入更有可能发生。目前使用最广泛的DNA测序平台是Illumina测序仪，它是基于图像处理和合成测序的概念，其本身存在10^-3～10^-4的错误率。另一种正在快速发展的DNA测序方法是纳米孔测序技术，例如第三代高通量测序：牛津纳米孔(OxfordNanopore Technology，ONT)测序中，利用DNA分子通过纳米孔引起孔内电流变化的电信号来确定DNA链的核苷酸序列，其测序读段长，速度快，无需PCR扩增过程且极具便携性，受到越来越广泛的重视，但其读段精度较低，导致测序后的数据中错误率较高，达到10％-15％。

一条含有4个天然核苷酸的DNA链上每个DNA字符最多可编码2比特。在最大码率下，信息中没有冗余，不能进行错误校正。然而，DNA合成和测序过程引入了错误，需要有效的错误纠正码来保证信息的可靠性恢复。纠错码降低了码率，但在将信息编码为DNA字符，以及读取端将DNA字符解码为信息位时，纠错码对于防止错误是必要的。

目前，数字通信领域的几种重要纠错码已经在DNA数据存储中进行了尝试。Grass等人在DNA序列设计方面提出使用里德-所罗门(Reed-Solomon,RS)码，利用逻辑冗余纠正碱基删除与随机错误。Elich等人提出了喷泉编码(纠删码)的概念，用于纠正寡核苷酸分子丢失造成的删除错误，其编码效率接近香农容量，同时在抵抗数据损坏方面具有很高的鲁棒性。最近一项研究中采用低密度奇偶校验(LDPC)码与RS码构成的乘积码纠正删除与随机错误。William等人采用哈希编码，贪婪穷举搜索解码的纠错码(HEDGES)作为内码，外码采用RS编码，HEDGES编码用于纠正插入删节错误，RS编码用于纠正残余错误。测试与仿真结果表明，从错误高达10％的受损DNA中也可以无错恢复千兆字节的数据。Organick等人同样采用RS编码用于纠错，并且译码之前采用聚类、一致性序列合并算法对测序读段进行合并纠错，消除了插入删节错误，实现数据可靠读出。在此基础上，该研究组进行了深入研究，将寡核苷酸池的短序列进行了组装，以支持三代纳米孔测序仪，设计并验证了一种采用寡核苷酸池存储的先组装后三代测序的策略，初步验证了纳米孔测序用于DNA数据存储读出的可行性。

插入/删除错误将导致发送序列与读出序列长度不一致，传统的针对替代错误的纠错码技术将不再适用。研究者提出了一系列性能优良的插入/删节纠错码，大致可分为四类：(1)可恢复同步的标记(Marker)码；(2)以Varshamov-Tenengolts(VT)码为代表的代数构造码；(3)基于卷积码的同步纠错码；(4)基于水印码的级联码。级联码是一种新型的插入/删除纠错码，它将检测同步的码与传统纠错码相结合，由同步码定位插入/删除，由纠错码纠正残留插入/删除及信道替代错误。针对插入/删除随机出现的情况，Davey和MacKay提出了一种具有纠正插入、删除能力的高效级联编码方案(DM构造)，DM构造将水印码与纠正替代错误的纠错码相结合，由水印码估计每个符号发生插入、删除的概率，并由替代错误纠错码纠正残留插入/删除及信道替代错误。Liu和Chen提出结合内外译码算法的迭代译码方案，可以提高系统的纠错能力。为降低水印译码算法的运算量，提出一种基于自适应删剪网格图的低复杂度水印译码算法。Kracht和Schober在DM级联码构造基础上修改和扩展了水印码的概念，将其用于DNA测序，证明水印概念在DNA测序应用中的可行性。Chen和Wang将分组纠错码与预先确定的伪随机序列相结合，生成一组用于标记不同样本的碱基序列，针对译码提出了一种软判决识别方法，可以有效地纠正标记序列中的多个错误，但该方法适用于短码长的条形码序列barcode，不适合三代测序数据中的长序列。

一项新的研究根据三代纳米孔测序的高错误率，重点考虑难以处理的碱基插入/删节错误，采用LDPC码叠加伪随机序列构造的水印码，设计了可纠正严重插入删节错误的高效编码方案，结合自主复制序列(ARS)，从头编码设计合成了一条长度为254,886bp专用于数据存储的酵母人工染色体。在读出方面，利用三代纳米孔测序器件实现了碱基的快速读出与无错恢复。碱基识别后的测序读段错误率高于10％，经过组装、定位后序列仍包含插入删节错误，为处理这些插入删节错误，设计了一个基于叠加水印的前向-后向硬判决纠错算法，该算法基于隐马尔科夫模型，利用前向与后向度量估计偏移路径，并通过最大概率估计每个碱基位置与初始位置的偏移，从而识别插入与删节错误，但硬判决未考虑发送端的稀疏码本，其纠错性能有限，不能达到最优效果。基于软判决的纠错算法被提出，软判决算法基于稀疏码本对每个符号对应的对数似然概率进行计算，相比于硬判决算法，增加了中间度量计算，由于采用了递归运算，复杂度较高。当所设计的插入/删除纠错码的码长较长，插入/删除信道环境变得非常恶劣时(插入删节错误数目不同，造成偏移路径较大)，为保证水印译码器的同步能力，前向-后向算法译码网格图中的状态数必须随之大量增加，与之对应的译码网格图的长度也需要增加。以上情况将导致水印译码算法的复杂度升高，进而增加了系统的译码延时。

发明内容

本发明提供了长DNA序列存储的插入删节分段识别方法，本发明主要解决由叠加水印的长DNA分子经过三代测序、组装得到的存在插入删节错误的长序列在采用硬判决纠错时性能存在损失，以及采用软判决算法时，状态转移网格图过大导致中间度量计算复杂度高的问题，详见下文描述：

长DNA序列存储的插入删节分段识别方法，具体包括以下步骤：

(1)针对由叠加水印的长DNA分子测序得到的长度为

的存在插入、删节与替代错误的读出序列r，基于水印序列和错误传输模型建立隐马尔科夫模型，确定隐藏状态和观测向量，建立状态转移网格图；

(2)改进观测向量的输出概率公式，基于网格图运行前向传递和后向传递，利用硬判决算法估计每个碱基位置的最大可能概率，推断每个碱基距离原始位置的偏移；

(3)对水印序列按照固定长度分段，给每个水印短片段的边界位置建立偏移索引，并基于索引计算与其对应的读出序列边界位置，将读出序列划分为若干长度不同的数据片段；

(4)对分段后的若干测序读出片段，依次递归计算中间度量，结合前向-后向软判决算法对向量d上的符号似然函数l进行估计，根据估计值作硬判决处理，然后与水印序列相异或，得到纠正后的片段，进而拼接为长DNA序列。

所述针对叠加水印的长DNA分子测序得到的长度为

的存在插入删节与替代错误的读出序列r，基于水印序列和错误传输模型建立隐马尔科夫模型，确定隐藏状态和观测向量，建立状态转移网格图，具体的步骤为：

(1.1)根据三代测序的错误特性与前期测序实验的统计数据的积累，估计读出序列的插入错误概率P_i，删除错误概率P_d，替代错误概率P_s以及传输概率P_t＝1-P_i-P_d；

(1.2)定义第i时刻的碱基偏移量x_i为从发送碱基t₀至待发送碱基t_i间存在的插入数目减去删节数目，将碱基偏移量x_i作为隐马尔科夫模型的隐藏状态，将碱基最大偏移量限定为x_max，以降低算法的计算量，x_max定义为

N为设计的长DNA序列的初始长度，x_i的取值范围为X＝{-x_max,...,-1,0,1,...,x_max}；

(1.3)定义测序读出长DNA序列r为r＝r ⁽⁰⁾ r ⁽¹⁾...r ^(N-1)，其子序列为

将r作为观测向量，子序列r ⁽ⁱ⁾的传输过程由x_i-1→x_i的状态转移决定；

(1.4)在状态x_i＝a转变为状态x_i+1＝b时，对应的输出序列被定义为u，根据映射规则{(00)→Α，(01)→T，(10)→G，(11)→C}将读出序列解映射为比特对u₁和u₂，双层水印序列分别为w_1i和w_2i，在水印序列上的噪声包含两部分：信息序列m经编码后，得到符号序列d，经过稀疏器后生成的稀疏码字s及过信道产生的替代错误P_s，故定义有效替代概率P_f为P_f＝f(1-P_s)+(1-f)P_s，其中f为长度为n的稀疏码的平均密度；

(1.5)根据观测向量的传输过程建立隐马尔科夫模型状态转移网格图，网格图中的每个节点对应第i个位置的偏移量x_i，将(x₀,x₁,…,x_i,…,x_N-1)表示为隐马尔科夫模型的状态序列，定义I为最大的连续插入错误数目，每个发送碱基对应一个长度从0到I+1的读出碱基序列，即当x_i＝a时，x_i+1＝b的取值范围是{a-1,…,a+I}。

所述改进观测向量的输出概率公式，基于网格图运行前向传递和后向传递，利用硬判决算法估计每个碱基位置的最大可能概率，推断每个碱基距离原始位置的偏移，具体的步骤为：

(2.1)由上一状态x_i＝a转移到下一状态x_i+1＝b存在两种情况，第一种，经过信道时，t_i前被插入(b-a+1)个碱基且发送碱基t_i丢失，第二种，经过信道时，t_i前被插入(b-a)个碱基且t_i被传输，基于该传输情形，定义状态转移概率P_a,_b＝P(x_i+1＝b|x_i＝a)，计算公式为，

其中，α_I等于1/(1-(P_i)^I)；

(2.2)对双层错误概率改进，发生插入错误时，插入碱基有4种选择(A/T/G/C)，每个碱基的插入概率为1/4，发生替代错误时，映射后的双层比特存在两种情况，(a)上下两层比特均不同，类型包括{A(00)→C(11)，T(01)→G(10)，C(11)→A(00)，G(10)→T(01)}，(b)只有1个比特不同，类型包括{A(00)→T(01)，G(10)→C(11)，T(01)→C(11)，A(00)→G(10)，T(01)→A(00)，C(11)→G(10)，C(11)→T(01)，G(10)→A(00)}，两种情况的错误类型占比1:2，定义有效替代错误概率分别为

和

定义改进的观测向量的输出概率

计算公式为，

其中，该式满足∑ _u Q_i,a,b(u)＝1，

是与水印比特w_1i相关联的读出比特，

是与水印比特w_2i相关联的读出比特；

(2.3)定义前向度量F_i(y)＝P(r₀,…,r_i-1+y,x_i＝y),y∈X，F_i(y)表示x_i＝y且接收到前i-1+y个碱基的概率，定义后向度量B_i(y)＝P(r_i+y,…|x_i＝y),y∈X，B_i(y)表示x_i＝y的条件下输出碱基串(r_i+y,…)的概率，使用分支度量P_a,y·Q_i-1,a,y递归运行前向传递与后向传递，估计每个碱基位置的最大可能概率F_i(y)B_i(y)，推断每个碱基距离原始位置的偏移，其中，第i时刻状态y的前向度量计算公式表示为，

同样的，第i时刻状态y的后向度量计算公式表示为，

所述对水印序列按照固定长度分段，给每个水印短片段的边界位置建立偏移索引，并基于索引计算与其对应的读出序列边界位置，将读出序列划分为若干长度不同的数据片段，具体的步骤为：

(3.1)将长度为N水印序列分为若干长度为N₁的定长片段，N₁满足可以被稀疏码长度n整除且为固定值，然后将长度为

的读出序列分为长度为N₂的若干片段，N₂长度可变，用t_i表示水印片段的起始节点，则水印片段的边界节点为

假定r_i表示读出序列片段的起始节点，读出序列的边界节点为r_border，其分段边界计算公式表示为，

其中，按照三代测序的特点，N与

不一定等长，则下一个水印分段片段的起始节点为

下一个读出序列分段的起始节点为r_border+1，依次计算，得出N/N₁个片段的起始节点和边界节点，与每个边界位置t_borde的最大概率

相对应的最大可能偏移量max_drift_border的计算公式为，

max_drift_border＝argmax(F_border(y)B_border(y))，

(3.2)初始化前向度量与后向度量，假定起始节点固定，初始化

边界位置的后向度量的初始化由水印片段长度N₁和读出序列片段可变长度N₂之间的大小关系判定，

i)当N₂小于N₁时，

ii)当N₂等于N₁时，

iii)当N₂大于N₁时，

所述对分段后的若干测序读出片段，依次递归计算中间度量，结合前向-后向软判决算法对向量d上的符号似然函数l进行估计，根据估计值作硬判决处理，然后与水印序列相异或，得到纠正后的片段，进而拼接为长DNA序列，具体的步骤为：

(4.1)向量d上的符号d_i与长度为n的稀疏向量(s_i,0,…,s_i,n-1)一一对应，固定d_i＝λ，执行

到

的前向传递，定义中间度量为，

其中符号d_i等概地选取GF(q)＝{0,1,…,λ,…,q-1}中的值，r ⁰表示接收比特

为符号d_i对应稀疏串起始位置的同步偏移状态，

为d_i对应稀疏码第k个位置的同步偏移状态；

(4.2)进一步，利用递归公式，第k时刻的中间度量计算公式为，

其中，N_L＝N/n，并初始化M₀(z,d_i)＝1；

(4.3)在硬判决算法的基础上增加中间度量

来获取向量d上的似然函数l，似然函数l输出概率公式可表示为，

和

分别为

时刻的前向度量和

时刻的后向度量，对映射后的上下两层比特单独进行纠正，水印序列和稀疏码表已知，所有的替代错误仅来自于信道，将w_j替换为

P_f替换为P_s，输出概率

可表示为，

(4.4)依次对符号d_i对应的输出概率进行计算，选取概率最大的符号λ所对应的稀疏码(s_i,0,…,s_i,n-1)与本地水印序列w_1i和w_2i异或，得到纠正后的片段，进而组合为长DNA序列。

本发明提供的技术方案的有益效果是：

本发明针对基于隐马尔科夫模型的前向-后向硬判决纠错算法的纠错性能存在损失，采用软判决算法时，因插入删节错误数目不同造成前向-后向算法的偏移状态累加，网格图过大导致纠错复杂度过高的问题，发明了长DNA序列存储的插入删节分段识别方法。具体而言，首先基于已知的水印序列与估计的错误概率建立隐马尔科夫模型，利用硬判决前向-后向算法计算每个碱基距离原始位置的偏移，给每个碱基建立偏移索引；然后对水印序列分段并按照分段后的水印序列边界的偏移索引将读出序列划分为若干长度不同的片段；最后依次对短序列片段采用软判决前向-后向算法估计符号似然概率，完成纠错。本发明提出的译码方案对于纠正DNA存储中插入、删节和替代错误具有明显性能增益，与传统软判决算法相比，分段式处理策略将长序列分为若干片段序列，使插入/删节错误平均分布在各个片段，降低纠错所需的偏移量，减少偏移网格图的大小降低了纠错复杂度，减少纠错时延。

附图说明

图1为长DNA序列存储的插入删节分段识别方法流程图；

图2为本发明建立的隐马尔科夫模型状态转移网格图；

图3为本发明提供的计算碱基偏移索引的流程图；

图4为本发明提供的读出序列分段的方法示意图；

图5为本发明提供的分段式软判决算法纠错流程图；

图6为本发明提供的读出序列原始错误分布示意图；

图7为本发明提供的软判决前后向纠错算法的最大可能路径示意图。

图8为本发明提供的6条长测序读出序列分别由软判决算法与硬判决算法纠错后的汉明距离误码率对比示意图。

图9为本发明提供的分段式软判决算法与传统软判决算法的复杂度对比示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面对本发明的实施方式作进一步地详细描述。

为了纠正DNA存储中的插入、删节和替代错误并降低纠错的复杂度，本发明提供了长DNA序列存储的插入删节分段识别方法，基本思想是：基于叠加的水印序列和插入删节错误传输模型建立隐马尔科夫模型来估计测序读出序列上的插入与删节错误位置；然后采用分段式策略对长DNA序列进行分段，以降低因插入删节错误数目不同造成前向-后向算法的网格状态偏移较大，导致纠错复杂度过高的问题，最后利用软判决前向-后向算法纠正错误，并输出码字估计值。

下面结合附图对本发明的实施方式做出详细说明。

本实施方式针对本发明提出的基于长序列的插入删节错误的分段式纠错方法进行详细描述，参见图1，具体包括以下步骤：

(1)针对由叠加水印的长DNA分子测序得到的长度为

的存在插入、删节与替代错误的读出序列r，基于水印序列和错误传输模型建立隐马尔科夫模型，确定隐藏状态和观测向量，建立状态转移网格图；

(2)改进观测向量的输出概率公式，基于网格图运行前向传递和后向传递，利用硬判决算法估计每个碱基位置的最大可能概率，推断每个碱基距离原始位置的偏移；

(3)对水印序列按照固定长度分段，给每个水印短片段的边界位置建立偏移索引，并基于索引计算与其对应的读出序列边界位置，将读出序列划分为若干长度不同的数据片段；

(4)对分段后的若干测序读出片段，依次递归计算中间度量，结合前向-后向软判决算法对向量d上的符号似然函数l进行估计，根据估计值作硬判决处理，然后与水印序列相异或，得到纠正后的片段，进而拼接为长DNA序列。

所述步骤(1)针对叠加水印的长DNA分子测序得到的长度为

的存在插入删节与替代错误的读出序列r，基于水印序列和错误传输模型建立隐马尔科夫模型，确定隐藏状态和观测向量，建立状态转移网格图，具体的操作为：

(1.1)根据三代测序的错误特性与前期测序实验的统计数据的积累，估计读出序列的插入错误概率P_i，删除错误概率P_d，替代错误概率P_s以及传输概率P_t＝1-P_i-P_d；

(1.2)定义第i时刻的碱基偏移量x_i为从发送碱基t₀至待发送碱基t_i间存在的插入数目减去删节数目，将碱基偏移量x_i作为隐马尔科夫模型的隐藏状态，将符号最大偏移量限定为x_max，以降低算法的计算量，x_max定义为

N为设计的长DNA序列的初始长度，x_i的取值范围为X＝{-x_max,...,-1,0,1,...,x_max}；

(1.3)定义测序读出长DNA序列r为r＝r ⁽⁰⁾ r ⁽¹⁾...r ^(N-1)，其子序列为

将r作为观测向量，子序列r ⁽ⁱ⁾的传输过程由x_i-1→x_i的状态转移决定；

(1.4)在状态x_i＝a转变为状态x_i+1＝b时，输出序列被定义为u，根据映射规则{(00)→Α，(01)→T，(10)→G，(11)→C}，读出序列解映射为比特对u₁和u₂，双层水印序列分别为w_1i和w_2i，在水印序列上的噪声包含两部分：信息序列m经编码后，得到符号序列d，经过稀疏器后生成的稀疏码字s及过信道产生的替代错误P_s，故定义有效替代概率P_f为P_f＝f(1-P_s)+(1-f)P_s，其中f为长度为n的稀疏码的平均密度；

(1.5)根据信道模型中存在的所有传输情形建立隐马尔科夫模型状态转移网格图，如图2所示，网格图中的每个节点对应第i个位置的偏移量x_i，将(x₀,x₁,…,x_i,…,x_N-1)表示为隐马尔科夫模型的状态序列，定义I为最大的连续插入错误数目，每个发送碱基对应一个长度从0到I+1的读出碱基序列，即当x_i＝a时，x_i+1＝b的取值范围是{a-1,…,a+I}，本发明中I设置为2。

其中，如图3所示，步骤(2)改进观测向量的输出概率公式，基于网格图运行前向传递和后向传递，利用硬判决算法估计每个碱基位置的最大可能概率，推断每个碱基距离原始位置的偏移，具体的操作为：

(2.1)由上一状态x_i＝a转移到下一状态x_i+1＝b存在两种情况，第一种，经过信道时，t_i前被插入(b-a+1)个碱基且发送碱基t_i丢失，第二种，经过信道时，t_i前被插入(b-a)个碱基且t_i被传输，基于该传输情形，定义状态转移概率P_a,b＝P(x_i+1＝b|x_i＝a)，计算公式为，

其中，α_I等于1/(1-(P_i)^I)；

(2.2)对双层错误概率改进，发生插入错误时，插入碱基有4种选择(A/T/G/C)，每个碱基的插入概率为1/4，发生替代错误时，映射后的双层比特存在两种情况，(a)上下两层比特均不同，类型包括{A(00)→C(11)，T(01)→G(10)，C(11)→A(00)，G(10)→T(01)}，(b)只有1个比特不同，类型包括{A(00)→T(01)，G(10)→C(11)，T(01)→C(11)，A(00)→G(10)，T(01)→A(00)，C(11)→G(10)，C(11)→T(01)，G(10)→A(00)}，两种情况的错误类型占比1:2，定义有效替代错误概率分别为

和

定义改进的观测向量的输出概率

计算公式为，

其中，该式满足∑ _u Q_i,a,b(u)＝1，

是与水印比特w_1i相关联的读出比特，

是与水印比特w_2i相关联的读出比特；

(2.3)定义前向度量F_i(y)＝P(r₀,…,r_i-1+y,x_i＝y),y∈X，F_i(y)表示x_i＝y且接收到前i-1+y个碱基的概率，定义后向度量B_i(y)＝P(r_i+y,…|x_i＝y),y∈X，B_i(y)表示x_i＝y的条件下输出碱基串(r_i+y,…)的概率，使用分支度量P_a,y·Q_i-1,a,y递归运行前向传递与后向传递，估计每个碱基位置的最大可能概率F_i(y)B_i(y)，推断每个碱基距离原始位置的偏移，其中，第i时刻状态y的前向度量计算公式表示为，

同样的，第i时刻状态y的后向度量计算公式表示为，

其中，如图4所示，步骤(3)对水印序列按照固定长度分段，给每个水印短片段的边界位置建立偏移索引，并基于索引计算与其对应的读出序列边界位置，将读出序列划分为若干长度不同的数据片段，具体的操作为：

(3.1)将长度为N水印序列分为若干长度为N₁的定长片段，N₁满足可以被稀疏码长度n整除且为固定值，然后将长度为

的读出序列分为长度为N₂的若干片段，N₂长度可变，用t_i表示水印片段的起始节点，则水印片段的边界节点为

假定r_i表示读出序列片段的起始节点，读出序列的边界节点为r_border，其分段边界计算公式表示为，

其中，按照三代测序的特点，N与

不一定等长，则下一个水印分段片段的起始节点为

下一个读出序列分段的起始节点为r_border+1，依次计算，得出N/N₁个片段的起始节点和边界节点，与每个边界位置t_borde的最大概率

相对应的最大可能偏移量max_drift_border的计算公式为，

max_drift_border＝argmax(F_border(y)B_border(y))，

(3.2)初始化前向度量与后向度量，假定起始节点固定，初始化

边界位置的后向度量的初始化由水印片段长度N₁和读出序列片段可变长度N₂之间的大小关系判定，

i)当N₂小于N₁时，

ii)当N₂等于N₁时，

iii)当N₂大于N₁时，

所述对分段后的若干测序读出片段，依次递归计算中间度量，结合前向-后向软判决算法对向量d上的符号似然函数l进行估计，根据估计值作硬判决处理，然后与水印序列相异或，得到纠正后的片段，进而拼接为长DNA序列，如图5所示，具体的步骤为：

(4.1)向量d上的符号d_i与长度为n的稀疏向量(s_i,0,…,s_i,n-1)一一对应，固定d_i＝λ，执行

到

的前向传递，定义中间度量为，

其中符号d_i等概地选取GF(q)＝{0,1,…,λ,…,q-1}中的值，r ⁰表示接收比特

i_-＝n×i，i₊＝n×(i+1)，

为符号d_i对应稀疏串起始位置的同步偏移状态，

为d_i对应稀疏码第k个位置的同步偏移状态；

(4.2)进一步，利用递归公式，第k时刻的中间度量计算公式为，

其中，N_L＝N/n，并初始化M₀(z,d_i)＝1；

(4.3)在硬判决算法的基础上增加中间度量

来获取向量d上的似然函数l，似然函数l输出概率公式可表示为，

和

分别为

时刻的前向度量和

时刻的后向度量，对映射后的上下两层比特单独进行纠正，水印序列和稀疏码表已知，所有的替代错误仅来自于信道，将w_j替换为

P_f替换为P_s，输出概率

可表示为，

(4.4)依次对符号d_i对应的输出概率进行计算，选取概率最大的符号λ所对应的稀疏码(s_i,0,…,s_i,n-1)与本地水印序列w_1i和w_2i异或，得到纠正后的片段，进而组合为长DNA序列。

具体实施例

下面给出具体的实施例，说明本发明给出的一种针对长DNA序列的插入删节错误的分段式纠错方法的可行性。

本发明对基于纳米孔测序的DNA储存方案文献中测序得到的长DNA序列进行纠错，验证本发明的可行性。通过比较在DNA序列映射成的双层比特码字在纠正译码后的汉明距离的误比特率，来验证发明的方法的纠错性能。其中，在编码端，数字信息经过LDPC编码，经过4比特转5比特的稀疏化后(即n＝5)，得到数据序列，进而与伪随机序列构成的水印码结合并映射为碱基序列。接收端，设计的长DNA序列经过三代纳米孔测序后，得到测序数据，进而经过组装、定位得到6条DNA序列。在本发明中使用的6条DNA序列分别命名为pic1，pic2，video1，video2，video3，video4，其中pic2初始设计长度为40320bp，其余序列全部设计为40500bp。首先根据测序得到的6条DNA序列的错误特性建立契合提出的分段式前向-后向软判决算法的隐马尔科夫模型。图6统计了6条接收的DNA序列的初始错误概率分布，6条序列呈现统一的错误分布特性：删除错误最多，替代错误其次，插入错误较少，且删除错误概率是插入错误概率的4倍左右，总的错误率低于千分之六。根据初始概率信息，设置隐马尔科夫模型中插入错误概率P_i＝0.001，删除错误概率P_d＝0.004，替代错误概率P_s＝0.001，并且设置最大偏移量x_max＝150。用于分段的水印固定长度N₁设置为1000，最大连续插入数目设置I＝2。

为了清晰的描述纠错算法所需最大漂移量x_max的选取原则和分段式处理策略在减少网格图大小，降低复杂度方面的贡献。本发明分析了用于纠正插入/删节的水印译码网格图的特点，首先记录传统纠错过程中每个时刻中前向/后向度量值相乘概率最大的状态，然后绘制最大可能路径。图7给出了传统水印译码算法的最大可能路径图，图中展示了长度为40500nt的序列的路径偏移状态。从图7中，可以看到由于接收序列发生了大量的删除错误，路径从第一个比特的偏移状态0偏移到了最终比特的状态-90，体现了偏移状态的累加效应，随着接收序列的增长，累加效应越明显，所对应的偏移状态越来越大，这将导致译码所需网格图过大，大大增加了译码复杂度。图7中的细节图更加清晰的展示了第2000个比特到第3000个比特的路径偏移状态，可以看出，在序列较短范围内，状态从-5转移到了-9，偏移状态的累加效应并不明显，起点与终点的偏差值很小，故译码所需网格图相应较小。由此可知，随着片段长度的增加，错误固然增多，则所需的最大偏移量x_max也随着最大可能路径的偏移而增加，提高运算复杂度。本发明提出的分段式策略，将长序列片段化，逐片段处理，减小所需的最大偏移量x_max，减少了网格大小，使偏移状态始终稳定在网格图中间部分，消除了偏移累加效应。

图8给出了6条测序序列在硬判决算法和分段式软判决算法纠正后在第一层和第二层的汉明距离误码率的比较，6条测序数据在第一层和第二层软判决纠错后的错误率都较硬判决纠错后的错误率发生了明显的降低，pic2序列由于初始错误比其他5条序列高，所以纠正后错误率要比其他序列错误率稍高，但分段式软判决算法纠正后错误率明显降低，比硬判决算法表现出了更好的纠错性能。其余5条序列纠正后的汉明距离错误率均低于我们预设的千分之六的初始错误率，充分验证了本发明所提出的基于隐马尔科夫模型的分段式前向-后向软判决纠错算法对于纠正碱基序列中发生的插入删除错误具有的优越的性能。

图9给出了本发明提供的分段式软判决算法与传统软判决算法的复杂度对比图，为验证分段式软判决算法在降低计算复杂度，提高运算速度，减少译码时延的贡献，我们对整体处理和分段式处理的计算量进行了精细化计算。为了便于计算，在固定位置j和偏移a下，将P_a,b·Q_j,a,b(u)作为一个整体定义为分支度量。假设计算一次分支度量所需的乘法次数和加法次数为单位1，设置整体运行状态数2x_max+1为300，分段处理状态数为16。图9描述了两种方案所需的乘法运算次数与加法运算次数的对比，前向/后向度量以及中间度量，分段式处理方法都有明显的运算数目的降低，显著降低了运算量。

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅为了描述，不代表实施例的优劣。

以上所述仅为本发明较佳的实施例，并不用于限制本发明，凡在本发明的精神原则内，所做的任何修改、等同替换以及改进等，均应包含在本发明的保护范围内。

Claims (5)

1.长DNA序列存储的插入删节分段识别方法，其特征在于，所述方法包括以下步骤：

(1)针对由叠加水印的长DNA分子测序得到的长度为

的存在插入、删节与替代错误的读出序列r，基于水印序列和错误传输模型建立隐马尔科夫模型，确定隐藏状态和观测向量，建立状态转移网格图；

(2)改进观测向量的输出概率公式，基于网格图运行前向传递和后向传递，利用硬判决算法估计每个碱基位置的最大可能概率，推断每个碱基距离原始位置的偏移；

(3)对水印序列按照固定长度分段，给每个水印短片段的边界位置建立偏移索引，并基于索引计算与其对应的读出序列边界位置，将读出序列划分为若干长度不同的数据片段；

(4)对分段后的若干测序读出片段，依次递归计算中间度量，结合前向-后向软判决算法对向量d上的符号似然函数l进行估计，根据估计值作硬判决处理，然后与水印序列相异或，得到纠正后的片段，进而拼接为长DNA序列。

2.根据权利要求1所述的长DNA序列存储的插入删节分段识别方法，其特征在于，所述针对由叠加水印的长DNA分子测序得到的长度为

的存在插入、删节与替代错误的读出序列r，基于水印序列和错误传输模型建立隐马尔科夫模型，确定隐藏状态和观测向量，建立状态转移网格图，具体步骤为：

(1.1)根据三代测序的错误特性与前期测序实验的统计数据的积累，估计读出序列的插入错误概率P_i，删除错误概率P_d，替代错误概率P_s以及传输概率P_t＝1-P_i-P_d；

(1.2)定义第i时刻的碱基偏移量x_i为从发送碱基t₀至待发送碱基t_i间存在的插入数目减去删节数目，将碱基偏移量x_i作为隐马尔科夫模型的隐藏状态，将碱基最大偏移量限定为x_max，以降低算法的计算量，x_max定义为

N为设计的长DNA序列的初始长度，x_i的取值范围为X＝{-x_max,...,-1,0,1,...,x_max}；

(1.3)定义测序读出长DNA序列r为r＝r ⁽⁰⁾ r ⁽¹⁾...r ^(N-1)，其子序列为

将r作为观测向量，子序列r ⁽ⁱ⁾的传输过程由x_i-1→x_i的状态转移决定；

(1.4)在状态x_i＝a转变为状态x_i+1＝b时，对应的输出序列被定义为u，根据映射规则{(00)→Α，(01)→T，(10)→G，(11)→C}将读出序列解映射为比特对u₁和u₂，双层水印序列分别为w_1i和w_2i，在水印序列上的噪声包含两部分：信息序列m经编码后，得到符号序列d，经过稀疏器后生成的稀疏码字s及过信道产生的替代错误P_s，故定义有效替代概率P_f为P_f＝f(1-P_s)+(1-f)P_s，其中f为长度为n的稀疏码的平均密度；

(1.5)根据观测向量的传输过程建立隐马尔科夫模型状态转移网格图，网格图中的每个节点对应第i个位置的偏移量x_i，将(x₀,x₁,…,x_i,…,x_N-1)表示为隐马尔科夫模型的状态序列，定义I为最大的连续插入错误数目，每个发送碱基对应一个长度从0到I+1的读出碱基序列，即当x_i＝a时，x_i+1＝b的取值范围是{a-1,…,a+I}。

3.根据权利要求1所述的长DNA序列存储的插入删节分段识别方法，其特征在于，所述改进观测向量的输出概率公式，基于网格图运行前向传递和后向传递，利用硬判决算法估计每个碱基位置的最大可能概率，推断每个碱基距离原始位置的偏移，具体步骤为：

(2.1)由上一状态x_i＝a转移到下一状态x_i+1＝b存在两种情况，第一种，经过信道时，t_i前被插入(b-a+1)个碱基且发送碱基t_i丢失，第二种，经过信道时，t_i前被插入(b-a)个碱基且t_i被传输，基于该传输情形，定义状态转移概率P_a,b＝P(x_i+1＝b|x_i＝a)，计算公式为，

其中，α_I等于1/(1-(P_i)^I)；

(2.2)对双层错误概率改进，发生插入错误时，插入碱基有4种选择(A/T/G/C)，每个碱基的插入概率为1/4，发生替代错误时，映射后的双层比特存在两种情况，(a)上下两层比特均不同，类型包括{A(00)→C(11)，T(01)→G(10)，C(11)→A(00)，G(10)→T(01)}，(b)只有1个比特不同，类型包括{A(00)→T(01)，G(10)→C(11)，T(01)→C(11)，A(00)→G(10)，T(01)→A(00)，C(11)→G(10)，C(11)→T(01)，G(10)→A(00)}，两种情况的错误类型占比1:2，定义有效替代错误概率分别为

和

定义改进的观测向量的输出概率Q_i,a,b(u)＝P(u|x_i＝a,x_i+1＝b,w_1iw_2i)，计算公式为，

其中，该式满足∑ _u Q_i,a,b(u)＝1，

是与水印比特w_1i相关联的读出比特，

是与水印比特w_2i相关联的读出比特；

(2.3)定义前向度量F_i(y)＝P(r₀,…,r_i-1+y,x_i＝y),y∈X，F_i(y)表示x_i＝y且接收到前i-1+y个碱基的概率，定义后向度量B_i(y)＝P(r_i+y,…|x_i＝y),y∈X，B_i(y)表示x_i＝y的条件下输出碱基串(r_i+y,…)的概率，使用分支度量P_a,y·Q_i-1,a,y递归运行前向传递与后向传递，估计每个碱基位置的最大可能概率F_i(y)B_i(y)，推断每个碱基距离原始位置的偏移，其中，第i时刻状态y的前向度量计算公式表示为，

同样的，第i时刻状态y的后向度量计算公式表示为，

4.根据权利要求1所述的长DNA序列存储的插入删节分段识别方法，其特征在于，所述对水印序列按照固定长度分段，给每个水印短片段的边界位置建立偏移索引，并基于索引计算与其对应的读出序列边界位置，将读出序列划分为若干长度不同的数据片段，具体步骤为：

(3.1)将长度为N水印序列分为若干长度为N₁的定长片段，N₁满足可以被稀疏码长度n整除且为固定值，然后将长度为

的读出序列分为长度为N₂的若干片段，N₂长度可变，用t_i表示水印片段的起始节点，则水印片段的边界节点为

假定r_i表示读出序列片段的起始节点，读出序列的边界节点为r_border，其分段边界计算公式表示为，

其中，按照三代测序的特点，N与

不一定等长，则下一个水印分段片段的起始节点为

下一个读出序列分段的起始节点为r_border+1，依次计算，得出N/N₁个片段的起始节点和边界节点，与每个边界位置t_borde的最大概率F_border(y)B_border(y)相对应的最大可能偏移量max_drift_border的计算公式为，

max_drift_border＝argmax(F_border(y)B_border(y))，

(3.2)初始化前向度量与后向度量，假定起始节点固定，初始化

边界位置的后向度量的初始化由水印片段长度N₁和读出序列片段可变长度N₂之间的大小关系判定，

i)当N₂小于N₁时，

ii)当N₂等于N₁时，

iii)当N₂大于N₁时，

5.根据权利要求1所述的长DNA序列存储的插入删节分段识别方法，其特征在于，所述对分段后的若干测序读出片段，依次递归计算中间度量，结合前向-后向软判决算法对向量d上的符号似然函数l进行估计，根据估计值作硬判决处理，然后与水印序列相异或，得到纠正后的片段，进而拼接为长DNA序列，具体步骤为：

(4.1)向量d上的符号d_i与长度为n的稀疏向量(s_i,0,…,s_i,n-1)一一对应，固定d_i＝λ，执行

到

的前向传递，定义中间度量为，

其中符号d_i等概地选取

中的值，r ⁰表示接收比特

为符号d_i对应稀疏串起始位置的同步偏移状态，

为d_i对应稀疏码第k个位置的同步偏移状态；

(4.2)进一步，利用递归公式，第k时刻的中间度量计算公式为，

其中，N_L＝N/n，并初始化M₀(z,d_i)＝1；

(4.3)在硬判决算法的基础上增加中间度量

来获取向量d上的似然函数l，似然函数l输出概率公式可表示为，

和

分别为

时刻的前向度量和

时刻的后向度量，对映射后的上下两层比特单独进行纠正，水印序列和稀疏码表已知，所有的替代错误仅来自于信道，将w_j替换为

P_f替换为P_s，输出概率

可表示为，

(4.4)依次对符号d_i对应的输出概率进行计算，选取概率最大的符号λ所对应的稀疏码(s_i,0,…,s_i,n-1)与本地水印序列w_1i和w_2i异或，得到纠正后的片段，进而组合为长DNA序列。

Images