本发明涉及基因测序技术领域，特别涉及一种基于bionano平台检测长串联重复序列的方法。
背景技术：
长串联重复序列是指dna序列中的多个核苷酸(单个重复单元大于1kb)前后首尾相连而构成的重复序列，重复单元数目的变化会对基因组结构造成重要影响。
bionano光学图谱是单个dna分子有序的全基因组限制性内切酶酶切位点图谱。利用内切酶对dna进行识别、酶切并标记荧光，再通过纳米级毛细管电泳来把dna分子拉直，将每个dna分子线性展开，进行超长单分子高分辨率荧光成像，生成酶切位点分布图。利用这些极长读长片段进行基因组比对克服了传统的使用小于重复单位的读长片段处理基因组重复区域的不可靠性。
saphyr是bionano第二代单分子基因组结构分析平台。对检测和分析基因组结构变异具有异常的敏感性和特异性，能够揭示很多基因组真实结构。saphyr将高速、高通量以及对结构变异出色敏感性相结合，使之成为人类和转化研究应用的理想解决方案。在许多领域皆可使用高分辨率的物理基因组图谱了解基因组结构，包括未确诊的遗传性疾病诊断、基因发现与治疗进展、癌症、细胞系的研究、选择育种、进化生物学、参考基因组组装。saphyr融合了专有的纳米通道和光学基因组图谱，使极长高分子量的dna在其原始状态下成像。这一技术对结构变化异常敏感，基因组组装接近于单独使用短读序列测序的100倍，并准确地纠正了基于序列的组装错误。解决由下一代测序系统(ngs)遗漏的大片段结构变异，大片段结构变异与多种疾病和病症的联系密不可分。
现有检测长串联重复序列的方法是基于bionano光学图谱技术的组装算法(参考https://bionanogenomics.com/wp-content/uploads/2014/10/bionano-poster-ashg2014-chan-long-repeats-cnv.pdf)，将reads*(读长)组装成contig，再比对到参考基因组序列，对串联重复序列进行肉眼计数。
基于组装算法检测长串联重复序列有以下几点缺陷：1、bionano数据中存在大量插入、缺失错误，易导致组装错误；2、计算时间长；3、消耗大量计算资源；4、不能准确检出嵌合体样本。
技术实现要素：
为了解决现有技术中存在的问题，本发明提供一种基于bionano平台检测长串联重复序列的方法。本发明的方法通过构建机器学习模型，对bionano数据进行过滤，去除插入、缺失位点的假阳性错误，并且基于比对算法实现长串联重复单元计数，减少运行时间、计算资源的消耗。所述方法还可以用于确定样本的基因类型，结合聚类分析算法和每条reads上的重复单元数目，确定样本基因型(纯合、杂合)或者嵌合体。
本发明提供的一种基于bionano平台检测长串联重复序列的方法，包括如下步骤：
(1)提取样本dna，采用内切酶对dna进行酶切、标记、修复、染色，用bionanosaphyr系统定量处理；
(2)构建基于内切酶酶切位点的bionano参考基因组；
(3)对原始数据进行信噪比过滤；
(4)将过滤后的数据比对到步骤(2)的参考基因组；
(5)对比对后的数据进行质量评估，如果质量不合格，则终止分析，如果质量合格，则进行步骤(6)；
(6)构建朴素贝叶斯分类器机器学习模型，利用朴素贝叶斯分类器对样本中指定区域的reads假阳性位点进行过滤；
(7)矩阵构建，根据reads在参考基因组上的比对情况结合reads上刻痕的距离信息，构建距离矩阵m，所述刻痕为酶切位点；对缺失位点以0填充，超过10％的reads在相同位置发生插入的位点按照真实插入位点处理，矩阵增加一列数据；
(8)聚类分析，采用欧氏距离计算步骤(7)中构建矩阵的各reads之间的距离，采用平均距离计算组间距离，进行层次聚类分析；
(9)重复单元计数，根据候选reads中刻痕的位置信息，识别重复单元；
(10)确定样本基因型，根据步骤(8)聚类分析中reads的距离关系以及每条reads上的重复单元数目判断样本基因型。
根据本发明的实施方式，在上述方法中，步骤(6)中的构建朴素贝叶斯分类器机器学习模型的步骤为：
(a)构建数据集
采用hx1数据构建的中国人参考基因组及其bionano光学图谱数据，将bionano数据比对到hx1参考基因组上，比对到参考基因组上的位点定为真阳性位点，未比对到的位点定为假阳性位点；分别随机选择1000个真阳性位点和1000个假阳性位点作为数据集；
(b)特征选择
针对bionano的数据特点，根据reads比对到参考基因组的置信度对位点的强度、信噪比、覆盖度进行加权；同时，结合位点上下游数据得到用于描述该位点的分类特征；
(c)构建模型
基于朴素贝叶斯分类器公式
并通过公式
确定分类结果，其中y表示分类标签，y为0表示假阳性位点，y为1表示真阳性位点，x1至xn表示步骤(b)中所述分类特征的值，n表示所述分类特征的编号。
根据本发明的实施方式，在上述构建朴素贝叶斯分类器机器学习模型的步骤中，步骤(b)对位点的强度、信噪比、覆盖度进行加权的公式为：
其中，n表示reads数目；c表示reads的置信度；w表示权重；d、d表示加权前后位点强度，r、r表示加权前后位点信噪比、v、v表示加权前后位点覆盖率。
根据本发明的实施方式，在上述构建朴素贝叶斯分类器机器学习模型的步骤中，步骤(b)中所述分类特征为：位点加权强度、位点加权信噪比、位点加权覆盖度、当前位点上游比对到参考基因组的位点数目、当前位点上游比对到参考基因组的位点的平均加权强度、当前位点上游比对到参考基因组的位点的平均加权信噪比、当前位点下游比对到参考基因组的位点数目、当前位点下游比对到参考基因组的位点的平均加权强度、以及当前位点下游比对到参考基因组的位点的平均加权信噪比。
根据本发明的实施方式，在上述方法中，步骤(10)中所述重复单元计数具体步骤为：
首先计算重复区域的位置间距di，di＝li+1-li，li表示reads上刻痕的位置坐标；
设定阈值α＝0.1，当相邻两个间距的距离差与两个间距的最小值的比值ti小于阈值α时，确定为一个真实重复，否则不进行计数，同时当连续3个比值ti均不小于α时，计数终止，ti的计算公式为：
重复单元数目公式为：
n表示该reads上的重复单元数目，n表示最终的重复单元数目。
根据本发明的实施方式，在上述方法中，步骤(10)确定样本基因型的方法是根据步骤(8)聚类分析中reads的距离关系，由近及远，排除包含reads数目小于总reads数目5％的类，随后做如下判断：
(i)纯合：reads聚为1类，该类别reads数目占总reads数目的80％～100％，且重复单元的数目仅有1类，则样本基因型为纯合；
(ii)杂合：reads聚为2类，每个类别reads数目占总reads数目的40％～60％，且重复单元的数目有2类，则样本基因型为杂合；
(iii)嵌合体：reads聚为3类，每个类别reads数目占总reads数目的20％～40％，且重复单元的数目有3类，则样本为嵌合体。
根据本发明的实施方式，在上述方法中，步骤(1)所述样本为人血液白细胞样本。
根据本发明的实施方式，在上述方法中，步骤(2)所述bionano参考基因组为bionanohg38参考基因组。
根据本发明的实施方式，在上述方法中，步骤(3)中的信噪比过滤根据基于histogram的过滤算法进行。
根据本发明的实施方式，在上述方法中，步骤(4)中所述比对所用的软件为bionanorefaligner。
根据本发明的实施方式，在上述方法中，步骤(9)中所述重复单元为人类4q35区域d4z4重复单元。
在本发明中，所述内切酶可以是但不限于bsssi酶。
本发明技术方案带来的有益效果：
1、构建机器学习模型的方法，根据具体分析要求选择分类特征，评估多种机器学习算法分类性能，采用最佳方案，可以最大限度的去除reads中可能存在的假阳性位点，提高后续处理的准确性，使得对于长串联重复序列的重复单元计数，更为准确。
2、通过聚类算法，相同基因型的reads可以很容易的聚在一起，可以很容易的区分基因型(纯合、杂合)以及嵌合体。
3、提出位点距离差算法，设置特定阈值，排除数据波动造成的影响，准确计算重复单元数目。
4、相比于bionano组装算法，运行速度更快、消耗内存资源更少。
附图说明
图1为基于bionano平台检测长串联重复序列的方法流程图。
图2为不同机器学习模型受试者工作特征曲线图。
图3为本发明方法与bionano组装算法、sb方法在基因型-重复单元数目检测中的比较图。
图4为本发明方法与bionano组装算法在运行时间和运行内存方面的比较图。
具体实施方式
以下通过实施例对本发明的方案进行详细说明。本领域的技术人员应该明白，下面的实施例仅用于解释说明本发明，而不是限定本发明的范围。
检测58例人类样本4q35区域d4z4重复单元数目(d4z4是一段长3.3kb左右的重复单元，4q35区域存在多个d4z4重复，不同的个体重复数目不同(参见lemmers,r.j.,dekievit,p.,sandkuijl,l.,padberg,g.w.,vanommen,g.j.b.,frants,r.r.,&vandermaarel,s.m.(2002).facioscapulohumeralmusculardystrophyisuniquelyassociatedwithoneofthetwovariantsofthe4qsubtelomere.naturegenetics,32(2),235-237.)，该实施例流程图如图1所示。
实施例一、构建机器学习模型
a.数据集
采用hx1数据(参见shi,l.,guo,y.,dong,c.,huddleston,j.,yang,h.,han,x.&lintner,k.e.(2016).long-readsequencinganddenovoassemblyofachinesegenome.naturecommunication.)构建的中国人参考基因组及其bionano光学图谱数据，将bionano数据比对到hx1参考基因组上，可以比对到参考基因组上的位点定为真阳性位点，未比对到的位点定为假阳性位点。据此，我们根据对两种情况分别随机选择1000个真阳性位点和1000个假阳性位点作为数据集。
b.特征选择
针对bionano的数据特点，根据reads比对到参考基因组的置信度(confidence)对位点的强度(intensity)、信噪比(snr)、覆盖度进行加权(公式1-4所示)。同时，结合位点上下游数据得到了用于描述该位点的9个分类特征(表1)。
注：n表示reads数目；c表示reads的置信度；w表示权重；d、d表示加权前后位点强度，r、r表示加权前后位点信噪比、v、v表示加权前后位点覆盖率(reads上有位点记为1，没有记为0)。
表1.构建模型所使用的分类特征
c.模型的构建
朴素贝叶斯分类器(nb)是一种基于贝叶斯理论的有监督的机器学习模型，在生物医学领域有着广泛的应用。本发明结合bionano数据特点，利用朴素贝叶斯分类器对假阳性位点进行过滤。
朴素贝叶斯分类器满足公式5的条件，y表示分类标签(y∈{0,1}，0表示假阳性位点，1表示真阳性位点)，xi表示表1.中提到的分类特征的值。
基于条件独立性假设，公式5可以变换为公式6。
又因为对于指定的输入p(x1,…,xn)为常量，所以公式7成立。
假定连续型特征变量p(xi|y)符合正态分布(公式8)。
随后，利用训练集对σy,i、μy,i进行估计。
对于测试集，通过判断两种分类情况的概率值(公式9)，确定分类结果。
d.模型训练和评估
采用十折交叉验证，将数据集分成训练集和测试集进行训练和评估。评估指标选用准确性(accuracy，公式10)、敏感性(sensitivity，公式11)特异性(specificity，公式12)等指标进行衡量。
注：tp表示真阳性位点的个数，tn表示真阴性位点的个数，fp表示假阳性位点的数目，fn表示假阴性位点的个数
评估后发现，针对bionano假阳性位点的评估，朴素贝叶斯分类器准确性为0.977，敏感性为0.976，特异性为0.978，达到较好的评估效能。
e.与其他分类器的比较
为了更好的说明朴素贝叶斯分类器对bionano假阳性位点的识别能力，同时对其他常用分类器，随机森林(rf)、决策树(dt)、支持向量机(svm)、k近邻(knn)、逻辑回归(lr)、人工神经网络(ann)等进行了评估，评估结果如图2、表2所示。
综上所述，本专利针对bionano平台的测序特点，筛选合适的分类特征，构建朴素贝叶斯分类算法，对假阳性位点有着很强的识别能力。识别、过滤假阳性位点，对后续长串联重复序列的基因型判断和重复单元计数有着重要意义。
表2.不同机器学习模型的准确性、敏感度、特异度以及曲线下面积
实施例二、长串联重复序列检测
1、实验方法
人血液红细胞裂解处理(1hour)
白细胞定量计数(5min)
对白细胞进行包埋处理(～1hour)
采用蛋白酶k进行消化
洗胶来固定dna
dna回收
dna透析和均质化
dna浓度定量(10μl,2小时30分钟)
采用bsssi酶对dna进行酶切(10μl,2小时30分钟)
标记(15μl,1小时15分钟)
修复(20μl,45分钟)
染色处理(60μl,16小时/过夜)
bionanosaphyr系统进行定量处理
实验细节参考：
https://bionanogenomics.com/wp-content/uploads/2017/03/30033-rev-c-bionano-prep-blood-dna-isolation-protocol.pdf；
https://bionanogenomics.com/wp-content/uploads/2017/07/30024-rev-j-bionano-prep-labeling-nlrs-protocol.pdf。
2、构建基于bsssi酶切位点的bionanohg38参考基因组。根据bsssi酶的特异性识别位点cacgag，对hg38参考基因组的fasta文件进行处理，得到bionano指定的cmap基因组文件，基因组信息表3所示。
表3.hg38基因组酶切位点统计
3、原始数据信噪比(snr)过滤，根据基于histogram(histogram-based)的过滤算法(参见：pedregosa,f.,varoquaux,g.,gramfort,a.,michel,v.,thirion,b.,grisel,o.,...&vanderplas,j.(2011).scikit-learn:machinelearninginpython.journalofmachinelearningresearch,12(oct),2825-2830.；wang,j.h.,liu,w.j.,&lin,l.d.(2002).histogram-basedfuzzyfilterforimagerestoration.ieeetransactionsonsystems,man,andcybernetics,partb(cybernetics),32(2),230-238)对原始数据进行处理，获取高质量的数据进行后续操作。
4、比对到参考基因组，采用bionano官方推荐的refaligner(version:6700.6902)软件将过滤后的数据比对到步骤2构建的参考基因组。
5、质量评估，根据bionano官方给出的质控标准(https://bionanogenomics.com/wp-content/uploads/2017/05/30175-rev-a-bionano-molecule-quality-report-guidelines.pdf、https://bionanogenomics.com/wp-content/uploads/2017/03/30110-rev-b-bionano-solve-theory-of-operation-structural-variant-calling.pdf)结合项目经验制定过滤标准(如表4所示)，58例样本均通过质控评估，可以进行后续分析。
表4.质量评估标准
6、根据实施例一构建的机器学习模型，对样本中指定区域的reads进行假阳性位点进行过滤。
7、矩阵构建，根据reads在参考基因组上的比对情况结合reads上label(刻痕)的距离信息，构建距离矩阵m。对缺失位点以0填充，超过10％的reads在相同位置发生插入的位点按照真实插入位点处理，即矩阵增加一列数据。
8、聚类分析，采用欧氏距离计算步骤7中构建矩阵的各reads之间的距离，采用平均距离计算组间距离，进行层次聚类分析。
9、d4z4重复单元计数，根据候选reads中刻痕的位置信息，识别d4z4重复单元。首先计算重复区域的位置间距di(公式13)，每个重复单元的间距di理论上应该是一致的，但是因为部分测序错误引入了误差，造成数据波动，为了消除误差影响，我们设定了阈值α＝0.1，当相邻两个间距的距离差与两个间距的最小值的比值(公式14)小于阈值α时(公式15)，确定为一个真实重复，否则不进行计数。同时当连续3个比值(ti)均不小于α时，计数终止。
di＝li+1-li公式13
注：l表示reads上label的位置坐标，n表示该reads上的重复单元数目，n表示最终的重复单元数目。
10、确定样本基因型，根据步骤8聚类分析中reads的距离关系，由近及远，排除包含reads数目小于总reads数目5％的类，随后做如下判断。
纯合：reads聚为1类(该类别reads数目占总reads数目的80％～100％)，且重复单元的数目仅有1类；
杂合：reads聚为2类(每个类别reads数目占总reads数目的40％～60％)，且重复单元的数目有2类；
嵌合体：reads聚为3类(每个类别reads数目占总reads数目的20％～40％)，且重复单元的数目有3类。
11、结果评估
采用上述方法对58例样本进行处理，获取其基因型与重复单元数目，同时与bionano组装算法和southern印迹杂交方法(southernblot，sb，该方法仅适用于d4z4重复单元数目的检测，不适用于其他长串联重复序列的检测)进行比较。结果发现(见表5、图3)，58例样本中，51例(87.93％)样本3种方法得出的结果一致。剩余的7例样本中，2例(3.45％)样本(s042、s057)southernblot未能找到合适的ecori/blni酶切位点，未得到合理结果实验失败，而本发明方法与bionano组装算法结果一致，具有很高的可信性；5例(8.62％)样本本发明方法与sb方法均检出为嵌合体，而bionano组装算法检出为杂合，怀疑是由于bionano组装算法固有缺陷(组装人类基因组只能组装二倍体)引起的错误。
在时间和内存的消耗上，我们也对本发明方法与bionano组装算法进行了比较(wilcoxon符号秩检验)，结果发现(表3，图4)，采用相同集群配置，与bionano组装算法比较，在运行时间上，本发明方法平均运行核时为5,786.67s，而bionano组装算法平均运行核时为52,911.58s,二者差异极显著(p-value<0.01)。在内存消耗上，本发明方法平均使用内存为248.09m，而bionano组装算法平均使用内存为1,005.95m，二者差异极显著(p-value<0.01)。
综上，与southernblot方法比较，本发明方法操作简单，适用广泛，不会因为特殊酶切位点的缺失而影响实验结果。与bionano组装算法相比较，本发明方法对嵌合体样本的处理准确性更高，同时运算速度更快、消耗资源更少。
表5.本发明方法与bionano组装算法、southernblot方法进行比较