一种确定靶向测序数据中代表性序列的方法、设备和介质与流程

本发明属于生物数据处理，具体地，涉及一种确定靶向测序数据中代表性序列的方法、设备和介质。

背景技术：

1、下一代测序（next-generation sequencing，ngs）又称为高通量测序（high-throughput sequencing），是基于pcr和基因芯片发展而来的边合成边测序技术。高通量测序技术的特点主要有：测序读长短，通量高，准确度高。高通量测序相比一代测序大幅降低了成本，同时保持了较高准确性，并且大幅降低了测序时间，目前高通量测序已经在全组学得到广泛应用。

2、高通量测序得到的原始图像数据经碱基识别（basecalling）转化为原始测序序列（sequenced reads），我们称之为raw data或raw reads，结果以fastq格式存储，其中包含测序序列（reads）信息以及其对应的测序质量信息。

3、fastq格式文件中每个read由四行描述，如下所示：

4、@hwi-st507:248:d29jdacxx:8:1101:1715:1919:1:1:0:acttga/1

5、ntaatattgggctagaaagtatctttgggattgcatgttttgatgcagaatcattgtgccgtagaatgc

6、+

7、bpyccaceceggghhfhhhhhhhhhhffhfhhgfahhchhhhhhfhbfghh_gfhhhhgghefffhhhh

8、其中第一行以“@”开头后跟illumina测序标识符，包括机器型号、上机次数、试剂型号、第几个lane、在flowcell上的坐标、barcode等；第二行是碱基序列；第三行以“+”开头后跟illumina测序标识符（为节省存储空间，部分fq文件会省略“＋”后的信息）；第四行是对应第二行碱基序列的质量值，是用来衡量测序准确度的，字符范围[b,h]，对应质量范围[2,40]。第四行每个字符对应的ascii值减去64，即为该碱基的测序质量值，例如h对应ascii值为104，104-64=40。质量值越高，测序错误率越低。

9、在ngs过程中可能会出现pcr重复，虽然重复似乎是单独的reads，但它们实际上是由于pcr和测序过程中的错误导致的技术噪音。分子标签（mtag）技术是在文库制备过程中与dna片段连接的随机短核苷酸序列。这些mtag序列充当唯一识别码，将每个reads标记为来自单个片段的扩增，为确定pcr重复提供了更准确的机制。但目前如何对针对pcr重复引入的分子标签技术进行去噪尚缺乏统一的标准和方法。

技术实现思路

1、为了解决上述技术问题中的至少一个，本发明采用的技术方案如下：

2、本发明第一方面提供一种基于mtag的靶向测序数据预处理的方法，包括以下步骤：

3、s1，mtag分类与修正：将测序reads基于mtag序列进行分类，包含reads数量小于第一阈值p1的mtag类别为第一类别，其余为第二类别，对于各第一类别中各read，基于其mtag序列中与各第二类别的mtag序列的差异碱基的质量值对其mtag序列进行修正并重新分类；

4、s2，mtag类别选择：选择包含reads数量大于或等于第二阈值p2的mtag类别；

5、s3，代表序列选择：从s2选择的各mtag类别中选择一条代表性序列，

6、其中，p1=3~5；p2=1~10。

7、在本发明中，靶向测序的测序reads包括接头序列、mtag序列及目标序列。在修正前，首先基于目标序列与参考基因组进行比对，过滤掉无法比对至目标区域的测序reads。

8、在本发明中，mtag即分子标签（molecule tag），也叫特异性分子标签（uniquemolecular indentifier，umi），是一段随机或固定或随机固定混合的核苷酸短序列，通常设计为完全随机的核苷酸序列（如nnnnnn），也可以设计为包含固定碱基和随机碱基的核苷酸序列（如nnncnnnngnnnn）。在靶向测序的文库构建过程中，mtag通过连接的方式导入，如同分子条形码一样特异性地标记每个模板。mtag通过给每一个原始dna片段加上一段特有的核苷酸标签序列，经过文库构建及pcr扩增过程之后，一起进行测序。根据不同的mtag序列区分不同来源的dna模板，可以分辨哪些是pcr扩增及测序过程中产生的随机错误造成的假阳性突变，哪些是患者真正携带的突变，从而提高检测的灵敏度和特异性。根据上述分子标签的原理可知，具有相同mtag的测序序列（read）来源于同一dna模板，因为进行分类时归为一个mtag类别。理论上，不同dna模板产生的reads数量不会相差太大，如果某个mtag类别包含的reads数量过低，例如低于上述第一阈值p1，有可能是由于pcr过程中产生的错误，导致mtag序列出错，因此需要修正。

9、在本发明中，第一阈值p1的选择与测序深度也有一定的关系，一般情况下，由于目的区域靶向测序更加经济高效，可以实现500~1000×，甚至更高的测序深度，本发明设定第一阈值p1=3~10，即包含reads数目低于所述第一阈值时，预设该mtag是由于pcr扩增引入的错误（噪声）或者由于测序错误引入。但是如果经过与第二类别的mtag对比，在综合考虑mtag序列差异碱基的质量值后无法进行修正，则相应的mtag类别可能真实来自某个模板，包含的reads数目少是由于pcr的偏好性导致的，该mtag类别需要被保留用于后续分析。

10、在本发明中，质量值表示测序质量值，用于衡量测序准确度，即质量值越高，表明测序的错误率越低。测序错误率用e表示，碱基质量值用q表示，则有下列关系：

11、q=-10log10(e)

12、若碱基的质量值为10，则错误率为10%；若碱基质量值为20，则错误率为1%；若碱基质量值为30，则错误率为0.1%。

13、在本发明的一些实施方案中，步骤s1中基于各第一分类中各read的mtag序列中与各第二类别的mtag序列的差异碱基的质量值对其mtag序列进行修正并重新分类的步骤具体包括：

14、s11，将待修正的第一类别的mtag序列分别与各第二类别的mtag序列进行比对，找相似性最高的第二类别，针对所述待修正的第一类别中的各read分别进行修正：

15、①如果该read的mtag序列与相似性最高的第二类别的mtag序列的差异碱基位于该第二类别的mtag序列的固定碱基位点，则无论该read的分子标签该差异碱基质量值如何，均修正为该第二类别的mtag序列相应位置上的碱基类型；对于一些捕获测序，引入的mtag序列可能全部是随机的序列，则无需根据本步骤进行修正；

16、②如果该read的mtag序列与相似性最高的第二类别的mtag序列的差异碱基位于该第二类别的mtag序列的随机碱基位点，并且该差异碱基的质量值小于第三阈值p3，则将该差异碱基修正为该第二类别的mtag序列相应位置上的碱基类型，否则不进行修正；

17、③若经步骤①~②进行修正后，该read的mtag序列与相似性最高的第二类别的mtag序列的差异碱基占mtag序列长度比例小于第四阈值p4，则无论差异碱基质量值如何，直接将差异碱基修正为该第二类别的mtag序列相应位置的碱基类型，否则该read的mtag序列所有碱基均不进行修正；

18、④对于经步骤①~③无法完成修正的reads，利用相似性第二高的第二类别重复步骤①~③，直至与所有第二类别进行上述修正步骤后均无法完成修正，

19、s12，基于修正后的mtag序列对reads重新进行分类，

20、其中，p3=10~13；p4=10%~30%。

21、在本发明的一些实施方案中，位于随机碱基位点的差异碱基，若其质量值小于第三阈值p3，则认为该差异碱基的质量值过低，可能是测序错误，依据所述相似性最高的第二类别的mtag序列相应位置的碱基进行修正。在本发明的一些具体实施方案中，p3设置为10，即错误率不高于10%，否则需要进行修正，但也不能苛求质量值过高，否则会导致分类错误，当质量值为13时，错误率仅为5%，质量值更高时，错误率更低，此时，再进行修正可能会导致分类错误。因此，设置p3=10~13。

22、在本发明的一些实施方案中，p4值设置过低，会导致一些错误的reads的mtag序列不能被正确修正，p4值设置过高，会导致一些reads的mtag序列被错误修正。发明人经过反复探索，发现p4设置10%~30%之间时，修正的准确率最高。

23、在本发明的一些具体实施方案中，所述mtag序列的长度为16，p4设置为10%。即若经步骤①~②进行修正后，该read的mtag序列与相似性最高的第二类别的mtag序列的差异碱基占mtag序列长度比例小于10%（即仅1个差异碱基），则无论差异碱基质量值如何，直接将差异碱基修正为该第二类别的mtag序列相应位置的碱基类型。

24、在本发明的一些实施方案中，所述第二阈值p2利用以下方法确定：

25、s21，选择具有至少n个reads的mtag类别进行分析时，得到在n个样本中位点的均方根误差，均方根误差的计算公式如下：

26、

27、其中，n为1~5的正整数；rmse为均方根误差，为第i个样本中位点的观察频率，为第i个样本中位点的期望频率；

28、s22，n取不同值时，得到不同的rmse值，将rmse值最小时对应的n值设置为第二阈值p2。

29、在本发明的一些具体实施方案中，若使得rmse值最小时对应的n值为2，则设置第二阈值p2=2，即选择包含reads数量大于或等于2的mtag类别进行后续分析。

30、在本发明的另一些实施方案中，所述第二阈值p2利用以下方法确定：

31、s21’，基于每种mtag类别对应的reads数和累计reads数绘制散点图；

32、s22’，从散点图中第一个点开始，计算任意两个相邻点连线的斜率，当首次出现某条连线的斜率低于前一条连线的斜率时，该连线的起点对应的横坐标值设置为所述第二阈值p2。

33、在本发明的一些具体实施方案中，当首次出现某条连线的斜率低于前一条连线的斜率时，表明该连线的终点对应的横坐标值（每种mtag类别对应的reads数）的reads数目增加开始变缓甚至不变（斜率0），则该连线的起点对应的横坐标值的mtag也需要被考虑纳入后续分析中。

34、在本发明的又一些实施方案中，所述第二阈值p2利用以下方法确定：s21’’，基于每种mtag类别对应的reads数和累计reads数绘制散点图；

35、s22’’，对散点图进行拟合得到曲线方程f（x），对曲线方程f(x)进行求导得到f'(x)，当f’(x)的斜率从小到大首次开始下降时，与对应的横坐标最接近的正整数设置为所述第二阈值p2。

36、在该实施方案中，与上述实施方案核心思想一致，只不过是本实施方案对散点图进行拟合，使得更加精准。在本发明的一些具体实施方案中，采用4阶函数进行曲线拟合，进一步拟合出a、b、c、d、e值，进一步对该进行求导f'(x)。

37、在本发明的一些实施方案中，步骤s3中，利用以下步骤选取代表性序列：

38、s31，对mtag类别中的序列根据目标序列进行计数，计算数目最多的reads在所述mtag类别中的占比，若占比低于第五阈值p5，则舍弃该mtag类别所有reads，否则选择数目最多的reads中的任意一条作为代表性序列，

39、其中，p5=75%~90%。

40、在本发明的一些具体实施方案中，p5是基于四分位法的一种优化方法。由四分位法可知，排在后1/4（即75%）位置上的数叫做第三四分位数，为了进一步提高数据分析质量，对临床基因诊断提供精准辅助，将p5设置为75%~90%，优选地，p5=80%。在本发明的一些实施方案中，若数目最多的reads在所述mtag类别中的占比超过第五阈值，表明其序列确实能够代表该mtag类别，其他reads（如有）可能是由于pcr扩增错误或测序错误引入的。也可能是mtag序列存在错误导致被错误分类至该mtag类别。

41、然而，由于第五阈值p5设置相对较高，可能会导致一些mtag类别由于一些目标序列在pcr过程中引入错误或测序引入错误导致整个mtag类别被舍弃，造成后续分析结果不准确。因此，在本发明的另一些实施方案中，步骤s3中，利用以下步骤选取代表性序列：

42、s31’，从目标序列第一个碱基位置开始，统计所述mtag类别中所有reads该位置的碱基类型及数目：

43、1）若该位置只有一种碱基类型，则无需修正，继续下一碱基的修正；

44、2）若该位置不止一种碱基类型，则对数目非第一多的碱基类型进行修正：

45、①若该碱基质量值不小于第六阈值p6，则无需修正；

46、②若该碱基质量值小于第六阈值p6，则将其修正为数目第一多的碱基类型；

47、③统计该mtag中所有reads该位置修正后的碱基类型及数目，若数目第一多的碱基类型占比不低于第七阈值p7，则将所有非第一多的碱基类型均修正为数目第一多的碱基类型，继续下一碱基的修正；若数目第一多的碱基类型占比低于第八阈值p8，则直接舍弃该mtag类别所有reads，

48、s32’，若按照步骤s31’完成所有碱基的修正，则选择其中任意一条reads作为代表性序列，

49、其中，p6=15~20；p7=75%~90%。

50、在本发明的一些具体实施方案中，由于靶向测序可能是为了寻找变异位点，因此，对碱基的质量值要求更高，不能接受错误率超过3%，则质量值不能低于15。当然，同样不能设置p6过高，否则会导致错误修正。当设置p6=15~20时，修正的正确率最为理想。

51、在本发明的一些具体实施方案中，设置p7的核心思想与设置p5的核心思想相同。

52、在本发明的一些实施方案中，在步骤s31’中，若数目第一多的碱基类型不止一种，则先选择一种碱基类型作为数目第一多的碱基类型，其余作为数目非第一多的碱基类型进行修正，若无法完成修正，则再选择另外一种碱基类型作为数目第一多的碱基类型进行修正，直至完成修正或者舍弃该mtag类别所有reads。

53、在本发明的一些具体实施方案中，若对于一个碱基，g=40%，c=40%，a=20%。此时，先将g作为数目第一多的碱基类型，将c和a均作为数目非第一多的碱基类型，按照上述方法进行修正；若无法完成修正（即修正后数目第一多的碱基类型占比仍低于第八阈值p8），则再针c作为数目第一多的碱基类型，将g和a作为数目非第一多的碱基类型，按照上述方法进行修正，若仍然无法完成修正（占比低于第八阈值p8），则直接舍弃该mtag类别所有reads。若能完成修正，则继续下一碱基的修正。

54、本发明第二方面提供一种计算机设备，包括：

55、存储器，用于存储计算机程序；

56、处理器，用于执行所述计算机程序时实现如本发明第一方面任一所述的一种基于mtag的靶向测序数据预处理的方法的步骤。

57、本发明第三方面提供一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时实现如本发明第一方面任一所述的一种基于mtag的靶向测序数据预处理的方法的步骤。

58、在本发明中，步骤s1、步骤s2和步骤s3可以作为单独的技术方案分别解决相应的技术问题。

59、本发明的有益效果

60、相对于现有技术，本发明取得了以下有益效果：

61、利用本发明的方法，可以对mtag进行修正，避免pcr过程中或测序过程中引入的错误，降低mtag被错误分类的可能性，可以分辨哪些是pcr扩增及测序过程中产生的随机错误造成的假阳性突变，哪些是患者真正携带的突变，从而提高检测的灵敏度和特异性。

62、利用本发明的方法，可以精准地选择合适的mtag类别进行后续分析，进一步提升分析的准确度。

63、利用本发明的方法，可以精准选择mtag类别中的代表性序列，数据利用率更高，分析准确性也得到大大提升。