一种可保留分子标签信息的在测序数据中模拟小变异的方法、装置及计算机可读存储介质

摘要

本发明公开了一种可保留分子标签信息的在测序数据中模拟小变异的方法、装置及计算机可读存储介质。实验证明，与现有技术模拟小变异工具(BAMSurgeon和VarBen)进行比较，本发明建立的模拟小变异的方法SafeMut能够更准确地模拟高通量测序数据中的小变异特征；通过SafeMut计算模拟得到的变异等位基因突变丰度和通过做生物学实验得到的变异等位基因丰度相比，有非常类似的平均值，方差和统计学分布。本发明的SafeMut可应用于建立更真实的生物信息变异检测使用的标准数据、通过更准确的变异模拟设定更准确的空白检测线(LOB)和样品检测线(LOD)等，可以弥补样本稀缺和实验消耗人力物力等问题。

权利要求

1.基于测序数据模拟变异的方法，其特征在于，所述方法包括将所述测序数据比对到参考基因组获得比对文件，将已知变异文件整合到所述比对文件得到整合比对文件的步骤；

所述整合通过对所述比对文件中的每条测序读长进行如下步骤实现：A1）所述测序数据为带分子标签测序数据，根据测序读长的位置信息和测序读长的分子标签生成伪随机数；

或，所述测序数据为不带分子标签的测序数据，根据测序读长的位置信息和BAM文件格式中定义的读长ID（QNAME）生成伪随机数；

A2）根据需要模拟变异的染色体号和基因组位置，生成对数正态分布随机变量；

A3）计算模拟变异后的测序读长的突变丰度；

A4）比较所述伪随机数和所述突变丰度，确定是否将模拟后的变异合并入测序读长。

2.根据权利要求1所述的方法，其特征在于：所述变异为单核苷酸变异、缺失变异和/或插入变异。

3.根据权利要求1或2所述的方法，其特征在于：A1）步骤中所述伪随机数为U，所述突变丰度为G；所述U小于所述G时，将所述模拟后的变异合并入所述测序读长，得到整合模拟变异后的测序读长。

4.根据权利要求3所述的方法，其特征在于：A4）步骤具体包括如下步骤：A4‑1）所述变异为单核苷酸变异，将所述变异单核苷酸的碱基测序质量值转换成错误概率，同时生成一个伪随机数，比较所述伪随机数和所述错误概率，确定合并入测序读长的核苷酸类型；

A4‑2）所述变异为缺失变异，删除测序读长中所述缺失变异对应的野生型核苷酸位置的核苷酸；

A4‑3）所述变异为插入变异，在测序读长的野生型核苷酸中所述插入变异对应位置插入所述插入变异的核苷酸。

5.基于测序数据模拟变异的装置，其特征在于，所述装置包括如下模块：B1）数据比对模块：用于将所述测序数据比对到参考基因组获得比对文件；

B2）模拟变异模块：用于将已知变异文件整合到所述比对文件得到整合比对文件；

所述整合通过包括如下步骤的方法实现：对所述比对文件中的每条测序读长进行如下步骤：C1）所述测序数据为带分子标签测序数据，根据测序读长的位置信息和测序读长的分子标签生成伪随机数；

或，所述测序数据为不带分子标签的测序数据，根据测序读长的位置信息和BAM文件格式中定义的读长ID（QNAME）生成伪随机数；

C2）根据需要模拟变异的染色体号和基因组位置，生成对数正态分布随机变量；

C3）计算模拟变异后的测序读长的突变丰度；

C4）比较所述伪随机数和所述突变丰度，确定是否将模拟后的变异合并入测序读长。

6.根据权利要求5所述的装置，其特征在于：所述变异为单核苷酸变异、缺失变异和/或插入变异。

7.根据权利要求5或6所述的装置，其特征在于：C1）步骤中所述伪随机数为U，所述突变丰度为G；所述U小于所述G时，将所述模拟后的变异合并入所述测序读长，得到整合模拟变异后的测序读长。

8.一种存储有计算机程序的计算机可读存储介质，所述计算机程序使计算机建立如权利要求1‑4中任一权利要求所述方法的步骤，或所述计算机程序使计算机建立如权利要求

5‑7中任一权利要求所述装置的模块。

9.权利要求1‑4中任一权利要求所述方法或权利要求5‑7中任一权利要求所述装置的下述任一应用：D1）在制备保留分子标签测序数据的生物信息学分析标准参考数据中的应用；

D2）在制备生物测序数据的生物信息学分析标准参考数据中的应用；

D3）在测试生物信息学分析软件或流程中的应用。

说明书

一种可保留分子标签信息的在测序数据中模拟小变异的方

法、装置及计算机可读存储介质

技术领域

[0001] 本发明属于生物信息技术领域，具体涉及一种可保留分子标签信息的在测序数据中模拟小变异的方法、装置及计算机可读存储介质。

背景技术

[0002] 随着测序技术的发展，高通量测序(high—throughput sequencing(HTS)又称next‑generation sequencing(NGS))的通量越来越大，因此产生的数据量也越来越大。与此同时，高通量测序的应用也越来越广泛并且越来越重要。

[0003] 2011年，Dr Bert Vogelstein发明了用分子标签(molecular barcode)给高通量测序数据纠错的方法(Kinde,I.,Wu,J.,Papadopoulos,N.,Kinzler,K.W.,&Vogelstein,B.(2011).Detection and quantification of rare mutations with massively parallel sequencing.Proceedings of the National Academy of Sciences,108(23),9530‑9535.)。分子标签的别名是UMI(unique molecular identifier，唯一分子标识符)。高通量测度的测序错误率在千分之一左右，但是通过使用分子标签校正测序错误，可以把错误率降到万分之一。因此目前通常使用分子标签去重纠错的测序检测cell‑free DNA里低丰度的肿瘤突变。分子标签的发明让液态活检行业得到了极大的发展，因为与做手术得到肿瘤组织样本相比，抽血得到肿瘤游离DNA要容易得多。

[0004] 分子标签(molecular barcode，别称unique molecular identifier(UMI))去重纠错的测序已经有广泛的应用，尤其在小变异(small variant)的检测纠错上有极大的应用，其中小变异指的是点突变(SNV)和插入缺失(insertion‑deletion,简称InDel，包括插入和缺失两种子类型的变异)。因此目前需要在分子标签去重纠错数据中模拟小变异的方法，但是目前并没有这方面的方法。虽然有研究人员尝试完全用in silico(只基于人参考基因组序列(也就是说不基于任何测序真实数据)的模拟)的方法模拟分子标签去重纠错的测序数据(Sater,V.,Viailly,P.‑J.,Lecroq,T.,Ruminy,P.,Bérard,C.,Prieur‑Gaston,&Jardin,F.(2020).UMI‑Gen:A UMI‑based read simulator for variant calling evaluation in paired‑end sequencing NGS libraries.Computational and structural biotechnology journal,18,2270‑2280.)，但是完全用in silico的方法无法模拟所有的测序噪音。除此之外，在高深度测序的情况下，突变丰度并不遵从二项分布(Muyas,F.,Bosio,M.,Puig,A.,Susak,H.,Domènech,L.,Escaramis,G.,...Rabionet,R.(2019).Allele balance bias identifies systematic genotyping errors and false disease associations.Human mutation,40(1),115‑126.)。研究人员曾令发现在每个指定位点，系统性错误的REF跟ALT的深度比值比遵循某个中位数是u，并且偏离到u÷2和u×2‑15/10的概率正好在中位数概率的10 的lognormal分布(Zhao,X.,Hu,A.C.,Wang,S.,&Wang,X.(2021).Calling small variants using universality with Bayes‑factor‑adjusted odds ratios.Briefings in Bioinformatics.)，因此系统性错误也可以被模拟。

[0005] 通过对样本做湿实验得到的数据是行业公认的金标准数据。但是做湿实验既耗费人力物力，也需要比较长的时间周期，并且含有稀有变异的样本是难以或得到的。因此湿实验数据经常不够用。所以必须需要对变异进行模拟，产生模拟变异的NGS数据，才可以弥补样本稀缺和实验消耗人力物力的问题。但是通过计算机模拟呈现的变异和通过对真实样本做湿实验呈现的变异有较大的区别，从而影响用于变异检测的数据的有效性。

[0006] 伪随机数(pseudo‑random number)：用确定性的算法计算得到的来自[0,1](包括0但是不包括1)均匀分布的随机小数序列，并不真正的随机，但具有类似于随机数的统计特征，如均匀性、独立性等。在计算伪随机数时，若使用的初值(种子)不变，那么伪随机数的数序也不变。伪随机数可以用计算机大量生成，在模拟研究中为了提高模拟效率，一般采用伪随机数代替真正的随机数。伪随机数生成器(pseudo‑random number generator)用于生成伪随机数。

[0007] 突变丰度(allele fraction，即等位基因突变丰度)：某个基因位点所有的等位基因中，突变的等位基因的所占相对比例(相对野生型等位基因)，即等于：

[0008] 突变型支持数目/(野生型支持数目+突变型支持数目)。

发明内容

[0009] 本发明所要解决的技术问题是如何模拟高通量测序数据中的小变异和/或如何在含有分子标签信息的高通量测序数据中模拟小变异和/或如何制备可保留分子标签信息的高通量测序数据的生物信息学分析的参考数据。

[0010] 为了解决上述技术问题，本发明首先提供了基于测序数据模拟变异的方法，所述方法可包括将所述测序数据比对到参考基因组获得比对文件，将已知变异文件整合到所述比对文件得到整合比对文件的步骤。

[0011] 所述整合可通过对所述比对文件中的每条测序读长进行如下步骤实现：

[0012] A1)所述测序数据为带分子标签测序数据，根据测序读长的位置信息和测序读长的分子标签生成伪随机数；

[0013] 或，所述测序数据为不带分子标签的测序数据，根据测序读长的位置信息和BAM文件格式中定义的读长ID(QNAME)生成伪随机数；

[0014] A2)根据需要模拟变异的染色体号和基因组位置，生成对数正态分布随机变量；

[0015] A3)计算模拟变异后的测序读长的突变丰度；

[0016] A4)比较所述伪随机数和所述突变丰度，确定是否将模拟后的变异合并入测序读长。

[0017] 上述方法中，所述变异可为单核苷酸变异、缺失变异和/或插入变异。

[0018] 上述方法中，A1)步骤中所述伪随机数为U，所述突变丰度为G；所述U小于所述G时，将所述模拟后的变异合并入所述测序读长，得到整合模拟变异后的测序读长。

[0019] 上述方法中，所述突变丰度G可通过如下式Ⅴ计算获得：

[0020]

[0021] 式Ⅴ中， NRV＝exp(N(0，σ^2))；其中保证σ能让NRV＝log(2)(‑15/10)的概率是NRV＝0的概率的10 。这里的σ的取值约等于0.186484。

[0022] 上述方法中，所述伪随机数U的计算公式可为如下式Ⅰ：

[0023] U＝Wang_hash(X31_hash_string(UMI序列)

[0024] xor Wang_hash(read终止位置xor Wang_hash(read起始位置xor Wang_hash(伪随机种子))))式Ⅰ

[0025] 式Ⅰ中的函数“Wang_hash”和“X31_hash_string”可来源于khash.h文件(网址https://github.com/attractivechaos/klib/blob/master/khash.h，“Wang_hash”在khash.h文件中名字为__ac_Wang_hash，“X31_hash_string”在khash.h文件中名字为__ac_X31_hash_string；式Ⅰ中的xor可代表是异或运算，伪随机种子可为默认值int2randint(13，1)。int2randint可为生成伪随机数种子的函数。

[0026] 所述待模拟小变异测序数据可为带有分子标签的测序数据或不带有分子标签的测序数据。在不带有分子标签的测序数据中，所述伪随机数U的计算公式中可使用BAM文件格式中定义的读长ID(QNAME)代替分子标签UMI序列计算获得所述伪随机数U。

[0027] 上述方法中，A4)步骤可具体包括如下步骤：A4‑1)所述变异为单核苷酸变异，将所述变异单核苷酸的碱基测序质量值转换成错误概率，同时生成一个伪随机数，比较所述伪随机数和所述错误概率，确定合并入测序读长的核苷酸类型；

[0028] A4‑2)所述变异为缺失变异，删除测序读长中所述缺失变异对应的野生型核苷酸位置的核苷酸；

[0029] A4‑3)所述变异为插入变异，在测序读长的野生型核苷酸中所述插入变异对应位置插入所述插入变异的核苷酸。

[0030] 上述方法中，A4‑1)所述伪随机数可为伪随机数x，所述伪随机数x的数值计算公式可为如下式Ⅵ：

[0031] Wang_hash(Wang_hash(Wang_hash(int2randint(13,4))xor X31_hash_string(read名字))xor(read起始位置))

[0032] 式Ⅵ。

[0033] 式Ⅵ中的函数“Wang_hash”和“X31_hash_string”可来源于khash.h文件(网址https://github.com/attractivechaos/klib/blob/master/khash.h，“Wang_hash”在khash.h文件中名字为__ac_Wang_hash，“X31_hash_string”在khash.h文件中名字为__ac_X31_hash_string；式Ⅰ中的xor代表是异或运算，伪随机种子的默认值是int2randint(13，1)。int2randint是一个生成伪随机数种子的函数。

[0034] 为了解决上述技术问题，本发明还提供了基于测序数据模拟变异的装置，所述装置可包括如下模块：

[0035] B1)数据比对模块：用于将所述测序数据比对到参考基因组获得比对文件；

[0036] B2)模拟变异模块：用于将已知变异文件整合到所述比对文件得到整合比对文件；

[0037] 所述整合可通过包括如下步骤的方法实现：对所述比对文件中的每条测序读长进行如下步骤：

[0038] C1)所述测序数据为带分子标签测序数据，根据测序读长的位置信息和测序读长的分子标签生成伪随机数；

[0039] 或，所述测序数据为不带分子标签的测序数据，根据测序读长的位置信息和BAM文件格式中定义的读长ID(QNAME)生成伪随机数；

[0040] C2)根据需要模拟变异的染色体号和基因组位置，生成对数正态分布随机变量；

[0041] C3)计算模拟变异后的测序读长的突变丰度；

[0042] C4)比较所述伪随机数和所述突变丰度，确定是否将模拟后的变异合并入测序读长。

[0043] 上述装置中，所述突变丰度G可通过如下式Ⅴ计算获得：

[0044]

[0045] 式Ⅴ中， NRV＝exp(N(0，σ^2))；其中保证σ能让NRV＝log(2)(‑15/10)的概率是NRV＝0的概率的10 。这里的σ的取值约等于0.186484。

[0046] 上述装置中，所述伪随机数U的计算公式可为如下式Ⅰ：

[0047] U＝Wang_hash(X31_hash_string(UMI序列)

[0048] xor Wang_hash(read终止位置xor Wang_hash(read起始位置xor Wang_hash(伪随机种子))))式Ⅰ

[0049] 式Ⅰ中的函数“Wang_hash”和“X31_hash_string”可来源于khash.h文件(网址https://github.com/attractivechaos/klib/blob/master/khash.h，“Wang_hash”在khash.h文件中名字为__ac_Wang_hash，“X31_hash_string”在khash.h文件中名字为__ac_X31_hash_string；式Ⅰ中的xor代表是异或运算，伪随机种子可为默认值int2randint(13，1)。int2randint可为生成伪随机数种子的函数。

[0050] 所述待模拟小变异测序数据可为不带有分子标签的测序数据。如果不带有分子标签，所述伪随机数U的计算公式中就使用BAM文件格式中定义的读长ID(QNAME)代替分子标签UMI序列计算获得所述伪随机数U。

[0051] 上述C4)步骤具体可包括如下步骤：

[0052] C4‑1)所述变异为单核苷酸变异，将所述变异单核苷酸的碱基测序质量值转换成错误概率，同时生成一个伪随机数，比较所述伪随机数和所述错误概率，确定合并入测序读长的核苷酸类型；

[0053] C4‑2)所述变异为缺失变异，删除测序读长中所述缺失变异对应的野生型核苷酸位置的核苷酸；

[0054] C4‑3)所述变异为插入变异，在测序读长的野生型核苷酸中所述插入变异对应位置插入所述插入变异的核苷酸。

[0055] C4‑1)所述伪随机数可为伪随机数x，所述伪随机数x的数值计算公式可为如下式Ⅵ：

[0056] Wang_hash(Wang_hash(Wang_hash(int2randint(13,4))xor X31_hash_string(read名字))xor(read起始位置))

[0057] 式Ⅵ。

[0058] 式Ⅵ中的函数“Wang_hash”和“X31_hash_string”可来源于khash.h文件(网址https://github.com/attractivechaos/klib/blob/master/khash.h，“Wang_hash”在khash.h文件中名字为__ac_Wang_hash，“X31_hash_string”在khash.h文件中名字为__ac_X31_hash_string；式Ⅰ中的xor代表是异或运算，伪随机种子的默认值是int2randint(13，1)。int2randint是一个生成伪随机数种子的函数。

[0059] 上述装置中，所述变异可为单核苷酸变异、缺失变异和/或插入变异。

[0060] 上述装置中，C1)步骤中所述伪随机数可为U，所述突变丰度可为G。所述U小于所述G时，将所述模拟后的变异合并入所述测序读长，得到整合模拟变异后的测序读长。

[0061] 为了解决上述技术问题，本发明还提供了一种存储有计算机程序的计算机可读存储介质。所述计算机程序可使计算机建立如上文任一所述方法的步骤。所述计算机程序也可使计算机建立如上文任一所述装置的模块。

[0062] 为了解决上述技术问题，本发明还提供了一种存储有计算机程序的计算机可读存储介质。所述计算机程序可使计算机运行如上文任一所述方法的步骤。所述计算机程序也可使计算机运行如上文任一所述装置的模块。

[0063] 上文所述方法或上文所述装置的下述任一应用也属于本发明的保护范围：

[0064] D1)在制备保留分子标签测序数据的生物信息学分析标准参考数据中的应用；

[0065] D2)在制备生物测序数据的生物信息学分析标准参考数据中的应用；

[0066] D3)在测试生物信息学分析软件或流程中的应用。

[0067] 上述应用中，所述测序数据可为全基因组重测序数据、外显子测序数据和/或Panel测序数据。所述测序数据可为带有分子标签的测序数据，所述测序数据也可为不带有分子标签的测序数据。

[0068] 本发明的目的是提供一种可保留分子标签信息的在高通量测序数据中模拟小变异(点突变(single‑nucleotide variant，简称SNV)和插入缺失(insertion‑deletion，也称之为InDel))的方法SafeMut。

[0069] 实验证明，使用本发明建立的SafeMut和现有技术模拟小变异工具(BAMSurgeon和VarBen)进行比较，能够更准确地模拟高通量测序数据中的小变异特征。通过使用本发明计算模拟得到的变异等位基因突变丰度，和通过做生物学实验得到的变异等位基因丰度相比，有非常类似的平均值，方差和统计学分布。因此本发明有很多重要的应用，例如建立更真实的生物信息变异检测使用的标准数据，通过更准确的变异模拟设定更准确的空白检测线(LOB)和样品检测线(LOD)等，可以弥补样本稀缺和实验消耗人力物力等问题。

附图说明

[0070] 图1为本专利技术方案的示意图。

[0071] 图2为使用SRA索引号是SRP162370的全外显子测序(whole exome sequencing,WES)数据，评估SafeMut软件模拟正常组织得到的(期望的)和在肿瘤组织实际观测到的等位基因的平均值相似性

[0072] 图3为使用SRA索引号是SRP162370的全外显子测序(whole exome sequencing,WES)数据，评估BAMSurgeon软件模拟正常组织得到的(期望的)和在肿瘤组织实际观测到的等位基因的平均值相似性。

[0073] 图4为使用SRA索引号是SRP162370的全外显子测序(whole exome sequencing,WES)数据，评估VarBen软件模拟正常组织得到的(期望的)和在肿瘤组织实际观测到的等位基因的平均值相似性。

[0074] 图5为使用SRA索引号是SRP162370的全外显子测序(whole exome sequencing,WES)数据，评估SafeMut软件模拟正常组织得到的(期望的)和在肿瘤组织实际观测到的等位基因的方差相似性。

[0075] 图6为使用SRA索引号是SRP162370的全外显子测序(whole exome sequencing,WES)数据，评估BAMSurgeon软件模拟正常组织得到的(期望的)和在肿瘤组织实际观测到的等位基因的方差相似性。

[0076] 图7为使用SRA索引号是SRP162370的全外显子测序(whole exome sequencing,WES)数据，评估VarBen软件模拟正常组织得到的(期望的)和在肿瘤组织实际观测到的等位基因的方差相似性。图8为使用SRA索引号是SRP162370的全外显子测序(whole exome sequencing,WES)数据，评估SafeMut软件模拟正常组织得到的(期望的)和在肿瘤组织实际观测到的等位基因丰度之间的分布形状相似性。或者更具体而言：使用所有相关样本，计算每个位点的等位基因丰度的平均值和方差，然后计算这些丰度的Z值并且以分位图的形式在图中展现这些Z值。

[0077] 图9为使用SRA索引号是SRP162370的全外显子测序(whole exome sequencing,WES)数据，评估BAMSurgeon软件模拟正常组织得到的(期望的)和在肿瘤组织实际观测到的等位基因丰度之间的分布形状相似性。或者更具体而言：使用所有相关样本，计算每个位点的等位基因丰度的平均值和方差，然后计算这些丰度的Z值并且以分位图的形式在图中展现这些Z值。

[0078] 图10为使用SRA索引号是SRP162370的全外显子测序(whole exome sequencing,WES)数据，评估VarBen软件模拟正常组织得到的(期望的)和在肿瘤组织实际观测到的等位基因丰度之间的分布形状相似性。或者更具体而言：使用所有相关样本，计算每个位点的等位基因丰度的平均值和方差，然后计算这些丰度的Z值并且以分位图的形式在图中展现这些Z值。

[0079] 图11为使用SRA索引号是SRP296025的分子标签测序数据，评估SafeMut软件模拟白细胞得到的(期望的)和在稀释的游离DNA(cell‑free DNA)实际观测到的等位基因的平均值相似性。或者更具体而言：使用所有相关样本，计算每个位点的等位基因丰度的平均值，并且在图中展现这些平均值。

[0080] 图12为使用SRA索引号是SRP296025的分子标签测序数据，评估VarBen软件模拟白细胞得到的(期望的)和在稀释的游离DNA(cell‑free DNA)实际观测到的等位基因的平均值相似性。或者更具体而言：使用所有相关样本，计算每个位点的等位基因丰度的平均值，并且在图中展现这些平均值。

[0081] 图13为使用SRA索引号是SRP296025的分子标签测序数据，评估SafeMut软件模拟白细胞得到的(期望的)和在稀释的游离DNA(cell‑free DNA)实际观测到的等位基因的方差相似性。或者更具体而言：使用所有相关样本，计算每个位点的等位基因丰度的方差，并且在图中展现这些方差。

[0082] 图14为使用SRA索引号是SRP296025的分子标签测序数据，评估VarBen软件模拟白细胞得到的(期望的)和在稀释的游离DNA(cell‑free DNA)实际观测到的等位基因的方差相似性。或者更具体而言：使用所有相关样本，计算每个位点的等位基因丰度的方差，并且在图中展现这些方差。

[0083] 图15为使用SRA索引号是SRP296025的分子标签测序数据，评估SafeMut软件模拟白细胞得到的(期望的)和在稀释的游离DNA(cell‑free DNA)实际观测到的等位基因的分布形状相似性。或者更具体而言：使用所有相关样本，计算每个位点的等位基因丰度的平均值和方差，然后在算出Z值，并且在图以分位图的形式展现这些Z值。

[0084] 图16为使用SRA索引号是SRP296025的分子标签测序数据，评估VarBen软件模拟白细胞得到的(期望的)和在稀释的游离DNA(cell‑free DNA)实际观测到的等位基因的分布形状相似性。或者更具体而言：使用所有相关样本，计算每个位点的等位基因丰度的平均值和方差，然后在算出Z值，并且在图以分位图的形式展现这些Z值。

具体实施方式

[0085] 下面结合具体实施方式对本发明进行进一步的详细描述，给出的实施例仅为了阐明本发明，而不是为了限制本发明的范围。以下提供的实施例可作为本技术领域普通技术人员进行进一步改进的指南，并不以任何方式构成对本发明的限制。

[0086] 下述实施例中的实验方法，如无特殊说明，均为常规方法，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。

[0087] 实施例1、基于高通量分子标签测序数据模拟小变异(基因点突变和插入缺失突变)的方法SafeMut

[0088] 1.测序数据获得(FASTQ格式文件)：通过高通量测序分别得到疾病样本和配对的正常样本的测序数据。

[0089] 2.测序数据比对：将测序数据比对到参考基因组。

[0090] 用BWA MEM算法(下载网址：https://github.com/lh3/bwa)将样本测序数据FASTQ文件中的测序读长序列(reads，即测序过程中生成的带有碱基质量值的序列)比对到参考基因组，输出BAM格式比对文件。

[0091] 3.使用SafeMut对正常样本测序数据模拟小变异(包括基因点突变和插入缺失突变变异)

[0092] 使用SafeMut模拟小变异的过程为：将已知变异文件整合到比对文件得到整合比对文件(BAM格式或FASTQ格式)。

[0093] 具体为输入步骤2得到的BAM格式的比对文件(名称为inBAM)和VCF格式的已知变异文件(名称为inVCF)进行整合，输出FASTQ格式的整合比对文件。其中整合过程是通过把inVCF中的已知变异放置到inBAM中的reads而实现的。

[0094] inVCF的具体下载地址在实施例2和实施例3中有详述。

[0095] 具体通过如下步骤实现：

[0096] 对于每条来源于inBAM的read和每个这条read覆盖的需要模拟的变异(来源于inVCF)，做以下操作：

[0097] 3.1生成伪随机数

[0098] 根据插入片段的位置信息(包括插入片段起始位置，插入片段终止位置)和插入片段的分子标签(UMI序列)生成一个从0到1之间的伪随机数U，其中插入片段指的是通过PCR扩增之后得到的序列，而read指的是通过测序仪测序得到的序列，插入片段长度可大于read的长度；伪随机数U的计算公式如下：

[0099] U＝Wang_hash(X31_hash_string(UMI序列)

[0100] xor Wang_hash(read终止位置xor Wang_hash(read起始位置xor Wang_hash(伪随机种子))))式Ⅰ

[0101] 式Ⅰ中的函数“Wang_hash”和“X31_hash_string”来源于khash.h文件(网址https://github.com/attractivechaos/klib/blob/master/khash.h，“Wang_hash”在khash.h文件中名字为__ac_Wang_hash，“X31_hash_string”在khash.h文件中名字为__ac_X31_hash_string；式Ⅰ中的xor指的是异或运算，伪随机种子的默认值是int2randint(13,1)。int2randint是一个生成伪随机数种子的函数。函数int2randint出自于以下标准：ISO/IEC 9899:201x函数srand和rand(int2randint(m,n)首先运行srand(m)，然后运行n次rand()，取最后rand()的返回值作为int2randint的返回值)。

[0102] 3.2根据需要模拟变异read的染色体号和基因组位置，生成对数正态分布随机变量。

[0103] 具体为根据read覆盖的对应于已知变异文件(inVCF)中相同位置的需要模拟的变异的染色体号和基因组位置，生成两个从0到1之间的伪随机数A和B，然后通过使用Box‑Muller transform(Box,1958)生成一个正态(normal)分布随机变量，然后再把这个正态分布随机变量转换成对数正态分布(lognormal distribution)随机变量。

[0104] X(伪随机种子)＝Wang_hash(参考基因组位置xor Wang_hash(染色体ID xor(Wang_hash(伪随机种子))))式Ⅱ。

[0105] A＝X(seedA)  式Ⅲ。

[0106] B＝X(seedB)  式Ⅳ。

[0107] 其中比率seedA和seedB默认值分别是比对文件inBAM的前50和500条reads的read名字的X31_hash_string函数数值相加的结果。例如read名字可以是：SRR7890831.1。

[0108] 3.3计算模拟变异后的read的突变丰度。

[0109] 如果模拟前read的突变丰度是F，那么模拟后read的突变丰度的计算公式如下：

[0110]

[0111] 式Ⅴ中， NRV＝exp(N(0,σ^2))；其中保证σ能让NRV＝log(2)(‑15/10)
的概率是NRV＝0的概率的10 。

[0112] 这里的N是通过使用前面所述的Box‑Muller transform生成的正态分布随机变量，而exp函数把这个正态分布转换成前面所述的对数正态分布，而且这里的σ的取值约等于0.18648(相关文献：Niu,D.,Li,L.,Yu,Y.,Zang,W.,Li,Z.,Zhou,L.,...Cheng,G.(2020).Evaluation of next generation sequencing for detecting HER2 copy number in breast and gastric cancers.Pathology&Oncology Research,26(4),2577‑2585；Zhao,X.,Hu,A.C.,Wang,S.,&Wang,X.(2021).Calling small variants using universality with Bayes‑factor‑adjusted odds  ratios.Briefings in Bioinformatics)。

[0113] 3.4比较伪随机数和模拟变异后的read的突变丰度，确定是否将模拟后的变异合并入read。

[0114] 如果伪随机数小于模拟变异后的read的突变丰度，即U

[0115] 在整合模拟变异的过程中，根据模拟变异的类型，修改read对应位置的碱基：

[0116] 3.4.1对于单核苷酸变异的修改

[0117] 如果需要模拟的变异是单核苷酸变异(SNV)，那么把变异核苷酸的碱基测序质量值(Phred score，Qphred)转换成一个0到1之间的错误概率值；同时生成一个0到1之间的伪随机数x。

[0118] 比较伪随机数x和变异核苷酸错误概率，如果伪随机数x小于错误概率，那就把read上的对应的核苷酸转换成{A，C，G，T}中随机挑选的核苷酸；否则(即伪随机数x大于或等于错误概率)就把read上的对应的核苷酸转换成已知变异文件中的对应的突变型的核苷酸。其中伪随机数x的数值计算公式是：

[0119] Wang_hash(Wang_hash(Wang_hash(int2randint(13,4))xor X31_hash_string(read名字))xor(read起始位置))式Ⅵ。

[0120] 式Ⅵ中的函数“Wang_hash”和“X31_hash_string”来源于khash.h文件(网址https://github.com/attractivechaos/klib/blob/master/khash.h，“Wang_hash”在khash.h文件中名字为__ac_Wang_hash，“X31_hash_string”在khash.h文件中名字为__ac_X31_hash_string；式Ⅰ中的xor指的是异或运算，伪随机种子的默认值是int2randint(13,1)。int2randint是一个生成伪随机数种子的函数。函数int2randint出自于以下标准：ISO/IEC 9899:201x函数srand和rand。

[0121] 3.4.2对于缺失变异的修改

[0122] 如果需要模拟的变异是缺失，那就在read的相应位置删除核苷酸，即根据已知变异文件中与read对应位置的野生型序列相比缺失的核苷酸，删除该read的野生型核苷酸，得到模拟后的缺失变异。

[0123] 3.4.3对于插入变异的修改

[0124] 如果需要模拟的变异是插入，那就在read的相应位置插入核苷酸，即根据已知变异文件中与read对应位置的野生型序列相比插入的核苷酸，在该read的野生型核苷酸中插入已知变异文件的插入核苷酸，得到模拟后的插入变异。

[0125] 插入的核苷酸默认的碱基质量值是30。

[0126] 实施例2、本发明SafeMut方法与现有技术对全外显子测序数据模拟小变异的结果对比

[0127] 下载NCBI数据库SRA索引号是SRP162370的测序数据(Fang,L.T.,Zhu,B.,Zhao,Y.,Chen,W.,Yang,Z.,Kerrigan,L.,...Idler,K.(2021).Establishing community reference samples,data and call sets for benchmarking cancer mutation detection using whole‑genome sequencing.Nature biotechnology,39(9),1151‑1160.)，分别使用本发明实施例1中建立的SafeMut和现有技术的模拟小变异工具BAMSurgeon及VarBen进行模拟变异后的结果比较。

[0128] SRA索引号是SRP162370的测序数据由美国食品药品管理局(FDA)领导的SEQC2联盟生成。SEQC2通过对HCC1395三阴乳腺癌肿瘤组织细胞系样本(T)和HCC1395BL非肿瘤配对样本(N)分别进行6次全外显子测序，得到了6组T‑N配对的测序数据集，其中有572个已知的真实变异(包含2个插入和26个缺失，其余的变异都是点突变)被所有12次(6次T测序和6次N测序)的外显子测序覆盖。

[0129] 该组测序数据来源于疾病组织，不需要用分子标签纠错，因此测序数据中没有分子标签信息。在使用本发明建立的SafeMut方法模拟测序数据小变异的过程中，对于该没有分子标签的测序数据，在计算伪随机数U的过程中，使用BAM文件格式中定义的读长ID(QNAME)代替分子标签进行伪随机数U的计算，然后再进行变异模拟(相关文献：Li,H.,Handsaker,B.,Wysoker,A.,Fennell,T.,Ruan,J.,Homer,N.,...Subgroup,G.P.D.P.(2009).The Sequence alignment/map(SAM)format and SAMtools.Bioinformatics,25(16),2078‑2079)。

[0130] 其中BAMSurgeon模拟小变异的方法参照如下文献进行：Ewing,A.D.,Houlahan,K.E.,Hu,Y.,et.,al.(2015).Combining tumor genome simulation with crowdsourcing to benchmark somatic single‑nucleotide‑variant detection.Nature methods,12(7),623‑630。

[0131] VarBen模拟小变异的方法参照如下文献进行：李子阳.(2019).体细胞基因突变高通量测序检测生物信息学分析参考物质的研究.北京协和医学院。

[0132] 本实施例所用的inVCF文件(已知变异文件)的具体下载地址为：https://ftp‑trace.ncbi.nlm.nih.gov/ReferenceSamples/seqc/Somatic_Mutation_WG/release/v1.2/，具体文件为：high‑confidence_sINDEL_in_HC_regions_v1.2.vcf.gz和high‑confidence_sSNV_in_HC_regions_v1.2.vcf.gz。

[0133] 1.对正常组织样本测序数据模拟小变异

[0134] 1.1测序数据获得(FASTQ格式文件)：包括疾病样本数据和配对的正常组织样本数据。

[0135] 1.2测序数据比对：将测序数据比对到参考基因组。

[0136] 用BWA MEM算法把所有样本(包括疾病样本数据和配对的正常组织样本数据)的FASTQ文件的reads分别比对到GRCh38参考基因组(Li,2013)(Li,H.(2013).Aligning sequence reads,clone sequences and assembly contigs with BWA‑MEM.arXiv preprint arXiv:1303.3997.)，输出BAM格式比对文件。BWA MEM算法对应的软件下载网址为：https://github.com/lh3/bwa。

[0137] 1.3分别用SafeMut、BAMSurgeon和VarBen这三种方法对正常组织样本测序数据的BAM格式比对文件进行小变异模拟得到相应的整合模拟变异后的数据文件。每款软件都会生成模拟后的BAM格式文件或者FASTQ格式文件。如果输出是FASTQ文件，就用上一步的方法通过BWA MEM序列比对把FASTQ文件转化成BAM文件。BAMSurgeon软件的下载网址为https://github.com/adamewing/bamsurgeon/；VarBen软件的下载网址为：https://github.com/nccl‑jmli/VarBen。模拟的变异在文献给出的链接中有详细介绍(Fang,L.T.,Zhu,B.,Zhao,Y.,Chen,W.,Yang,Z.,Kerrigan,L.,...Idler,K.(2021).Establishing community reference samples,data and call sets for benchmarking cancer mutation detection using whole‑genome sequencing.Nature biotechnology,39(9),1151‑1160.)。此文献成果出自于与美国FDA领导的变异检测标准化项目。研究人员通过对HCC1395/HCC1395BL肿瘤/正常组织配对样本在不同时间不同地点用不同测序技术进行大量测序，得到肿瘤金标准体系变异。

[0138] 2.不同软件模拟结果比较

[0139] 对于每一种方法，对于数据中的每一批次，计算每个模拟的变异在(1)肿瘤组织样本的BAM格式比对文件数据中和(2)配对的正常组织样本经过模拟变异后得到的BAM格式整合比对文件数据中的等位基因突变丰度(简称突变丰度)相关信息，包括：

[0140] 数据中所有批次的等位基因突变丰度的平均值(mean)(所有样本突变丰度的相加的总和除以样本数目)、方差(variance)和用Z值归一化后的Z值分位图(z‑score quantile‑quantile plot)。突变型测序深度是零的位点没有相关信息。

[0141] 2.1比较三种软件模拟小变异后得到的(期望的)正常组织样本的等位基因突变丰度的平均值和在肿瘤组织实际观测到的该等位基因突变丰度的平均值的相似性。

[0142] 具体为使用所有相关样本，计算每个模拟用的位点的等位基因突变丰度的平均值，并且在图中展现这些平均值和比较两种样本平均值的相似性。SafeMut方法的结果如图2所示，BAMSurgeon方法的结果如图3所示，VarBen方法的结果如图4所示。通过比较可以看出，使用本发明实施例1建立的模拟小变异的方法SafeMut模拟小变异得到的模拟位点的等位基因突变丰度的平均值和实际实验操作得到的等位基因突变丰度(简称突变丰度)的平均值高度相似，大幅度超过了BAMSurgeon和VarBen模拟结果得到的相似度，因此本发明的方法SafeMut对小变异的模拟在三种方法中性能是最好的。

[0143] 突变丰度(allele fraction)：某个基因位点所有的等位基因中，突变的等位基因所占相对比例(相对野生型等位基因)，即等于：

[0144] 突变型支持数目/(野生型支持数目+突变型支持数目)。

[0145] 2.2比较三种软件模拟变异后得到的(期望的)正常组织样本的等位基因突变丰度的方差和在肿瘤组织实际观测到的等位基因的方差相似性。

[0146] 具体为使用所有相关样本，计算每个模拟位点的等位基因突变丰度的方差，并且在图中展现这些方差和比较两种样本方差的相似性。SafeMut方法的结果如图5所示，BAMSurgeon方法的结果如图6所示，VarBen方法的结果如图7所示。通过比较可以看出，使用本发明实施例1建立的模拟小变异的方法SafeMut模拟小变异得到模拟位点的等位基因丰度的方差和实际实验操作得到的等位基因丰度的方差高度相似，大幅度超过了BAMSurgeon和VarBen模拟结果得到的相似度，因此本发明方法SafeMut对小变异的模拟在三种方法中性能是最好的。

[0147] 2.3比较三种软件模拟变异后得到的(期望的)正常组织样本的等位基因丰度之间的分布形状和在肿瘤组织实际观测到的等位基因丰度之间的分布形状的相似性。

[0148] 具体为使用所有相关样本，计算每个模拟位点的等位基因突变丰度的平均值和方差，然后计算这些丰度的Z值并且以分位图的形式在图中展现这些Z值，并比较两种样本方差的相似性。SafeMut方法的结果如图8所示，BAMSurgeon方法的结果如图9所示，VarBen方法的结果如图10所示。通过比较可以看出，使用本发明实施例1建立的模拟小变异的方法SafeMut模拟得到的等位基因突变丰度的统计学分布和实际实验操作得到的等位基因突变丰度的统计学分布高度相似，大幅度超过了BAMSurgeon和VarBen模拟结果得到的相似度，因此本发明方法SafeMut的对小变异的模拟性能是最好的。

[0149] 因此，使用本发明实施例1中建立的SafeMut和现有技术模拟小变异工具进行比较，能够更准确地模拟变异特征。通过使用本发明计算模拟得到的变异等位基因突变丰度，和通过做生物学实验得到的变异等位基因丰度相比，有非常类似的平均值，方差和统计学分布。因此本发明有很多重要的应用，例如建立更真实的生物信息变异检测使用的标准数据，通过更准确的变异模拟设定更准确的空白检测线(LOB)和样品检测线(LOD)等。

[0150] 实施例3、本发明SafeMut方法与现有技术对分子标签测序数据模拟小变异的结果对比

[0151] 下载NCBI数据库SRA索引号是SRP296025的测序数据(相关文献：Deveson,I.W.,Gong,B.,Lai,K.,LoCoco,J.S.,Richmond,T.A.,Schageman,J.,...Jones,W.(2021).Evaluating the analytical validity of circulating tumor DNA sequencing assays for precision oncology.Nature biotechnology,1‑14.)，分别使用本发明实施例1中建立的SafeMut和现有技术的模拟小变异工具BAMSurgeon及VarBen进行模拟小变异后的结果比较。该数据为Illumina和IDT公司上传的FASTQ格式的分子标签数据。原始数据来源于不同批次的样本，每个批次有两个样本，分别是一个游离DNA(cell‑free DNA)样本和配对的正常配对样本(gDNA)，其中游离DNA含有肿瘤游离DNA(circulating tumor DNA)。然后提取分子标签。

[0152] SRA索引号是SRP296025的测序数据由美国食品药品管理局(FDA)领导的SEQC2联盟生成。SEQC2通过对安捷伦UHRR细胞系(T，(Novoradovskaya et al.,2004))和安捷伦已知阴性细胞系(N，安捷伦产品号#:5190‑8848)进行物理混合，得到物理混合之后的细胞系。然后SEQC2把物理混合得到的样本当作肿瘤游离DNA(T),把已知阴性细胞系当作正常配对样本(N)，对T和N分别进行了12次测序，得到了12组T‑N配对的测序数据集，其中有425个已知的真实小变异(包含2个插入和8个缺失，其余的变异都是点突变)被24次测序(12次T测序和12次N测序)覆盖。

[0153] 1.对正常配对样本测序数据模拟小变异

[0154] 1.1分子标签测序数据获得(FASTQ格式文件)：包括疾病样本数据和配对的正常样本数据。

[0155] 1.2分子标签测序数据比对：将测序数据比对到参考基因组。

[0156] 用BWA MEM算法把所有FASTQ文件的reads比对到GRCh37参考基因组(相关文献：Li,H.(2013).Aligning sequence reads,clone sequences and assembly contigs with BWA‑MEM.arXiv preprint arXiv:1303.3997)，输出BAM格式文件。

[0157] 1.3分别用用SafeMut、BAMSurgeon和VarBen这三种方法对BAM格式比对文件的正常配对样本测序数据进行变异模拟得到相应的整合模拟变异后的数据文件。每款软件都会生成模拟后的BAM格式文件或者FASTQ格式文件，如果输出是FASTQ文件，就用上一步的方法通过BWA MEM序列比对把FASTQ文件转化成BAM文件。

[0158] 本实施例所用的inVCF文件(已知变异文件)的下载地址为：https://figshare.com/articles/dataset/Consensus_Target_Region/13511829(KnownPositives_h g19.vcf.gz)

[0159] 2.不同软件模拟结果比较

[0160] 分别对于每一种方法，对于每一批次数据，按照实施例2中的方法计算每个模拟变异在(1)肿瘤游离DNA样本的BAM比对文件数据和(2)配对的非肿瘤白细胞经过变异模拟后的BAM格式整合比对文件数据中的丰度相关信息，包括：所有批次的丰度的平均值(mean)，方差(variance)和用Z值归一化后的Z值分位图(z‑score quantile‑quantile plot)。

[0161] 2.1比较三种软件模拟小变异后得到的(期望的)正常的非肿瘤白细胞样本的等位基因平均值和在肿瘤游离DNA样本中实际观测到的该等位基因的平均值相似性。

[0162] 具体为使用所有相关样本，计算每个模拟位点的等位基因突变丰度的平均值，并且在图中展现这些平均值和比较两种样本平均值的相似性。SafeMut方法的结果如图11所示，VarBen方法的结果如图12所示，BAMSurgeon因为不适用于高深度分子标签测序数据所以运行出错，没有相关结果。通过比较可以看出使用本发明实施例1建立的模拟小变异的方法SafeMut模拟得到的等位基因突变丰度的平均值和实际实验操作得到的等位基因丰度的平均值高度相似，大幅度超过了BAMSurgeon和VarBen模拟结果得到的相似度，因此本发明方法SafeMut对小变异的模拟性能是最好的)。

[0163] 突变丰度(allele fraction)：某个基因位点所有的等位基因中，突变的等位基因的所占相对比例(相对野生型等位基因)，即等于：

[0164] 突变型支持数目/(野生型支持数目+突变型支持数目)。

[0165] 2.2比较三种软件模拟小变异后得到的(期望的)正常组织样本的等位基因的方差和在肿瘤组织实际观测到的等位基因的方差相似性。

[0166] 具体为使用所有相关样本，计算每个模拟位点的等位基因突变丰度的方差，并且在图中展现这些方差和比较两种样本方差的相似性。SafeMut方法的结果如图13所示，VarBen方法的结果如图14所示，BAMSurgeon因为不适用于高深度分子标签测序数据所以运行出错，没有相关结果。通过比较可以看出使用本发明实施例1建立的模拟小变异的方法SafeMut模拟得到的等位基因突变丰度的方差和实际实验操作得到的等位基因丰度的方差高度相似，大幅度超过了BAMSurgeon和VarBen模拟结果得到的相似度，因此本发明方法SafeMut对小变异的模拟性能是最好的。

[0167] 2.3比较三种软件模拟小变异后得到的(期望的)正常样本的等位基因突变丰度之间的分布形状和在肿瘤组织实际观测到的等位基因丰度之间的分布形状的相似性。

[0168] 具体为使用所有相关样本，计算每个模拟位点的等位基因突变丰度的平均值和方差，然后计算这些丰度的Z值并且以分位图的形式在图中展现这些Z值，并比较两种样本方差的相似性。SafeMut方法的结果如图15所示，VarBen方法的结果如图16所示，BAMSurgeon运行出错，没有相关结果。通过比较可以看出使用本发明实施例1建立的模拟小变异的方法SafeMut模拟得到的等位基因突变丰度的统计学分布和实际实验操作得到的等位基因突变丰度的统计学分布高度相似，大幅度超过了BAMSurgeon和VarBen模拟结果得到的相似度，因此本发明方法SafeMut对小变异的模拟性能是最好的。

[0169] 因此，使用本发明实施例1中建立的SafeMut和现有技术模拟小变异工具进行比较，能够更准确地模拟变异特征。通过使用本发明计算模拟得到的变异等位基因突变丰度，和通过做生物学实验得到的变异等位基因丰度相比，有非常类似的平均值，方差和统计学分布。因此本发明有很多重要的应用，例如建立更真实的生物信息变异检测使用的标准数据，通过更准确的变异模拟设定更准确的空白检测线(LOB)和样品检测线(LOD)等。

[0170] 本发明建立的模拟小变异的方法SafeMut能更真实地模拟小变异，让通过计算机模拟得到的结果与做湿实验得到的结果高度相似。可以弥补样本稀缺和实验消耗人力物力的问题。

[0171] 以上对本发明进行了详述。对于本领域技术人员来说，在不脱离本发明的宗旨和范围，以及无需进行不必要的实验情况下，可在等同参数、浓度和条件下，在较宽范围内实施本发明。虽然本发明给出了特殊的实施例，应该理解为，可以对本发明作进一步的改进。总之，按本发明的原理，本申请欲包括任何变更、用途或对本发明的改进，包括脱离了本申请中已公开范围，而用本领域已知的常规技术进行的改变。按以下附带的权利要求的范围，可以进行一些基本特征的应用。