超低频突变分子标签聚类分群算法
【专利摘要】本发明公开了一种对测序读段进行聚类的方法，所述测序读段携带标签序列，该方法包括：(1)将多个测序读段与参考序列进行比对，并确定各测序读段的两端位置，将两端位置一致的测序读段归类至相同的一级群；(2)对属于同一个一级群的测序读段根据其标签序列进一步分二级群，将分子标签序列相似的测序读段分为同一个二级群。通过该方法能准确有效地对测序读段进行聚类分群，为后期通过各个群的一致性序列来精确检测低频突变奠定坚实的基础。
【专利说明】
超低频突变分子标签聚类分群算法
技术领域
[0001] 本发明设及测序技术领域，特别是超低频突变分子标签聚类分群算法，具体地，本 发明设及对测序读段进行聚类的方法。
【背景技术】
[0002] 随着二代测序的迅速发展，测序费用的降低，二代测序在各个方面的检测研究中 得到了越来越广泛的应用。而相对于全基因组测序，目标区间测序能大幅度降低测序成本 和数据的复杂性，使我们感兴趣的目标区间在较低的成本的同时达到很高的测序覆盖度， 运使得检测癌症突变中的低频突变成为了可能。
[0003] 目标区间测序方法中，采用特异性引物对目标区间进行PCR扩增的方法由于其操 作简单、快速，且只需少量DNA等优点，已被人们广泛应用。然而，特异性引物扩增测序中，不 可避免会存在严重的扩增偏好性，同时也存在扩增测序引入的各种错误。运些问题一方面 直接影响定量的准确性，因为测序数据中的数量已不能代表原始DNA片段的数量;另一方面 会影响分析结果的准确性，引入大量的假阳性。而在肿瘤突变研究中，由于肿瘤的高异质 性，存在大量的低频突变，使得运些问题尤为突出。
[0004] 因而，目前的特异性引物扩增测序仍有待改进。
【发明内容】
[0005] 本发明旨在至少解决现有技术中存在的技术问题之一。为此，本发明的一个目的 在于提出一种对测序读段进行聚类的方法，从而实现对DNA分子精确的定量，同时为后期利 用一致性序列进行精确的超低频突变检测奠定坚实的基础。
[0006] 需要说明的是，本发明是基于发明人的下列工作而完成的：
[0007] 现阶段，针对特异性引物扩增测序的上述问题，研究者引入了分子标签，在原始 DNA分子上连接一段能代表该DNA分子的unique标签序列。不同的DNA分子连接不同的分子 标签，通过分子标签序列可W准确的识另化NA分子。分子标签的引入，可W对DNA分子和突变 进行准确的定量，同时也可W降低甚至消除由扩增和测序等造成的错误。
[000引针对添加分子标签的二代测序数据，在数据处理时，需要根据其分子标签将reads 进行分群，将reads起止位置一样，且分子标签也一样的reads分为一群，认为运是由同一个 DNA分子片段通过PCR扩增生成的多个复本。然后针对每个群，找到其最终的一致性序列(在 本文中，有时也将"一致性序列"称为"共有序列"），即是该群所对应的原始DNA分子的序列。 最后，再利用运些一致性序列进行后续的突变检测等分析。
[0009]然而，由于实验中对添加分子标签后的分子模板进行PCR扩增，同一个分子模板会 产生一群一模一样的子分子;但在实验测序过程中，又不可避免引入一些错误，最后得到一 些含有少量错误的分子模板被多次重复测序的化S化数据。本发明即是针对运种情况，致力 于根据分子标签和read(测序读段）的自身序列（与基因组的比对位置），在考虑测序错误的 前提下，把来源于同一个分子模板的reads进行聚类分群，W便后续分析。
[0010] 进而，在本发明的第一方面，本发明提供了一种对测序读段进行聚类的方法，所述 测序读段携带标签序列。根据本发明的实施例，所述方法包括：
[0011] (1)将多个测序读段与参考序列进行比对，并确定各测序读段两端的位置，将两端 位置一致的测序读段归类至相同的一级群；
[0012] (2)对属于同一个一级群的测序读段根据其标签序列进一步分二级群，将分子标 签序列相似的测序读段分为同一个二级群。
[0013] 根据本发明的实施例，所述步骤(2)的详细步骤包括：
[0014] (a)确定所述一级群内的各标签的深度；
[0015] (b)将所述各标签按深度从高到低进行排序；
[0016] (C)针对深度由高至低的标签依次实施下列步骤：
[0017] 如果所述标签与已有的种子标签序列的错配不超过指定错配数，则将具有所述标 签的测序读段分配至所述种子标签子群中；
[0018] 如果所述标签与已有的种子标签序列的错配超过指定错配数，则选择所述标签为 新的种子标签，并将具有所述标签的测序读段分配至相应的种子标签子群中；
[0019] 经过上述二级群处理后，所有的测序读段都分成了若干个二级群，运些二级群即 最后的分群结果。
[0020] 发明人惊奇地发现，通过该方法能准确有效地对测序读段进行聚类分群，为后期 通过各个群的一致性序列来精确检测低频突变奠定坚实的基础。
[0021] 根据本发明的实施例，（C)中所述种子标签是指该二级群的深度最高的标签序列， 可W认为是该群的真实的标签序列，同时该群中存在一些深度较低的含有错误的标签序 列。由此，测序读段的聚类分群结果可靠，后续测序分析结果准确。
[0022] 根据本发明的实施例，在(C)中，依据所采用的测序平台确定指定错配数，其中，当 采用Il Iumina测序平台时，由于Illumina测序平台主要^mismatch (错配数)为主要的测序 错误，所WSbp的分子标签容1个mismatch,也即所述指定错配数为1。由此，聚类分群结果可 靠，后续测序分析结果准确。
[0023] 本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变 得明显，或通过本发明的实践了解到。
【附图说明】
[0024] 本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得 明显和容易理解，其中：
[0025] 图1显示了根据本发明实施例的对测序读段进行聚类的方法的流程示意图。
【具体实施方式】
[0026] 下面详细描述本发明的实施例。下面描述的实施例是示例性的，仅用于解释本发 明，而不能理解为对本发明的限制。
[0027] 在本发明的第一方面，本发明提供了一种对测序读段进行聚类的方法，所述测序 读段携带标签序列。
[0028] 根据本发明的实施例，参照图1，所述方法包括：
[0029] (I)将多个测序读段与参考序列进行比对，并确定各测序读段的两端位置，将两端 位置一致的测序读段归类至相同的一级群；
[0030] (2)对属于同一个一级群的测序读段根据其标签序列进一步分二级群，将分子标 签序列相似的测序读段分为同一个二级群，其具体步骤为：
[0031] (a)确定所述一级群内的各标签的深度；
[0032] (b)将所述各标签按深度从高到低进行排序；
[0033] (C)针对深度由高至低的标签依次实施下列步骤：
[0034] 如果所述标签与已有的种子标签序列的错配不超过指定错配数，则将具有所述标 签的测序读段分配至所述种子标签子群中；
[0035] 如果所述标签与已有的种子标签序列的错配超过指定错配数，则选择所述标签为 新的种子标签，并将具有所述标签的测序读段分配至相应的种子标签子群中；
[0036] 经过上述二级群处理后，所有的测序读段都分成了若干个二级群，运些二级群即 最后的分群结果。
[0037] 根据本发明的实施例，（C)中所述种子标签是指该二级群的深度最高的标签序列， 认为是该群的真实的标签序列，同时该群中存在一些深度较低的含有错误的标签序列。由 此，测序读段的聚类分群结果可靠，后续测序分析结果准确。
[0038] 根据本发明的实施例，在(C)中，依据所采用的测序平台确定指定错配数，其中，当 采用Il Iumina测序平台时，由于Illumina测序平台主要^mismatch (错配数)为主要的测序 错误，所WSbp的分子标签容1个mismatch,也即所述指定错配数为1。由此，聚类分群结果可 靠，后续测序分析结果准确。
[0039] 由此，测序读段的聚类分群结果可靠，从而实现对DNA分子精确的定量，同时为后 期利用一致性序列进行精确的超低频突变检测奠定坚实的基础。
[0040] 下面将结合实施例对本发明的方案进行解释。本领域技术人员将会理解，下面的 实施例仅用于说明本发明，而不应视为限定本发明的范围。实施例中未注明具体技术或条 件的，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。所用试剂或仪 器未注明生产厂商者，均为可W通过市购获得的常规产品，例如可W采购自Ilhimina公司。 [0041 ] 实施例1:
[0042] 本实施例针对已知8个突变位点（如下表1所示）的突变频率为1%的样本(人类）， 采用Sbp随机分子标签对DNA分子进行标记，然后，采用Ampli化q(;ol(峽360Master Mix对样 本进行针对各已知突变位点的特异性引物扩增，最后利用11 Iumina NS500测序平台对各扩 增产物进行75PE测序。
[0043] 然后，根据本发明的对测序读段进行聚类的方法，参照图1，按照W下步骤对获得 的测序读段进行聚类，获得测序读段群：
[0044] (1)将5,475,216个测序读段与UCSC数据库中人类参考基因组化gl9)进行比对，并 确定各测序读段的两端位置，将两端位置一致的测序读段归类至相同的一级群，得到共 25540个一级群。
[0045] (2)对属于同一个一级群的测序读段根据其标签序列进一步分二级群，将分子标 签序列相似的测序读段分为同一个二级群，具体步骤如下：
[0046] (a)确定所述一级群内的各标签的深度；
[0047] (b)将所述各标签按深度从高到低进行排序；
[0048] (C)针对深度由高至低的标签依次实施下列步骤：
[0049] 如果所述分子标签(Sbp)与已有的种子标签序列的错配不超过1个，则将具有所述 标签的测序读段分配至所述种子标签子群中；
[0050] 如果所述标签与已有的种子标签序列的错配超过1个，则选择所述标签为新的种 子标签，并将具有所述标签的测序读段分配至相应的种子标签子群中；
[0051] 经过上述二级群处理后，所有的测序读段都分成了 71187个二级群，运些二级群即 最后的分群结果。
[0052] 其中，本实施例中所用参考基因组来源于UCSC数据库的人类基因组化gl9版本）， 网址：http: / Agdownload. cse. UCSC. edu/goldenPath/hgl9/chromosomes/。
[0053] 本实施例中分子标签序列为8bp，标签与种子标签序列的错配数选为1，也可W根 据情况调整。
[0054] 获得71187个二级群后，对运些二级群进行过滤和处理，W分别确定各测序读段群 的共有序列的方法，具体步骤如下：
[0化5] 1、过滤；
[0056] 测序读段(read)聚类分群得到测序读段群(reads groups)后，对运些测序读段群 按照W下条件进行过滤：
[0057] a)对双端比到不同染色体的read groups进行过滤；
[005引 b)对插入片段大小<30,或MOO的read groups进行过滤；
[0化9]由于CfDNA的片段大小主要在16化P和330bp左右，所W插入片段大小最大不应超 过4(K)bp;而扩增引物的长度一般为20多bp，故插入片段大小最小不应小于30bp。
[0060] C)对read的起始位置不在扩增引物起始位置的read groups进行过滤；
[0061] 由于是扩增引物的扩增产物，read的起始位置应该是引物的起始位置。
[0062] 2、确定共有序列(有时也称为乂 onsensus序列"）
[0063] 基本原理：
[0064] 每个测序读段群中的reads是同一个分子模板产生的，所W原则上同一个group中 的reads应该序列一样，且barcode-样;但是由于在实验和测序过程中，不可避免存在一些 错误，group中的reads会有一些错误。而确定Consensus序列的过程，就是排除运些错误，得 到分子模板的真实序列。
[00化]处理步骤：
[0066] a)针对read各个位置，进行W下操作：
[0067] i.统计ATCG 4种碱基各自的深度；
[006引 ii .对ATCG4种碱基的深度从高到低排序，得到max、sec、third、fou;rth
[0069] iii .计算系数C= (max-sec)/max,若该系数0 = 0.65,贝认为max深度的碱基即为 该位置Consensus碱基，而该Consensus碱基的质量为Q = 20+(ma巧C~ 2)/2，当QMO时，取40; 若C<0.65，贝认为read运个位置的碱基不确定，Consensus序列该位置为N，相应质量值Q = 2。
[0070] 对read各个碱基进行运些操作后，得到该group的Consensus序列W及对应的质量 值;但Consensus序列中可能有一些碱基不确定，为N。
[0071] b)若整个read中不确定的碱基数超过5，则过滤该group;若不超过5，则进行下一 步(C)判断；
[0072] C)统计该gro叩中barcode (即分子标签）的深度，同上方法，判断该gro叩中 barcode是否能确定；若不确定，则过滤该group;若确定，该group保留，且最终的Consensus 序列、相应质量值，W及其barcode序列都已获得。
[0073] 由此，最终得到10970条一致性序列。
[0074] 然后利用获得的一致性序列进行突变检测，检测结果如下表1:
[0075] 表 1 [00761
[0077] 表1中第一列是染色体编号，第二列是突变位点在染色体上的位置，第=列是基因 名，第四列是基因在染色体上的方向，第五列是具体的CDS和蛋白突变信息，第六列是突变 频率，第屯列是该实验的检测结果(YES是检测到，NO是未检测到）。
[0078] 由上述结果可知，本实施案例采用添加分子标签的技术，结合分子标签聚类分群， 在仅约5M reads的测序情况下，成功精确地检测到了所有的突变频率仅1 %的突变。
[0079] 实施例2:
[0080] 本实施例针对已知8个突变位点（如下表2所示）的突变频率为0.1 %的样本（人 类），采用Sbp随机分子标签对DNA分子进行标记，然后，采用Amplihq.如Id壞360Master Mix 对样本进行针对各已知突变位点的特异性引物扩增，最后利用11 Iumina NS500测序平台对 各扩增产物进行75PE测序。
[0081] 然后，根据本发明的对测序读段进行聚类的方法，参照图1，按照W下步骤对获得 的测序读段进行聚类，获得测序读段群：
[0082] (1)将5,328,887个测序读段与UCSC数据库中人类参考基因组化gl9)进行比对，并 确定各测序读段的两端位置，将两端位置一致的测序读段归类至相同的一级群，得到共 25634个一级群。
[0083] (2)对属于同一个一级群的测序读段根据其标签序列进一步分二级群，将分子标 签序列相似的测序读段分为同一个二级群，具体步骤如下：
[0084] (a)确定所述一级群内的各标签的深度；
[0085] (b)将所述各标签按深度从高到低进行排序；
[0086] (C)针对深度由高至低的标签依次实施下列步骤：
[0087] 如果所述分子标签(8bp)与已有的种子标签序列的错配不超过I个，则将具有所述 标签的测序读段分配至所述种子标签子群中；
[0088] 如果所述标签与已有的种子标签序列的错配超过1个，则选择所述标签为新的种 子标签，并将具有所述标签的测序读段分配至相应的种子标签子群中；
[0089] 经过上述二级群处理后，所有的测序读段都分成了 61557个二级群。
[0090] 其中，本实施例中所用参考基因组来源于UCSC数据库的人类基因组化gl9版本）， 网址：http: / Agdownload. cse. UCSC. edu/goldenPath/hgl9/chromosomes/。
[0091] 本实施例中分子标签序列为8bp，标签与种子标签序列的错配数选为1，也可W根 据情况调整。
[0092] 获得61557个二级群后，对运些二级群进行过滤和处理，W确定各测序读段群的共 有序列，具体方法步骤如实施例1。由此，最终得到10584条一致性序列。
[0093] 然后利用获得的一致性序列进行突变检测，检测结果如下表2:
[0094] 表 2
[0095]
[0
[0097] 表2中第一列是染色体编号，第二列是突变位点在染色体上的位置，第=列是基因 名，第四列是基因在染色体上的方向，第五列是具体的CDS和蛋白突变信息，第六列是突变 频率，第屯列是该实验的检测结果(YES是检测到，NO是未检测到）。
[0098] 由上述结果可知，本实施案例采用添加分子标签的技术，结合分子标签聚类分群， 在仅约5化eads的测序数据量下，成功精确地检测到了6个突变频率低至0.1 %的突变，另2 个突变在提高测序数据量的情况下也能检测到。
[0099] 目前检测低频突变的技术，例如ARMS和Digi化1 PCR等技术才能检测到低至0.1% 的突变，但运些技术存在通量低，成本高，且只能检测已知突变位点的缺点，而普通的二代 测序技术只能检测2%的突变频率。而由上述实施例的结果可知，本发明在添加分子标签的 技术基础上，结合分子标签聚类分群方法，对测序数据进行分析，即克服了 ARMS和Digital PCR等技术的缺点，同时又成功检测到了突变频率低至0.1 %的突变。
[0100] 在本说明书的描述中，参考术语"一个实施例"、"一些实施例"、"示例"、"具体示 例"、或"一些示例"等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特 点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不 一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可W在任何 的一个或多个实施例或示例中W合适的方式结合。
[0101]尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可W理解:在不 脱离本发明的原理和宗旨的情况下可W对运些实施例进行多种变化、修改、替换和变型，本 发明的范围由权利要求及其等同物限定。
【主权项】
1. 一种对测序读段进行聚类的方法，所述测序读段携带标签序列，其特征在于，所述方 法包括： (1) 将多个测序读段与参考序列进行比对，并确定各测序读段两端的位置，将两端位置 一致的测序读段归类至相同的一级群； (2) 对属于同一个一级群的测序读段根据其标签序列进一步分二级群，将分子标签序 列相似的测序读段分为同一个二级群。2. 根据权利要求1所述的方法，其特征在于，所述步骤(2)的详细步骤包括： (a) 确定所述一级群内的各标签的深度； (b) 将所述各标签按深度从高到低进行排序； (c) 针对深度由高至低的标签依次实施下列步骤： 如果所述标签与已有的种子标签序列的错配不超过指定错配数，则将具有所述标签的 测序读段分配至所述种子标签子群中； 如果所述标签与已有的种子标签序列的错配超过指定错配数，则选择所述标签为新的 种子标签，并将具有所述标签的测序读段分配至相应的种子标签子群中； 经过上述二级群处理后，所有的测序读段都分成了若干个二级群，这些二级群即最后 的分群结果。3. 根据权利要求2所述的方法，其特征在于，（c)中所述种子标签是指该二级群的深度 最高的标签序列，认为是该群的真实的标签序列，同时该群中存在一些深度较低的含有错 误的标签序列。
【文档编号】G06F19/24GK106021987SQ201610350317
【公开日】2016年10月12日
【申请日】2016年5月24日
【发明人】曾华萍, 宋卓, 袁梦兮
【申请人】人和未来生物科技（长沙）有限公司