一种SNP序列的填充方法与流程

一种snp序列的填充方法
技术领域
1.本发明涉及一种snp序列的填充方法，属于生物技术领域。


背景技术：

2.单核苷酸多态性（snp）主要是指在基因组水平上由单个核苷酸的变异所引起的dna序列多态性。它是人类可遗传的变异中最常见的一种，占所有已知多态性的90%以上。snp在人类基因组中广泛存在，平均每1000个碱基对中就有1个，其总数可达300万个甚至更多。snp是一种二态的标记，由单个碱基的转换或颠换所引起，也可由碱基的插入或缺失所致。snp对群体遗传学、制药业、法医学、癌症及遗传性疾病甚至进化的研究都将产生不可估量的影响。
3.snp(单核苷酸多态性标记)芯片测序的过程中导致的基因数据的丢失给全基因组关联分析研究带来很大的挑战，基因型数据的丢失分为遗传性丢失和检测性丢失。我们在基因型缺失的分析过程中，一般讨论的是技术性缺失，而不是人为的缺失，主要有下列原因导致：全基因组重测序导致的缺失、简化基因测序导致的缺失、外显子测序以及目标区域捕获测序导致的缺失以及snp芯片导致的缺失等。
4.目前普遍通过带有缺失值的基因序列拟合一个参数，学习缺失数据的总体特征，然后根据特征对缺失值进行填充，这种方式需要数据缺失值对数据整体的分布产生一个比较小的影响，但是当下的基因样本数量还不足以支持如此大的数据量，导致填充效率低下且预测得到的基因填充值错误率高。


技术实现要素：

5.本发明所要解决的技术问题是，克服现有技术的缺点，提供一种snp序列的填充方法，该方法简单易行，填充具有较高的准确率。
6.为了解决以上技术问题，本发明提供一种snp序列的填充方法，具体包括以下步骤：步骤一：获取已有基因序列原始数据，原始数据包括版本号为38的人类全基因组参考序列和所有染色体上的snvindels文件；步骤二：对获取的数据预处理；步骤三：根据数据构建填充神经网络模型；步骤四：将含有缺失值的snp序列输入到填充神经网络模型，实现snp序列缺失数据的填充。
7.本发明进一步限定的技术方案是：进一步的，前述snp序列的填充方法中，步骤二中数据预处理具体为：1）确定原始数据人类全基因组参考序列文件中存储的碱基位置索引；2）提取snp序列数据；3）确定高频率的snp序列数据，提取出不含缺失值的序列数据作为标签序列；
4）生成训练样本数据，将其作为神经网络模型的输入序列，通过不断优化神经网络模型的输出结果与标签序列的误差，实现神经网络模型中各参数的确定，然后利用该神经网络模型，实现对含有缺失值的序列数据进行预测。前述snp序列的填充方法中，步骤三中填充神经网络模型由卷积神经网络cnn、循环神经网络rnn及连通时序分类器ctc组成。
8.前述snp序列的填充方法中，步骤四的具体操作为：1）输入一维含缺失值的snp序列数据，首先进行卷积神经网络cnn层的卷积运算；2）将步骤1）的结果输出到具有双向lstm结构的rnn层；3）对rnn层的输出结果进行concat后，输入到全连接层中；4）对全连接层的输出结果进行ctc解码，得到与输入序列长度完全相同的预测结果即填充好的snp序列。
9.前述snp序列的填充方法中，rnn层网络结构采用4个双向lstm层串联构造。
10.本发明的有益效果是：现有的方法都是基于传统意义上的统计分析方法，如时间序列分析方法，主成分分析方法等，本发明基于cnn、rnn和ctc组合构架的神经网络模型，该模型对人类基因组参考序列中的常染色体上的snp序列的填充具有较高的准确率，特别是当snp位点缺失率0.2≤r≤0.8时，其填充准确率在80%左右。
附图说明
11.图1为本发明实施例snp序列的填充方法的流程图；图2为图1中原始数据预处理的具体流程图；图3为本发明实施例snp序列的填充方法中神经网络模型的示意图；图4为本发明实施例snp序列的填充方法中神经网络模型的另一种示意图。
具体实施方式
12.实施例1本实施例提供的一种snp序列的填充方法，流程如图1-2所示，具体包括以下步骤：步骤一：获取已有基因序列原始数据；基于现有ncbi或1000genome网站：到https://www.ncbi.nlm.nih.gov/或https://www.internationalgenome.org/网站上获取版本号为38的人类全基因组参考序列和所有染色体上的snvindels数据文件；人类全基因组参考序列文件的文件名为：gca_000001405.15_grch38_full_analysis_set.fna.gz，简称g文件；对应染色体编号为i(i=1,2,...,22,x,y)的snvindel文件名为：all.chri.shapeit2_integrated_snvindels_v2a_27022019.grch38.phased.vcf.gz，简称si文件。
13.si文件中记录了编号为i的染色体所有的单核苷酸变异（snp）、基因组结构性变异（sv）和在基因组的某个位置上发生的小片段序列的插入或删除（indel是insertion和deletion的简称）的所有数据。例如，sx文件，即是记录x号染色体的svnindel文件名为all.chri.shapeit2_integrated_snvindels_v2a_27022019.grch38.phased.vcf.gz。
14.染色体文件中包含26个区域的2548个样本数据，https://www.internationalgenome.org/网站首页给出了这些地区和样本的详细说明；步骤二：对获取的数据预处理；1）确定原始数据全基因组序列文件中存储的碱基位置索引，确定snp位点对应的碱基在全基因组序列中所在的位置；解压g文件，得到24个染色体的全基因组参考序列数据，分别是chr1.fa，chr2.fa，......，chr22.fa，chrx.fa和chry.fa文件，对于编号为i的染色体全基因组文件chri.fa，简称ci文件；以x号染色体为例，读取cx文件，同时解压并读取sx文件，发现cx中若碱基的位置索引从1开始编号，则第12568位置的碱基为c，第13587位置的碱基为t，正好对应于sx文件中的第一个和第三个snp位置上对应的ref碱基，即染色体si文件中存储的碱基位置对应的是其全基因组序列文件ci中存储的碱基位置从1开始编号的数据；2）提取snp序列数据；解压并打开步骤一中的si文件，遍历文件中每一行的有效数据，若ref和alt条目下对应的数据都是一个碱基(这种情形称为snp，即本实施例方法研究的只是单核苷酸多态性，不考虑结构性变异sv和插入删除indel对应的数据)，则需要存储当前碱基位置索引对应的所有样本等位基因数据（因为任一样本数据对应的等位基因的存储格式都为“a|b”的形式(a=0或1，b=0或1)，故需要确定是保存a还是保存b），不失一般性，保存左边的a数据，则对于所有的样本数据，在确定snp序列中某一个碱基位置索引对应的数据时，只需保存a数据即可，当a=0时，表示对应的是ref碱基；当a=1时，表示对应的是alt碱基；3）确定高频率的snp序列数据；统计每一个snp位点的频数，即对于所有的2548个样本数据，若某一个snp位点的a=1的个数f大于2548/2=1274，则其频数为2548-f，否则为f；重新生成一个新的高频率的snp序列数据，按照snp位点索引从小到大的顺序，逐一判断当前snp位点的频数f，若当前位点的频数f不小于指定的频数阈值th（这里取th=500），则将当前位点对应的所有样本数据提取出来，写到要保存的高频率的snp序列数据文件中，令n表示最终找出的snp位点频数不小于th的个数，则生成的高频率snp序列数据为n行2548列的矩阵，矩阵中的元素值为0或1，0表示与当前行对应的ref碱基相同，1表示与当前对应的alt碱基相同，记当前矩阵为l，l可表示当前所有样本在不同snp高频位点的标签值，这样就提取出不含缺失值的序列数据作为标签序列；4）生成训练样本数据，将其作为神经网络模型的输入序列，通过不断优化神经网络模型的输出结果与标签序列的误差，实现神经网络模型中各参数的确定，然后利用该神经网络模型，实现对含有缺失值的序列数据进行预测；定义snp位点缺失率r，计算出丢失数据的位点总个数为n=int(n * r)，然后从[1,2,...,n]序列中随机生成n个不重复的数，并将这n个不重复的数按照从小到大的顺序排列，组成一个表征丢失位点索引的序列d=[i1,i2,...,in]，最后生成一个表征训练样本数据矩阵t作为输入序列数据，t的大小与l相同，t矩阵中的数据在i1，i2，...，in行的数据为0.5，其余行的数据与l相同；训练时：输入神经网络模型的数据是t，因为t中含有0.5，神经网络模型的目的是将t矩阵中的0.5变为0或1，神经网络模型输出值也叫预测值，l中的数值叫真实值，也叫做
标签序列，目的是通过训练神经网络模型，使得预测值与真实值之间的误差（在线性等权最小二乘意义上的误差）最小；步骤三：根据数据构建填充神经网络模型；填充神经网络模型由卷积神经网络、循环神经网络及连通时序分类器组成，如图3-4所示；按图3方式cnn与rnn连接，concat指的是将正反向的rnn结果进行拼接，全连接层（fully connected layer）的输出结果是与ctc decoder 连接的；步骤四：将含有缺失值的snp序列输入到填充神经网络模型，实现snp序列缺失数据的填充，具体操作为：1）由图3可看出，输入一维含缺失值的snp序列数据（每次输入就是矩阵t中每一列的数据）于训练好的经网络模型，首先进行卷积神经网络cnn层的卷积运算；2）将步骤1）的结果输出到具有双向lstm结构的rnn层；3）对rnn层的输出结果进行concat后，输入到全连接层中；4）对全连接层的输出结果进行ctc解码，得到与输入序列长度完全相同的预测结果即填充好的snp序列；图4是对图3的详细描述，它详细给出了本发明使用的cnn模型是基于残差结构的残差网络模型，每个残差块b可表示为b=f(b1+b2)，其中b1是将输入当前残差块的序列依次进行1x1卷积、1x3卷积和1x1卷积后得到的结果，b2是输入当前残差块的序列进行1x1卷积的结果，f表示对b1与b2的累加结果进行relu激活函数运算，本发明中的rnn网络结构是采用4个双向lstm层串联的构造。
[0015]
本实施例基于cnn、rnn和ctc组合构架的神经网络模型，该模型对人类基因组参考序列中的常染色体上的snp序列的填充具有较高的准确率，本实施例在不同缺失率下对应的填充准确率见表1所示；表1不同缺失率情形下的填充准确率缺失率20%30%40%50%60%70%80%填充准确率82.3%81.7%82.1%79.7%77.6%77.5%78.2%如表1所示，从表1中可以看出，随着缺失率的增加，填充准确率整体呈现下降趋势，但并不是呈现单调下降的趋势，这也体现了人类全基因组的snp序列具有随机性的一面，当snp位点缺失率0.2《r《0.8时，其填充准确率在80%左右。
[0016]
除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。