一种基于神经网络的真实纳米孔测序信号滤波方法及装置
阅读说明：本技术 一种基于神经网络的真实纳米孔测序信号滤波方法及装置 (Real nanopore sequencing signal filtering method and device based on neural network ) 是由 陈为刚 张鹏 韩昌彩 赵毅强 于 2020-12-28 设计创作，主要内容包括：本发明公开了一种基于神经网络的真实纳米孔测序信号滤波方法及装置,所述方法包括以下步骤：在纳米孔测序K-mer孔模型中输入核苷酸序列将其转换为与其对应的预期测序信号序列,将预期测序信号序列中的每个信号值重复多次生成待测真实测序信号序列；构建基于双向门控循环单元神经网络的真实测序信号处理模型；构建信号处理模型的损失函数,进行信号处理模型参数的初始化,通过自适应优化器最小化损失函数实现模型参数的训练；将待测真实测序信号序列输入到完成参数训练的信号处理模型中实现对待测真实测序信号序列的滤波处理；本发明的优点在于：能够准确地滤掉待测真实测序信号序列中与实际测序信号序列无关的高频分量并保留有用的高频分量。(The invention discloses a real nanopore sequencing signal filtering method and a device based on a neural network, wherein the method comprises the following steps: inputting a nucleotide sequence into a nanopore sequencing K-mer pore model, converting the nucleotide sequence into an expected sequencing signal sequence corresponding to the nucleotide sequence, and repeating each signal value in the expected sequencing signal sequence for multiple times to generate a real sequencing signal sequence to be detected; constructing a real sequencing signal processing model based on a bidirectional gating cycle unit neural network; constructing a loss function of the signal processing model, initializing parameters of the signal processing model, and realizing the training of the model parameters by minimizing the loss function through a self-adaptive optimizer; inputting the real sequencing signal sequence to be tested into a signal processing model for completing parameter training to realize filtering processing on the real sequencing signal sequence to be tested; the invention has the advantages that: the method can accurately filter out high-frequency components irrelevant to the actual sequencing signal sequence in the real sequencing signal sequence to be detected and reserve useful high-frequency components.)
技术领域
本发明涉及纳米孔测序信号处理领域，更具体涉及一种基于神经网络的真实纳米孔测序信号滤波方法及装置。
背景技术
新一代测序技术的提出使研究人员能够以高通量的方式对DNA和RNA进行测序，这促进了在基因组学、转录组学和表观基因组学等领域的众多突破。目前市场上最受欢迎的新一代测序技术主要包括Illumina、PacBio和Nanopore等测序平台。与其他测序技术不同的是，纳米孔测序技术主要是通过嵌入在分隔两个电解液室的薄膜上的纳米孔来直接检测碱基序列。通过在薄膜上施加电压差，将单链核苷酸序列通过纳米孔穿过薄膜，孔中的碱基会影响孔的电阻，因此可以由随时间变化的电流信号检测到通过孔的碱基序列。与目前常用的短读段测序技术(例如Illumina MiSeq测序平台)相比，纳米孔测序技术具有多个优势。首先，纳米孔测序技术可以从单个核苷酸分子中实时产生测序读段，结合快速的文库制备工作，有效地缩短了从测序生物样本预处理到测序数据分析之间的时间。此外，纳米孔测序技术还可以直接用于RNA测序，而无需提前进行逆转录或扩增工作。任意长度的DNA或RNA分子都可以使用纳米孔测序技术来进行测序，长读段是非常有价值的，它们简化了基因组组装和结构变异检测等研究工作。与此同时，纳米孔测序平台(如MinION测序仪)要比当前的短读段测序平台更加轻便，这使得它能够在传统实验室环境之外进行测序工作。
尽管纳米孔测序技术的优势众多，但复杂的测序环境会导致低信噪比的测序信号序列，这对于进一步的测序数据分析提出了挑战。在纳米孔测序过程中，待测序碱基序列移动通过纳米孔时输出的电流信号测量值存储为16位整数值。对于当前的MinION孔化学，单个DNA链以平均450bp/s的速度穿过纳米孔，而测序电流信号的采样频率为4kHz，这意味着平均每个k-mer有9个离散的测量信号值。尽管由于运动蛋白易位速度的波动，这些值可能是不同的。为了将输出的原始测序电流信号序列转换成碱基序列，需要复杂的碱基识别软件来实现。在纳米孔测序过程中，有以下几个因素会导致低信噪比的原始电流信号序列：通过纳米孔的DNA或RNA的四种碱基结构较为相似，导致不同碱基产生的原始电信号差异较小；原始电流信号主要受到同时占据纳米孔的5或6个碱基的影响，所以一个电信号测量值对应45或46个可能的k-mers；牵引DNA或RNA序列移动的运动蛋白的随机性导致碱基序列的移速不均匀；均聚物通过纳米孔时的电信号没有发生变化，从而导致多个相同碱基的长链被错误解释。上述几个因素使得纳米孔测序技术产生的测序读段精度较低，从而限制了纳米孔测序技术在实际中的应用。
随着纳米孔测序技术的快速发展，针对纳米孔测序数据的分析方法和工具也迅速涌现。例如，用于纳米孔测序技术的基因序列比对工具Graphmap，Minimap2和MashMap2，用于组装纳米孔测序长读段的基因组组装工具Canu和Racon，用于纳米孔测序信号的可视化工具BulkVis和SquiggleKit。可以预见，在不久的将来研究人员将开发出更多的用于纳米孔测序的数据分析方法和工具。研究人员通常使用实验测序数据或仿真测序数据对这些新算法和新工具进行性能测试，与实验测序数据相比，仿真测序数据可以大大节约研发成本，降低数据分析难度，提高研发效率，因此对纳米孔测序数据进行更准确的仿真有助于为纳米孔测序技术开发出更多性能更优的方法和工具。
纳米孔测序数据的仿真算法可以分为纳米孔测序读段仿真和纳米孔测序信号仿真两种。ReadSim、SiLiCO和NanoSim三个纳米孔测序读段仿真软件利用输入的核苷酸序列和配置文件生成仿真测序读段，其中配置文件包含一组预先设置的参数，例如插入率、删节率、替代率、读段长度和质量分数等参数。这三个纳米孔测序读段仿真软件的不同之处在于ReadSim使用固定的配置文件，SiLiCO使用用户提供的配置文件，NanoSim使用用户提供的实际测序数据来学习将在仿真阶段使用的配置文件。尽管纳米孔测序读段仿真软件可以生成高质量的仿真读段，但是纳米孔测序电信号才是纳米孔测序技术的本质，这些纳米孔测序读段仿真软件无法生成纳米孔测序仿真信号，从而限制了这些测序读段仿真软件的应用。
DeepSimulator是第一个完整地模拟纳米孔测序技术整个流程的仿真软件，它可以同时生成纳米孔测序仿真信号和仿真读段。纳米孔测序技术整个流程主要包括三个阶段：首先通过样本制备过程产生用于测序实验的核苷酸样品；然后使用纳米孔测序设备(如MinION)测量核苷酸序列通过纳米孔时输出的电流信号序列，采集到的测序电流信号序列通常存储在fast5文件中；最后使用碱基识别软件将测序电流信号序列转换为碱基序列。相应地，DeepSimulator的主要工作框架包括序列生成模块，信号生成模块和碱基识别模块。首先，序列生成模块随机选择输入参考基因组序列上的起始位置，生成满足实际测序读段长度分布的较短核苷酸序列。然后，信号生成模块根据已知的纳米孔测序6-mer孔模型生成序列生成模块输出的核苷酸序列所对应的纳米孔测序仿真信号序列。最后，碱基识别模块使用碱基识别软件(Albacore、Guppy)将仿真测序信号序列转换为仿真测序读段。在DeepSimulator主要工作框架中的三个模块中，信号生成模块是DeepSimulator软件的核心模块，其首先根据已知的纳米孔测序6-mer孔模型和实际测序信号序列的重复次数分布生成输入核苷酸序列所对应的真实测序信号序列，然后使用低通滤波器滤掉嵌入在真实测序信号序列中与真实测序信号序列无关的高频分量，最后在滤波信号序列上添加高斯噪声输出最终的仿真测序信号序列。
为了生成高质量的纳米孔仿真测序信号序列，DeepSimulator信号生成模块中的一个关键步骤是使用低通滤波器对真实测序信号序列进行滤波处理。真实测序信号序列是由一系列方波组成，其频谱是无限正弦波的组合。为了更加准确地对纳米孔测序信号序列进行仿真，必须对嵌入在真实测序信号序列方波中的高频分量进行滤波处理。DeepSimulator通过将windowed-sinc函数和真实测序信号序列进行卷积来实现对真实测序信号序列的低通滤波处理。考虑到单链核苷酸序列通过纳米孔的速度约为450bp/s，因此低通滤波器的截止频率应大于450Hz。当低通滤波器的截止频率设为950Hz时，可以生成与实际测序信号序列最相似的仿真测序信号序列。低通滤波器是一种可以通过频率低于所选截止频率的信号并衰减频率高于截止频率的信号的滤波器。由于低通滤波器会衰减真实测序信号序列中所有高于截止频率的高频分量，无法保留其中与实际测序信号序列相关的高频分量，这对于某些输入核苷酸序列来说可能会造成仿真测序信号序列与真实测序信号序列之间存在着较大的差异，这给关心仿真测序信号序列输出的用户带来不便。
发明内容
本发明所要解决的技术问题在于传统低通滤波器无法准确地滤掉真实测序信号序列中与实际测序信号序列无关的高频分量的问题。
本发明通过以下技术手段实现解决上述技术问题的：一种基于神经网络的真实纳米孔测序信号滤波方法，所述方法包括以下步骤：
步骤一：在纳米孔测序K-mer孔模型中输入核苷酸序列将其转换为与其对应的预期测序信号序列，将预期测序信号序列中的每个信号值重复多次生成待测真实测序信号序列；
步骤二：构建基于双向门控循环单元神经网络的真实测序信号处理模型，该信号处理模型包括三层双向门控循环单元神经网络和一层全连接层，其中信号处理模型的输入为输入到第一层双向门控循环单元神经网络的经中值归一化的真实测序信号序列，信号处理模型的输出为全连接层输出的与输入的待测真实测序信号序列等长的滤波信号序列；
步骤三：获取纳米孔测序平台输出的待训练实际测序信号序列并根据待训练实际测序信号序列计算每个待训练实际测序信号序列所对应的待训练真实测序信号序列，将待训练实际测序信号序列和与其对应的待训练真实测序信号序列作为信号处理模型参数训练所需的监督训练数据，构建信号处理模型的损失函数，进行信号处理模型参数的初始化，通过自适应优化器最小化损失函数实现模型参数的训练；
步骤四：将待测真实测序信号序列输入到完成参数训练的信号处理模型中实现对待测真实测序信号序列的滤波处理。
本发明基于已知的纳米孔测序K-mer孔模型将输入核苷酸序列转换为预期测序信号序列，将预期测序信号序列中的每个信号值重复多次生成待测真实测序信号序列；然后建立真实测序信号序列处理模型，使用监督训练数据完成信号处理模型参数的训练；最后使用完成训练的信号处理模型对待测真实测序信号序列进行滤波处理，利用神经网络学习到实际测序信号序列的时频特性，进而能够准确地滤掉待测真实测序信号序列中与实际测序信号序列无关的高频分量并保留有用的高频分量。
进一步地，所述步骤一包括：
步骤101：在纳米孔测序K-mer孔模型中输入一条具有T个碱基的核苷酸序列X＝x1,x2,...,xT，xT表示第T个碱基，从输入核苷酸序列的第一个碱基开始，每次移动一个碱基取出一个K-mer，直到移动到倒数第K个碱基位置结束，T个碱基共得到T-K+1个K-mer；
步骤102：纳米孔测序K-mer孔模型包含每个K-mer所对应的预期电流信号值，对照纳米孔测序K-mer孔模型依次查找T-K+1个K-mer所对应的预期电流信号值，生成输入核苷酸序列X对应的预期测序信号序列Y＝y1,y2,...,yT-5，其中yi表示从X中的位置i开始的K-mer所对应的预期电流信号值；
步骤103：根据实际测序信号序列中的信号值重复次数分布，将预期测序信号序列中的每个信号值重复多次，生成与实际测序信号序列之间具有相似长度分布的待测真实测序信号序列。
进一步地，所述步骤二包括：
步骤201：将中值归一化后的待测真实测序信号序列I＝i1,i2,...,iT作为信号处理模型中第一层双向门控循环单元神经网络的输入向量，其中双向门控循环单元神经网络的基本单元由一个前向传播的门控循环单元和一个后向传播的门控循环单元组成，首先计算t时刻的前向输出向量计算公式为：
其中，it为t时刻的输入信号值，为t-1时刻的前向输出向量，σ为sigmoid函数且σ(z)＝1/(1+e-z)，表示第一中间变量，表示第二中间变量，表示第三中间变量，表示两个等维向量的对应元素相乘，为第一权重矩阵，为第二权重矩阵，Wf为第三权重矩阵，为第一偏置向量，为第二偏置向量，然后计算t时刻的后向输出向量计算公式为：
其中，为t+1时刻的后向输出向量，表示第四中间变量，表示第五中间变量，表示第六中间变量，为第四权重矩阵，为第五权重矩阵，Wr为第六权重矩阵，为第三偏置向量，为第四偏置向量，最后t时刻的输出向量ht是由t时刻的前向输出向量和t时刻的后向输出向量连接得到，即其中||表示向量的连接符号；
步骤202：将第一层双向门控循环单元神经网络的输出向量作为第二层双向门控循环单元神经网络的输入向量，第二层双向门控循环单元神经网络的输出向量作为第三层双向门控循环单元神经网络的输入向量，其中第二层双向门控循环单元神经网络、第三层双向门控循环单元神经网络的计算过程与第一层双向门控循环单元神经网络的计算过程相同，不同层之间的权重矩阵与偏置向量参数不同；
步骤203：将三层双向门控循环单元神经网络中最后一层的输出向量作为全连接层的输入向量进行处理，全连接层的输出向量为与输入的待测真实测序信号序列等长的滤波信号序列O＝o1,o2,...,oT。
进一步地，所述步骤三包括：
步骤301：信号处理模型的监督训练数据包含待训练实际测序信号序列和与其对应的待训练真实测序信号序列，首先读取纳米孔测序平台输出的fast5文件得到待训练实际测序信号序列，然后使用连续小波动态时间规整算法计算每个待训练实际测序信号序列所对应的待训练真实测序信号序列；
步骤302：构建信号处理模型的损失函数为其中，cosh表示双曲函数，其中，oi表示信号处理模型的第i个输出信号，ri表示实际测序信号序列中第i个信号；
步骤303：将每层双向门控循环单元神经网络中权重矩阵以及偏置向量初始化为区间内的均匀分布，其中ni表示每层双向门控循环单元神经网络的输入向量维度，no表示每层双向门控循环单元神经网络的输出向量维度；
步骤304：使用Adam自适应优化器迭代计算得到最小化损失函数的模型参数，完成信号处理模型的参数训练过程。
更进一步地，所述步骤301包括：
步骤3011：读取纳米孔测序平台输出的fast5测序文件得到待训练实际测序信号序列，其中fast5文件为HDF5文件格式；
步骤3012：对待训练实际测序信号序列进行碱基识别得到测序读段，采用基因序列比对算法将测序读段比对到参考基因组序列，根据比对结果得到与待训练实际测序信号序列对应的参考基因组序列片段，然后使用参考基因组序列片段和纳米孔测序K-mer孔模计算出待训练预期测序信号序列；
步骤3013：使用连续小波动态时间规整算法完成待训练实际测序信号序列与待训练预期测序信号序列之间的点到点映射，根据映射结果，求得待训练实际测序信号序列中每个信号值所对应的预期信号值，进而计算出待训练实际测序信号序列所对应的待训练真实测序信号序列。
本发明还提供一种基于神经网络的真实纳米孔测序信号滤波装置，所述装置包括以下步骤：
待测真实测序信号序列生成模块，用于在纳米孔测序K-mer孔模型中输入核苷酸序列将其转换为与其对应的预期测序信号序列，将预期测序信号序列中的每个信号值重复多次生成待测真实测序信号序列；
信号处理模型构建模块，用于构建基于双向门控循环单元神经网络的真实测序信号处理模型，该信号处理模型包括三层双向门控循环单元神经网络和一层全连接层，其中信号处理模型的输入为输入到第一层双向门控循环单元神经网络的经中值归一化的真实测序信号序列，信号处理模型的输出为全连接层输出的与输入的待测真实测序信号序列等长的滤波信号序列；
信号处理模型训练模块，用于获取纳米孔测序平台输出的待训练实际测序信号序列并根据待训练实际测序信号序列计算每个待训练实际测序信号序列所对应的待训练真实测序信号序列，将待训练实际测序信号序列和与其对应的待训练真实测序信号序列作为信号处理模型参数训练所需的监督训练数据，构建信号处理模型的损失函数，进行信号处理模型参数的初始化，通过自适应优化器最小化损失函数实现模型参数的训练；
滤波处理模块，用于将待测真实测序信号序列输入到完成参数训练的信号处理模型中实现对待测真实测序信号序列的滤波处理。
进一步地，所述待测真实测序信号序列生成模块还用于：
步骤101：在纳米孔测序K-mer孔模型中输入一条具有T个碱基的核苷酸序列X＝x1,x2,...,xT，xT表示第T个碱基，从输入核苷酸序列的第一个碱基开始，每次移动一个碱基取出一个K-mer，直到移动到倒数第K个碱基位置结束，T个碱基共得到T-K+1个K-mer；
步骤102：纳米孔测序K-mer孔模型包含每个K-mer所对应的预期电流信号值，对照纳米孔测序K-mer孔模型依次查找T-K+1个K-mer所对应的预期电流信号值，生成输入核苷酸序列X对应的预期测序信号序列Y＝y1,y2,...,yT-5，其中yi表示从X中的位置i开始的K-mer所对应的预期电流信号值；
步骤103：根据实际测序信号序列中的信号值重复次数分布，将预期测序信号序列中的每个信号值重复多次，生成与实际测序信号序列之间具有相似长度分布的待测真实测序信号序列。
进一步地，所述信号处理模型构建模块还用于：
步骤201：将中值归一化后的待测真实测序信号序列I＝i1,i2,...,iT作为信号处理模型中第一层双向门控循环单元神经网络的输入向量，其中双向门控循环单元神经网络的基本单元由一个前向传播的门控循环单元和一个后向传播的门控循环单元组成，首先计算t时刻的前向输出向量计算公式为：
其中，it为t时刻的输入信号值，为t-1时刻的前向输出向量，σ为sigmoid函数且σ(z)＝1/(1+e-z)，表示第一中间变量，表示第二中间变量，表示第三中间变量，表示两个等维向量的对应元素相乘，为第一权重矩阵，为第二权重矩阵，Wf为第三权重矩阵，为第一偏置向量，为第二偏置向量，然后计算t时刻的后向输出向量计算公式为：
其中，为t+1时刻的后向输出向量，表示第四中间变量，表示第五中间变量，表示第六中间变量，为第四权重矩阵，为第五权重矩阵，Wr为第六权重矩阵，为第三偏置向量，为第四偏置向量，最后t时刻的输出向量ht是由t时刻的前向输出向量和t时刻的后向输出向量连接得到，即其中||表示向量的连接符号；
步骤202：将第一层双向门控循环单元神经网络的输出向量作为第二层双向门控循环单元神经网络的输入向量，第二层双向门控循环单元神经网络的输出向量作为第三层双向门控循环单元神经网络的输入向量，其中第二层双向门控循环单元神经网络、第三层双向门控循环单元神经网络的计算过程与第一层双向门控循环单元神经网络的计算过程相同，不同层之间的权重矩阵与偏置向量参数不同；
步骤203：将三层双向门控循环单元神经网络中最后一层的输出向量作为全连接层的输入向量进行处理，全连接层的输出向量为与输入的待测真实测序信号序列等长的滤波信号序列O＝o1,o2,...,oT。
进一步地，所述信号处理模型训练模块还用于：
步骤301：信号处理模型的监督训练数据包含待训练实际测序信号序列和与其对应的待训练真实测序信号序列，首先读取纳米孔测序平台输出的fast5文件得到待训练实际测序信号序列，然后使用连续小波动态时间规整算法计算每个待训练实际测序信号序列所对应的待训练真实测序信号序列；
步骤302：构建信号处理模型的损失函数为其中，cosh表示双曲函数，其中，oi表示信号处理模型的第i个输出信号，ri表示实际测序信号序列中第i个信号；
步骤303：将每层双向门控循环单元神经网络中权重矩阵以及偏置向量初始化为区间内的均匀分布，其中ni表示每层双向门控循环单元神经网络的输入向量维度，no表示每层双向门控循环单元神经网络的输出向量维度；
步骤304：使用Adam自适应优化器迭代计算得到最小化损失函数的模型参数，完成信号处理模型的参数训练过程。
更进一步地，所述步骤301包括：
步骤3011：读取纳米孔测序平台输出的fast5测序文件得到待训练实际测序信号序列，其中fast5文件为HDF5文件格式；
步骤3012：对待训练实际测序信号序列进行碱基识别得到测序读段，采用基因序列比对算法将测序读段比对到参考基因组序列，根据比对结果得到与待训练实际测序信号序列对应的参考基因组序列片段，然后使用参考基因组序列片段和纳米孔测序K-mer孔模计算出待训练预期测序信号序列；
步骤3013：使用连续小波动态时间规整算法完成待训练实际测序信号序列与待训练预期测序信号序列之间的点到点映射，根据映射结果，求得待训练实际测序信号序列中每个信号值所对应的预期信号值，进而计算出待训练实际测序信号序列所对应的待训练真实测序信号序列。
本发明的优点在于：本发明基于已知的纳米孔测序K-mer孔模型将输入核苷酸序列转换为预期测序信号序列，将预期测序信号序列中的每个信号值重复多次生成待测真实测序信号序列；然后建立真实测序信号序列处理模型，使用监督训练数据完成信号处理模型参数的训练；最后使用完成训练的信号处理模型对待测真实测序信号序列进行滤波处理，利用神经网络学习到实际测序信号序列的时频特性，进而能够准确地滤掉待测真实测序信号序列中与实际测序信号序列无关的高频分量并保留有用的高频分量。
附图说明
图1为本发明实施例所公开的一种基于神经网络的真实纳米孔测序信号滤波方法的示意图；
图2为本发明实施例所公开的一种基于神经网络的真实纳米孔测序信号滤波方法的实际测序信号序列中信号值重复次数分布示意图；
图3为本发明实施例所公开的一种基于神经网络的真实纳米孔测序信号滤波方法中门控循环单元神经网络结构示意图；
图4为本发明实施例所公开的一种基于神经网络的真实纳米孔测序信号滤波方法中三层双向门控循环单元神经网络结构示意图；
图5为本发明实施例所公开的一种基于神经网络的真实纳米孔测序信号滤波方法中信号处理模型参数训练过程中损失函数下降示意图；
图6为本发明实施例所公开的一种基于神经网络的真实纳米孔测序信号滤波方法中信号处理模型的监督训练数据示意图；
图7为现有技术的低通滤波器和本发明信号处理模型输出滤波信号波形对比示意图，图7中(A)为低通滤波器输出信号示意图，图7中(B)为本发明信号处理模型输出信号示意图；
图8为真实测序信号、DeepSimulator仿真测序信号、本发明仿真测序信号波形对比示意图，图8中(A)为真实测序信号波形示意图，图8中(B)为DeepSimulator仿真测序信号波形示意图，图8中(C)为本发明仿真测序信号波形示意图；
图9为DeepSimulator仿真测序信号、本发明仿真测序信号与真实测序信号的DTW距离对比示意图；
图10为真实测序信号、DeepSimulator仿真测序信号和本发明仿真测序信号连续小波变换时频图，图10中(A)为真实测序信号时频图，图10中(B)为DeepSimulator仿真测序信号时频图，图10中(C)为本发明仿真测序信号时频图；
图11为DeepSimulator仿真测序信号、本发明仿真测序信号与真实测序信号连续小波变换时频图之间的PCC比较示意图；
图12为真实测序读段、DeepSimulator仿真测序读段、本发明仿真测序读段的错误特性对比示意图；
图13为真实测序读段、DeepSimulator仿真测序读段、本发明仿真测序读段和肺炎克雷伯菌定制模型仿真测序读段的错误特性对比示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，一种基于神经网络的真实纳米孔测序信号滤波方法，主要以牛津纳米孔技术公司的MinION测序平台(R9.4测序芯片)测序得到的4000条Lambda病毒测序样本、4000条大肠杆菌测序样本、4467条不动杆菌测序样本、15178条肺炎克雷伯菌测序样本、16742条粘质沙雷氏菌测序样本和11047条金黄色葡萄球菌测序样本为例，将Lambda病毒测序样本和大肠杆菌测序样本作为训练数据集也即本发明所述的待训练真实测序信号序列，剩余测序数据作为测试数据集也即本发明所述的待测真实测序信号序列，所述方法包括以下步骤：
步骤S1：在纳米孔测序K-mer孔模型中输入核苷酸序列将其转换为与其对应的预期测序信号序列，将预期测序信号序列中的每个信号值重复多次生成待测真实测序信号序列；具体过程为：
步骤101：在纳米孔测序K-mer孔模型中输入一条具有T个碱基的核苷酸序列X＝x1,x2,...,xT，xT表示第T个碱基，从输入核苷酸序列的第一个碱基开始，每次移动一个碱基取出一个K-mer，直到移动到倒数第K个碱基位置结束，T个碱基共得到T-K+1个K-mer；需要说明的是，给定一个整数值K，从输入核苷酸序列的的第一个位置开始，取一连续K个碱基的短串，称之为K-mer，本发明实施例中，K取6，因此，本发明采用纳米孔测序6-mer孔模型，每次移动一个碱基取出一个6-mer。
步骤102：纳米孔测序K-mer孔模型包含每个K-mer所对应的预期电流信号值，对照纳米孔测序K-mer孔模型依次查找T-K+1个K-mer所对应的预期电流信号值，生成输入核苷酸序列X对应的预期测序信号序列Y＝y1,y2,...,yT-5，其中yi表示从X中的位置i开始的K-mer所对应的预期电流信号值；
步骤103：如图2所示，根据实际测序信号序列中的信号值重复次数分布，将预期测序信号序列中的每个信号值重复多次，生成与实际测序信号序列之间具有相似长度分布的待测真实测序信号序列。
步骤S2：构建基于双向门控循环单元神经网络的真实测序信号处理模型，该信号处理模型包括三层双向门控循环单元神经网络和一层全连接层，其中信号处理模型的输入为输入到第一层双向门控循环单元神经网络的经中值归一化的真实测序信号序列，信号处理模型的输出为全连接层输出的与输入的待测真实测序信号序列等长的滤波信号序列；具体过程为：
步骤201：将中值归一化后的待测真实测序信号序列I＝i1,i2,...,iT作为信号处理模型中第一层双向门控循环单元神经网络的输入向量，其中双向门控循环单元神经网络的基本单元由一个前向传播的门控循环单元和一个后向传播的门控循环单元组成，首先计算t时刻的前向输出向量如图3所示，计算公式为：
其中，it为t时刻的输入信号值，为t-1时刻的前向输出向量，σ为sigmoid函数且σ(z)＝1/(1+e-z)，表示第一中间变量，表示第二中间变量，表示第三中间变量，表示两个等维向量的对应元素相乘，为第一权重矩阵，为第二权重矩阵，Wf为第三权重矩阵，为第一偏置向量，为第二偏置向量，然后计算t时刻的后向输出向量计算公式为：
其中，为t+1时刻的后向输出向量，表示第四中间变量，表示第五中间变量，表示第六中间变量，为第四权重矩阵，为第五权重矩阵，Wr为第六权重矩阵，为第三偏置向量，为第四偏置向量，最后t时刻的输出向量ht是由t时刻的前向输出向量和t时刻的后向输出向量连接得到，即其中||表示向量的连接符号；
步骤202：如图4所示，将第一层双向门控循环单元神经网络的输出向量作为第二层双向门控循环单元神经网络的输入向量，第二层双向门控循环单元神经网络的输出向量作为第三层双向门控循环单元神经网络的输入向量，其中第二层双向门控循环单元神经网络、第三层双向门控循环单元神经网络的计算过程与第一层双向门控循环单元神经网络的计算过程相同，不同层之间的权重矩阵与偏置向量参数不同；
步骤203：将三层双向门控循环单元神经网络中最后一层的输出向量作为全连接层的输入向量进行处理，全连接层的输出向量为与输入的待测真实测序信号序列等长的滤波信号序列O＝o1,o2,...,oT。
步骤S3：获取纳米孔测序平台输出的待训练实际测序信号序列并根据待训练实际测序信号序列计算每个待训练实际测序信号序列所对应的待训练真实测序信号序列，将待训练实际测序信号序列和与其对应的待训练真实测序信号序列作为信号处理模型参数训练所需的监督训练数据，构建信号处理模型的损失函数，进行信号处理模型参数的初始化，通过自适应优化器最小化损失函数实现模型参数的训练；具体过程为：
步骤301：信号处理模型的监督训练数据包含待训练实际测序信号序列和与其对应的待训练真实测序信号序列，首先读取纳米孔测序平台输出的fast5文件得到待训练实际测序信号序列，然后使用连续小波动态时间规整算法计算每个待训练实际测序信号序列所对应的待训练真实测序信号序列；
步骤302：构建信号处理模型的损失函数为其中，cosh表示双曲函数，其中，oi表示信号处理模型的第i个输出信号，ri表示实际测序信号序列中第i个信号；
步骤303：将每层双向门控循环单元神经网络中权重矩阵以及偏置向量初始化为区间内的均匀分布，其中ni表示每层双向门控循环单元神经网络的输入向量维度，no表示每层双向门控循环单元神经网络的输出向量维度；
步骤304：使用学习率为0.0001的Adam自适应优化器迭代计算得到最小化损失函数的模型参数，完成信号处理模型的参数训练过程，训练期间将batch_size设置为64，迭代次数设置为1000，参数训练中损失函数下降过程参见图5。
其中，所述步骤301包括：
步骤3011：读取纳米孔测序平台输出的fast5测序文件得到待训练实际测序信号序列，其中fast5文件为HDF5文件格式；
步骤3012：对待训练实际测序信号序列进行碱基识别得到测序读段，采用基因序列比对算法将测序读段比对到参考基因组序列，根据比对结果得到与待训练实际测序信号序列对应的参考基因组序列片段，然后使用参考基因组序列片段和纳米孔测序K-mer孔模计算出待训练预期测序信号序列；
步骤3013：使用连续小波动态时间规整算法完成待训练实际测序信号序列与待训练预期测序信号序列之间的点到点映射，参见图6，根据映射结果，求得待训练实际测序信号序列中每个信号值所对应的预期信号值，进而计算出待训练实际测序信号序列所对应的待训练真实测序信号序列。
步骤S4：将待测真实测序信号序列输入到完成参数训练的信号处理模型中实现对待测真实测序信号序列的滤波处理。
为了验证本发明提出的基于双向门控循环单元神经网络的信号处理模型可以高效且准确地滤掉标准测序信号序列中无用的高频分量，分别使用DeepSimulator软件信号生成模块的低通滤波器和本发明构建的信号处理模型对标准测序信号序列进行了滤波处理，输出的滤波信号波形对比结果参见图7。根据滤波信号波形对比结果可以看到，与传统的低通滤波器相比本发明提出的基于双向门控循环单元神经网络的信号处理模型只对标准测序信号序列中6-mer变化处的信号进行滤波处理，这更符合实际纳米孔测序的原理。本发明所提出的信号处理模型利用神经网络可以准确地学习到标准测序信号序列与真实测序信号之间的内在关系。
通过在输出滤波信号上添加高斯噪声可以生成最终的仿真测序信号，参见图8，DeepSimulator软件与本发明都可以生成与真实测序信号波形相似的仿真测序信号。为了定量分析两种仿真测序信号与真实测序信号之间的相似度，使用动态时间规整算法进行分析，这是比较两个信号之间差异的标准方法。在测试数据集中随机选取1000条测序数据进行了测试，计算结果显示DeepSimulator软件生成的仿真测序信号与真实测序信号的平均归一化DTW距离为0.132，而本发明生成的仿真测序信号与真实测序信号的平均归一化DTW距离为0.121，比DeepSimulator软件降低了7.5％，这表明本发明可以生成更接近于真实测序信号的仿真测序信号。图9显示了1000条测序样本的实验测试结果，图中每一个点代表一条测序数据，对角线以上的点表示本发明仿真效果更好，而对角线以下的点表示DeepSimulator软件更好，根据统计结果，大约89.1％的测试样本表明本发明方法更好。
除了比较仿真测序信号与真实测序信号之间的DTW距离，同时采用连续小波变换算法对三类测序信号进行了时频分析，得到的连续小波变换时频图如图10所示。总体而言，两种仿真测序信号的时频图与真实测序信号的时频图都非常相似，但本发明生成的仿真测序信号时频图更接近于真实测序信号时频图，特别是高频部分。为了比较两种仿真测序信号时频图与真实测序信号时频图之间的相似度，分别计算了两种仿真测序信号时频图与真实测序信号时频图之间的皮尔逊相关系数(Pearson correlation coefficient，PCC)，并且进一步计算了信号时频图低频部分和高频部分处的PCC。使用测试数据集中4000条测序数据计算的平均值如图11所示。整体而言，本发明生成的仿真测序信号时频图与真实测序信号时频图之间的PCC较DeepSimulator提升了9.08％。对于时频图的高频部分，本发明较DeepSimulator提升了18.88％，这进一步说明了本发明提出的基于双向门控循环单元神经网络的信号处理模型可以更准确地滤掉标准测序信号序列中无用的高频分量，保留有用的高频分量，从而生成时域和频域上更接近于真实测序信号的仿真测序信号。
评估仿真测序信号质量的一个重要指标是比较由仿真测序信号生成的仿真读段是否与实际实验生成的测序读段具有相似的错误特性。使用牛津纳米孔技术公司最新的碱基识别软件Guppy分别对两种仿真测序信号和真实测序信号进行碱基识别，然后使用edlib软件计算三种测序读段的错误特性，包括插入率、删节率以及替代率。使用测试数据集提供的四种生物测序样品获得的统计结果如图12所示。由于实际测序读段的准确性与复杂的测序环境相关，因此仿真测序读段的准确性与实际测序读段的准确性有很大不同。尽管通过不同的测序信号仿真方法生成的两种仿真测序读段的准确性相似，但是在仿真读段中三类错误所占的比例并不一致。为了详细比较两种仿真测序读段和四种真实测序读段的错误特性，进一步计算了六种测序读段中三类错误所占的比例。除了肺炎克雷伯菌测序样本外，本发明所产生的仿真测序读段与实际测序读段具有相似的错误分布。然而，对于所有的生物测序样本而言，由DeepSimulator生成的仿真测序读段与实际测序读段之间的差异较大。针对仿真测序读段的分析结果表明，本发明提出的纳米孔测序信号仿真方法大大提高了仿真测序信号的质量。
对于肺炎克雷伯菌测序样本，由DeepSimulator和本发明提出的仿真方法生成的两种仿真测序读段都与实际测序读段之间有很大的差异。为了进一步探索本发明提出的测序信号仿真方法的可扩展性，我们从肺炎克雷伯菌测序样本中随机选择了4000个测序读段来训练提出的基于双向门控循环单元神经网络的信号处理模型，其余的测序读段用作测试数据集以验证肺炎克雷伯菌定制模型的性能。如图13所示，由肺炎克雷伯菌定制模型产生的仿真测序读段具有与真实测序读段相似的插入率和替代率，而仿真测序读段的删节率小于真实测序读段的删节率。与其他两类仿真测序读段相比，肺炎克雷伯菌定制模型产生的仿真测序读段的错误特性更接近于真实测序读段，这表明所提出的信号仿真方法可以有效地学习到真实测序信号的时频特性，并且比DeepSimulator具有更好的可扩展性。
通过以上技术方案，本发明基于已知的纳米孔测序K-mer孔模型将输入核苷酸序列转换为预期测序信号序列，将预期测序信号序列中的每个信号值重复多次生成待测真实测序信号序列；然后建立真实测序信号序列处理模型，使用监督训练数据完成信号处理模型参数的训练；最后使用完成训练的信号处理模型对待测真实测序信号序列进行滤波处理，利用神经网络学习到实际测序信号序列的时频特性，进而能够准确地滤掉待测真实测序信号序列中与实际测序信号序列无关的高频分量并保留有用的高频分量。
实施例2
与本发明实施例1相对应的，本发明实施例2还提供一种基于神经网络的真实纳米孔测序信号滤波装置，所述装置包括以下步骤：
待测真实测序信号序列生成模块，用于在纳米孔测序K-mer孔模型中输入核苷酸序列将其转换为与其对应的预期测序信号序列，将预期测序信号序列中的每个信号值重复多次生成待测真实测序信号序列；
信号处理模型构建模块，用于构建基于双向门控循环单元神经网络的真实测序信号处理模型，该信号处理模型包括三层双向门控循环单元神经网络和一层全连接层，其中信号处理模型的输入为输入到第一层双向门控循环单元神经网络的经中值归一化的真实测序信号序列，信号处理模型的输出为全连接层输出的与输入的待测真实测序信号序列等长的滤波信号序列；
信号处理模型训练模块，用于获取纳米孔测序平台输出的待训练实际测序信号序列并根据待训练实际测序信号序列计算每个待训练实际测序信号序列所对应的待训练真实测序信号序列，将待训练实际测序信号序列和与其对应的待训练真实测序信号序列作为信号处理模型参数训练所需的监督训练数据，构建信号处理模型的损失函数，进行信号处理模型参数的初始化，通过自适应优化器最小化损失函数实现模型参数的训练；
滤波处理模块，用于将待测真实测序信号序列输入到完成参数训练的信号处理模型中实现对待测真实测序信号序列的滤波处理。
具体的，所述待测真实测序信号序列生成模块还用于：
步骤101：在纳米孔测序K-mer孔模型中输入一条具有T个碱基的核苷酸序列X＝x1,x2,...,xT，xT表示第T个碱基，从输入核苷酸序列的第一个碱基开始，每次移动一个碱基取出一个K-mer，直到移动到倒数第K个碱基位置结束，T个碱基共得到T-K+1个K-mer；
步骤102：纳米孔测序K-mer孔模型包含每个K-mer所对应的预期电流信号值，对照纳米孔测序K-mer孔模型依次查找T-K+1个K-mer所对应的预期电流信号值，生成输入核苷酸序列X对应的预期测序信号序列Y＝y1,y2,...,yT-5，其中yi表示从X中的位置i开始的K-mer所对应的预期电流信号值；
步骤103：根据实际测序信号序列中的信号值重复次数分布，将预期测序信号序列中的每个信号值重复多次，生成与实际测序信号序列之间具有相似长度分布的待测真实测序信号序列。
具体的，所述信号处理模型构建模块还用于：
步骤201：将中值归一化后的待测真实测序信号序列I＝i1,i2,...,iT作为信号处理模型中第一层双向门控循环单元神经网络的输入向量，其中双向门控循环单元神经网络的基本单元由一个前向传播的门控循环单元和一个后向传播的门控循环单元组成，首先计算t时刻的前向输出向量计算公式为：
其中，it为t时刻的输入信号值，为t-1时刻的前向输出向量，σ为sigmoid函数且σ(z)＝1/(1+e-z)，表示第一中间变量，表示第二中间变量，表示第三中间变量，表示两个等维向量的对应元素相乘，为第一权重矩阵，为第二权重矩阵，Wf为第三权重矩阵，为第一偏置向量，为第二偏置向量，然后计算t时刻的后向输出向量计算公式为：
其中，为t+1时刻的后向输出向量，表示第四中间变量，表示第五中间变量，表示第六中间变量，为第四权重矩阵，为第五权重矩阵，Wr为第六权重矩阵，为第三偏置向量，为第四偏置向量，最后t时刻的输出向量ht是由t时刻的前向输出向量和t时刻的后向输出向量连接得到，即其中||表示向量的连接符号；
步骤202：将第一层双向门控循环单元神经网络的输出向量作为第二层双向门控循环单元神经网络的输入向量，第二层双向门控循环单元神经网络的输出向量作为第三层双向门控循环单元神经网络的输入向量，其中第二层双向门控循环单元神经网络、第三层双向门控循环单元神经网络的计算过程与第一层双向门控循环单元神经网络的计算过程相同，不同层之间的权重矩阵与偏置向量参数不同；
步骤203：将三层双向门控循环单元神经网络中最后一层的输出向量作为全连接层的输入向量进行处理，全连接层的输出向量为与输入的待测真实测序信号序列等长的滤波信号序列O＝o1,o2,...,oT。
具体的，所述信号处理模型训练模块还用于：
步骤301：信号处理模型的监督训练数据包含待训练实际测序信号序列和与其对应的待训练真实测序信号序列，首先读取纳米孔测序平台输出的fast5文件得到待训练实际测序信号序列，然后使用连续小波动态时间规整算法计算每个待训练实际测序信号序列所对应的待训练真实测序信号序列；
步骤302：构建信号处理模型的损失函数为其中，cosh表示双曲函数，其中，oi表示信号处理模型的第i个输出信号，ri表示实际测序信号序列中第i个信号；
步骤303：将每层双向门控循环单元神经网络中权重矩阵以及偏置向量初始化为区间内的均匀分布，其中ni表示每层双向门控循环单元神经网络的输入向量维度，no表示每层双向门控循环单元神经网络的输出向量维度；
步骤304：使用Adam自适应优化器迭代计算得到最小化损失函数的模型参数，完成信号处理模型的参数训练过程。
更具体的，所述步骤301包括：
步骤3011：读取纳米孔测序平台输出的fast5测序文件得到待训练实际测序信号序列，其中fast5文件为HDF5文件格式；
步骤3012：对待训练实际测序信号序列进行碱基识别得到测序读段，采用基因序列比对算法将测序读段比对到参考基因组序列，根据比对结果得到与待训练实际测序信号序列对应的参考基因组序列片段，然后使用参考基因组序列片段和纳米孔测序K-mer孔模计算出待训练预期测序信号序列；
步骤3013：使用连续小波动态时间规整算法完成待训练实际测序信号序列与待训练预期测序信号序列之间的点到点映射，根据映射结果，求得待训练实际测序信号序列中每个信号值所对应的预期信号值，进而计算出待训练实际测序信号序列所对应的待训练真实测序信号序列。
以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。
- 上一篇：一种医用注射器针头装配设备
- 下一篇：一种基于分治法的mRNA序列优化的方法与装置