CN110491441A - 一种模拟人群背景信息的基因测序数据仿真系统及方法 - Google Patents

Abstract

本发明公开了一种模拟人群背景信息的基因测序数据仿真系统及方法，输入目标捕获区域文件、测序深度、七种变异类型、变异发生频率和变异在参考基因组上的坐标，根据测序深度确定模板数量，使用接受拒绝算法产生相应模板长度的概率分布，遍历每一条模板并对当前已经遍历完的模板数进行判断，当算法尚未遍历所有模板时，分别对取出的长度模板进行拷贝数变异模拟、单核苷酸位点变异模拟、基因融合模拟、串联重复模拟、倒位变异模拟、插入片段模拟和缺失片段模拟，生成读段并写入测序文件；当已经遍历完所有模板时，测序文件生成完毕；进行读段比对输出模拟的测序文件及其比对文件，仿真结束。本发明可以简便和快速的获得包含特定变异样本。

Description

一种模拟人群背景信息的基因测序数据仿真系统及方法

技术领域

本发明属于以精准医学为应用背景的数据科学技术领域，具体涉及一种模拟人群背景信息的基因测序数据仿真系统及方法。

背景技术

精准诊疗是现代医学发展的主流方向。精准诊疗的基础是基因大数据分析。近年来，随着世界各国的基因大数据计划纷纷落地，大数据快速积累，各种数据分析软件也应运而生。这其中，变异检测是大数据分析的基础，主流的变异检测软件有Samtools、GATK、Pindel、Delly等数十种。但是，在临床应用中，这些变异检测软件的准确率仍然有待提高，面对多种多样的数据特点和临床任务，假阳性和假阴性的检测结果十分常见。理论研究和实践普遍认为，这些变异检测软件存在如下共性问题：

(1)在参数调试阶段，受到样本提取成本高、突变稀有等客观因素的制约，变异检测软件的开发者缺乏类型多样的样本对变异检测软件进行调试，特别是参数调试。

(2)在软件测试阶段，由于拿到的测试样本中所包含突变的真实情况无从知晓，因而对变异检测软件的准确性无法做出全面和准确的衡量。这在很大程度上影响了下游的基因变异分析软件的选取和参数调试。

因此，迫切需要高通量测序数据模拟软件，来模拟生成大量的测试数据，解决上述两个共性问题。由于客观需求旺盛，目前常见的测序数据模拟软件有：454sim、ART、BEAR、SinC、wgsim、dwgSIM、Grinder、Mason、SeqMaker、GemSIM、simseq、Flowsim、pbsim、genfrag、metasim、pIRS、bamsurgeon等。然而，上述软件大多数仅针对特定场景。与本发明面向的场景相同的软件的特点及其主要不足如下：bamsurgeon的使用需要输入一个比对文件，通过对比对文件的直接修改来生成含有特定变异的比对文件，但是其参数设置不够灵活，局限性较大。GemSIM只支持单核苷酸位点变异的模拟，功能单一。dwgsim支持单核苷酸位点变异、小片段插入删除(英文名称：insert and deletion，英文缩写：indel)、染色体倒位变异和基因融合变异的模拟，但是不支持基因拷贝数变异和串联重复变异的模拟。SinC支持基因拷贝数变异的模拟，以及基于此的单核苷酸位点变异和小片段插入删除的模拟，但是不支持基因融合变异、染色体倒位变异以及串联重复变异的模拟。SeqMaker支持单核苷酸位点变异、小片段插入删除、基因融合变异、拷贝数变异以及倒位变异的模拟，但是不支持大片段的插入、复杂结构变异(英文名称：Complex structural variant，英文缩写：CSV)、串联重复变异的模拟。

基于文献检索的结果，目前没有软件能够全面支持已知的所有主要变异类型，也不能模拟出人群多态，不支持用户按照自身需要训练模板长度分布、相邻位点深度分布、整体深度分布以及质量值分布等几个主要的数据特征；模板是指将参考基因组随机打断成长度为几十到几百个碱基对不等的碱基序列片段。此外，面对海量数据需求，现有软件不具备特定目标精度下样本批量生成、验证变异检测软件的特异度和敏感度指标的功能。目标精度是拷贝数变异的拷贝倍数或者单核苷酸位点变异、插入突变、删除突变、基因融合变异、串联重复变异、倒位变异的频率值。

发明内容

本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种模拟人群背景信息的基因测序数据仿真系统及方法，用于模拟生成已知突变结果的测序数据并验证变异检测软件的技术性能指标。

本发明采用以下技术方案：

一种模拟人群背景信息的基因测序数据仿真方法，输入目标捕获区域文件、测序深度、单核苷酸位点变异、插入突变、删除突变、拷贝数变异、倒位变异、基因融合变异、串联重复变异、变异发生频率和变异在参考基因组上的坐标，加载参考基因组和目标捕获区域文件；目标捕获区域文件中记录参考基因组上用户关注的每个目标区域的起始坐标和终止坐标；根据测序深度确定模板数量；通过线性回归法构建模型，使用接受拒绝算法产生相应模板长度的概率分布，遍历每一条模板并对当前已经遍历完的模板数进行判断，当算法尚未遍历所有模板时，先根据模板长度的概率分布确定当前模板的长度，然后从参考基因组取出该长度模板，分别对取出的长度模板进行拷贝数变异模拟、单核苷酸位点变异模拟、基因融合模拟、串联重复模拟、倒位变异模拟、插入片段模拟和缺失片段模拟，最后生成读段，并将读段写入测序文件；当已经遍历完所有模板时，测序文件生成完毕；最后，对测序文件进行读段比对并生成比对文件，比对结束之后输出模拟的测序文件及其比对文件，仿真结束。

具体的，当尚未遍历所有模板时，首先在目标区域生成随机断点，然后将生成的断点定位到参考基因组的对应位置，采用接受拒绝算法在该位置前后1000碱基对范围内生成拟合相应混合高斯模型的新断点，再使用建索引的方式实现从参考基因组的该位置处快速取出模板，建索引的方式为：先创建字典记录参考基因组每行序列所在相对位置，然后在读取时使用文件指针在查阅字典后直接跳到要读取的模板处。

具体的，当尚未遍历当前模板的所有碱基时，以每个碱基在参考基因组中的坐标作为随机种子，使用梅森旋转算法生成随机数，比对千人基因组数据库相应位点，当该位置随机数小于数据库的人群发生率，则发生变异；当该位置随机数大于人群发生率，则不对该碱基做任何改变；当已经遍历完当前模板的所有碱基时，人群背景信息的模拟结束。

4、根据权利要求1所述的模拟人群背景信息的基因测序数据仿真方法，其特征在于，已经遍历当前模板的所有碱基时，结束人群背景信息模拟并依次进行拷贝数变异模拟、单核苷酸位点变异模拟基因融合模拟、串联重复模拟、倒位变异模拟、插入片段模拟和缺失片段模拟，然后采用内维尔算法生成测序碱基质量并写入对应读段，最后将该读段写入以fastq标准格式定义的输出文件，重新返回判断是否遍历所有模板。

一种模拟人群背景信息的基因测序数据仿真系统，包括单样本变异模拟模块、质量分布统计模块、相邻位点深度分布训练模块、整体深度分布统计模块、模板长度分布训练模块和特定目标精度下批量样本生成与变异检测软件验证模块；单样本变异模拟模块用于单核苷酸位点变异模拟、插入片段模拟、删除片段模拟、拷贝数变异模拟、基因融合变异模拟、串联重复变异模拟和倒位变异模拟；特定目标精度下批量样本生成与变异检测软件验证模块用于批量生成正常样本和变异样本，并对变异检测软件进行验证；通过将变异写在文件中，然后输入文件路径并启动，实现同时模拟多种变异。

具体的，单核苷酸位点变异模拟由用户指定变异类型为单核苷酸位点变异，然后依次输入要模拟的单核苷酸位点变异所在参考基因组上的坐标、发生突变之前的碱基和发生突变之后的碱基，以及变异发生率；

插入片段模拟由用户在命令行参数中输入插入片段长度或插入片段的碱基序列，并指定变异类型为插入，然后依次输入变异在参考基因组上的坐标和变异的发生率；

删除片段模拟由用户在命令行参数中指定变异类型删除，然后依次输入需要删除片段在参考基因组上的坐标、变异发生率以及发生删除后是否随机插入片段，若是，在命令行参数中指定随机插入片段的长度和该片段的碱基序列；

拷贝数变异模拟由用户在命令行参数中指定变异类型为拷贝数变异，然后依次输入要模拟的拷贝数变异在参考基因组上的坐标和拷贝数；当要模拟拷贝数增加，则输入一个大于1的整数或浮点数，当模拟拷贝数减少，则输入一个小于1的浮点数；或者在命令行参数中输入原始比对文件的路径以及发生的拷贝数变异在参考基因组上的坐标和拷贝数，实现直接对真实样本的比对文件进行修改来模拟拷贝数变异的功能；

基因融合变异模拟先指定变异类型为融合，然后先依次输入发生融合的左边基因的参考基因组坐标和正负链向，再依次输入发生融合的右边基因的参考基因组坐标和正负链向；最后再输入基因融合变异发生率以及发生融合时是否在融合位置随机插入片段，在命令行参数中指定随机插入片段的长度和该片段的碱基序列；

串联重复变异模拟由用户先在命令行参数中指定变异类型为串联重复；然后依次输入要模拟的串联重复变异在参考基因组上的坐标、重复次数以及串联重复变异发生率；

倒位变异模拟由用户在命令行参数中指定变异参数为倒位，再依次输入要模拟的倒位变异在染色体上的坐标和倒位变异发生率。

具体的，特定目标精度下批量样本生成模块用于为变异检测软件分梯度提供批量测试样本，产出阴性标准集和阳性标准集；并生成样本信息记录文件；目标精度是拷贝数变异的拷贝倍数或者单核苷酸位点变异、插入突变、删除突变、基因融合变异、串联重复变异、倒位变异的频率值；梯度是目标精度梯度变化下呈现出的不同的频率或不同的拷贝倍数；然后，测试样本经变异检测软件检测后，将变异检测软件的检测结果填入样本信息记录文件，输入文件路径并调用变异检测软件验证功能，显示变异检测软件在对应目标精度下的敏感度和特异度。

具体的，仿真系统还包括人群多态模拟、扩增错误模拟、测序错误模拟、单端测序和双端测序模拟、接头模拟和单分子索引模拟；

人群多态模拟以千人基因组数据库为参照，预置突变随机发生实现，以每个碱基在参考基因组中的坐标做为随机种子并使用梅森旋转算法产生高质量的随机数，每次模拟样本产生一个随机数，使每个样本具有惟一的随机性，再加上每个碱基在参考基因组中具有固定的坐标值，保证每个样本呈现出种群基因组学中的人群多态性；同一位点在首次访问数据库后便写入缓存，此后查找均优先在缓存中查找，未能在缓存中找到时再访问千人基因组数据库，用户在使用时设置命令行参数的值为整体、东亚、美洲、非洲、欧洲或南亚，分别代表模拟对应的人群信息，包括等位基因的频率、东亚人群等位基因的频率、美洲人群等位基因的频率、非洲人群等位基因的频率、欧洲人群等位基因的频率、南亚人群等位基因的频率。

扩增错误模拟由用户在命令行参数中指定扩增错误率实现；

测序错误模拟由用户在命令行参数中指定测序错误率实现；

单端测序和双端测序模拟由用户在命令行通过对测序方式进行设置来实现单端或双端测序数据的输出，设置测序方式值为1表示双端测序，设置测序方式值为0表示双端测序，同时用户设置单端或双端测序时读段的读长；

接头模拟由用户在命令行参数指定两个接头的碱基序列，以及循环周期数实现；

索引序列模拟由用户在命令行参数中设置索引序列1和索引序列2的值为1表示进行索引序列模拟，设置为0时表示不加索引实现。

具体的，质量分布统计模块用于统计真实样本的质量值，生成参考质量相对分布，然后再使用内维尔插值算法来模拟不同读长时各个位点的碱基质量；

相邻位点深度分布训练模块采用最大期望算法，先根据已有的观测数据估计出相邻位点深度的概率模型的参数值；然后依据上一次估计出的参数值估计缺失数据的值，再根据估计出的缺失数据结合之前已有的观测数据重新对参数值进行估计；反复迭代，直至似然函数值在相邻两次循环中的变化小于用户指定的阈值，迭代结束；

整体深度分布统计模块用于统计一个真实样本整体深度的相对分布，用户输入一个测序样本的比对文件和与真实样本对应的目标捕获区域文件，然后算法将该样本整体深度分布的统计数据生成在户指定或者自动生成的文件夹中。

具体的，模板长度分布训练模块使用平面曲线离散点集拐点的快速查找算法和线性回归法构建数学模型；先通过平面曲线离散点集拐点快速查找算法找到曲线拐点，然后按照拐点将曲线分割，对每一段分割出来的曲线分别使用线性回归法进行拟合；在模板长度生成时，采用接受拒绝算法产生相应分布的随机数，实现了模板长度分布的模拟。

与现有技术相比，本发明至少具有以下有益效果：

本发明一种模拟人群背景信息的基因测序数据仿真方法，只需要设置参数并运行本系统就可以简便和快速的获得包含特定变异样本、符合数据标准的测序文件和比对文件。

进一步的，可以更加真实的模拟测序样本中邻近位点的测序深度的分布，且在从参考基因组读取序列时，本发明的设计方案可以显著提高算法的运行效率，节省用户获取样本的计算耗时。

进一步的，通过对人群背景信息进行模拟，使得变异检测软件在使用该方法产生的模拟样本进行测试时，样本范围能够体现不同人群的基因组特征，使得软件测试具有通用性，有助于变异检测软件在常见人群数据中都能够产出相对更加可靠的结果。

进一步的，通过对测序数据的碱基质量这一数据特征进行拟合使得模拟生成的测序数据在测序质量值的分布上与真实测序数据更为接近。

本发明还公开了一种模拟人群背景信息的基因测序数据仿真系统，不仅支持包含特定变异样本的生成与样本数据特征的拟合，还支持变异检测软件的验证，从而方便了变异检测软件的开发者快速获取大量样本并定位自己所开发软件的技术性能指标，进而指导其软件进行设计完善。

进一步的，实现了更丰富的变异模拟功能，能够模拟目前已知的所有常见变异类型和一些肿瘤基因组中发现的罕见变异类型。

进一步的，可以产生各种类型的变异在不同目标精度下的批量样本，进而在分析变异检测软件对批量样本的检测结果后反馈其性能指标。

进一步的，通过对测序数据的各项数据特征进行拟合，使得模拟数据的各项数据特征指标更加贴近真实测序数据。

进一步的，可以方便用户灵活拟合数据特征，用户根据自身需要，训练特定样本的深度分布模型，然后根据该模型进行样本模拟。

进一步的，可以方便用户灵活拟合数据特征，用户根据自身需要，训练特定样本的模板长度分布模型，然后根据该模型进行样本模拟。

综上所述，本发明可以支持单核苷酸位点变异、插入突变、删除突变和删除片段时插入片段变异、拷贝数变异、基因融合变异和基因融合时插入片段变异、串联重复变异以及倒位变异等目前已知的所有主要变异类型和一些罕见变异类型的模拟。该系统新增了人群背景信息模拟、扩增错误模拟、测序错误模拟、单双端测序模拟、接头模拟、单分子索引模拟、质量分布拟合、相邻位点深度分布拟合、整体深度分布拟合和模板长度分布拟合以及质量分布、相邻位点深度分布、整体深度分布和模板长度分布训练模块，使得该系统模拟生成的测序数据相比于其他任何已有软件更加贴近真实的测序数据。最后，该系统添加了特定目标精度下批量样本生成和变异检测软件功能验证的功能来批量生成变异样本并辅助变异检测软件变异检测软件在调参时快速验证其软件性能。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为本发明流程图；

图2为GSDcreator系统功能结构图；

图3为变异配置文件；

图4为正常样本信息记录文件；

图5为单核苷酸位点变异样本信息记录文件；

图6为插入片段样本信息记录文件；

图7为删除片段样本信息记录文件；

图8为拷贝数变异样本信息记录文件；

图9为基因融合变异信息记录文件；

图10为串联重复变异信息记录文件；

图11为倒位变异信息记录文件；

图12为干扰变异信息记录文件；

图13为人群多态模拟结果；

图14为单核苷酸位点变异模拟结果；

图15为插入片段模拟结果；

图16为删除片段模拟结果；

图17为删除片段时随机插入片段模拟结果；

图18为拷贝数变异模拟结果；

图19为基因融合变异模拟结果，其中，(a)为基因融合变异模拟中融合位点左侧的序列，(b)为基因融合变异模拟中融合位点右侧的序列；

图20为基因融合时随机插入片段变异模拟结果，其中，(a)为基因融合时随机插入片段变异模拟中融合位点左侧的序列，(b)为基因融合时随机插入片段变异模拟中融合位点右侧的序列；

图21为串联重复变异模拟结果；

图22为倒位变异模拟结果；

图23为真实样本GC含量图；

图24为模拟样本GC含量图；

图25为真实样本临近位点深度分布图；

图26为模拟样本临近位点深度分布图；

图27为真实样本质量分布图；

图28为模拟样本质量分图；

图29为真实样本模板长度分布图；

图30为模拟样本模板长度分布图。

具体实施方式

请参阅图1，本发明一种模拟人群背景信息的基因测序数据仿真方法，先进行参考基因组文件和目标捕获区域文件的加载；目标捕获区域文件中记录了参考基因组上用户所关注的每个目标区域的起始坐标和终止坐标；文件加载完成后本系统开始进入各个变异模拟模块，根据用户设置的测序深度、七种变异类型(分别是单核苷酸位点变异、插入突变、删除突变、拷贝数变异、倒位变异、基因融合变异、串联重复变异)、变异发生频率、变异在参考基因组上的坐标完成相应模拟；此时程序的打印台显示一个百分数，表示程序的运行进度；测序文件模拟完成，本系统自动调用比对软件BWA生成比对文件；本系统运行结束之后在一个以本系统运行时的日期和时刻命名的文件夹中生成相应测序文件和比对文件；为了加快程序的运行效率，本系统支持多进程运行，用户可以在命令行指定开启的进程数。

具体步骤如下：

S1、输入变异参数、目标捕获区域文件。其中变异参数包括：变异类型及其梯度信息、变异发生的染色体(包括1-22号常染色体和X、Y性染色体)、变异在参考基因组中的坐标、测序深度、读长、索引序列(命令行参数1表示添加索引序列、0表示不添加索引序列)、循环周期数、接头序列；变异类型包括：单核苷酸位点变异、插入突变、删除突变、删除出入叠加变异、拷贝数变异、融合变异、融合插入叠加变异、串联重复变异、倒位变异；梯度信息包括：单核苷酸位点变异、插入突变、删除突变、删除插入叠加变异、融合变异、融合插入叠加变异和倒位变异的频率梯度(一个大于0小于1的值)、拷贝数变异的拷贝倍数梯度(一个非负数值)、串联重复变异的重复数梯度(一个非负整数值)；测序深度是测序得到的碱基总量与基因组大小的比值，默认值是1000X，其中X是测序深度的度量单位；读长是读段的长度，默认值是75bp，其中bp是读长的度量单位；循环周期数是测序时一次通过的最长序列长度；加载参考基因组文件和目标捕获区域文件，创建输出文件夹，根据测序深度确定模板数量，测序深度是测序得到的碱基总量与基因组大小的比值，模板是参考基因组随机打断后生成的长度为几十至几百个碱基对不等的碱基序列片段，当用户没有指定目标区域则自动生成目标区域；

S2、通过线性回归法构建模型，使用接受拒绝算法产生相应分布模板长度；模板是参考基因组随机打断后生成的长度为几十至几百个碱基对不等的碱基序列片段；对当前已经遍历的模板数目进行判断，当程序已经遍历完所有模板时，测序文件生成完毕，而后本系统对测序文件进行读段比对，比对结束之后输出符合BAM格式标准的比对文件，结束；

S3、当尚未遍历当前所有模板时，根据相应混合高斯模型在目标区域生成断点，然后将生成的断点定位到参考基因组的对应位置，接着采用接受拒绝算法在该位置前后1000bp范围内生成拟合相应混合高斯模型的新断点，再使用建索引的方式实现从参考基因组的该位置处快速取出模板。建索引的方式是，先创建字典记录参考基因组每行序列所在相对位置，然后在读取时使用文件指针在查阅字典后直接跳到要读取的模板处，将模板取出，然后对模板进行人群背景拟合；

S4、当尚未遍历当前模板的所有碱基时，以每个碱基在参考基因组中的坐标作为随机种子，基于此使用梅森旋转算法生成随机数，而后比对千人基因组数据库相应位点；若该位置随机数小于数据库的人群发生率，则发生变异；若该位置随机数大于人群发生率，则不对该碱基做任何改变；当已经遍历完当前模板的所有碱基时，人群背景信息的模拟结束；

S5、当已经遍历当前模板所有碱基时，人群背景模拟结束并依次进行拷贝数变异模拟、单核苷酸位点变异模拟、基因融合模拟、串联重复模拟、倒位变异模拟、插入片段模拟和缺失片段模拟，然后采用内维尔算法生成测序碱基质量并生成测序序列，然后将测序序列写入测序文件，并重新返回步骤S2判断是否遍历所有模板。

请参阅图2，本发明一种模拟人群背景信息的基因测序数据仿真系统包括以下模块：单样本变异模拟模块、质量分布统计模块、相邻位点深度分布训练模块、整体深度分布统计模块、模板长度分布训练模块和特定目标精度下批量样本生成与变异检测软件验证模块，特定目标精度下批量样本生成与变异检测软件验证模块用于批量生成正常样本或变异样本，并对变异检测软件进行验证。

单样本变异模拟模块用于单核苷酸位点变异模拟、插入片段模拟、删除片段模拟、拷贝数变异模拟、基因融合变异模拟、串联重复变异模拟和倒位变异模拟。

单核苷酸位点变异模拟

用户在命令行参数中指定变异类型为单核苷酸位点变异，然后依次输入要模拟的单核苷酸位点变异在参考基因组中的坐标、发生突变之前的碱基和发生突变之后的碱基，以及变异发生率。

插入片段模拟

用户先在命令行参数中输入插入片段长度或输入插入片段的碱基序列，并指定变异类型为插入，然后依次输入变异在参考基因组中的坐标和变异的发生率即可。

删除片段模拟

用户在命令行参数中指定变异类型删除，然后依次输入需要删除片段在参考基因组中的坐标、变异发生率以及发生删除后是否随机插入片段，可以在命令行中指定插入片段的长度和插入片段碱基序列。

拷贝数变异模拟

用户在命令行参数中指定变异类型为拷贝数变异，然后依次输入要模拟的拷贝数变异在参考基因组中的坐标和拷贝数即可；当要模拟拷贝数增加，则输入一个大于1的整数或浮点数，当要模拟拷贝数减少，则输入一个小于1的浮点数；或者也可以输入原始比对文件的路径以及要模拟的拷贝数变异在参考基因组中的坐标和拷贝数，则可以实现直接对真实样本的比对文件进行修改来模拟拷贝数变异的功能。

基因融合变异模拟

基因融合变异模拟的设置需要先指定变异类型为融合，然后先依次输入发生融合的左边基因在参考基因组中的坐标和正负链向，再依次输入发生融合的右边基因在参考基因组中的坐标和正负链向；最后再输入基因融合变异发生率以及发生融合时是否在融合位置随机插入片段。同样，可以在命令行中指定插入片段的长度和插入片段碱基序列。

串联重复变异模拟

用户先在命令行参数中指定变异类型为串联重复；然后依次输入要模拟的串联重复变异在参考基因组中的坐标、重复次数以及串联重复变异发生率。

倒位变异模拟

倒位变异的模拟只需用户在命令行参数中指定变异参数为倒位，再依次输入要模拟的倒位变异在参考基因组中的坐标以及倒位变异发生率即可。

同时模拟多种变异

如果需要同时模拟多种变异，除了采用依次在命令行中将每种变异列出的方式外。还可以将变异写在文件中，然后使用直接传入文件路径进行模拟。文件截图如图3所示。

人群多态模拟

人群多态模拟的实现主要以千人基因组计划的千人基因组数据库为参照，预置突变随机发生。以每个碱基在参考基因组中的坐标做为随机种子并使用梅森旋转算法产生高质量的随机数，每次模拟样本产生一个随机数，使每个样本具有惟一的随机性，再加上每个碱基在参考基因组中具有固定的坐标值，保证每个样本呈现出人群多态性。由于模拟人群多态时遍历每一个碱基位点都要使用基因组数据处理包pysam访问一次千人基因组数据库，影响了算法的执行效率，所以考虑到大多数参考基因组坐标都会重复访问，采用了类似缓存的机制，同一参考基因组坐标在第一次访问数据库后便写入缓存，下一次查找时就会优先去缓存中查找，找不到时再访问千人基因组数据库，大大加快了算法的运行效率。用户在使用时，可以设置人群多态参数的值为整体、东亚、美洲、非洲、欧洲或南亚分别代表模拟整体人群等位基因的频率、东亚人群等位基因的频率、美洲人群等位基因的频率、非洲人群等位基因的频率、欧洲人群等位基因的频率、南亚人群等位基因的频率。

扩增错误模拟

高通量测序中的模板都要进行扩增，在扩增过程中就避免不了发生扩增偏好性问题，有的片段扩增的比较多，而有的片段的则扩增的较少。因此，为了使模拟数据更加拟合真实测序数据，加入了扩增错误模拟功能。用户只需在命令行参数中指定扩增错误率即可。

测序错误模拟

在高通量测序中，测序错误是不可避免的。因此在测序数据分析流程中必须考虑对测序错误的处理。为了验证数据指控流程中对测序错误处理的有效性，在模拟数据中加入了随机的测序错误。用户可以在命令行参数中指定测序错误率。

单端测序和双端测序模拟

高通量测序主要有单端测序和双端测序两种测序方式，单端测序只将引物序列连接到模板的一端，而双端测序是在两端的接头上都加上测序引物结合位点，测得其两端的序列。用户可以在命令中设置来单端和双端测序。设置测序方式值为1表示双端测序，设置测序方式值为0表示双端测序。同时可以设置单端或双端测序时候的读长。

接头模拟

在高通量测序中，文库构建过程需要将接头连接到待测片段的两端。在测序数据分析中，如果待测片段的长度小于测序周期的时候，就会将接头序列包含在测序结果中。所以数据分析流程的第一步就是要去除接头序列。为了验证数据分析时是否能准确去除接头，设计了接头模拟功能。用户可以使用命令行参数指定接头1和接头2序列以及循环周期数。

索引条码模拟

在高通量测序中，一个泳道可以同时测多个样品，为了避免混淆，每种样品的序列由一种索引修饰，这样实际得到的测序数据中的读段数据都是由索引标记的。用户可以在命令行参数中设置索引序列1和索引序列2的值为1表示进行索引条码模拟，设置为0时表示不加索引。

特定目标精度下批量样本生成

特定目标精度下批量样本生成模块用来为变异检测软件分梯度(梯度是指不同的目标精度，即不同的频率或不同的拷贝倍数)提供批量测试样本，产出阴性标准集(指不含目标变异的样本集)和阳性标准集(指含有目标变异的样本集)。用户在使用时可以指定想要模拟的变异类型、目标精度、每个样本所包含的变异数、生成的样本数以及目标区域和输出路径。如果没有指定目标区域，则默认在全基因组范围内产生随机变异。如果没有指定变异类型，则按照默认样本数生成一批正常样本。如果指定了变异类型但没有指定频率梯度，则会在随机频率梯度下，按照默认样本数和默认单个样本包含变异数生成一批样本，并输出在一个以变异类型和梯度值命名的文件夹中。在输出文件夹中，除包含生成的模拟样本外，还有一个记录该批次样本信息的文件。和一个记录干扰变异的文件。记录样本信息的文件中最后一列用来记录变异检测软件检测结果，检测到就填入1，没有检测到就填入0。图4、图5、图6、图7、图8、图9、图10、图11分别展示了正常样本、单核苷酸位点变异样本、插入片段样本、删除片段样本、拷贝数变异样本、基因融合变异样本、串联重复变异样本、倒位变异样本中记录样本信息的文件截图。图12展示了一个记录干扰变异的文件截图。

变异检测软件验证

将变异检测软件的检测结果填入记录样本信息的文件之后，输入文件路径并调用变异检测软件验证功能，会显示出该变异检测软件在这一目标精度下的敏感度和特异度。

质量分布统计模块

该功能模块通过统计真实样本的质量值，生成测序质量相对分布，然后再使用内维尔插值算法来模拟不同读长时各个位点的碱基质量。用户使用该功能模块可以生成任意真实样本中读段数据的质量值相对分布，只需要输入真实样本的质量值数据文件即可。关于质量值相对分布的统计数据会生成在用户指定或者自动生成的文件夹中。用户在调用样本模拟功能模块时，只需传入该文件夹的路径，就会模拟生成符合该质量分布的测序数据。同时，支持用户使用参数正常碱基质量最大值和正常碱基质量最小值设置正常碱基质量的最大值和最小值范围，使用参数错误碱基质量最大值和错误碱基质量最小值设置测序错误碱基质量的最大值和最小值范围。

相邻位点深度分布训练模块

通过实验观察发现，实际测序数据中相邻位点的深度分布呈混合高斯分布。在该训练模块中采用混合高斯模型的最大期望算法(expectation maximization algorithm,EM)，其基本思想是：先根据已有的观测数据，估计出模型参数的值；然后再依据上一步估计出的参数值估计缺失数据的值，再根据估计出的缺失数据加上之前已经观测到的数据重新再对参数值进行估计，然后反复迭代，直至对数似然函数值不再有明显的变化，迭代结束。用户在使用时只需要输入真实样本的比对文件以及要训练的区域在参考基因组中的坐标即可。训练完成之后，会将训练结果生成在用户指定或者自动生成的文件夹中。用户在调用样本模拟功能模块时，只需传入该文件夹的路径，就会模拟生成符合该相邻位点深度分布的测序数据。

整体深度分布统计模块

该功能模块用来统计一个真实样本整体深度的相对分布，需要用户输出一个真实样本的比对文件以及与真实样本相对应的目标捕获区域文件，运行本系统之后，会将该样本整体深度分布的统计数据生成在户指定或者自动生成的文件夹中。用户在调用样本模拟功能模块时，只需传入该文件夹的路径，就会模拟生成符合该整体深度分布的测序数据。同时支持用户使用参数测序深度对模拟样本的平均深度进行设置。

模板长度分布训练模块

该功能模块使用平面曲线离散点集拐点的快速查找算法和机器学习的线性回归法构建模型。为了提高曲线的拟合效果，采用了将平面曲线离散点集拐点查找算法和线性回归法相结合的方法。先通过平面曲线离散点集拐点快速查找算法找到曲线拐点，然后按照拐点将曲线分割，对每一段分割出来的曲线分别使用线性回归法进行拟合，从而达到了很好的拟合效果。在模板长度生成时，采用接受拒绝算法产生相应分布的随机数，进而实现了模板长度分布的模拟。用户在使用时，只需输入记录样本模板长度的文件即可，运行本系统后，会将训练好的模型生成在用户指定或者自动生成的文件夹中。用户在调用样本模拟功能模块时，只需传入该文件夹的路径，就会模拟生成符合该模板长度分布的测序数据。同时，支持用户使用模板长度最大值和模板长度最小值参数设置模板长度的最大值和最小值范围。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参阅图13，该图为人群背景模拟结果使用可视化工具IGV进行展示的图。图中，方框所标注位置的第一行所显示碱基为参考基因组中该位点处记录的碱基(如第二个方框中的碱基C)，下面的碱基为按照千人数据库替换为特定人群该位置的碱基(如第二个方框中的碱基T)。说明成功模拟了特定人群中该位置处出现的替换碱基。

请参阅图14，该图为单核苷酸位点变异的模拟结果使用可视化工具IGV进行展示的图。图中，显示出了与参考基因组不匹配的碱基C，说明该位点按照预设模拟发生了单核苷酸位点变异。

请参阅图15，该图为插入突变的模拟结果使用可视化工具IGV进行展示的图。从图中可以看出，许多读段数据出现了与参考基因组匹配不一致的情况。说明此处按照预设模拟发生了插入突变。

请参阅图16，该图为删除突变的模拟结果使用可视化工具IGV进行展示的图。从图中可以看出，许多读段数据发生了一长串的空缺。说明此处按照预设模拟发生了缺失突变。

请参阅图17，该图为删除插入叠加变异的模拟结果使用可视化工具IGV进行展示的图。从图中可以看出，在和图16所示变异的相同位置处重新插入了碱基序列，并出现了匹配不一致的情况。说明此处按照预设模拟发生了删除插入叠加变异。

请参阅图18，该图为拷贝数变异的模拟结果使用可视化工具IGV进行展示的图。从图中可以看出，在正常样本中本来应该呈现正态分布的曲线，在图中出现了明显的凹陷。说明此处按照预设模拟发生了拷贝数目减少的拷贝数变异。

请参阅图19，这两幅图为融合变异的模拟结果使用可视化工具IGV进行展示的图。其中，19(a)图为融合位点左边序列匹配到参考基因上的图示，19(b)图为融合位点右边序列匹配到参考基因上的图示。所以，融合位点左右的序列分别匹配到参考基因组上两处相距较远的位置，说明此处按照预设模拟发生了融合变异。

请参阅图20，该图为融合插入叠加变异的模拟结果使用可视化工具IGV进行展示的图。其中，20(a)图为融合位点左边序列匹配到参考基因上的图示，20(b)图为融合位点右边序列匹配到参考基因上的图示。所以，融合位点左右的序列分别匹配到参考基因组上两处相距较远的位置，并且还出现了明显的与参考基因组匹配不一致的情况。说明此处按照预设发生了融合插入叠加变异。

请参阅图21，该图为串联重复变异的模拟结果使用可视化工具IGV进行展示的图。从图中可以看出，曲线中间有明显凸起，深度显著加深，表示此处有较多重复的碱基序列。说明按照预设发生了串联重复变异。

请参阅图22，该图为倒位变异的模拟结果使用可视化工具IGV进行展示的图。从图中可以看出，双末端读对的链向信息两正同时出现，而正常的双末端读对的链向信息为一正一负。所以，双末端读对出现了不一致比对，链向信息异常。说明按照预设发生了倒位变异。

请参阅图23和24，图23和图24分别为一个真实样本的GC含量图和拟合该真实样本得到的模拟样本的GC含量图。可以看到模拟样本较好的拟合出了真实样本的GC含量分布。

请参阅图25和26，图25和图26分别为一个真实样本的邻近位点深度分布图和拟合该真实样本得到的模拟样本的邻近位点深度分布图。可以看到两个分布图非常相似。

请参阅图27和28，图27和图28分别为一个真实样本的质量分布图和拟合该真实样本得到的模拟样本的质量分布图。可以看到真实样本和对应模拟样本的质量分布基本一致。

请参阅图29和30，图29和图30分别为一个真实样本的模板长度分布图和拟合该真实样本得到的模拟样本的模板长度分布图。可以看到在模板长度分布的模拟中，模拟样本也几乎与真实样本是一致的。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。

Claims (10)

1.一种模拟人群背景信息的基因测序数据仿真方法，其特征在于，输入目标捕获区域文件、测序深度、单核苷酸位点变异、插入突变、删除突变、拷贝数变异、倒位变异、基因融合变异、串联重复变异、变异发生频率和变异在参考基因组上的坐标，加载参考基因组和目标捕获区域文件；目标捕获区域文件中记录参考基因组上用户关注的每个目标区域的起始坐标和终止坐标；根据测序深度确定模板数量；通过线性回归法构建模型，使用接受拒绝算法产生相应模板长度的概率分布，遍历每一条模板并对当前已经遍历完的模板数进行判断，当算法尚未遍历所有模板时，先根据模板长度的概率分布确定当前模板的长度，然后从参考基因组取出该长度模板，分别对取出的长度模板进行拷贝数变异模拟、单核苷酸位点变异模拟、基因融合模拟、串联重复模拟、倒位变异模拟、插入片段模拟和缺失片段模拟，最后生成读段，并将读段写入测序文件；当已经遍历完所有模板时，测序文件生成完毕；最后，对测序文件进行读段比对并生成比对文件，比对结束之后输出模拟的测序文件及其比对文件，仿真结束。

2.根据权利要求1所述的模拟人群背景信息的基因测序数据仿真方法，其特征在于，当尚未遍历所有模板时，首先在目标区域生成随机断点，然后将生成的断点定位到参考基因组的对应位置，采用接受拒绝算法在该位置前后1000碱基对范围内生成拟合相应混合高斯模型的新断点，再使用建索引的方式实现从参考基因组的该位置处快速取出模板，建索引的方式为：先创建字典记录参考基因组每行序列所在相对位置，然后在读取时使用文件指针在查阅字典后直接跳到要读取的模板处。

3.根据权利要求1所述的模拟人群背景信息的基因测序数据仿真方法，其特征在于，当尚未遍历当前模板的所有碱基时，以每个碱基在参考基因组中的坐标作为随机种子，使用梅森旋转算法生成随机数，比对千人基因组数据库相应位点，当该位置随机数小于数据库的人群发生率，则发生变异；当该位置随机数大于人群发生率，则不对该碱基做任何改变；当已经遍历完当前模板的所有碱基时，人群背景信息的模拟结束。

4.根据权利要求1所述的模拟人群背景信息的基因测序数据仿真方法，其特征在于，已经遍历当前模板的所有碱基时，结束人群背景信息模拟并依次进行拷贝数变异模拟、单核苷酸位点变异模拟基因融合模拟、串联重复模拟、倒位变异模拟、插入片段模拟和缺失片段模拟，然后采用内维尔算法生成测序碱基质量并写入对应读段，最后将该读段写入以fastq标准格式定义的输出文件，重新返回判断是否遍历所有模板。

5.一种根据权利要求1至4中任一项所述仿真方法的模拟人群背景信息的基因测序数据仿真系统，其特征在于，包括单样本变异模拟模块、质量分布统计模块、相邻位点深度分布训练模块、整体深度分布统计模块、模板长度分布训练模块和特定目标精度下批量样本生成与变异检测软件验证模块；单样本变异模拟模块用于单核苷酸位点变异模拟、插入片段模拟、删除片段模拟、拷贝数变异模拟、基因融合变异模拟、串联重复变异模拟和倒位变异模拟；特定目标精度下批量样本生成与变异检测软件验证模块用于批量生成正常样本和变异样本，并对变异检测软件进行验证；通过将变异写在文件中，然后输入文件路径并启动，实现同时模拟多种变异。

6.根据权利要求5所述的模拟人群背景信息的基因测序数据仿真系统，其特征在于，单核苷酸位点变异模拟由用户指定变异类型为单核苷酸位点变异，然后依次输入要模拟的单核苷酸位点变异所在参考基因组上的坐标、发生突变之前的碱基和发生突变之后的碱基，以及变异发生率；

插入片段模拟由用户在命令行参数中输入插入片段长度或插入片段的碱基序列，并指定变异类型为插入，然后依次输入变异在参考基因组上的坐标和变异的发生率；

删除片段模拟由用户在命令行参数中指定变异类型删除，然后依次输入需要删除片段在参考基因组上的坐标、变异发生率以及发生删除后是否随机插入片段，若是，在命令行参数中指定随机插入片段的长度和该片段的碱基序列；

拷贝数变异模拟由用户在命令行参数中指定变异类型为拷贝数变异，然后依次输入要模拟的拷贝数变异在参考基因组上的坐标和拷贝数；当要模拟拷贝数增加，则输入一个大于1的整数或浮点数，当模拟拷贝数减少，则输入一个小于1的浮点数；或者在命令行参数中输入原始比对文件的路径以及发生的拷贝数变异在参考基因组上的坐标和拷贝数，实现直接对真实样本的比对文件进行修改来模拟拷贝数变异的功能；

基因融合变异模拟先指定变异类型为融合，然后先依次输入发生融合的左边基因的参考基因组坐标和正负链向，再依次输入发生融合的右边基因的参考基因组坐标和正负链向；最后再输入基因融合变异发生率以及发生融合时是否在融合位置随机插入片段，在命令行参数中指定随机插入片段的长度和该片段的碱基序列；

串联重复变异模拟由用户先在命令行参数中指定变异类型为串联重复；然后依次输入要模拟的串联重复变异在参考基因组上的坐标、重复次数以及串联重复变异发生率；

倒位变异模拟由用户在命令行参数中指定变异参数为倒位，再依次输入要模拟的倒位变异在染色体上的坐标和倒位变异发生率。

7.根据权利要求5所述的模拟人群背景信息的基因测序数据仿真系统，其特征在于，特定目标精度下批量样本生成模块用于为变异检测软件分梯度提供批量测试样本，产出阴性标准集和阳性标准集；并生成样本信息记录文件；目标精度是拷贝数变异的拷贝倍数或者单核苷酸位点变异、插入突变、删除突变、基因融合变异、串联重复变异、倒位变异的频率值；梯度是目标精度梯度变化下呈现出的不同的频率或不同的拷贝倍数；然后，测试样本经变异检测软件检测后，将变异检测软件的检测结果填入样本信息记录文件，输入文件路径并调用变异检测软件验证功能，显示变异检测软件在对应目标精度下的敏感度和特异度。

8.根据权利要求5所述的模拟人群背景信息的基因测序数据仿真系统，其特征在于，仿真系统还包括人群多态模拟、扩增错误模拟、测序错误模拟、单端测序和双端测序模拟、接头模拟和单分子索引模拟；

人群多态模拟以千人基因组数据库为参照，预置突变随机发生实现，以每个碱基在参考基因组中的坐标做为随机种子并使用梅森旋转算法产生高质量的随机数，每次模拟样本产生一个随机数，使每个样本具有惟一的随机性，再加上每个碱基在参考基因组中具有固定的坐标值，保证每个样本呈现出种群基因组学中的人群多态性；同一位点在首次访问数据库后便写入缓存，此后查找均优先在缓存中查找，未能在缓存中找到时再访问千人基因组数据库，用户在使用时设置命令行参数的值为整体、东亚、美洲、非洲、欧洲或南亚，分别代表模拟对应的人群信息，包括等位基因的频率、东亚人群等位基因的频率、美洲人群等位基因的频率、非洲人群等位基因的频率、欧洲人群等位基因的频率、南亚人群等位基因的频率；

扩增错误模拟由用户在命令行参数中指定扩增错误率实现；

测序错误模拟由用户在命令行参数中指定测序错误率实现；

单端测序和双端测序模拟由用户在命令行通过对测序方式进行设置来实现单端或双端测序数据的输出，设置测序方式值为1表示双端测序，设置测序方式值为0表示双端测序，同时用户设置单端或双端测序时读段的读长；

接头模拟由用户在命令行参数指定两个接头的碱基序列，以及循环周期数实现；

索引序列模拟由用户在命令行参数中设置索引序列1和索引序列2的值为1表示进行索引序列模拟，设置为0时表示不加索引实现。

9.根据权利要求5所述的模拟人群背景信息的基因测序数据仿真系统，其特征在于，质量分布统计模块用于统计真实样本的质量值，生成参考质量相对分布，然后再使用内维尔插值算法来模拟不同读长时各个位点的碱基质量；

相邻位点深度分布训练模块采用最大期望算法，先根据已有的观测数据估计出相邻位点深度的概率模型的参数值；然后依据上一次估计出的参数值估计缺失数据的值，再根据估计出的缺失数据结合之前已有的观测数据重新对参数值进行估计；反复迭代，直至似然函数值在相邻两次循环中的变化小于用户指定的阈值，迭代结束；

整体深度分布统计模块用于统计一个真实样本整体深度的相对分布，用户输入一个测序样本的比对文件和与真实样本对应的目标捕获区域文件，然后算法将该样本整体深度分布的统计数据生成在户指定或者自动生成的文件夹中。

10.根据权利要求5所述的模拟人群背景信息的基因测序数据仿真系统，其特征在于，模板长度分布训练模块使用平面曲线离散点集拐点的快速查找算法和线性回归法构建数学模型；先通过平面曲线离散点集拐点快速查找算法找到曲线拐点，然后按照拐点将曲线分割，对每一段分割出来的曲线分别使用线性回归法进行拟合；在模板长度生成时，采用接受拒绝算法产生相应分布的随机数，实现了模板长度分布的模拟。