CN114507723A - 一种测序信号归一化的方法 - Google Patents
一种测序信号归一化的方法 Download PDFInfo
- Publication number
- CN114507723A CN114507723A CN202210103665.XA CN202210103665A CN114507723A CN 114507723 A CN114507723 A CN 114507723A CN 202210103665 A CN202210103665 A CN 202210103665A CN 114507723 A CN114507723 A CN 114507723A
- Authority
- CN
- China
- Prior art keywords
- sequencing
- signal
- reaction
- sequencing reaction
- correction
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000012163 sequencing technique Methods 0.000 title claims abstract description 278
- 238000010606 normalization Methods 0.000 title claims abstract description 21
- 238000012937 correction Methods 0.000 claims abstract description 62
- 238000000034 method Methods 0.000 claims abstract description 48
- 238000006243 chemical reaction Methods 0.000 claims description 160
- 239000002773 nucleotide Substances 0.000 claims description 50
- 125000003729 nucleotide group Chemical group 0.000 claims description 50
- 239000003153 chemical reaction reagent Substances 0.000 claims description 37
- 239000000178 monomer Substances 0.000 claims description 29
- 108090000623 proteins and genes Proteins 0.000 claims description 27
- 238000010348 incorporation Methods 0.000 claims description 6
- 229920000388 Polyphosphate Polymers 0.000 claims description 3
- 239000001205 polyphosphate Substances 0.000 claims description 3
- 235000011176 polyphosphates Nutrition 0.000 claims description 3
- 230000000875 corresponding effect Effects 0.000 description 27
- 238000004364 calculation method Methods 0.000 description 11
- 239000011159 matrix material Substances 0.000 description 8
- 239000000758 substrate Substances 0.000 description 6
- 108020004707 nucleic acids Proteins 0.000 description 4
- 150000007523 nucleic acids Chemical class 0.000 description 4
- 102000039446 nucleic acids Human genes 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 3
- 150000002500 ions Chemical class 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 238000003786 synthesis reaction Methods 0.000 description 3
- 230000002159 abnormal effect Effects 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 238000010438 heat treatment Methods 0.000 description 2
- 238000005286 illumination Methods 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 239000004005 microsphere Substances 0.000 description 2
- 239000000203 mixture Substances 0.000 description 2
- 229920000642 polymer Polymers 0.000 description 2
- 239000000376 reactant Substances 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 238000005406 washing Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- XPPKVPWEQAFLFU-UHFFFAOYSA-J diphosphate(4-) Chemical compound [O-]P([O-])(=O)OP([O-])([O-])=O XPPKVPWEQAFLFU-UHFFFAOYSA-J 0.000 description 1
- 235000011180 diphosphates Nutrition 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000007274 generation of a signal involved in cell-cell signaling Effects 0.000 description 1
- 239000000017 hydrogel Substances 0.000 description 1
- 238000002347 injection Methods 0.000 description 1
- 239000007924 injection Substances 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Images
Classifications
-
- C—CHEMISTRY; METALLURGY
- C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
- C12Q—MEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
- C12Q1/00—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions
- C12Q1/68—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions involving nucleic acids
- C12Q1/6869—Methods for sequencing
- C12Q1/6874—Methods for sequencing involving nucleic acid arrays, e.g. sequencing by hybridisation
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
- G16B20/30—Detection of binding sites or motifs
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Chemical & Material Sciences (AREA)
- Health & Medical Sciences (AREA)
- Engineering & Computer Science (AREA)
- Organic Chemistry (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Physics & Mathematics (AREA)
- Biotechnology (AREA)
- Analytical Chemistry (AREA)
- Genetics & Genomics (AREA)
- Zoology (AREA)
- Biophysics (AREA)
- Molecular Biology (AREA)
- General Health & Medical Sciences (AREA)
- Wood Science & Technology (AREA)
- Evolutionary Biology (AREA)
- Theoretical Computer Science (AREA)
- Immunology (AREA)
- Microbiology (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Medical Informatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biochemistry (AREA)
- General Engineering & Computer Science (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
- Investigating, Analyzing Materials By Fluorescence Or Luminescence (AREA)
Abstract
本发明提供了一种稳定高效的测序信号归一化方法，独创地利用了随机序列区的群体信号均值进行衰减拟合和均值校正，利用校正序列区的已知信号，可以获得测序的单位信号。相对于单点归一化方法，本发明的群体归一化方法具有更高的可靠性和准确性，并能够更加高效而有针对性地处理测序数据。
Description
技术领域
本发明涉及一种测序信号归一化的方法，属于基因测序领域。
背景技术
在所有的边合成边测序的测序体系中，测序信号的强度都与作为反应物的模板DNA分子量成正比，随着测序反应轮数的增加，模板DNA分子会随着加热、光照、通电、试剂冲刷等各种物理化学扰动逐渐损失，其损伤的机理在不同测序仪中可能各不相同，但造成的结果都是相应的测序信号逐渐下降。测序信号的下降若不能被正确地估计和校正则会影响对测序信息的读出，可能导致测得碱基个数产生偏差，或者对测得碱基种类发生误判。因此对测序信号衰减的校正是所有测序仪技术中必不可少的一个关键步骤。
Ion Torrent测序仪采用了基于单个测序点的一倍信号分段拟合衰减的方法，单个数据点能提供的数据量是十分有限的，如果要提高信息量需要尽量提高读长增加测序轮数，然而测序轮数大的信号本身质量也会越来越差，总体而言这种方法的信噪比十分受限，对衰减的校正准确度不高。在单位信号的计算上Ion Torrent采用了校正序列区信号的前4轮对ATCG四碱基各测量一个一倍信号，分别用这4轮信号标定单位信号，此方法的设计受限于其电化学测序的信号生成机理，4碱基信号不一致性较大，必须分别进行标定，效率偏低。
因此，需要开发一种可以用于简并测序，且更加准确、效率更高、抗噪性更强的衰减校正方法。
发明内容
本发明公开一种测序信号归一化的方法，其特征在于，包括以下步骤：
A.利用基因测序芯片对由待测基因构建的文库进行多轮测序反应，得到多幅原始测序图像；所述文库中校正序列区紧邻待测基因，所述校正序列区的核苷酸序列已知；一轮测序反应指的是，提供一种测序试剂，将该种测序试剂中具有可检测标记的核苷酸单体掺入所述文库之后检测所述可检测标记生成的荧光信号的过程；
B.在每轮测序反应所得的所述原始测序图像中，对图像中的亮点提取获得初步信号，计算初步信号的平均值作为该轮测序反应的初始群体信号强度；所述亮点对应于芯片上发生核苷酸单体掺入的微反应室；
C.利用除去所述校正序列区之后的每一轮测序反应对应的所述初始群体信号强度进行指数衰减拟合，得到拟合曲线和相关参数；
D.根据所述拟合曲线和相关参数，计算得到所述校正序列区的调整后的信号强度，根据调整后的部分校正序列区信号总和及其对应的核苷酸总数，得到所述亮点处的单位信号；
E.结合单位信号和所述拟合曲线得到所述亮点在每一轮测序反应的归一化信号值。
根据优选的实施方式，所述测序为3’端不封闭的测序反应，在测序反应中，包括两种不同的测序试剂：第一测序试剂和第二测序试剂；两种测序试剂循环加入；其中所述第一测序试剂包含具有可检测标记的两种不同的核苷酸单体；所述第二测序试剂包含具有可检测标记的两种核苷酸单体，且所述两种核苷酸单体不同于所述第一测序试剂中存在的所述核苷酸单体，并且其中所述第二测序试剂是在提供了所述第一测序试剂随后提供的，将所述核苷酸单体掺入所述文库之后检测所述可检测标记生成的荧光信号。
根据优选的实施方式，所述核苷酸单体为5’端多磷酸修饰有荧光切换性质的荧光团的核苷酸，所述的荧光切换性质指的是测序反应后荧光信号相比测序反应前有明显改变。
根据优选的实施方式，所述文库包括：至少一个测序引物结合位点、0个或至少一个索引序列、校正序列区、待测基因。
根据优选的实施方式，所述校正序列区的核苷酸序列已知，组成所述校正序列区的核苷酸个数为6-10个，所述校正序列区对应的测序反应轮数大于2。
根据优选的实施方式，所述提取获得初步信号包括：扣除芯片背景，识别微反应室位置，通过计算每个微反应室的灰度值总和得到微反应室强度，区分有核苷酸单体掺入的亮室和无核苷酸单体掺入的暗室，通过暗室强度计算荧光本底，亮室强度扣除所述荧光本底得到初步信号。
根据优选的实施方式，所述拟合曲线为：
其中c是测序反应轮数，
为第c轮测序反应的初始群体信号强度，
为扣除h的第1轮测序反应的初始群体信号强度，d是拟合所得的群体衰减常数，h是零倍信号校正常数，反映了实际信号与理想线性信号的少量截距偏差，lkey是校正序列区对应的测序反应轮数。
根据优选的实施方式，按照以下公式对校正序列区的初始群体信号强度进行调整：
其中c是测序反应轮数，sa(c)为第c轮测序反应的调整后信号强度，s(c)为第c轮测序反应的初始群体信号强度，n(c)为第c轮测序反应的校正序列区对应的碱基个数，
为扣除h的第1轮测序反应的群体信号强度，d是拟合所得的群体衰减常数，h是零倍信号校正常数，反映了实际信号与理想线性信号的少量截距偏差，lkey是校正序列区对应的测序反应轮数。
根据优选的实施方式，步骤D中所述的部分校正序列区中不包括校正序列区对应的最后一轮测序反应的碱基。
根据优选的实施方式，所述归一化信号值可通过以下校正公式进行求解：
其中c是测序反应轮数，sa(c)为第c轮测序反应的调整后信号强度，d是拟合所得的群体衰减常数，c是测序反应轮数，h是零倍信号校正常数，u是对应的单位信号。
根据优选的实施方式，根据步骤E所得的归一化信号值得到所述待测基因的简并碱基序列。
根据优选的实施方式，所述方法进一步包括，在得到所述归一化信号值后，对所述归一化信号值进行失相校正，去除超前反应和滞后反应造成的信号偏差。
本发明的有益之处
本发明为简并测序体系提供了一种稳定高效的信号归一化方法，独创地利用了随机序列的群体均值进行衰减拟合和均值校正，具有以下优点：
(1)相对于单点归一化，本方法具有更高的可靠性和准确性，并能够更加高效而有针对性地处理简并测序数据。
(2)利用随机序列的群体均值拟合衰减曲线，整体校正信号衰减，抗噪性能优秀。
(3)利用随机序列群体信息校准校正序列区信号均值，在单位信号计算上对校正序列区序列设计无依赖性，归一化计算性能稳定，鲁棒性好。
附图说明
图1.归一化信号获取流程图。
图2.简并测序不同测序反应轮的初始群体信号强度图，其中横坐标表示测序轮次，纵坐标表示初始群体信号强度。
图3.初始群体信号强度的拟合曲线图。图中叉号表示的是校正序列区对应的信号，圆点表示的是待测基因对应的测序信号。
图4.芯片上某一测序反应位点在不同测序轮次的归一化信号值，图中，横坐标表示的是测序反应轮次，纵坐标表示的是归一化信号值，对于大多数测序反应轮，归一化信号值可近似为整数。
具体实施方式
除非另有定义，否则本文使用的所有科技术语具有与本领域普通技术人员通常理解的含义相同的含义。为了更好地公开本发明的方法和内容，在此对于本发明中较为关键的术语做详细的解释说明。
2+2简并测序
所述简并测序为3’端不封闭的测序反应，2+2简并测序中，包括两种不同的测序试剂：第一测序试剂和第二测序试剂；两种测序试剂循环加入；其中所述第一测序试剂包含具有可检测标记的两种不同的核苷酸单体；所述第二测序试剂包含具有可检测标记的两种核苷酸单体，且所述两种核苷酸单体不同于所述第一测序试剂中存在的所述核苷酸单体，并且其中所述第二测序试剂是在提供了所述第一测序试剂随后提供的，将所述核苷酸单体掺入所述文库之后检测所述可检测标记生成的荧光信号。
一轮测序反应
即一个测序cycle，或称一个测序反应轮，一轮测序反应指的是，提供一种测序试剂，将该种测序试剂中具有可检测标记的核苷酸单体掺入所述文库之后检测所述可检测标记生成的荧光信号的过程；对于2+2简并测序，以MK测序流程为例，当待测序列为AGTAGG时，每一轮测序反应所延伸的碱基分别为(A,GT,A,GG)，即进行了4轮测序反应。
测序信号归一化
本发明中，测序信号归一化主要指的是从原始测序信号中经过衰减校正和单位信号计算两步，将测序信号强度转化为延伸的碱基数目的过程。
衰减(decay)
在所有边合成边测序的测序体系中，测序信号的强度都与作为反应物的模板DNA分子量成正比，随着测序反应轮数的增加，模板DNA分子会随着加热、光照、试剂冲刷等各种物理化学扰动逐渐损失，也伴随着测序延伸链脱落，模板损伤造成无法继续延伸等，使得相应的测序信号逐渐下降，上述过程即为衰减。测序信号的下降若不能被正确地估计和校正则会影响对测序信息的读出，可能导致测得碱基个数产生偏差，或者对测得碱基种类发生误判。因此对测序信号衰减的校正是测序技术中必不可少的一个关键步骤。
单位信号
单位信号指的是每个测序位点上延伸一个碱基的信号强度，是一个正比于该测序点模板DNA分子数的物理量，测序信号的强度即等于单位信号乘延伸反应的碱基个数。无论对任何测序技术来说，单位1的测量精度都直接影响到测序结果的准确度。因此单位信号的计算方法对于任何测序方法都意义重大。
2倍信号均值曲线
进行2+2简并测序反应，当碱基分布完全随机时，对于每一轮测序反应，测序信号的分布符合如下数学模型：1倍信号的发生概率为1/2，2倍信号的发生概率为1/22,3倍信号的发生概率为1/23，n倍信号的发生概率为1/2n。对一个GC比均衡的用随机打断建库方式构建的测序文库进行测序时，只要选取一个足够大的群体，序列数目足够多，则其碱基分布就可以视为完全随机，即符合上述信号分布规律。根据此分布，信号倍数的数学期望为：
因此信号均值等价于2倍单位信号的均值。
校正序列区
校正序列区或者叫做key区，紧邻待测基因，且位于其3’端，是在测序时先于待测基因序列测得的序列，组成key区的核苷酸数量及序列是已知的，组成所述校正序列区的核苷酸个数为6-10个，所述校正序列区对应的测序反应轮数大于2。例如，校正序列可以是AGTAGG，以MK测序流程为例，每一轮反应所延伸的碱基分别为(A,GT,A,GG)，即进行4轮测序反应，满足反应轮数大于2；若序列为ACTTTG，则以MK流程测序每一轮反应所延伸的碱基分别为(AC,TTTG)，即进行了2轮测序反应，不满足反应轮数大于2，所以该序列不能作为校正序列。
简并多聚物长度(DPL)
2+2简并测序为多碱基测序，区别于单碱基测序每轮反应只延伸一个核苷酸分子，多碱基测序每轮反应延伸的核苷酸可能是多个，测序反应释放的荧光信号强度与释放的荧光基团数目成正相关，在没有衰减和失相的理想条件下，每轮反应的荧光信号反映了该轮延伸的碱基数，被称为简并多聚物长度(degenerate polymer length,DPL)。以序列AAGCTGTCCAGG的MK流程为例，每一轮反应所延伸的碱基分别为(AA,G,C,TGT,CCA,GG)，因此DPL为(2,1,1,3,3,2)。实际测得的简并多聚物长度往往并不是整数，此值叫做DPM(简并多聚物长度测量值)。
具体的，本发明公开了一种利用群体信息和指数衰减模型对测序信号归一化的方法，其特征在于，包括以下步骤：
A.利用基因测序芯片对由待测基因构建的文库进行多轮测序反应，得到多幅原始测序图像；所述文库中校正序列区紧邻待测基因，所述校正序列区的核苷酸序列已知；一轮测序反应指的是，提供一种测序试剂，将该种测序试剂中具有可检测标记的核苷酸单体掺入所述文库之后检测所述可检测标记生成的荧光信号的过程；
B.在每轮测序反应所得的所述原始测序图像中，对图像中的亮点提取获得初步信号，计算初步信号的平均值作为该轮测序反应的初始群体信号强度；所述亮点对应于芯片上发生核苷酸单体掺入的微反应室；
C.利用除去所述校正序列区之后的每一轮测序反应对应的所述初始群体信号强度进行指数衰减拟合，得到拟合曲线和相关参数；
D.根据所述拟合曲线和相关参数，计算得到所述校正序列区的调整后的信号强度，根据调整后的部分校正序列区信号总和及其对应的核苷酸总数，得到所述亮点处的单位信号；
E.结合单位信号和所述拟合曲线得到所述亮点在每一轮测序反应的归一化信号值。
本发明公开的方法的总体流程可参见图1，在利用基因测序芯片对待测文库测序得到原始测序数据(荧光图像)后，首先判断此数据是否合规，例如不符合对焦标准、旋转角度过大的图片为不合规的图片，将其标记为异常，这些数据不进行后续的校正流程；对于合规的数据，首先对其进行衰减校正，拟合衰减曲线，根据衰减曲线调整校正序列区的信号值，据此计算单位信号，结合单位信号和衰减曲线最终得到每个测序位点在每一轮测序反应中的归一化信号值，并输出。
在优选的实施方式中，本发明的测序方法为2+2简并测序，2+2简并测序是边合成边测序技术的一种新形式，是一种3’端不封闭的基因测序方法，除了第一轮进样外每轮测序反应可以有1至n个碱基反应，测序中没有空轮次。与其他测序方法不同，一次2+2简并测序并不直接得到精确的由4种碱基组成的序列，而是得到由2种简并碱基组成的简并碱基序列。简并序列按照标准简并碱基(degenerate nucleotide)标识，写作M/K，R/Y，W/S，M为A/C，K为G/T；R为A/G，Y为C/T，W为A/T，S为G/C。简并碱基序列中所包含的信息略小于4碱基序列，但测序速度得到大幅提高。其获得一种简并序列所需的反应轮数是逐碱基测序的一半，是Ion Torrent的1*4测序的1/6。
根据优选的实施方式，所述核苷酸单体为5’端多磷酸修饰有荧光切换性质的荧光团的核苷酸，所述的荧光切换性质指的是测序反应后荧光信号相比测序反应前有明显改变。优选的，测序反应后荧光信号相比测序反应前有明显增强。
本发明中，按照IUPAC符号命名规则(Nucleic acid notation)，使用下面表1的字母表示简并碱基，例如字母M表示A和/或C。
表1
| 字母 | 所代表的碱基 |
| M | A/C |
| K | G/T |
| R | A/G |
| Y | C/T |
| W | A/T |
| S | C/G |
| B | C/G/T |
| D | A/G/T |
| H | A/C/T |
| V | A/C/G |
根据优选的实施方式，本发明使用的基因测序芯片表面具有彼此分隔的、阵列排布或非阵列排布的微反应室，或称为微坑、微井或微阱。对于绝大多数微坑，其内部都加载了待测核酸文库中的一条序列，少数微坑没有加载序列或加载超过一条待测序列。核酸分子与微坑之间优选的通过载体进行连接，例如，水凝胶微球或磁性微球等。在进行测序反应的过程中，通过油封将微坑封住，保证每轮测序反应每个微坑产生的荧光分子封闭在对应的微坑内，而不会游离出来。在进行测序信号采集的过程中，发生延伸反应的微反应室在采集到的原始测序图像上显示为一个个亮点，而未发生反应的微反应室则显示为暗点或暗区域。
提取获得初步信号：在得到原始测序图像后，首先判断图像是否符合标准，将不符合对焦标准、旋转角度过大或偏移量过大的图片去除，保留符合上述标准的图片。再对上述图片进行初步信号提取，主要包括：扣除芯片背景，识别微反应室位置，通过计算每个微反应室的灰度值总和得到微反应室强度，区分有测序反应的亮室和无测序反应的暗室，通过暗室强度计算荧光本底，亮室强度扣除所述荧光本底得到初步信号，一个亮室对应一个亮点。
对于衰减校正，现有的技术方案中，以Ion Torrent为例，其采用的是基于单个测序点的一倍信号分段拟合衰减的方法，这种方法的弊端是比较明显的，因为单个数据点能提供的数据量是十分有限的，而且一旦出现噪音，产生的影响是巨大的，总体而言这种方法的信噪比十分受限，会严重降低校正的准确度。与基于单个测序点分段拟合的校正方法相比，本发明采用了基于群体信号(即很多个测序点)进行衰减校正的方法，以测序图像中发生延伸反应的初步信号作为一个群体计算其信号的平均值，作为该轮测序反应的初始群体信号强度，能够有效抵抗噪音的干扰，信噪比高。此群体中包含成千上万个反应信号，每轮反应中延伸的碱基数量总体来看是随机的，对于理想的随机序列群体，在进行2+2简并测序时，群体中的1/2延伸1个碱基，1/4延伸2个碱基，1/8延伸3个碱基，1/2n延伸n个碱基。每个单独序列每轮延伸的碱基数是随机的，但整体的每轮反应平均延伸碱基数是2个。
测序中每个测序点上DNA模板分子数的损失速率是基本恒定的，将每轮衰减的比例记为1-d，即第n+1轮剩余的有效模板分子数是第n轮的d倍。在整个测序过程中模板分子数的变化符合指数衰减，第n轮剩余的有效模板分子数是第1轮的dn-1倍。去除校正序列区对应的测序反应的初始群体信号强度，将其余每一轮测序反应对应的初始群体信号强度进行指数衰减拟合，得到的拟合曲线为：
其中c是测序反应轮数，
为第c轮测序反应的初始群体信号强度，
为扣除h的第1轮测序反应的初始群体信号强度，d是拟合所得的群体衰减常数，h是零倍信号校正常数，反映了实际信号与理想线性信号的少量截距偏差，lkey是校正序列区对应的测序反应轮数。
单位信号指每个测序点上延伸一个碱基的信号强度，是一个正比于该测序点模板DNA分子数的物理量，测序信号的强度即等于单位信号乘延伸反应的碱基个数。无论对任何测序技术来说，单位1的测量精度都直接影响到测序结果的准确度。因此，单位信号的计算方法对于任何测序方法都意义重大。
归一化是一种简化计算的方式，还可以保证程序运行时收敛加快。归一化方法有两种形式，一种是把数变为(0，1)之间的小数，一种是把有量纲表达式变为无量纲表达式。本发明中的归一化指的是将每轮测序反应释放的荧光信号强度值变为延伸的碱基数量，即把有量纲的表达式变为无量纲表达式。在得到原始测序图像后，对图像中的亮点提取获得初步信号，例如，一个反应信号在某个测序循环中的亮度值为1000，在下一个测序循环的亮度值为2000(可以理解的，计算过程中需要排除亮度值异常的点)，亮度的绝对值并不是我们关注的重点，亮度之间的倍数关系才是重点，比如，若亮度为1000时表示延伸了1个碱基，则亮度为2000时表示延伸了2个碱基。因此通过归一化处理便可以得到每轮测序反应延伸的碱基个数。因此，归一化处理的关键是得到单位信号。
根据优选的实施方式，本发明中的测序文库的组成有特殊要求，所述文库包括：至少一个测序引物结合位点、0个或至少一个索引序列、校正序列区、待测基因。在测序时，先于待测基因序列测得的几个碱基是校正序列区，其碱基组成和数量是已知的，例如，校正序列区可以是AGTAGG，以MK流程为例，每一轮反应所延伸的碱基分别为(A,GT,A,GG)，即进行4轮测序反应，对应的DPL为(1,2,1,2)。由于校正序列区至少有2轮的反应碱基数是设计好的已知数，因此用各反应位点在这几轮的校正后信号值除以反应的碱基数即可计算得到各反应位点的单位信号。
根据优选的实施方式，本发明中求得单位信号的关键是校正序列区，校正序列区是在待测基因序列之前测得的一小段设计好的已知序列，由于其序列是已知的，我们可以得到其归一化信号的理论值(即DPL)，而实际测得的DPM与理论值之间的差异需要被校正，以去除衰减的影响。利用2+2简并测序信号的特点，随机序列的信号均值等于2倍信号的均值，因此所述拟合曲线也是2倍信号均值曲线，对序列已知的校正序列区信号可以进行如下平移变换将其均值调整到符合2倍信号曲线的相应位置：
其中c是测序反应轮数，sa(c)为第c轮测序反应的调整后信号强度，s(c)为第c轮测序反应的初始群体信号强度，n(c)为第c轮测序反应的校正序列区对应的碱基个数，
为扣除h的第1轮测序反应的群体信号强度，d是拟合所得的群体衰减常数，h是零倍信号校正常数，反映了实际信号与理想线性信号的少量截距偏差，lkey是校正序列区对应的测序反应轮数。
根据优选的实施方式，得到调整后的校正序列区信号以后，需要进一步得到单位信号。即算出调整后的校正序列区信号的总和，再除以校正序列区碱基总数，则可以得到每个点的单位信号u，即：
需要注意的是，由于校正序列区最后一轮反应的碱基可能还包含了若干待测基因的碱基，因此为了方便计算，参与最后一轮测序反应的校正序列区的碱基不参与上述单位信号u的计算，即步骤D中所述的部分校正序列区中不包括校正序列区对应的最后一轮测序反应的碱基。此外，第一轮测序反应有可能是空反应，例如，校正序列为GTAGCC，第一轮通入的底物为M，则第一轮为空反应，此时，可以将第一轮反应的数据舍弃。
根据优选的实施方式，上述计算得到了亮点对应的单位信号，可以据此计算出每一轮测序反应中每个亮点的归一化信号值，即每轮反应每个位点延伸了多少个核苷酸，通过以下校正公式计算得到每一轮的每个反应信号的归一化信号值：
其中c是测序反应轮数，sa(c)为第c轮测序反应的调整后信号强度，d是拟合所得的群体衰减常数，c是测序反应轮数，h是零倍信号校正常数，u是对应的单位信号。
根据优选的实施方式，在得到每一轮测序反应每个反应信号的归一化信号值后，以单个测序位点为单位按照测序反应轮数由小到大得到每个测序位点的核酸链的简并碱基序列，得到整个视野所有反应信号的简并序列数据。
根据优选的实施方式，所有的原始测序图像都需要进行本发明所述的衰减校正和归一化处理，如此便可以得到所述待测基因的简并碱基序列数据。
根据优选的实施方式，利用上述归一化方法得到归一化信号值后，进一步包括，对所述归一化信号值进行失相校正，去除超前反应和滞后反应造成的信号偏差。
需要注意的是，凡是能够满足2+2简并测序原理的测序反应都可以适用本发明的方法，例如，对于3’端不封闭的焦磷酸测序，进行简单的调整，即可满足本发明的方法，即在一轮测序反应中加入两种核苷酸底物，在相邻的反应轮中加入另外两种核苷酸底物，交替进行。
实施例1
对待测基因文库进行2+2简并测序反应，使用MK测序流程，即奇数轮测序反应通入简并底物M(包括A和C)，偶数轮反应通入简并底物K(包括G和T)，每一轮测序反应中都对芯片进行拍照，得到多张原始测序图像，进行基本的图像处理和初步信号提取，在测序图像上，显示有未发生碱基延伸的暗坑，和发生碱基延伸的亮坑，对每个测序视野可以得到一个微坑强度值乘以测序轮数的原始信号强度矩阵。筛选出无反应的暗坑信号，取均值，从剩下所有信号值中扣除暗坑本底强度，得到亮微坑的信号强度矩阵，如表2所示，每一行为一个亮坑(亮点)的信号，每一列为一轮。
表2
| 70 | 886 | 498 | 535 | 470 | 494 | 495 | 528 | 489 | … |
| 80 | 897 | 492 | 519 | 463 | 1730 | 484 | 2938 | 888 | … |
| 283 | 738 | 565 | 459 | 443 | 466 | 1381 | 992 | 1226 | … |
| 61 | 824 | 461 | 513 | 455 | 473 | 453 | 815 | 444 | … |
| 95 | 926 | 534 | 623 | 513 | 546 | 950 | 1773 | 920 | … |
| … | … | … | … | … | … | … | … | … | … |
根据校正序列区TTCGCT和测序的底物进样顺序MK，可知校正序列区的测序信号为前6轮反应，信号倍数依次为[0,2,1,1,1,1+]。从第7轮起为待测随机序列信号。
对每一轮信号求均值，得到均值信号，如图2所示。前6轮信号均值变化趋势符合校正序列区信号特点，第7轮到最后的测序信号均值变化趋势符合信号衰减的指数模型。
使用除校正序列区外的轮次的信号均值进行单指数拟合，得到拟合曲线方程：
y＝849.1×0.996067x-1-32.6
拟合曲线如图3所示，图中离散分布的叉号表示的是校正序列区6轮测序反应对应的信号均值，黑点为待测基因序列信号均值。
在校正序列区，将拟合曲线上的值视为2倍信号均值，并据此调整校正序列区的信号。表达为以下公式：
其中，其中
为第c轮的信号均值，
为扣除h的第1轮群体信号的均值，d是拟合所得的群体衰减常数，c是信号的轮数，h是零倍信号校正常数，反映了实际信号与理想线性信号的少量截距偏差，lkey为校正序列区的测序轮数，sa(c)为第c轮的调整后信号值，s(c)为第c轮原始信号值，n(c)为第c轮反应的校正序列碱基个数。
例如第2轮反应，c＝2。反应2个碱基为2倍信号，即n(c)＝2。该轮信号均值即上述矩阵的第2列均值
由上述曲线拟合结果得到其余参数，
d＝0.996067，h＝-32.6。全部代入计算可得
sa(c)＝s(c)+64.2
即调整后的第2轮信号都加64.2。
校正序列区校正完毕得到新的矩阵：
| -80 | 845 | 377 | 392 | 387 | 494 | 495 | 528 | 489 | … |
| -70 | 856 | 371 | 376 | 380 | 1730 | 484 | 2938 | 888 | … |
| 133 | 697 | 444 | 316 | 360 | 466 | 1381 | 992 | 1226 | … |
| -89 | 783 | 340 | 370 | 372 | 473 | 453 | 815 | 444 | … |
| -55 | 885 | 413 | 480 | 430 | 546 | 950 | 1773 | 920 | … |
| … | … | … | … | … | … | … | … | … | … |
根据此校正后矩阵计算每个亮点的单位信号，最开始不反应的第1轮不参与计算，第6轮包含了校正序列区的最后一个碱基和若干待测碱基也不参与计算。即每行第2-5轮信号减h求和，除以校正序列这几轮反应的碱基总数5，得到每行一个单位信号。例如第一行的单位信号u为
最后根据公式
计算归一化信号值。其中，sa(c)为上述矩阵中第c轮的信号，u为该行的单位信号，计算得到最终的归一化信号矩阵：
| -0.11 | 2.07 | 0.97 | 1.01 | 1 | 1.26 | 1.27 | 1.35 | 1.26 | … |
| -0.09 | 2.11 | 0.96 | 0.98 | 0.99 | 4.25 | 1.25 | 7.22 | 2.25 | … |
| 0.43 | 1.88 | 1.23 | 0.91 | 1.02 | 1.31 | 3.72 | 2.7 | 3.33 | … |
| -0.14 | 2.05 | 0.94 | 1.02 | 1.03 | 1.29 | 1.25 | 2.18 | 1.23 | … |
| -0.05 | 1.97 | 0.96 | 1.11 | 1 | 1.26 | 2.15 | 3.97 | 2.1 | … |
| … | … | … | … | … | … | … | … | … | … |
此矩阵中的每一行数据即为校正了单位信号和衰减的一个亮坑的每一轮反应的信号。图4所示为一个亮坑每一轮测序反应的归一化信号图，横坐标为测序轮数，纵坐标为归一化信号值，可以看出在经过衰减校正和归一化处理后，大多数归一化信号值非常接近整数，有利于准确读出2+2测序每一轮反应延伸的碱基个数。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。
Claims (11)
1.一种测序信号归一化的方法，其特征在于，包括以下步骤：
A.利用基因测序芯片对由待测基因构建的文库进行多轮测序反应，得到多幅原始测序图像；所述文库中校正序列区紧邻待测基因，所述校正序列区的核苷酸序列已知；一轮测序反应指的是，提供一种测序试剂，将该种测序试剂中具有可检测标记的核苷酸单体掺入所述文库之后检测所述可检测标记生成的荧光信号的过程；
B.在每轮测序反应所得的所述原始测序图像中，对图像中的亮点提取获得初步信号，计算初步信号的平均值作为该轮测序反应的初始群体信号强度；所述亮点对应于芯片上发生核苷酸单体掺入的微反应室；
C.利用除去所述校正序列区之后的每一轮测序反应对应的所述初始群体信号强度进行指数衰减拟合，得到拟合曲线和相关参数；
D.根据所述拟合曲线和相关参数，计算得到所述校正序列区的调整后的信号强度，根据调整后的部分校正序列区信号总和及其对应的核苷酸总数，得到所述亮点处的单位信号；
E.结合单位信号和所述拟合曲线得到所述亮点在每一轮测序反应的归一化信号值。
2.根据权利要求1所述的方法，其特征在于，所述测序为3’端不封闭的测序反应，在测序反应中，包括两种不同的测序试剂：第一测序试剂和第二测序试剂；两种测序试剂循环加入；其中所述第一测序试剂包含具有可检测标记的两种不同的核苷酸单体；所述第二测序试剂包含具有可检测标记的两种核苷酸单体，且所述两种核苷酸单体不同于所述第一测序试剂中存在的所述核苷酸单体，并且其中所述第二测序试剂是在提供了所述第一测序试剂随后提供的，将所述核苷酸单体掺入所述文库之后检测所述可检测标记生成的荧光信号。
3.根据权利要求1或2所述的方法，其特征在于，所述核苷酸单体为5’端多磷酸修饰有荧光切换性质的荧光团的核苷酸，所述的荧光切换性质指的是测序反应后荧光信号相比测序反应前有明显改变。
4.根据权利要求1所述的方法，其特征在于，所述文库包括：至少一个测序引物结合位点、0个或至少一个索引序列、校正序列区、待测基因。
5.根据权利要求1或4所述的方法，其特征在于，所述校正序列区的核苷酸序列已知，组成所述校正序列区的核苷酸个数为6-10个，所述校正序列区对应的测序反应轮数大于2。
6.根据权利要求1所述的方法，其特征在于，所述提取获得初步信号包括：扣除芯片背景，识别微反应室位置，通过计算每个微反应室的灰度值总和得到微反应室强度，区分有核苷酸单体掺入的亮室和无核苷酸单体掺入的暗室，通过暗室强度计算荧光本底，亮室强度扣除所述荧光本底得到初步信号。
7.根据权利要求1所述的方法，其特征在于，所述拟合曲线为：
其中c是测序反应轮数，
为第c轮测序反应的初始群体信号强度，
为扣除h的第1轮测序反应的初始群体信号强度，d是拟合所得的群体衰减常数，h是零倍信号校正常数，反映了实际信号与理想线性信号的少量截距偏差，lkey是校正序列区对应的测序反应轮数。
8.根据权利要求1所述的方法，其特征在于，按照以下公式对校正序列区的初始群体信号强度进行调整：
其中c是测序反应轮数，sa(c)为第c轮测序反应的调整后信号强度，s(c)为第c轮测序反应的初始群体信号强度，n(c)为第c轮测序反应的校正序列区对应的碱基个数，
为扣除h的第1轮测序反应的群体信号强度，d是拟合所得的群体衰减常数，h是零倍信号校正常数，反映了实际信号与理想线性信号的少量截距偏差，lkey是校正序列区对应的测序反应轮数。
9.根据权利要求1所述的方法，其特征在于，步骤D中所述的部分校正序列区中不包括校正序列区对应的最后一轮测序反应的碱基。
10.根据权利要求1所述的方法，其特征在于，所述归一化信号值可通过以下校正公式进行求解：
其中c是测序反应轮数，sa(c)为第c轮测序反应的调整后信号强度，d是拟合所得的群体衰减常数，c是测序反应轮数，h是零倍信号校正常数，u是对应的单位信号。
11.根据权利要求1-10任一项所述的方法，其特征在于，所述方法进一步包括，在得到所述归一化信号值后，对所述归一化信号值进行失相校正，去除超前反应和滞后反应造成的信号偏差。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210103665.XA CN114507723B (zh) | 2022-01-28 | 2022-01-28 | 一种测序信号归一化的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210103665.XA CN114507723B (zh) | 2022-01-28 | 2022-01-28 | 一种测序信号归一化的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114507723A true CN114507723A (zh) | 2022-05-17 |
| CN114507723B CN114507723B (zh) | 2024-07-23 |
Family
ID=81548996
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210103665.XA Active CN114507723B (zh) | 2022-01-28 | 2022-01-28 | 一种测序信号归一化的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114507723B (zh) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080020378A1 (en) * | 2004-09-25 | 2008-01-24 | Yaodong Chen | Fluorescent base analogues' usage in the characterization of nucleic acid molecules and their interactions |
| CN107958138A (zh) * | 2016-10-14 | 2018-04-24 | 北京大学 | 一种从高通量dna测序的原始信号中读取序列信息的方法 |
| CN108699599A (zh) * | 2015-11-19 | 2018-10-23 | 北京大学 | 获得和校正生物序列信息的方法 |
| CN113257351A (zh) * | 2020-02-12 | 2021-08-13 | 赛纳生物科技(北京)有限公司 | 一种用于多碱基基因测序的基因文库及其构建方法 |
-
2022
- 2022-01-28 CN CN202210103665.XA patent/CN114507723B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080020378A1 (en) * | 2004-09-25 | 2008-01-24 | Yaodong Chen | Fluorescent base analogues' usage in the characterization of nucleic acid molecules and their interactions |
| CN108699599A (zh) * | 2015-11-19 | 2018-10-23 | 北京大学 | 获得和校正生物序列信息的方法 |
| CN107958138A (zh) * | 2016-10-14 | 2018-04-24 | 北京大学 | 一种从高通量dna测序的原始信号中读取序列信息的方法 |
| CN113257351A (zh) * | 2020-02-12 | 2021-08-13 | 赛纳生物科技(北京)有限公司 | 一种用于多碱基基因测序的基因文库及其构建方法 |
Non-Patent Citations (1)
| Title |
|---|
| 黄国亮等: "微阵列芯片的荧光光漂白特性", 《发光学报》, vol. 27, no. 2, pages 259 - 264 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114507723B (zh) | 2024-07-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CA2689626C (en) | Methods and processes for calling bases in sequence by incorporation methods | |
| US7133782B2 (en) | Sequencing method and apparatus | |
| CN116083547A (zh) | 一种校正测序期间超前量的方法 | |
| CN101633961B (zh) | 循环“连接－延伸”基因组测序法 | |
| CN113012757B (zh) | 识别核酸中的碱基的方法和系统 | |
| CN108504649B (zh) | 编码pcr二代测序建库方法、试剂盒及检测方法 | |
| US20090105961A1 (en) | Methods of nucleic acid identification in large-scale sequencing | |
| US11866778B2 (en) | Methods and systems for evaluating microsatellite instability status | |
| US20220389495A1 (en) | Methods and systems for multiplex analysis | |
| CN114507723A (zh) | 一种测序信号归一化的方法 | |
| EP1415279A2 (en) | Image metrics in the statistical analysis of dna microarray data | |
| US10614571B2 (en) | Object classification in digital images | |
| CN116246703A (zh) | 一种核酸测序数据的质量评估方法 | |
| CN114420214A (zh) | 核酸测序数据的质量评估方法和筛选方法 | |
| WO2006119996A1 (en) | Method of normalizing gene expression data | |
| WO2018235938A1 (ja) | 核酸をシークエンシングする方法および解析する方法 | |
| CN116426623A (zh) | 一种校正核酸测序产生的序列信息误差的方法 | |
| US20240177807A1 (en) | Cluster segmentation and conditional base calling | |
| US20240212791A1 (en) | Context-dependent base calling | |
| CN116825189A (zh) | 一种2+2测序的信号校正方法 | |
| CN116814753A (zh) | 一种测序信号归一化的方法 | |
| CN110241209B (zh) | 一种引物、试剂盒及用途 | |
| CN117912550A (zh) | 确定读段质量的方法和测序方法 | |
| CN117987522A (zh) | 测序方法 | |
| JP2006191844A (ja) | 標的核酸含有量推定方法、プログラム、並びにシステム |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant |