一种基于二代测序数据的hla分型方法
技术领域
本发明涉及一种基于二代测序数据的HLA分型方法，应用于HLA二代测序数据分型。
背景技术
基于二代测序数据的HLA分型由于其通量高、分型速度快等特点，日后必将逐渐替代一代测序成为HLA分型的主流测序手段。研发人员开发了多种基于二代测序技术的分型方法，比如：中国发明专利CN103221551A公开了HLA基因型别-SNP连锁数据库、其构建方法、以及HLA分型方法，是一种根据不同型别的SNP连锁关系来进行HLA分型的方法。中国发明专利申请CN109477143A公开了一种人类白细胞抗原分型方法，是一种使用已知HLA等位基因参考序列的多重序列比对(MSA)来确定一个或多个额外的HLA 等位基因参考序列的方法，最终选择出与所述个体匹配所接近的参考序列，从而完成HLA分型。
虽然现有技术存在多种基于二代测序数据的HLA分型方法，但仍存在以下几方面的问题：1.很多方案中不支持高分辨率的HLA分型。2.HLA基因多态性较高，重复序列多，比对难度大，出现错误比对的情况较多，难以保证分型结果的准确性。3、采用设定阈值方法判断是否为杂合碱基，方法不够灵活，在复杂情况下容易出错。
发明内容
针对上述现有技术，为了解决基于二代测序技术的HLA分型准确性的问题，本发明提供了一种基于二代测序数据的HLA分型方法。
为了解决所述技术问题，本发明采样的技术方案是：一种基于二代测序数据的HLA分型方法，包括以下步骤：
S01）、构建HLA参考基因序列数据库，具体为：
S11)、下载HLA型别数据库，将所有HLA型别对应的序列进行处理，标准化为长度一致的序列；
S12）、生成每一个型别序列的坐标映射向量，坐标映射向量是一条序列的每个碱基在数据库中的标准坐标位置的映射数值向量；
S13）、构建HLA参考基因序列数据库，在数据库中添加型别序列，每次增加一条型别序列，并确保该型别序列与数据库中的其他序列的差异值均大于预设阈值T1；
S02）、进行HLA二代测序数据的分析，具体为：
S21）、对原始HLA二代测序数据进行质控与过滤，去除低质量数据，所述低质量数据是指测序质量不合格和长度小于最低长度的reads；
S22）、进行序列比对，将经过步骤S21的原始HLA二代测序数据比对至构建好的HLA参考基因序列数据库，得到初次比对结果；
S23）、进行二次比对，将初次比对后的结果进行统计，计算出数据库中每个型别序列的最佳匹配次数，按照从高到低排序，将匹配次数排在前N个的型别序列重新组建成新的HLA参考基因序列数据库，缩小比对范围，重新进行二次比对；
S24）、统计参考基因序列每一个坐标上的比对结果信息；
S03）、通过聚类算法判断杂合碱基位置，对基因上同一个外显子或内含子区域的比对结果进行聚类算法分类，分离出杂合碱基位置与纯合碱基位置，纯合碱基位置为某一坐标位置处的序列比对结果为单一碱基，杂合碱基位置为某一坐标位置处的序列比对结果为两个碱基同时出现；
S04）、获取单倍体序列，提取所有至少覆盖两个杂合碱基位置的原始read序列，进行局部序列重比对，获取单倍体序列结果；
S05）、进行序列分别，将获取的单倍体序列比对至HLA型别数据库中，获取正确的分型结果。
进一步的，通过聚类算法判断杂合碱基位置的过程为：
S31）、计算出每个位置的碱基频率，由高到低进行排序，将反映杂合情况的第二高的频率值作为聚类算法的选取对象，随机选取2个对象作为初始的聚类中心，然后计算每个对象与各个种子聚类中心之间的距离，把每个对象分配给距离它最近的聚类中心；
S32）、聚类中心以及分配给它们的对象就代表一个聚类，每分配一个样本，聚类的聚类中心根据聚类中现有的对象被重新计算；
S33）、重复步骤S32，直到没有对象被重新分配给不同的聚类；
S34）、最终确定的两个聚类中心的频率值如果相差小于阈值T2，则判断整个区域均为纯合碱基，如果两个聚类中心的频率值相差大于或者等于阈值T2，则该测序区域包含杂合碱基，分配给频率值高的聚类中心的位置判定为杂合碱基位置，分配给频率值低的聚类中心的位置判定为纯合碱基位置。
进一步的，步骤S22中，采用bowtie2算法进行序列比对， Bowtie2使用FM索引对基因组进行索引，将测序reads与HLA参考基因序列数据库中的长参考序列进行比对，经过序列比对之后的文件包括比对到参考基因组上的位置以及比对的插入缺失信息。
进一步的，获取单倍体序列的过程为：对相邻且距离较近的杂合位置进行reads重新比对，统计在同一条read上的突变组合个数，个数多的组合判定两个突变位于同一条染色体上，进而最终确定单倍型序列。
进一步的，步骤S11中，登录IMGT官网，下载最新HLA型别数据库。
进一步的，步骤S11中，对不同型别序列中的插入和缺失在正确位置用符号表示。
本发明的有益效果：本发明的HLA二代测序数据分型方法，可对复杂的HLA区域进行准确的高分辨率分型，且通过软件的形式，支持界面操作，能在普通的计算机上运行复杂的二代测序数据分析，降低了对计算资源的消耗，并且解决了基于二代测序技术的HLA分型准确率不高的难题，具有分型准确率高、易用性强的特点。本发明的HLA二代测序数据分型方法，已在申请人所在公司HLA部门投入使用，在数据质控较好的情况下可实现较高的HLA分型准确率，整体分型准确率在99%以上，大大提高了HLA数据分析人员的工作效率，解决了以往数据分析难、依赖国外收费分析工具的难题。
附图说明
图1为实施例1所述方法的流程图；
图2为IMGT/HLA数据库收录的序列信息示意图；
图3为原始HLA二代测序数据的示意图；
图4为原始HLA二代测序数据过滤示意图；
图5为经过序列比对之后的示意图；
图6为统计序列比对结果的示意图；
图7为纯合碱基与杂合碱基位置的示意图；
图8为局部重比对的原理示意图。
具体实施方式
下面结合实施例对本发明作进一步的说明。然而，本发明的范围并不限于下述实施例。本领域的专业人员能够理解，在不背离本发明的精神和范围的前提下，可以对本发明进行各种变化和修饰。
实施例1
本实施例公开一种基于二代测序数据的HLA分型方法，如图1所示，包括以下步骤：
S01）、构建HLA参考基因序列数据库，包括以下具体步骤：
S11）、登录IMGT官网（http://www.imgt.org/），下载最新HLA型别数据库，将所有HLA型别对应的序列进行处理，对不同型别序列中的插入和缺失在正确位置用符号表示，标准化为长度一致的序列；
其中HLA型别数据库来自于IMGT（the international ImMunoGeneTicsinformation system）/HLA数据库，该数据库收录了HLA不同Allel 的序列信息，数据信息如图2所示。
S12）、基于HLA型别数据库的序列坐标标记方法，生成每一个型别序列的坐标映射向量；
其中的坐标映射向量是一条序列的每个碱基在数据库中的标准坐标位置的映射数值向量，由于不同型别之间的序列差异性很大，所以记录每个型别的坐标映射向量可以更准确地判断发生在HLA基因上的碱基插入与缺失。
S13）、构建HLA参考基因序列数据库，在数据库中添加型别序列，每次增加一条型别序列，并确保该型别序列与数据库中的其他序列的差异值均大于阈值T1。
本方法所构建的HLA参考基因序列数据库能够准备的表达发生在HLA基因上的碱基插入与缺失，并且添加的型别序列差异值均大于阈值T1，即添加的型别序列之间差异不会很大，方便后续的分析。
S02）、进行HLA二代测序数据的分析；
S21）、对一批原始下机数据（原始HLA二代测序数据）进行质控与过滤，去除低质量数据；
原始下机数据为FastQ格式的文件，其内容如图3所示，FastQ格式的序列一般都包含有四行，第一行由'@'开始，后面跟着序列的描述信息，第二行是序列。第三行由'+'开始，后面也可以跟着序列的描述信息。第四行是第二行序列的质量评价。
所述低质量数据是指测序质量不合格和长度小于最低长度的reads。通过设置阈值去除测序质量值不合格的reads，去除长度小于最低长度的reads，数据过滤的情况如图4所示。
S22）、进行序列比对，将经过步骤S21的数据比对至步骤S01构建好的HLA参考基因序列数据库。
所述序列比对采用bowtie2算法，Bowtie2 是将测序reads与长参考序列比对的工具。适用于将长度大约为50到100或1000字符的reads与相对较长的基因组进行比对。Bowtie2使用FM索引（基于Burrows-Wheeler Transform 或 BWT）对基因组进行索引，以此来保持其占用较小内存。经过序列比对之后的文件如图5所示，其中包含比对到参考基因组上的位置，以及比对的插入缺失信息等。
S23）、进行二次比对，将初次比对后的结果进行统计，计算出数据库中每个型别的最佳匹配次数，按照从高到低排序，将匹配上次数较多的前N个型别重新组建成新的HLA参考基因序列数据库，缩小比对范围，重新进行二次比对。
S24）、统计参考基因序列每一个坐标上的比对结果信息，统计信息如图6所示。
S03）、采用聚类算法判断杂合碱基位置，对基因上同一个外显子或内含子区域的比对结果信息进行聚类算法分类，分离出杂合碱基位置与纯合碱基位置，分离结果如图7所示。
本步骤中，通过聚类算法判断杂合碱基位置的过程为：
S31）、计算出每个位置的碱基频率，由高到低进行排序，将反映杂合情况的第二高的频率值作为聚类算法的选取对象，随机选取2个对象作为初始的聚类中心，然后计算每个对象与各个种子聚类中心之间的距离，把每个对象分配给距离它最近的聚类中心；
S32）、聚类中心以及分配给它们的对象就代表一个聚类，每分配一个样本，聚类的聚类中心根据聚类中现有的对象被重新计算；
S33）、重复步骤S32，直到没有对象被重新分配给不同的聚类；
S34）、最终确定的两个聚类中心的频率值如果相差小于阈值T2，则判断整个区域均为纯合碱基，如果两个聚类中心的频率值相差大于或者等于阈值T2，则该测序区域包含杂合碱基，分配给频率值高的聚类中心的位置判定为杂合碱基位置，分配给频率值低的聚类中心的位置判定为纯合碱基位置。
所述纯合碱基位置为某一坐标位置处的序列比对结果为单一碱基，该碱基频率统计理论上为100%（实际上由于存在测序错误与比对误差等因素，频率会大概率在90%到100%之间波动）。杂合碱基位置为某一坐标位置处的序列比对结果为两个碱基同时出现，每个碱基频率统计理论上为50%（实际上由于存在测序错误与比对误差，以及两条DNA序列测序产生的不平衡等因素，两个碱基的统计频率会大概率在10%到90%之间波动）。
S04）、采用单倍体序列获取算法，提取所有至少覆盖两个杂合碱基位置的原始read序列，进行局部序列重比对，获取单倍体序列结果，局部重比对原理如图8所示。
其中“单倍体序列获取算法”，对距离较近的杂合位置进行reads重新比对，统计在同一条read上的突变组合个数，个数多的组合可以判定两个突变位于同一条染色体上，进而最终确定单倍型序列。
本实施例采用二代测序数据，使用的是双端测序，即对一个Insert片段从两端分别进行测序，进而生成一对测序reads序列，如果两个杂合位置能在一对测序reads的范围之内，就可以统计突变组合的支持reads数，这种就算距离较近。
S05）、进行序列分型，将获取的单倍体序列比对至HLA型别数据库中，获取正确的分型结果。
为了验证该系统的分型能力，对一批HLA二代测序下机数据进行统计，所有数据均通过一代测序确定准确分型结果。本发明的HLA二代测序数据分型方法对一共2310个位点进行了HLA分型，最终，本发明的HLA二代测序数据分型方法准确分型了2297个位点，准确率达到99.4%，存在13个位点因数据质量差未能分型，没有分型错误的位点，错误率为0%。采用市场上认可度较高的GenDX NGSengine商业分析软件对该批数据进行分析，分型错误位点为32个，错误率1.4%，大部分错误产生在DRB1位点与DQB1位点，考虑错误原因主要为设置阈值分型的方法对多态性较高、较难进行序列比对的区域分型能力不足，本发明中的分析方法可以更好的解决这种问题，因此准确率更高。
本发明构建HLA参考基因数据库、设计适用于HLA序列特点的数据分析流程、设计用于分离出杂合碱基位置和纯合碱基位置的聚类算法，以及设计单倍体序列获取算法。可以很好的解决基于二代测序技术的HLA分型准确性的问题。
给本领域技术人员提供上述实施例，以完全公开和描述如何实施和使用所主张的实施方案，而不是用于限制本文公开的范围。对于本领域技术人员而言显而易见的修饰将在所附权利要求的范围内。
