基于核心家系的串联重复变异分型检测方法及其应用与流程

[0001][0002]
本发明涉及生物信息技术领域，具体地，涉及一种基于核心家系的串联重复变异分型检测方法及其应用。


背景技术：

[0003][0004]
人类基因组中的变异和人类演化、疾病风险等方面都有着密切的联系。随着二代短读长高通量测序技术（ngs）的发展，研究者们开发了一系列检测基因组结构变异（structure variations，简称svs）的方法。在不考虑组装的情况下，基于二代测序数据检测变异的方法主要有三种策略：(1) read pair (rp)，即根据pair reads在参考基因组上的映射距离和方向，将每个读取对象分类为normal或sv，然后识别出符合sv类别的读取次数较多的区域，并分配一个置信度得分；(2) split read (sr)，即两条pair-end的reads中只有一条可以比对到参考基因组，另一条却不可以，产生这种分裂reads的位点往往是有结构变异的发生；(3) read depth (rd)，即reads覆盖深度，主要用于序列丢失或者重复的检测。特别地，针对串联重复区域，研究者开发出以下几种软件用于检测其变异信息，包括：
[0005]
(1) tssv(https://pypi.org/project/tssv/)；
[0006]
(2) hipstr（https://github.com/tfwillems/hipstr）；
[0007]
(3) strscan（http://darwin.informatics.indiana.edu/str/）等。
[0008]
然而，由于二代测序读长短，很难覆盖到较长串联重复的整个区域，主要利用断点、深度信息等通过算法进行推断。例如，hipstr使用repeat-aware隐式马尔科夫模型来比对reads，训练一个locus-specific的stutter models，同时利用flanking的snp进行分型。此外，ngs测序在扩增时会有序列偏好性，也会影响tr检测。
[0009]
三代测序将reads长度从二代的150~200 bp增加到15 kb~4 mb，这种长读长reads可以直接覆盖绝大部分的tr区域，提供了准确获取拷贝数信息的测序数据基础，使得三代测序在串联重复变异检测方面比ngs技术有明显优势。不过目前基于三代测序数据开发的tr检测方法并不多，主要有：(1) straglr，不基于已知的tr区域进行检测，主要检测基因组中新发的tr扩增；(2) tricolor，采取无参模式，对全基因组中扩增区域进行检测。这两种方法检测到的变异位点数较少，相较于参考基因组上约一百万个串联重复区域，这两个方法仅能得到几千到几万个区域的tr扩增位点；(3) trgt，对已知的tr区域进行检测。该方法能对参考基因组上所有的串联重复区域进行检测。目前已有的方法主要通过聚类算法得到tr的分型信息，其准确性仍有待评估。此外，还有针对nanopore测序数据开发的tr检测软件，如deeprepeat，nanostr，但仍具有一定的局限性，并非针对pacbio ccs高精准三代测序数据开发的。
[0010]
综上，串联重复区域的变异检测是ngs技术的难点，因其短读长无法跨过重复区域，而且比对时，在串联重复区域，很难确定断点位置和深度信息。长读长测序技术的发展
为串联重复区域的变异检测带来技术上的支持。然而，目前的大多数算法要么存在有效位点数较低，要么无法得到准确的单倍型变异信息的问题。


技术实现要素：

[0011][0012]
针对上述现有技术存在的问题，本发明提供一种基于核心家系的串联重复变异分型检测方法及应用。本发明针对三代测序数据，尤其是pacbio ccs高精准三代测序数据，开发了高准确度、全基因组覆盖的串联重复区域变异检测方法；另一方面，结合分型技术，提出一种基于核心家系的串联重复变异分析技术，能得到遗传过程中串联重复发生变异的单倍型信息。
[0013]
具体来说，本发明涉及如下基于核心家系的串联重复变异分型检测方法和应用。
[0014]
1.一种基于核心家系的串联重复变异分型检测方法，包括：
[0015]
步骤1、基于父本的测序数据集和母本的测序数据集，对子代的测序数据集进行第一分型和组装，得到两套子代单倍型基因组；
[0016]
步骤2、以得到的子代单倍型基因组为参考，分别对父本和母本基因组进行单核苷酸位点变异检测，并对父本的测序数据集和母本的测序数据集分别进行第二分型，得到两套父本单倍型基因组、两套母本单倍型基因组，并明确由父本遗传到子代的单倍型基因组、由母本遗传到子代的单倍型基因组；
[0017]
步骤3、对父本遗传到子代的单倍型基因组、母本遗传到子代的单倍型基因组、两套子代单倍型基因组分别进行串联重复变异检测。
[0018]
2、根据上述的方法，所述测序数据集包括基因组的长reads；优选地，所述长reads为三代测序方法得到的长reads；更优选地，所述三代测序方法选自pacbio和nanopore中的至少一种；更优选地，所述测序数据集为hifi reads。
[0019]
3、根据上述的方法，所述第一分型的方法基于trio_binning方法。
[0020]
4、根据上述的方法，所述第一分型的方法包括基于trio_binning方法的find-unique-kmers分析，得到父本和母本特异的kmer序列。
[0021]
5、根据上述的方法，基于classify_by_kmers方法，并根据父本和母本特意的kmer序列判断子代测序数据集为属于父本的测序数据集、属于母本的测序数据集，或者为未分型的测序数据集。
[0022]
6、根据上述的方法，所述组装的方法选自hifiasm、wtdbg2、canu、nextdenovo中的至少一种。
[0023]
7、根据上述的方法，所述第二分型的方法基于whatshap方法。
[0024]
8、根据上述的方法，所述串联重复变异检测的方法为改进的tandem repeat finder方法。
[0025]
9、根据上述的方法，所述改进包括motif循环特点、检测区段的特征。
[0026]
10、根据上述的方法，所述方法还包括：分析遗传过程中子代的两套单倍型基因组的串联重复变异信息。
[0027]
11、上述的基于核心家系的串联重复变异分型检测方法在生物学中的应用。
[0028]
12、根据上述的应用，所述应用为在疾病病理研究、亲子鉴定、犯罪侦查、种群遗传
学中的应用。
附图说明
[0029][0030]
图1为本发明实施例1提供的基于核心家系的串联重复变异分型检测方法的流程图。
[0031]
图2为本发明实施例1提供的motif长度分布图。
[0032]
图3为实施例1和对比例2的motif长度和频数分布图。
具体实施方式
[0033][0034]
下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。
[0035]
术语定义
[0036]
（1）串联重复（tandem repeats, 简称tr）：dna中的一个或多个核苷酸前后相连接的重复。当重复单元的核苷酸数低于10时，称为短串联重复（str），当重复单元的核苷酸数在10到60之间时，被称为小卫星。当重复单元的拷贝数在人群中存在不同时，被称为可变数目串联重复序列（vntr）。
[0037]
（2）三代测序：三代测序技术，单分子测序技术，不需要经过pcr扩增，实现对每一条dna分子进行单独测序的技术，无gc偏好性，有着更快的数据读取速度。
[0038]
（3）核心家系：由一个孩子（子代）及其父母双亲组成的家系。
[0039]
（4）基因组分型：通过生物检定法检测某一个体的dna序列，并对比参照其他个体的基因型或序列的过程，可用于显示该个体等位基因从其父母遗传而来的情况。
[0040]
（5）基因组组装：基因组组装（genome assembly）是指使用测序方法将待测物种的基因组生成序列片段（即read），并根据reads 之间的重叠区域对片段进行拼接，先拼接成较长的连续序列（contig），再将contigs 拼接成更长的允许包含空白序列（gap）的scaffolds，通过消除scaffolds 的错误和gaps，将这些scaffolds 定位到染色体上，从而得到高质量的全基因组序列。
[0041]
本发明第一方面提供了一种基于核心家系的串联重复变异分型检测方法，包括：
[0042]
步骤1、基于父本的测序数据集和母本的测序数据集，对子代的测序数据集进行第一分型和组装，得到两套子代单倍型基因组；
[0043]
步骤2、以得到的子代单倍型基因组为参考，分别对父本和母本基因组进行单核苷酸位点变异检测，并对父本的测序数据集和母本的测序数据集分别进行第二分型，得到两套父本单倍型基因组、两套母本单倍型基因组，并明确由父本遗传到子代的单倍型基因组、由母本遗传到子代的单倍型基因组；
[0044]
步骤3、对父本遗传到子代的单倍型基因组、母本遗传到子代的单倍型基因组、两套子代单倍型基因组分别进行串联重复变异检测。
[0045]
根据本发明的分型检测方法，优选地，所述测序数据集包括基因组的长reads；优
选地，所述测序数据集包括基因组的三代测序方法得到的长reads；更优选地，所述三代测序方法选自pacbio和nanopore中的至少一种；更优选地，所述测序数据集为hifi reads。根据本发明的方法，平均测序深度优选为15
×
以上。
[0046]
根据本发明的分型检测方法，优选地，所述第一分型的方法基于trio_binning方法。trio_binning首先使用来自两个亲本基因组的高精度短读长数据将子代的长读长序列划分为单倍型特异性的集合，然后每个单倍型独立组装，形成一个完整的二倍体重建。
[0047]
根据本发明的分型检测方法，优选地，所述第一分型的方法包括基于trio_binning方法的find-unique-kmers分析，得到父本和母本特异的kmer序列。更优选地，基于classify_by_kmers方法，并根据父本和母本特意的kmer序列判断子代测序数据集为属于父本的测序数据集、属于母本的测序数据集，或者为未分型的测序数据集。优选地，对于未分型的测序数据集不进行组装。
[0048]
根据本发明的分型检测方法，优选地，所述组装的方法选自hifiasm、wtdbg2、canu、nextdenovo中的至少一种。
[0049]
根据本发明的分型检测方法，优选地，所述第二分型的方法基于whatshap方法。
[0050]
根据本发明的分型检测方法，优选地，所述串联重复变异检测的方法为改进的tandem repeat finder方法。
[0051]
根据本发明的分型检测方法，优选地，所述改进包括motif循环特点、检测区段的特征。其中，motif循环特点包括motif碱基序列等。例如但不限于在考虑atc的motif碱基序列的同时，也会考虑tca和cat的motif循环等。检测区段的特征是指在检测某个区间时，考虑该区段起止点之间整个区段的拷贝数。
[0052]
根据本发明的分型检测方法，优选地，所述方法还包括：分析遗传过程中子代的两套单倍型基因组的串联重复变异信息。计算得到遗传过程中串联重复区域拷贝数的变化，然后进行分析遗传过程中子代的两套单倍型基因组的串联重复变异信息。
[0053]
本发明第二方面提供了上述的基于核心家系的串联重复变异分型检测方法在生物学中的应用。
[0054]
根据本发明的应用，优选地，所述应用为在疾病病理研究、亲子鉴定、犯罪侦查、种群遗传学中的应用。
[0055]
本发明具有以下优点：
[0056]
（1）本发明的方法将组装、分型技术结合到变异检测中，首先从测序数据集中准确地分出了子代的两套单倍型基因组，以及父母亲本遗传给子代的基因组，这一点相较已有算法中通过聚类方法得到的两个单倍型基因型更准确。
[0057]
（2）本发明方法采用改进后的trf算法，综合考量到motif的循环特点、区段不同部分特征不同等情况，检测到的拷贝数更准确。
[0058]
（3）本发明方法针对参考基因组中的所有串联重复区域进行了全面检索，得到了上百万个位点的拷贝数信息，相较已有算法的无参模式中得到的几千到几万个位点，本发明方法更全面地挖掘了全基因组上的串联重复变异信息。
[0059]
本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
[0060]
下面参考具体实施例，对本发明进行说明，需要说明的是，这些实施例仅仅是说明
性的，而不能理解为对本发明的限制。
[0061]
【实施例1】
[0062]
一种基于核心家系的串联重复变异分型检测方法，该方法原理示意图如图1所示。图1中，p1代表子代的测序数据集；long reads代表长读长；p代表来自父本的子代单倍型基因组；m代表来自母本的子代单倍型基因组；fa代表父本测序数据集；mo代表母本测序数据集；hap1为第一套单倍型基因组，也是由亲本遗传给子代的单倍型基因组；hap2为第二套单倍型基因组。
[0063]
步骤1、准备一个核心家系，对父本、母本、子代进行取样，对样本进行三代hifi测序，测序数据量170g，平均测序深度为18.4
×
，分别得到父本的长reads、母本的长reads、子代的长reads（以下简写长reads为reads）；基于父本的reads和母本的reads，使用trio_binning方法的find-unique-kmers分析对子代的reads（图1中的p1）进行分型，得到父本和母本特异的kmer序列，基于classify_by_kmers方法，并根据父本和母本特异的kmer序列判断子代测序reads为属于父本的reads、属于母本的reads，或者为未分型的reads（不组装）。使用hifiasm进行组装，得到两套子代单倍型基因组（图1中的p和m）；
[0064]
步骤2、以得到的子代单倍型基因组为参考，对父本和母本基因组进行单核苷酸位点变异检测，并基于whatshap方法对父本的reads和母本的reads分别进行分型，得到两套父本单倍型基因组（fa的hap1、hap2）、两套母本单倍型基因组（mo的hap1、hap2），并确定由父本遗传到子代的单倍型基因组（例如fa的hap1）、由母本遗传到子代的单倍型基因组（例如mo的hap1）；结果如下表1。
[0065]
表1
[0066]
sampleread_numbase_num(gb)phased_read_rate(%)hap1hap2fa443407362.4498.7227431641634016mo363216053.0699.8022430351381715p1325125254.5686.5814205271394308
[0067]
表1中，sample指样品；read-num指read数；base-num指碱基数；phased_read_rate指分型效率；hap1指第一套单倍型基因组reads数；hap2指第二套单倍型基因组reads数。
[0068]
步骤3、使用改进的tandem repeat finder方法，对父本遗传到子代的单倍型基因组、母本遗传到子代的单倍型基因组、两套子代单倍型基因组分别进行串联重复变异检测。其中，改进的tandem repeat finder方法需考虑motif循环特点（motif碱基序列）、检测区段的特征（在检测某个区间时，考虑该区段起止点之间整个区段的拷贝数），从而得到参考基因组上所有已知重复区间在具体样本中的拷贝数信息。结果共得到1025347个串联重复（tr）位点拷贝信息，这些位点的motif长度分布如图2所示。统计得到每个样本的串联重复分型拷贝数信息，以三个区间为示例，结果如下表2。
[0069]
表2
[0070]
chrstartendmotifmotif lengthcopy numberfathermotherp1
02447-3公开的方法进行串联重复变异分型检测。
[0081]
根据此方法，父本、母本、子代的三个样本中检测到的有效位点个数如表4所示。同时，将本发明实施例检测到的有效位点个数做对比，结果如表4所示。
[0082]
表4. 有效位点个数
[0083]
fa（个）mo（个）p1（个）对比例1295929072962实施例1102534710253471025347
[0084]
将采用本发明方法的实施例1与对比例1进行比较，能够看出，本发明方法能够产生更多有效位点的拷贝数信息，远远多于straglr方法所能提供的信息。
[0085]
【对比例2】
[0086]
一种基于trgt方法的串联重复变异分型检测方法。
[0087]
该方法为pacbio官方推荐的一个方法，其基于已知tr区域进行检测，提供两个分型的拷贝数信息。
[0088]
对p1样本的两款软件的结果进行了比较，统计不同拷贝数差值区间的位点数，结果如表5和图3所示。
[0089]
表5
[0090]
拷贝数差值分布（个）计数（频数）~-1008339-100~-50412-50~-201358-20~-102779-10~-56900-5~01608790~51019205~10555510~20418020~50181750~100595100~5528
[0091]
表5中，分布指trgt方法与实施例1的方法得到的拷贝数值差值的分布；计数是指该分布下的频数。
[0092]
通过表5和图3能够看出，trgt方法与实施例1的方法得到的拷贝数值差值的分布主要集中在-5~0和0~5之间，发现87.5%的位点拷贝数差异在5个拷贝以内，说明本发明的方法的准确性比较高。
[0093]
相较于trgt，本方法提供的两套分型结果是基于家系的分型策略得到的两套基因组，而非简单聚类，所以本方法最后得到的遗传过程中tr的变异信息更加准确。
[0094]
此外，通过对局部位点验证，见表6，发现部分位点上，trgt会得到不存在的异常扩增拷贝数，即假阳性，而在这些位点上，本发明的方法的准确性更优于trgt。
[0095]
表6
[0096]
chrstartendmotif参考基因组拷贝数trgt检测拷贝数本方法检测拷贝数判定chr51759980817599840tttat6.44127.66.4trgt假阳性chr10134640370134640421ggag12.8333412.8trgt假阳性chr8101801252101801289tttta7.6558.617.2trgt假阳性chrx153748481153748523tttta8.4399.69.4trgt假阳性chr2155739815574006ctcc6.2391.213.2trgt假阳性chr194059099840591023aag8.3369.38.3trgt假阳性
[0097]
通过实施例1、对比例1和对比例2能够看出，本发明得到了1025347个串联重复（tr）位点从亲本遗传到子代的拷贝数变异信息（表3的diff_p1_fa，diff_p1_mo），该信息可用于遗传病的病因分析、新生突变分析等。然而对比例1检测到的有效位点数远远小于本发明的方法（表4）。对比例2在部分位点上存在较高的假阳性（表6）的问题。因此本发明的方法具有更高的准确性和全面性。
[0098]
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、
ꢀ“
示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
[0099]
尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。

技术特征：
1.一种基于核心家系的串联重复变异分型检测方法，包括：步骤1、基于父本的测序数据集和母本的测序数据集，对子代的测序数据集进行第一分型和组装，得到两套子代单倍型基因组；步骤2、以得到的子代单倍型基因组为参考，分别对父本和母本基因组进行单核苷酸位点变异检测，并对父本的测序数据集和母本的测序数据集分别进行第二分型，得到两套父本单倍型基因组、两套母本单倍型基因组，并明确由父本遗传到子代的单倍型基因组、由母本遗传到子代的单倍型基因组；步骤3、对父本遗传到子代的单倍型基因组、母本遗传到子代的单倍型基因组、两套子代单倍型基因组分别进行串联重复变异检测。2.根据权利要求1所述的方法，其特征在于，所述测序数据集包括基因组的三代测序方法得到的长reads。3.根据权利要求1所述的方法，其特征在于，所述第一分型的方法基于trio_binning方法。4.根据权利要求3所述的方法，其特征在于，所述第一分型的方法包括基于trio_binning方法的find-unique-kmers分析，得到父本和母本特异的kmer序列。5.根据权利要求4所述的方法，其特征在于，基于classify_by_kmers方法，并根据父本和母本特异的kmer序列判断子代测序数据集为属于父本的测序数据集、属于母本的测序数据集，或者为未分型的测序数据集。6.根据权利要求1所述的方法，其特征在于，所述组装的方法选自hifiasm、wtdbg2、canu、nextdenovo中的至少一种。7.根据权利要求1所述的方法，其特征在于，所述第二分型的方法基于whatshap方法。8.根据权利要求1-7中任一项所述的方法，其特征在于，所述串联重复变异检测的方法为改进的tandem repeat finder方法。9.根据权利要求8所述的方法，其特征在于，所述改进包括motif循环特点、检测区段的特征。10.根据权利要求1-7中任一项所述的方法，其特征在于，所述方法还包括：分析遗传过程中子代的两套单倍型基因组的串联重复变异信息。11.权利要求1-10中任一项所述的基于核心家系的串联重复变异分型检测方法在生物学中的应用。12.根据权利要求11所述的应用，其特征在于，所述应用为在疾病病理研究、亲子鉴定、犯罪侦查、种群遗传学中的应用。

技术总结
本发明公开了一种基于核心家系的串联重复变异分型检测方法及其应用。方法包括：基于父本的测序数据集和母本的测序数据集，对子代的测序数据集进行第一分型和组装；以得到的子代单倍型基因组为参考，分别对父本和母本基因组进行单核苷酸位点变异检测，并对父本的测序数据集和母本的测序数据集分别进行第二分型，得到两套父本单倍型基因组、两套母本单倍型基因组，并明确由父本遗传到子代的单倍型基因组、以及由母本遗传到子代的单倍型基因组；对单倍型基因组分别进行串联重复变异检测。本发明提供一种基于三代长读长测序技术，高准确度、全基因组覆盖的串联重复区域变异检测方法；而且本发明方法能得到遗传过程中串联重复发生变异的单倍型信息。发生变异的单倍型信息。发生变异的单倍型信息。


技术研发人员：权慧 杨伟飞 任雪 刘涛 车天栋 朱慧萍 王娟 李志民
受保护的技术使用者：北京安诺优达医学检验实验室有限公司 浙江安诺优达生物科技有限公司
技术研发日：2023.08.31
技术公布日：2023/10/11