CN110929542A - 基于分组纠错码的测序条形码构造与软判决识别方法 - Google Patents

Abstract

本发明公开了一种基于分组纠错码的测序条形码构造与软判决识别方法，所述方法包括：将分组纠错码表示为不同的样本，与预定伪随机序列组合生成测序条形码；内译码器接收测序条形码并进行解映射，建立估计碱基序列插入/删节的隐马尔科夫模型，通过对前向‑后向算法进行修正，输出分组码各个比特的软判决信息；内译码器将码字比特概率信息输入外译码器，作为先验信息进行译码，得到信息序列估计值，据此进行样本识别。本发明能够有效纠正DNA合成、建库与测序中的插入、删节和替代错误。

Description

基于分组纠错码的测序条形码构造与软判决识别方法

技术领域

本发明涉及基因组学的高通量测序领域，尤其涉及一种基于分组纠错码的测序条形码构造与软判决识别方法。

背景技术

下一代测序技术正在以其高测序通量和低测序成本改变传统基因组学研究范式。例如，Illumina测序平台能够在2天产生1.8Tb的数据量，华大基因最新推出的超高通量基因测序仪MGISEQ-T7日产出数据高达6Tb。但是，单个测序样本往往并不需要如此高的通量，可以通过多个样本复用的方法同时进行测序，可以充分利用测序仪不断增长的测序通量。使用多路复用技术，需利用测序条形码来区分混合物中的不同样品。测序后，通过识别测序条形码将测序读段进行分离，分配给不同样品。但是，在该过程中，由于引物合成、连接、样品扩增以及高通量测序中的错误，测序条形码上可能发生各类错误。更为严重的是，对于某些典型的三代测序技术，例如PacBio单分子实时(Single Molecule Real-Time,SMRT)测序技术和牛津纳米孔技术(Oxford Nanopore Technology，ONT)，读取的DNA测序读段(read)很长，并且插入、删节和替代错误率很高。例如，SMRT测序中的测序错误率甚至会达到11％，并且插入/删节错误占主要部分，只有1％左右的替代错误。ONT测序中有15％左右的测序错误率，其中插入/删节为主要的错误类型(10％)，其次是替代错误(约占5％)。这些错误在解复用过程中导致样本错误分配，严重影响测序效率。因此，需要通过设计纠错DNA测序用条形码来抵抗DNA合成和测序中引入的插入/删节和替代错误，从而增加正确识别的样本量，提高整体测序效率。

目前，已经提出了许多基于纠错码的测序条形码构造与识别方法。针对可能存在的插入、删节错误，研究者提出了基于Levenshtein距离的长度为6bp，最小距离为3的测序条形码集，因此能够纠正单个插入/删节/替代错误。又有研究者在此基础上，将传统的Levenshtein距离度量进行改进，提出了“Sequence-Levenshtein”码，将其作为DNA测序条形码，这种新的距离度量考虑了DNA文本的干扰以及由此产生的测序条形码之间的距离。该方法能够检测DNA文本中测序条形码的边界，并且纠正比传统的Levenshtein更多的随机突变；但是，这种方法对于选择符合纠错条件的测序条形码集具有更高的难度和复杂度。

研究者也提出利用Davey和Mackay提出的水印码作为DNA测序条形码，并对用于测序的信道模型和译码器进行适当的修改，使其适用于DNA文本。这种识别方法能够检测并恢复测序条形码在DNA测序的长读段中的位置，而且识别错误率较低。但是这种方法的译码复杂度较高，短码的性能并不突出。

针对上述问题，发明人发明了一种基于一般的分组纠错码结合已知的伪随机序列的高通量复用测序用条形码，一方面可生成不同长度、不同容错能力的条形码；另一方面，设计条形码具有很强的纠错能力，尤其是对抗三代测序中的插入与删节错误的能力，非常适合三代测序应用。进一步，本发明提出了针对上述条形码的检测方法，通过建立序列的隐马尔科夫模型结合前后向算法以及纠错纠删译码，能高效识别被错误破坏的读段恢复原始测序条形码，具有很高的鲁棒性。总之，本发明提供一种针对未来三代高通量、高错误率、长读段条件下复用测序的条形码的灵活设计方法与识别方法，具有很好的应用前景。

发明内容

本发明提供了一种基于分组纠错码的测序条形码构造与软判决识别方法，本发明能够有效纠正DNA合成、建库与测序中的插入、删节和替代错误，详见下文描述：

一种基于分组纠错码的测序条形码构造与软判决识别方法，所述方法包括：

将分组纠错码表示为不同的样本，与预定伪随机序列组合生成测序条形码；

内译码器接收测序条形码并进行解映射，建立估计碱基序列插入/删节的隐马尔科夫模型，通过对前向-后向算法进行修正，输出分组码各个比特的软判决信息；

内译码器将码字比特概率信息输入外译码器，作为先验信息进行译码，得到信息序列估计值，据此进行样本识别。

其中，所述将分组纠错码表示为不同的样本，与预定伪随机序列组合生成测序条形码具体为：

将表示某一特定样本序号的k比特的信息向量，经分组码编码器产生长度为n的分组码码字；将预定伪随机序列与分组码码字对应位置的比特组合成比特对，按照映射规则将比特对转换成碱基，生成长度为n的测序条形码。

进一步地，所述建立估计碱基序列插入/删节的隐马尔科夫模型具体为：

将碱基漂移量x_i作为HMM的隐状态，接收测序条形码估计值作为HMM的观测向量；当发送第i个碱基t_i时，发生从x_i-1到x_i的状态转移，接收测序条形码子序列r⁽ⁱ⁾输出。

其中，所述通过对前向-后向算法进行修正，输出分组码各个比特的软判决信息具体为：

计算由当前时刻的漂移状态a转移到下一时刻的漂移状态b的转移概率P_a,b；

对由x_i＝a转移到x_i+1＝b时，计算输出碱基序列的条件概率；

对0≤i＜n，-x_max≤a≤x_max，a-1≤b≤a+I且b∈X，

计算发送碱基为t_i＝λ时，从状态x_i＝a转移到x_i+1＝b得到接收碱基序列的中间度量值；

对0≤i＜n，-x_max≤a≤x_max，初始化第0时刻的前向度量值，利用转移概率P_c,a和输出概率，计算第i时刻状态x_i＝a的前向度量值；

对0＜i≤n，-x_max≤b≤x_max，初始化第n时刻的后向度量值，利用转移概率P_b,c和输出概率，递归计算第i时刻状态x_i＝b的后向度量值；

对0≤i＜n，

计算发送碱基t_i＝λ时的核苷酸级软信息；

对0≤i＜n，基于核苷酸级软信息，根据发送碱基与码字比特的映射关系，计算码字比特d_i取0和1的概率值，组成比特级软信息。

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

1、本发明将分组纠错码表示不同的样本，与已知伪随机序列组合构造高可靠的测序用条形码；

2、本发明建立估计碱基序列的隐马尔科夫模型(HMM)，采用前向-后向算法，依据已知的伪随机序列估计插入/删节，生成分组码各个比特的软判决信息；

3、本发明执行软判决译码，纠正测序条形码上的各类错误，这种测序条形码识别方案对于测序中的插入、删节与替代错误，具有高鲁棒性。

附图说明

图1为本发明提出的测序用条形码软判决识别方法框图；

图2为本发明提出的测序用条形码软判决识别方法中译码流程图；

图3为本发明提出的核苷酸级HMM模型以及HMM中输出概率Q_i,a,b的计算示意图；

图4为DNA测序用条形码在多路测序中应用的系统框图以及两种不同的测序错误场景图；

图5为本发明用于测序用条形码(BCH(24,9,3))时，在两种错误场景下，其识别错误率比较示意图；

图6为本发明用于具有不同k值的长度为24碱基的测序用条形码时，其识别错误率比较示意图。

具体实施方式

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

本发明提出一种高鲁棒性的基于分组纠错码的测序条形码构造与软判决识别方法。具体而言，将一般的分组纠错码表示不同的样本，与预定伪随机序列组合构造测序用条形码。进一步，建立估计碱基序列的隐马尔科夫模型(HMM)，采用前向-后向算法，依据已知的伪随机序列估计插入/删节，生成分组码各个比特的软判决信息。进一步执行软判决译码，有效地纠正测序条形码上的各类错误。这种识别方法能够抵抗DNA合成、建库与测序中的插入、删节和替代错误，具有很高鲁棒性。

参见图1和图2，一种基于分组纠错码的测序条形码构造与软判决识别方法的具体步骤为：

(1)用长度为k比特的不同信息序列分别表示样本序号{0,1,...,2^k-1}，则DNA测序用条形码最多能支持2^k个不同样本并行测序；将表示某一特定样本序号i(i∈{0,1,...,2^k-1})的k比特的信息向量m_i，经过分组码编码器(n,k,t)，产生长度为n的分组码码字d(参数t表示分组码码字的纠错能力，可根据具体测序质量进行选择)；将预定伪随机序列w＝(w₀,...,w_n-1)与分组码码字d＝(d₀,...,d_n-1)对应位置的比特组合成比特对(w₀d₀,w₁d₁,...,w_n-1d_n-1)，按照映射规则{00→A，01→T，10→G，11→C}，将比特对转换成碱基，生成长度为n的测序条形码t；经过样本建库、扩增与测序，得到包含各类错误的接收测序条形码估计值

其中，上述对样本进行建库、扩增与测序的操作为本领域技术人员所公知，本发明实施例对此不做赘述。

(2)对接收测序条形码估计值

进行解映射，按照解映射规则{A→00，T→01，G→10，C→11}，将碱基序列(r₀,...,r_n-1)转换为比特对(u₀c₀,u₁c₁,...,u_n-1c_n-1)，得到相应的接收序列估计值u＝(u₀,...,u_n-1)和接收码字估计值c＝(c₀,...,c_n-1)；

(3)参见图3(a)，建立估计碱基序列插入/删节的隐马尔科夫模型(HMM)，将碱基漂移量x_i作为HMM的隐状态，接收测序条形码估计值

作为HMM的观测向量，其中，

为接收测序条形码

的子序列；当发送第i个碱基t_i时，发生从x_i-1到x_i的状态转移，从而导致接收测序条形码子序列r⁽ⁱ⁾的输出；这里，第i个位置碱基的漂移量x_i(i＝0,1,...,n)表示从发送第一个碱基t₀到待发送碱基t_i之间存在的插入数目减去删节数目，漂移量x_i的状态集为X＝{-x_max,...,-1,0,1,...,x_max}，共有X＝2x_max+1个取值，x_max表示最大漂移量，设置为

(P_d为碱基删节概率)；

(4)计算由当前时刻的漂移状态a(a∈X)转移到下一时刻的漂移状态b(b∈{a-1,...,a+I}且b∈X)的转移概率P_a,b，具体公式为，

其中，P_i、P_d分别表示碱基插入、删节概率，根据测序后的碱基插入、删节数目估计得到；I表示最大连续插入碱基数，设置为大于1的常数；连续插入碱基序列的长度呈几何分布，当限定连续插入长度不超过I的条件下，归一化常数

(5)参见图3(b)，对0≤i＜n，-x_max≤a≤x_max，a-1≤b≤a+I且b∈X，计算由x_i＝a转移到x_i+1＝b时，输出碱基序列(r_i+a,...,r_i+b)的条件概率，即输出概率Q_i,a,b(r_i+a,...,r_i+b)；以w_i＝0，t_i＝A为例，当u^*＝w_i(u^*＝u_i+b是接收序列u中与伪随机比特w_i相关联的接收比特)，除发生插入/删节外，碱基可能有两种传输情形：(i)正确传输A(00)→A(00)，正确传输概率为(1-P_s)(P_s表示碱基替代概率)；(ii)替代A(00)→T(01)，该替代情形占w_i＝0对应的所有可能的碱基替代情形(A(00)→T(01)，A(00)→G(10)，A(00)→C(11))的1/3，因此，概率为1/3P_s；当u^*≠w_i时，除发生插入/删节外，碱基只有一种传输情形，即替代情形A(00)→G(10)，A(00)→C(11)，占所有可能符号替代情形的2/3，因此，概率为2/3P_s；w_i＝1的情况类似；因此，输出概率Q_i,a,b(r_i+a,...,r_i+b)修正后的计算公式为：

其中，u^*＝u_i+b是接收序列u中与伪随机比特w_i相关联的接收比特。

(6)对0≤i＜n，-x_max≤a≤x_max，a-1≤b≤a+I且b∈X，

计算发送碱基为t_i＝λ时，从状态x_i＝a转移到x_i+1＝b的中间度量值

中间度量包含码字信息(d_i＝0，1)，则发送符号t_i是已知的。

因此，直接通过比较发送符号t_i和对应的接收符号r^*，计算中间度量值；即若t_i＝r^*，除发生插入/删节外，碱基正确传输，概率为1-P_s；若t_i≠r^*，则碱基发生替代，概率为P_s；中间度量值

修正后的计算公式为：

其中，

r^*＝r_i+a是接收碱基序列

中与发送碱基t_i相关联的接收碱基；

表示伪随机比特w_i对应的碱基符号子集，由比特对到碱基的映射规则可知，当w_i＝0，

当w_i＝1，

(7)对0≤i＜n，-x_max≤a≤x_max，初始化第0时刻的前向度量值

利用转移概率P_c,a和输出概率Q_i,c,a(r_i+c,...,r_i+a)，递归计算第i时刻状态x_i＝a的前向度量值F_i(a)＝P(r₀,...,r_i-1+a,x_i＝a|w)，

其中，a,c表示同步漂移状态，c∈{a-I,...,a+1}且c∈X；F_i(a)表示x_i＝a且接收到前i+a个碱基的概率；

(8)对0＜i≤n，-x_max≤b≤x_max，初始化第n时刻的后向度量值

利用转移概率P_b,c和输出概率Q_i,b,c(r_i+b,...,r_i+c)，递归计算第i时刻状态x_i＝b的后向度量值B_i(b)＝P(r_i+b,...,r_n|x_i＝b,w)，

其中，b,c表示同步漂移状态，c∈{b-1,...,b+I}且c∈X；B_i(b)表示x_i＝b的条件下输出碱基序列(r_i+b,...)的概率；

(9)对0≤i＜n，

计算发送碱基t_i＝λ时的核苷酸级软信息P(r|t_i＝λ,w)：

(10)对于0≤i＜n，基于核苷酸级软信息

同时根据发送碱基

与码字比特d_i∈{0,1}的映射关系，计算码字比特d_i取0和1的概率值(l_i0,l_i1)，组成比特级软信息l＝((l₀₀,l₀₁),(l₁₀,l₁₁),...,(l_n-1,0,l_n-1,1))；

具体地，对于w_i＝0,χ⁰＝{A,T}，

对于w_i＝1,χ¹＝{G,C}，

(11)根据分组码各个比特对应的软信息l，进行软判决译码，输出译码码字

得到信息向量的估计值

从而识别不同测序样本。

下面给出一个具体的实施例，说明本发明给出的基于分组纠错码的测序条形码构造与软判决识别方法的具体执行步骤以及对各种非理想因素的高鲁棒性。

如图4所示，在本发明的实施例中，分组纠错码具体采用Bose-ChaudhuriHocquenghem(BCH)码和定义在高阶伽罗华域(Galois Field，GF)上的非二进制单奇偶校验(Non-Binary Single Parity Check，NB-SPC)码构造DNA测序用条形码。然后，经过DNA样本建库、扩增和测序，基于测序读上的引物位置，得到包含各类错误的接收测序条形码。最后，对接收测序条形码进行基于软信息的识别。

在步骤(1)中，分别用BCH(24,14,2)，BCH(24,9,3)和定义在GF(2⁶)的NB-SPC(4,3)生成编码码字d。这里，定义在GF(2⁶)的NB-SPC(4,3)可以视为n＝24,k＝18的二进制单奇偶校验码。然后将设计的长度为24的伪随机序列w与码字d对应位置的比特组合成比特对，通过映射{00→A，01→T，10→G，11→C}，生成三种不同的长度为24碱基的测序条形码t。

考虑到实际测序场景中，barcode需要嵌入到DNA文本中，由于没有“分隔符”将barcode与DNA文本分开，因而检测不到由于插入或删节而导致的barcode长度变化。为了模拟实际测序错误情形，按照以下两种不同的错误场景对barcode加入一定量的插入、删节和替代错误，如图4所示。(i)Barcode上随机发生相同数目的插入和删节，以及一定数目的替代，barcode长度不变。(ii)Barcode上发生插入或者删节，导致barcode长度发生变化。若barcode长度由于删节而缩短，则在barcode尾部补入来自后续样本序列的碱基。若barcode长度由于插入而延长，则删除barcode尾部多余碱基。然后，得到包含各类错误的长度固定的接收测序条形码估计值

其中，N_i＝N_d，则t与

的长度相等。

在步骤(2)中，首先，按照规则{A→00，T→01，G→10，C→11}，对

进行解映射，得到接收序列估计值u。然后，将

w和u输入内译码器，采用修正的基于HMM的前向-后向算法，输出码字比特软信息。各个参数设置如下：根据碱基插入、删节、替代数目N_i、N_d、N_s确定插入、删节、替代概率分别为：P_i＝N_i/n,P_d＝N_d/n,P_s＝N_s/n(n＝24)，则可得测序条形码上每碱基突变概率P_mut＝P_i+P_d+P_s；信道中最大连续插入数目设为I＝2，最大漂移量设为

每个时刻的状态数X＝2x_max+1。最后输出比特概率信息l＝((l₀₀,l₀₁),(l₁₀,l₁₁),...,(l_23,0,l_23,1))。最后，将比特概率信息l输入分组码译码器得到译码码字

其中，BCH(n,k,t)译码器采用Chase软判决译码算法进行译码，最小不可靠位数设为p＝8。NB-SPC(n,k)译码器采用基于快速傅里叶变换的置信传播(belief propagationdecoding algorithm based on fast Fourier transform，FFT-BP)译码算法。

图5给出了在两种错误场景下，由BCH(24,9,3)构造的24-nt的测序条形码的识别错误率与条形码上碱基突变概率P_mut∈[0.10,0.25]的关系曲线。本次仿真考虑了两种插入/删节概率情形P_i＝P_d＝0.033和P_i＝P_d＝0.066。仿真结果表明对于固定的碱基突变概率P_mut，测序条形码在插入/删节概率为P_i＝P_d＝0.066情形下的识别错误率低于在P_i＝P_d＝0.033情形下，表明测序条形码上不同错误类型的比例会对测序条形码的识别错误率产生影响。此外，对于P_mut∈[0.10,0.25]，错误场景(ii)下的识别错误率与错误场景(i)下的识别错误率差异很小，当P_mut≤0.10，识别错误率均低于10^-5，具有高鲁棒性。说明测序条形码在两种测序错误场景中均具有良好的纠错性能，这为测序条形码在DNA实际测序中的应用提供了可靠保证。

图6给出了在错误场景(ii)中，构造的三种具有不同k值的长度为24碱基的测序条形码的识别错误率与碱基突变概率P_mut的关系图，其中P_mut∈[0.10,0.22](P_i＝P_d＝0.033和P_i＝P_d＝0.066)。仿真结果表明，具有不同k值的测序条形码具有不同的纠错性能，k值越小，其纠错能力越强，但其能支持的并行样本数越少。因此，可以根据不同的测序质量和试验要求构造具有不同长度、不同纠错能力的测序条形码。

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

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

Claims (4)

1.一种基于分组纠错码的测序条形码构造与软判决识别方法，其特征在于，所述方法包括：

将分组纠错码表示为不同的样本，与预定伪随机序列组合生成测序条形码；

内译码器接收测序条形码并进行解映射，建立估计碱基序列插入/删节的隐马尔科夫模型，通过对前向-后向算法进行修正，输出分组码各个比特的软判决信息；

内译码器将码字比特概率信息输入外译码器，作为先验信息进行译码，得到信息序列估计值，据此进行样本识别。

2.根据权利要求1所述的一种基于分组纠错码的测序条形码构造与软判决识别方法，其特征在于，所述将分组纠错码表示为不同的样本，与预定伪随机序列组合生成测序条形码具体为：

将表示某一特定样本序号的k比特的信息向量，经分组码编码器产生长度为n的分组码码字；将预定伪随机序列与分组码码字对应位置的比特组合成比特对，按照映射规则将比特对转换成碱基，生成长度为n的测序条形码，所述测序条形码用于样本建库、扩增与测序。

3.根据权利要求1所述的一种基于分组纠错码的测序条形码构造与软判决识别方法，其特征在于，所述建立估计碱基序列插入/删节的隐马尔科夫模型具体为：

将碱基漂移量x_i作为HMM的隐状态，接收测序条形码估计值作为HMM的观测向量；当发送第i个碱基t_i时，发生从x_i-1到x_i的状态转移，接收测序条形码子序列r⁽ⁱ⁾输出。

4.根据权利要求1所述的一种基于分组纠错码的测序条形码构造与软判决识别方法，其特征在于，所述通过对前向-后向算法进行修正，输出分组码各个比特的软判决信息具体为：

计算由当前时刻的漂移状态a转移到下一时刻的漂移状态b的转移概率P_a,b；

对由x_i＝a转移到x_i+1＝b时，计算输出碱基序列的条件概率；

对0≤i＜n，-x_max≤a≤x_max，a-1≤b≤a+I且b∈X，

计算发送碱基为t_i＝λ时，从状态x_i＝a转移到x_i+1＝b得到接收碱基序列的中间度量值；

对0≤i＜n，-x_max≤a≤x_max，初始化第0时刻的前向度量值，利用转移概率P_c,a和输出概率，计算第i时刻状态x_i＝a的前向度量值；

对0＜i≤n，-x_max≤b≤x_max，初始化第n时刻的后向度量值，利用转移概率P_b,c和输出概率，递归计算第i时刻状态x_i＝b的后向度量值；

对0≤i＜n，

计算发送碱基t_i＝λ时的核苷酸级软信息；

对0≤i＜n，基于核苷酸级软信息，根据发送碱基与码字比特的映射关系，计算码字比特d_i取0和1的概率值，组成比特级软信息。

Images