本发明属于拷贝数变异技术领域，尤其涉及一种基于状态转移模型的新一代测序拷贝数变异仿真方法。
背景技术：
拷贝数变异是由基因组发生重排而导致的，一般指长度为1kb以上的基因组大片段的拷贝数增加或者减少，主要表现为亚显微水平的缺失(deletion)和重复(insertion)。拷贝数变异是基因组结构变异(Structural variation，SV)的重要组成部分，它的位点突变率远高于SNP(Single nucleotide polymorphism)，是人类疾病的重要致病因素之一。按照发生场所的不同，可以将拷贝数变异分为生殖细胞拷贝数变异(Copy number variation，CNV)和体细胞拷贝数变异(Copy number alternation，CNA)，顾名思义就是他们的发生场所分别为生殖细胞和体细胞。CNV具有遗传效应，CNA没有遗传效应，这是由它们所处细胞的机制所决定的。拷贝数变异的仿真就是设计仿真算法，用程序对拷贝数变异的过程进行模拟。近年来，基因组测序领域发展迅速，这非常有助于对许多生物系统的理解。在过去的五年中，计算机生物学家和生物信息学专家针对发现、分析和解释不同的基因组变异的高通量测序数据，提出了新的、更好的和更有效的检测拷贝数变异的工具。在使用检测工具时，可靠的模拟数据集是必不可少的，模拟数据的获得是测试新开发检测工具的第一步。虽然目前已经有很多可用的拷贝数变异仿真工具，但是这些工具的功能都不是很全面，要么就是只可以模拟CNV和CNA中的一个功能，要么就是没有一个可信的状态转移模型，要么就是只有序列生成部分。因此，开发一个有效的关于CNV和CNA的模拟器和序列生成器是必要的，它要能够模拟拷贝数变异且考虑到真实生物样品的错误率。不同的下一代测序仪所生成的reads的length和error profile也不同，目前最流行的测序数据是从Illumina测序平台所产生的，它采用了化学方法来进行序列合成并生成reads，要开发的高效模拟器正是基于Illumina测序平台的。正是由于Illumina平台产生的数据受欢迎且应用广泛这个特点，任何其它的测序平台通过提供一个特定的错误配置文件就可以对其进行使用。目前可用的仿真软件可以生成基于特定平台的相关错误配置文件的reads，也可以跨平台生成reads。已经存在的一些仿真软件都有各自的优点，但是同时它们也存在着一些缺陷。现有的仿真软件最大的缺陷是不能能同时仿真CNV和CNA，下面针对一些仿真软件的性质和功能分别加以说明。SInc是用C语言开发的一款仿真软件，它是开源的，拥有CLI接口，也有自己的error model，但是它存在的问题没有仿真CNA变异；MetaSim是用JAVA语言开发的一款仿真软件，它拥有CLI和GUI接口，不是开源的，可以仿真pair-end数据，它的缺陷是只有序列生成部分而没有变异仿真部分，即没有将quality value赋值给reads；FlowSim是用Haskell语言编写的，它拥有CLI接口，是一款开源仿真软件，它的缺陷是没有变异仿真部分且不能仿真Illumina平台的数据，即不能仿真paie-end数据；GenFrag仿真软件是开源的，有CLI接口，它的缺点是没有变异仿真部分且它的erroe model过分简单；DwgSim仿真软件是由变异仿真和序列生成两个部分构成的，它有CLI接口且是开源的，它的缺陷是不能模拟真实数据。当然，这些仿真软件共有的一个缺陷是不能仿真CNA变异。
现有拷贝数变异仿真软件的实现方法存在存在以下问题：没有将生殖细胞和体细胞的两种拷贝数变异集成在一起，即只可以仿真单个的CNV，不能同时仿真CNA,这就导致了仿真软件的功能不全面，仿真出来的数据比较片面，不是完整的拷贝数变异后的数据，限制了用户的使用；没有合适的model来确定拷贝数变异的各个状态，因为拷贝数变异中各个状态之间的转换是符合某种转换机制的，要是没有相应的model的话，仿真出来的数据和真实数据的差异较大，仿真结果的可信度将受到影响；没有将拷贝数变异的两种形式变异仿真和序列生成集成在一块，一般存在这种问题的仿真软件大都只有序列生成的部分，没有变异仿真的部分，即它能生成最终的fq文件，但是将生物变异的部分省略，这样显然是不符合实际情况的，因为拿到的真实样本不一定是完全没有发生拷贝数变异的，相反发生变异的比例还很大，所以加上变异仿真这一步是相当有必要的。
技术实现要素：
本发明的目的在于提供一种基于状态转移模型的新一代测序拷贝数变异仿真方法，旨在解决为拷贝数变异检测提供合适的模拟数据的问题。
本发明是这样实现的，一种基于状态转移模型的新一代测序拷贝数变异仿真方法，所述基于状态转移模型的新一代测序拷贝数变异仿真方法采用拷贝数变异仿真算法；在仿真算法中增加状态转移模型和序列生成部分；
所述拷贝数变异包括CNV和CNA；
基于Illumina测序平台的Profile文件的生成，核心步骤是将fq文件的reads说明部分的ASCii码转换成碱基的quality value，相应方法是对应字符的ASCii码减去33；
将仿变异真后的fa文件和生成的profile文件作为输入，设置read length，利用多线程和序列生成算法，生成并输出最终的fq文件。
进一步，所述CNV仿真算法和包括：
(a)确定发生CNV变异的位置、尺寸、类型；
(b)根据a中确定的CNV变异的参数执行CNV变异，并打印变异参数的记录文件和变异后的fa文件。
进一步，所述CNV状态转移模型为：
Normal:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pa＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pa)*Pnn
Pn＝1-Pa-Pd
Insertion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pd＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pd)*Paa
Pa＝1-Pn-Pd
Deletion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pn＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pn)*Pdd
Pd＝1-Pa-Pn。
进一步，所述CNA仿真算法包括：
(a)确定发生CNA变异的位置、尺寸、类型；
(b)根据a中确定的CNA变异的参数，执行CNA变异，并打印变异参数的记录文件和变异后的fa文件。
进一步，所述CNA状态转移模型为：
Normal:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pa＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pa)*Pnn
Pn＝1-Pa-Pd
Insertion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pd＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pd)*Paa
Pa＝1-Pn-Pd
Deletion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pn＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pn)*Pdd
Pd＝1-Pa-Pn。
本发明的另一目的在于提供一种应用所述基于状态转移模型的新一代测序拷贝数变异仿真方法的CNV和CNA模拟器。
本发明的另一目的在于提供一种应用所述基于状态转移模型的新一代测序拷贝数变异仿真方法的CNV和CNA序列生成器。
本发明提供的基于状态转移模型的新一代测序拷贝数变异仿真方法，能够解决现有拷贝数变异的仿真程序只考虑了CNV或CNA的情况，将CNA和CNV集成在同一个仿真软件中，并拥有自己独特的拷贝数变异(包括CNV和CNA)仿真算法；在CNV和CNA仿真算法的基础上还增加了状态转移模型，有了这个model之后，整个仿真的过程就更加具有可信度。没有model的时候，拷贝数变异类型是由程序设定按一定比例产生的，一般是发生缺失的比例与发生插入的比例为4:1：有状态转移模型的时候，下一状态发生变化的类型和上一状态有关，至于有什么关系，主要取决于model的设定。下面分别是没有状态转移模型和有状态转移模型时候拷贝数变异记录文件的对比图，由对比图可以发现，没有状态转移模型的时候，变异类型状态之间的转换服从稳定的4:1的比例，但是加上状态转移模型之后，变异状态之间的比例不一定服从4:1，这更符合实际数据，因为真实数据发生变异的情况受环境等多种因素影响，不可能服从一个稳定的变化比率，这需要反复训练真实数据来获得。在CNA和CNV及其model的基础之上，增加了序列生成部分，使得这个仿真软件不仅拥有变异仿真功能，还拥有序列生成功能。
本发明采用状态转移模型，将状态转移模型加在拷贝数变异的仿真算法中，在此基础上，还应加上序列生成的功能；在实现仿真拷贝数变异的过程中，不仅仿真了生殖细胞的拷贝数变异，还仿真了体细胞的变异。
本发明在仿真CNV和CNA的基础上，还加入了序列生成的部分，模拟了序列生成部分中可能发生的错误机制，增加了数据的真实性。综合了CNV和CNA于一个仿真软件中，方便了用户的使用和研究，并且仿真算法中加入了状态转移模型，使得生物变异的仿真更加具有可信度。除此之外，本发明还加入了序列生成部分，考虑了真实测序中可能发生的错误，生成错误配置文件，进而生成最终的reads。总而言之，和其它仿真软件相比，本拷贝数变异仿真软件功能更加完整，数据更加贴近真实数据。下面三幅图是本发明最终的成果图，由图可以看到本仿真软件有仿真CNV和CNA的功能，有生成错误配置文件的功能，有生成最终的reads的功能及一些其它的附加功能，功能完善齐全，还有完整的GUI(图形用户界面)，使用起来方便快捷。
附图说明
图1是本发明实施例提供的基于状态转移模型的新一代测序拷贝数变异仿真方法流程图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
下面结合附图对本发明的应用原理作详细的描述。
如图1所示，本发明实施例提供的基于状态转移模型的新一代测序拷贝数变异仿真方法包括以下步骤：
S101:CNV仿真算法和状态转移模型的设计；
S102:CNA仿真算法和状态转移模型的设计；
S103：基于Illumina测序平台的Profile文件的生成，核心步骤是将fq文件的reads说明部分的ASCii码转换成碱基的quality value，相应方法是对应字符的ASCii码减去33；
S104:将仿变异真后的fa文件和生成的profile文件作为输入，设置合适的read length，利用多线程和序列生成算法，生成并输出最终的fq文件。
下面结合具体实施例对本发明的应用原理作进一步的描述。
本发明的实施例在拷贝数变异生物特性及仿真算法的基础之上，建立状态转移模型，经过反复对仿真数据进行训练，设置合适的状态转移模型，对仿真算法进行改进。
本发明实施例的技术方案
(1)CNV仿真算法和状态转移模型的设计
CNV仿真算法:
(a)确定发生CNV变异的位置、尺寸、类型；
(b)根据a中确定的CNV变异的参数执行CNV变异，并打印变异参数的记录文件和变异后的fa文件。
CNV状态转移模型：
Normal:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pa＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pa)*Pnn
Pn＝1-Pa-Pd
Insertion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pd＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pd)*Paa
Pa＝1-Pn-Pd
Deletion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pn＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pn)*Pdd
Pd＝1-Pa-Pn
(2)CNA仿真算法和状态转移模型的设计
CNA仿真算法:
(a)确定发生CNA变异的位置、尺寸、类型；
(b)根据a中确定的CNA变异的参数，执行CNA变异，并打印变异参数的记录文件和变异后的fa文件。
CNA状态转移模型：
Normal:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pa＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pa)*Pnn
Pn＝1-Pa-Pd
Insertion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pd＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pd)*Paa
Pa＝1-Pn-Pd
Deletion:
Paa＝Pa Pnn＝Pn Pdd＝Pd
Pn＝Paa*Pnn*Pdd/(2-Paa*Pnn*Pdd)
Pd＝(1-Pn)*Pdd
Pd＝1-Pa-Pn
(3)Profile文件的生成
本发明的测序数据是基于Illumina测序平台，所以生成profile文件时所用的fq文件也是由Illumina测序平台产生。Profile文件其实是统计序列中某个碱基出现的次数，是根据它的quality value来判断它出现的次数的，所以在生成profile文件之前应该将fq文件中reads的序列说明部分的ASCii码转换成碱基的quality value，相应的方法是对应字符的ASCii码减去33。
(4)Reads的生成
将仿变异真后的fa文件和生成的profile文件作为输入，设置合适的read length，利用多线程和序列生成算法，生成并输出最终的文件。
(5)算法的性能评估
对设计的仿真算法和状态转移模型进行评估和改进，对序列生成过程进行优化，形成最终的仿真软件。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。