一种快速获取罕见病致病位点的方法和系统与流程

本发明涉及罕见病基因检测领域，具体涉及一种从基因一代测序结果中快速检测变异，并确定致病位点的方法。

背景技术：

罕见病主要是基因变异导致的遗传病，需要依靠基因检测报告来确诊病因和指导用药。如血友病、地中海贫血、苯丙酮尿症、法布雷病等致病基因明确的单基因遗传罕见病，用二代测序等技术成本高时间长，用一代测序可以低成本高速度地获得该基因的序列，除此之外一代测序还是验证变异位点的金标准技术，在基因检测领域应用广泛。

但现有的一代基因测序结果分析是通过人工判读色谱图分析是否存在变异，无法快速准确地识别一代序列中的变异和基因型，特别是未知的变异。一个基因需要多个一代才能覆盖全长，多个结果的变异分析对人力技术和时间成本都要求较高，而且流程无法质控，不同的技术人员识别的变异结果可能不同，而且存在漏检的风险；检出变异后还要交给报告解读人员查找收录变异的数据库和文献等信息，筛选出致病/可能致病的位点进行致病评级，才能出具变异检测结果报告。对于明确的单基因遗传罕见病，工作重复性高，整体效率低下，需要一款自动化的能够快速全面地获取罕见病致病位点的方法和系统，提高罕见病一代基因检测的报告效率和变异检测全面性和准确性，提高流程可控性。

技术实现要素：

针对现有技术的不足，本发明提供一种快速获取罕见病致病位点的方法和系统，能够从一代测序结果中快速检测并筛选出致病/可能致病变异，并直接输出结果，让检测更全面高效。

为实现上述目的，本发明提供一种快速获取罕见病致病位点的方法和系统，包括如下步骤。

步骤1：提取一代序列并比对到人类参考基因组上，构建比对配对坐标矩阵。其特征在于，从原始测序文件中按照信号阈值提取初级峰序列和次级峰序列两种序列，提取双序列才能检测杂合变异。进一步地，将两种序列同时比对到人类参考基因组上，记录连续正整数坐标上一一配对的参考序列和基因组序列。进一步地，按照比对的染色体位置将初级峰和次级峰矩阵分成同一组，按照比对得分从高到低对分组结果进行排序，构建比对配对坐标矩阵。

步骤2：遍历每组比对配对坐标矩阵，根据差异比对配对的坐标映射关系检测变异。其特征在于，原始序列坐标出现单核苷酸改变时，其参考基因组坐标与原始序列坐标一对一映射。原始序列坐标出现插入变异时，其一个参考基因组坐标映射多个原始序列坐标。原始序列坐标出现缺失变异时，其多个参考基因组坐标映射一个原始序列坐标。

步骤3：根据变异所在的序列进行基因分型。其特征在于，初级峰序列或者次级峰序列坐标出现单一峰检出变异时，输出为杂合变异。初级峰序列和次级峰序列坐标同时检出同一个变异时，输出为纯合变异。初级峰序列和次级峰序列坐标出现变异位点相同，但是变异碱基不同时，输出为多等位基因变异。进一步地，对变异进行左端对齐，最后合并变异结果。

步骤4：变异检测结果用人类罕见病相关数据库和文献进行注释，过滤低测序质量区域和非目标基因比对结果后根据注释结果进行打分，从高到低输出大于0的变异结果。

进一步地，步骤4所述的人类罕见病相关数据库包括但不限于：人类参考基因组功能注释数据库，千人全基因组、人类外显子组亚洲人群变异频率数据库，人类疾病相关变异位点数据库，蛋白质功能预测数据库，已评级变异数据库。其中已评级位点的变异评级数据库，收录变异评级、数据库是否收录、报道变异的文献等信息。

进一步地，步骤4致病/可能致病变异评分规则为：人群频率小于0.01或不存在于人群数据库的位点（1分）；发生在编码区或外显子剪接区的非同义突变（1分）；内含子上蛋白质功能预测分数不少于0.6（1分）；被人类疾病位点数据库或评级数据库注释为致病或可能致病的位点（5分）。

步骤5：用样本性别升级x性染色体上的纯合致病变异排名。若样本为男性，则校正x染色体上纯合变异基因型为半合子。优先报告致病/可能致病的半合子和纯合子的结果。

本发明的有益效果是，与现有技术相比，本方法能快速全面地检测罕见病基因测序中的所有变异，并判定致病/可能致病的变异位点，不仅可以用于已知变异的验证，还可以用于检测未知的变异。此方法避免人工干涉一代序列分析，可快速分析同一个样本的多个文件，同时致病位点的判定标准统一可控，避免人工判读致病位点的误差，极大提升罕见病致病位点报告和解读的效率，并保证整个分析流程可控可追溯。

附图说明

图1是本发明实施例提供的一种快速获取罕见病致病位点的方法和系统的流程图。

图2是本发明实施例中初级峰序列和比对结果示意图。

图3是本发明实施例中次级峰序列和比对结果示意图。

图4是本发明实施例中初级峰比对配对坐标矩阵示意图。

图5是本发明实施例中次级峰比对配对坐标矩阵示意图。

图6是本发明实施例中变异检测结果示意图。

图7是本发明实施例中变异结构vcf格式示意图。

图8是本发明实施例中变异评级数据库示意图。

图9是本发明实施例中致病/可能致病变异打分和筛选结果示意图。

图10是本发明实施例中致病/可能致病变异报告结果示意图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步的说明。以下实施例和附图用于说明本发明，但不用来限制本发明的范围。

本实施例以一个gla基因的一代测序序列（该序列仅用于验证本发明方法的效果及有效性）为例，通过识别文件名可获得全部8个一代测序abi文件，每个文件都利用本发明所述方法进行处理，即比对一代序列并进行变异检测，合并检测结果后注释并进行致病/可能致病变异的打分，最后通过样本基本信息校正变异结果，输出为致病/可能致病结果表格。本实施例提供的实现过程主要以检出致病/可能致病的一代结果文件fabry001.gla.6.f.ab1文件为例，其他文件处理过程同该文件，下文不赘述。

如图1，本发明提供的一种快速获取罕见病致病位点的方法和系统包括以下步骤。

步骤1：比对一代结果的初级峰和次级峰序列到人类参考基因组上，计算比对序列中每个比对配对在参考基因组和原始序列的位置，转化为连续正整数坐标，构建比对配对坐标矩阵。将比对到相同染色体位置的初级峰和次级峰矩阵分成同一组，按照比对得分从高到低对分组结果进行排序。

1.1本实施例中，读取fabry001.gla.6.f.ab1二进制文件，ab1文件包含header标题、directory目录和data数据三部分，提取data中data字段的数据为atcg四种碱基的荧光信号矩阵，b1pt字段为荧光矩阵峰值坐标。提取b1pt峰值坐标对应的荧光信号最强的碱基序列为初级峰序列，按照信号阈值提取次于最强荧光信号的碱基序列为次级峰序列。本实施例中信号阈值取值为0.33，优选地，若需要提高检测灵敏度则需要降低该信号阈值，但是会降低检测特异性。如图2和3，本实施例获得初级峰和次级峰序列均为353bp。

1.2本实施例中，将初级峰和次级峰比对到人类参考基因组上，参考基因组使用homosapiens(human)genomeassemblygrch37版本。如图2，初级峰的第8至351号碱基（共344个）唯一比对到x染色体100653303-100653648位置上，由于存在插入缺失，比对配对共有346对。如图3，次级峰第7-351号碱基（共345个）唯一比对到x染色体100653303-100653649上，比对配对共347个。分别计算比对序列中每个比对配对在参考基因组和原始序列的位置坐标，构建如图4和5的比对配对坐标矩阵。

1.3本实施例中，将比对到相同染色体位置的初级峰和次级峰比对配对坐标矩阵分成同一组，按照比对得分从高到低对分组结果进行排序。本实施例中，初级峰和次级峰均唯一比对到参考基因组x:100653303-100653649区域上，分为同一组，按比对得分排序为第一组，命名为rank1。

步骤2：遍历每组比对配对坐标矩阵，根据差异比对配对的坐标映射关系检测变异，当参考基因组坐标与原始序列坐标一对一映射，检测为单核苷酸改变；当一个参考基因组坐标映射多个原始序列坐标，则检测为插入，合并该坐标上的原始序列作为插入变异序列；当多个参考基因组坐标映射一个原始序列坐标，则检测为缺失，合并该坐标上的参考基因组序列作为缺失变异序列。

2.1本实施例中，如图4，遍历rank1组的初级峰矩阵检测变异，图4中reference列为参考基因组碱基序列，primary.seq为初级峰上与之对应的碱基序列。遍历至初级峰序列第14个碱基，即primary.seq.loc列值为14，图4中用实线黑色方框圈出，并标注为b.var1，第14个碱基映射了x:100653640-100653642三个参考基因组坐标，则判定第14个碱基处存在缺失变异；遍历至第184位碱基，图4中用实线黑框框出并标注为b.var2，识别其唯一映射了100653470位置的参考基因组坐标，则判定为第184位碱基处存在单核苷酸改变；第343位碱基同理不赘述，图4中黑框标注为b.var3。

2.2本实施例中，如图5遍历rank1组的次级峰矩阵检测变异，同初级峰，次级峰也检出了b.var1，b.var2，b.var3这三个变异，除此之外，次级峰上还检测到s.var1和s.var2两个单核苷酸改变，这两个变异不存在于初级峰中，图5中用虚线黑框框出。

2.3本实施例中，每个矩阵遍历结果均需要记录变异的主要特征，记录染色体名称chr，参考基因组坐标pos，一代序列的坐标位置seq.loc，参考基因组序列碱基reference，一代序列碱基seq，一代序列的比对方向strand，1为正向，-1为反向，变异位点离一代序列左右两端的距离trim5和trim3。其中，插入或者缺失变异记录的参考基因组序列和一代序列，需要将连续坐标和序列进行合并，如b.var1变异应位置合并成x:100653640-100653641，参考基因组序列合并为“tt”，一代序列为“-”，变异来源。每个矩阵的每条变异输出一行记录。

步骤3：根据变异所在的序列进行基因分型。当仅单一峰检出变异判定为杂合；当两种峰同时检出同一个变异，判定为纯合变异；当初级峰和次级峰变异位点相同，但是变异碱基不同，判定为多等位基因变异位点杂合变异。输出结果对变异进行左端对齐，然后合并变异结果。

3.1本实施例中，遍历每个分组两种峰的变异结果，根据变异所在的峰进行基因分型。rank1组经过步骤2得到初级峰和次级峰的变异检测结果，如图6，两种峰均检测到b.var1、b.var2这两个变异，则基因型gt一列记录为hom即纯合型；b.var3在两种峰中均检测到，但是碱基序列不同，如图4，在初级峰中b.var3一代序列primary.seq为g，但是在次级峰中b.var3一代序列secondary.seq为c，则b.var3为杂合型，gt一行记为het，同时多等位基因位点需要合并变异结果，如图6的b.var3行，合并alt为g,c，也可输出为2条结果并在合并vcf结果的时候再合并多等位基因位点。此外，还检出了s.var1和s.var2两个只存在于次级峰的单核苷酸改变，因此基因型为杂合。输出结果如图6，除了步骤2.3中记录的变异信息外，还记录了变异的基因型、序列来源和比对分组的排名，以便后续的筛选。

3.2本实施例中，需要把3.1输出的结果再转化为变异的标准记录格式，才能进行后续的数据库注释。如图7，变异标准记录格式为vcf文件格式，由文件题头和变异信息组成，变异信息为10列，即染色体名称chrom，参考基因组位置pos，变异名称id，参考基因组序列ref，变异序列alt，变异质量qual，变异过滤信息filter，变异详细信息info，变异结构化信息索引format和样本具体信息fabry001，其中样本信息一列存储格式化的变异信息，便于后续进行筛选计算，通常基因型gt除了存放在info列外，还会格式化后存在format和样本信息列，便于提取，按照规范gt为0/1表示杂合型，1/1表示纯合型。

3.3本实施例中，步骤3.1的变异结果还需要对齐才能输出为3.2中的变异标准格式，其中单核苷酸改变对齐位置为该变异所在的参考基因组的位置，插入和缺失改变对齐该变异发生的前一位，如图6中的缺失变异b.var1，在参考基因组x:100653640-100653641，应对齐到前一位见图7的x:100653639，同时ref也需要对齐到该位置上，由tt变为ttt，alt则变为t，即ttt缺失变为t，按照标准格式ref和alt不能出现横杠”-“。对齐步骤所需要的序列信息从参考基因组文件中获取。其他变异的对齐不再赘述。

3.4本实施例中，经过上述所有步骤，每一个一代测序的文件输出一个变异结果vcf文件，若无变异则vcf文件只有题头信息。根据vcf的样本具体信息名称，合并同一个样本检出的所有标准的vcf格式输出结果，多等位基因位点会合并为同一个位置，最终结果同一个样本输出一个vcf文件并按照参考基因组的位置进行排序。

步骤4：变异检测结果用人类罕见病相关数据库和文献进行注释，过滤低测序质量区域和非目标基因比对结果后根据注释结果进行打分，从高到低输出大于0的变异结果。

4.1本实施例中，用人类罕见病相关数据库注释步骤3所得的vcf文件，根据参考基因组的位置和变异序列进行注释。本实施例中数据库包括：人类参考基因组功能注释数据库，千人全基因组、人类外显子组亚洲人群变异频率数据库，人类疾病相关变异位点数据库，蛋白质功能预测数据库，已评级变异数据库。除了已评级数据库为自建，其他数据库均为公共数据库，可以到数据库官方网站进行下载。如图8，已评级位点的变异评级数据库为自建数据库，收录变异评级、数据库是否收录、报道变异的文献等信息，图8示例了本实施例中检测到的变异x:100653470a>g，它已被数据库收录，有2篇文献报道，历往acmg标准评级为可能致病突变。acmg变异评级标准为行业内共识的变异评级标准，此处不赘述。

4.2本实施例中，根据步骤4.1中注释结果对变异进行打分，致病/可能致病变异评分规则为：a人群频率小于0.01或不存在于人群数据库的位点（1分）；b发生在编码区或外显子剪接区的非同义突变（1分）；c内含子上蛋白质功能预测分数不少于0.6（1分）；d被人类疾病位点数据库或评级数据库注释为致病或可能致病的位点（5分）。如图9，本实施例中变异经打分后，x:100653470a>g得分为7分，得分最高。

4.3本实施例中，过滤步骤4.2种低测序质量区域和非目标基因比对的变异。本实施例中取trim5和trim3列小于30的变异为低测序质量变异，根据此标准过滤了图9中x:100653311-100653311和x:100653640-100653641位点上的变异，另外本实施例中所有变异均比对到gla基因上，比对排名rank均为1，无非目标基因比对变异。结果按照图9中打分结果score从大到小输出得分大于0的结果。

步骤5：用样本性别升级x性染色体上的纯合致病变异排名。若样本为男性，则校正x染色体上纯合变异基因型为半合子。优先报告致病/可能致病的半合子和纯合子的结果。

5.1本实施例中，样本为男性，x染色体为单个拷贝，因此将步骤4中x染色体上纯合变异结果校正为半合子。同时根据有限报告半合子和纯合子结果的规则，输出得分最高的半合子和纯合子结果，即x:100653470a>g，结果保存成过滤后的变异结果文件。

5.2本实施例中，读取步骤5.1得到的结果文件和原始一代测序原始文件，根据荧光信号矩阵绘制该变异周围的一代测序色谱图，报告变异检测结果，输出为pdf格式检测报告，如图10为去掉敏感信息后的报告结果示意图。

上述是本发明的优选实施方式，通过本发明所述方法能够快速获取罕见病致病位点，本实施例中示例文件检测变异在5s中完成，本实施例样本在3分钟内可完成所有文件分析和生成报告，大大提高了检测效率和准确性。

本发明所公开的实施例和上述说明，用于使本领域专业技术人员能够实现或使用本发明，熟悉本领域的专业技术人员在不违背本发明精神的前提下对本发明做的等同变形或替换的内容，均包含在本申请权利要求所限定的范围内。

sequencelisting

<110>厦门基源医疗科技有限公司

<120>一种快速获取罕见病致病位点的方法和系统

<130>2020.12.9

<160>6

<170>patentinversion3.3

<210>1

<211>353

<212>dna

<213>homosapiens

<400>1

cagacgagctgtggagtggtttctccatatgggtcatctaggtaactttaagaatgtttc60

ctcctctcttgtttgaattatttcattctttttctcagttagtgattggcaactttggcc120

tcagctggaatcagcaagtaactcagatggccctctgggctatcatggctgctcctttat180

tcacgtctaatgacctccgacacatcagccctcaagccaaagctctccttcaggataagg240

acgtaattgccatcaatcaggaccccttgggcaagcaagggtaccagcttagacaggtaa300

ataagagtatatattttaagatggctttatatacccaatacccactttgtcaa353

<210>2

<211>353

<212>dna

<213>homosapiens

<400>2

ctaacgtgctgtggagtggtttctccatatgggtcatctaggtaactttaagaatgtttc60

ctcctctcctgtttgaattatttcattctttttctcagttagtgattggcacctttggcc120

tcagctggaatcagcaagtaactcagatggccctctgggctatcatggctgctcctttat180

tcacgtctaatgacctccgacacatcagccctcaagccaaagctctccttcaggataagg240

acgtaattgccatcaatcaggaccccttgggcaagcaagggtaccagcttagacaggtaa300

ataagagtatatattttaagatggctttatatacccaataccgactttgtcaa353

<210>3

<211>346

<212>dna

<213>homosapiens

<400>3

gctgtggaaagtggtttctccatatgggtcatctaggtaactttaagaatgtttcctcct60

ctcttgtttgaattatttcattctttttctcagttagtgattggcaactttggcctcagc120

tggaatcagcaagtaactcagatggccctctgggctatcatggctgctcctttattcatg180

tctaatgacctccgacacatcagccctcaagccaaagctctccttcaggataaggacgta240

attgccatcaatcaggaccccttgggcaagcaagggtaccagcttagacaggtaaataag300

agtatatattttaagatggctttatatacccaataccaactttgtc346

<210>4

<211>347

<212>dna

<213>homosapiens

<400>4

tgctgtggaaagtggtttctccatatgggtcatctaggtaactttaagaatgtttcctcc60

tctcttgtttgaattatttcattctttttctcagttagtgattggcaactttggcctcag120

ctggaatcagcaagtaactcagatggccctctgggctatcatggctgctcctttattcat180

gtctaatgacctccgacacatcagccctcaagccaaagctctccttcaggataaggacgt240

aattgccatcaatcaggaccccttgggcaagcaagggtaccagcttagacaggtaaataa300

gagtatatattttaagatggctttatatacccaataccaactttgtc347

<210>5

<211>21

<212>dna

<213>homosapiens

<400>5

cctttattcatgtctaatgac21

<210>6

<211>21

<212>dna

<213>homosapiens

<400>6

cctttattcacgtctaatgac21