CN116030892A - 一种鉴定染色体相互易位断点位置的系统和方法 - Google Patents

Abstract

本发明涉及一种鉴定染色体相互易位断点位置的系统和方法，该系统包括单细胞扩增及测序单元、测序数据清洗比对单元、计算单元、归一化单元、样本筛选单元和断点位置确定单元，单细胞扩增及测序单元用于获得测序数据；测序数据清洗比对单元用于对获得的测序数据去除二代测序接头序列和低质量碱基；计算单元用于计算基因组每500kb的窗口中比对的测序数据数目；归一化单元用于归一化测序深度，计算每个窗口矫正后的测序数据数目；样本筛选单元用于计算每个样本的离散程度，以挑选离散程度满足要求的单倍体样本，进一步筛选出具有染色体部分缺失的单倍体样本用于下一步计算；断点位置确定单元用于动态划分窗口并确定相互易位染色体断点位置。

Description

一种鉴定染色体相互易位断点位置的系统和方法

技术领域

本发明涉及染色体相互易位技术领域，具体涉及一种使用单倍体低深度测序鉴定染色体相互易位断点位置的系统和方法。

背景技术

在染色体相互易位的研究中，确定相互易位染色体的断点位置是进行后续相关研究的首要任务。通过单细胞基因组扩增及二代测序，可以得到单细胞的基因组测序数据，继而对基因组测序数据进行计算，得到单细胞染色体拷贝数变异（CNV），进一步对CNV分析可以初步确定相互易位染色体的断点区域。然而，目前通过分析CNV仅能初步确定相互易位染色体断点的大致范围，其结果不够精确；且多是直接对二倍体细胞进行测序，得到的数据离散程度较大，基于二倍体数据的CNV确定染色体易位断点的结果不够准确。

在生殖遗传的科学研究及临床应用领域，染色体相互易位是一种常见的染色体结构异常，指非同源染色体发生断裂后相互交换，形成两条新的衍生染色体的结构变化。在此过程中，染色体断裂的位置称为染色体易位断点。据统计，染色体相互易位在新生儿中的患病率约为0.14%，在不育人群中的患病率约为0.6%。

染色体相互易位携带者由于基因总数未丢失，因此其外貌等表型通常正常，但也有一部分存在微缺失、重复等遗传学变异，从而导致自闭症、先天畸形等疾病。对于生育下一代来说，携带者的异常配子与配偶的正常配子结合会导致流产、死胎、畸形或染色体异常胎儿等。即使染色体数量平衡的配子与配偶的正常配子结合，后代同为染色体相互易位携带者的概率也为50%。因此，能够鉴定出胚胎中的染色体相互易位携带胚胎，有利于阻断相互易位染色体在家族中的遗传，有利于优生优育，而相互易位染色体断点位置的确定是该过程的重要步骤。

目前用于鉴定染色体相互易位的主要技术有：荧光原位杂交（FISH）技术，比较基因组杂交（CGH）技术，单核苷酸多态性芯片（SNP array），二代测序（NGS）技术等。

FISH技术与其他的检测技术相比，检测周期长，操作过程比较复杂，检测位点少，无法检测微缺失和微重复。由于每种相互易位都需要设计专门的荧光探针来区分，在大多数情况下，仅用于检测13，16，18，21，22，X和Y染色体的非整倍体，无法对全基因组染色体进行筛查。此外，受操作经验、杂交效率和探针质量的影响，存在较高的假阳性率和假阴性率，且FISH无法精确检测易位断点。

目前CGH技术可以对染色体易位导致的非整倍体进行筛查，但不能检测单倍体、多倍体，仅能检测染色体的非平衡易位：即染色体片段的缺失和重复。此外，存在非均一扩增现象。SNP array的方法较CGH检测位点多，覆盖度好，但SNP分布不会绝对均一，在进行分析诊断前，需要挑选突变位点周围的SNP位点进行扩增并检测，得到每个SNP位点的基因型信息，如果挑选的SNP位点在胚胎活检细胞DNA扩增产物中位于低扩增区域，则后续实验难以实现SNP位点的有效鉴定，则会给检测过程造成很大的难度且会极大地增加工作量。因此，CGH和SNP array应用于易位携带染色体诊断具有一定局限性，这两种方法仅能确定易位断点的大致所在范围。

NGS是目前具有高度潜力和发展空间的可检测染色体易位的技术。NGS技术最主要的特征是高通量测序，测序时间和成本都显著降低，需要的样本少，具有高灵敏度，大通量，高自动化的特点，能够检测包括点突变，基因拷贝数变化和染色体易位等在内的多种基因改变。

目前报道过多种定位相互易位染色体断点的方法，BAC-FISH+染色体步移的技术需要大量BAC和反复FISH寻找被打断的BAC，染色体步移费时费力，且如果断点附近有较多重复区域则无法精确定位。Flow sorting+NGS技术需要制备中期染色体，需要特殊流式细胞仪，测序成本高，难于推广和临床应用。BAC-FISH+区域捕获+NGS技术很难将区域缩小至较小的范围，需要设计特殊芯片捕获，仍然难以大规模应用。环化建库+NGS技术需要特殊建库方法，建库难度高，测序量大，在获得具体的断裂点位置时还要经过大量的PCR后采用一代测序验证，费用高、耗时长。染色体显微切割+NGS技术，该方法操作需要超高倍染色体显微切割仪器，显微切割的片段还要进行纯化后才可进行后续扩增实验，该法技术难度较大，过程较为繁琐。

现有技术中，专利申请文件CN105039569A公布了一种相互易位染色体的断点分析方法。该方法利用染色体显微切割技术获得易位染色体断点附近的染色质，先进行预扩增处理、PCR扩增处理，而后对扩增产物进行PCR荧光标记，并在有丝分裂中期核分裂相中进行FISH，以验证显微切割是否成功；后将单细胞扩增产物和多重PCR扩增产物一同混合建库，进行后续NGS测序分析，在获得相互易位染色体的断点位置后设计引物，对样本的gDNA进行PCR扩增，并对PCR产物进行Sanger测序验证，从而获取断点位置。

专利申请文件CN106650310A公布了一种鉴别染色体相互易位携带胚胎和正常胚胎的方法。该方法首先对染色体相互易位携带者夫妇双方、至少一名携带者亲属和染色体相互易位携带者夫妇体外受精胚胎进行大规模SNP基因型检测，再确定有效SNPs位点，通过家系连锁分析得到染色体单体型并判断染色体易位断裂点区域是否发生了同源重组，从而判断胚胎是否携带易位染色体。

专利申请文件CN106834490A公布了一种鉴定胚胎相互易位断裂点和易位携带状态的方法。该方法先对待测样本和父母DNA进行扩增，构建文库后测序，将参考基因组分成N个区域片段，其中每个区域片段为一个窗口，计算每个窗口的拷贝数，确定正常拷贝数的阈值范围，和异常染色体拷贝数范围，二者变化的窗口区域的起始和终止的位置即为拷贝数变异的起始和终止断点；选取多个拷贝数异常的胚胎，计算每个胚胎的两条染色体相互易位的断裂点，分别计算两条染色体易位断裂点的三均值，即为胚胎的易位断裂点。

专利文件CN106929595B公布了一种鉴定胚胎平衡易位携带状态的系统和方法。该方法首先对单精子进行扩增，构建文库后测序，对测序结果进行染色体拷贝数分析，仅能以500kb的固定窗口大小初步确定染色体易位位置的大致范围，再通过SNP分析比对，确定胚胎平衡易位携带状态。

专利申请文件CN105039569A涉及染色体显微切割+NGS技术，该方法需要制备有丝分裂中期的细胞样本，且需要超高倍显微切割仪器，显微切割的片段还要进行多次扩增实验。该方法技术难度较大，检测周期长，仪器要求高，不适用于大规模推广。

专利申请文件CN106650310A需要对携带者夫妇双方、至少一名携带者亲属进行大规模SNP检测，耗费时间长，花费成本高，仅能通过SNP位点来粗略判定染色体断点位置，为后续检测带来可能的技术风险。

专利申请文件CN106834490A将测序区段划分为长度50kb的窗口，导致其测量结果仅能定位到长度为100kb的窗口，不够精确，且该方法是对胚胎及父母的二倍体细胞进行测序，测序范围较大。

专利文件CN106929595B采取单倍体精子细胞扩增测序，通过染色体拷贝数分析仅能初步确定相互易位染色体断点所在500kb的范围，其结果不够精确，分辨率较低，后续仍需要进行SNP分析对比，过程较为繁琐。

已有的检测技术都不是直接通过对单倍体细胞测序得到相互易位染色体的断点位置，且已有技术得到的断点位置都不够精确，多是完成测序后进行snp分析，过程较为复杂。因此，本领域迫切需要开发一种能够更高效精确定位相互易位染色体断点位置的方法。

本申请中部分缩略语和关键术语的定义如下：

PGT: Pre-implantation genetic testing，胚胎植入前遗传学检测。

CNV: Copy number variation，拷贝数变异。

FISH: Fluorescence in situ hybridization，荧光原位杂交。

CGH: Comparative genomic hybridization，比较基因组杂交。

SNP array: Single nucleotide polymorphism array，单核苷酸多态性微阵列芯片。

NGS: Next generation sequencing，二代测序。

WGA: Whole genome amplification，全基因组扩增。

MALBAC: Multiple annealing and looping-based amplification cycles，多次退火环状循环扩增技术。

PCR: Polymerase chain reaction，聚合酶链反应。

发明内容

本发明旨在提供一种使用单倍体低深度测序鉴定染色体相互易位断点位置的系统和方法，以解决如何精确地定位相互易位染色体的断点位置的问题。

本发明的目的是解决现有技术的不足，提供一种使用单倍体低深度测序鉴定染色体相互易位断点位置的系统，包括单细胞扩增及测序单元、测序数据清洗比对单元、计算单元、归一化单元、样本筛选单元和断点位置确定单元，所述的单细胞扩增及测序单元用于获得相互易位携带者的精子细胞，进行单细胞DNA扩增，扩增后对DNA进行纯化，继而构建二代测序文库，获得测序数据（reads）；所述的测序数据清洗比对单元用于对获得的测序数据去除二代测序接头序列和低质量碱基，保留序列长度大于36 bp的测序数据；将处理后的测序数据比对至UCSC hg38人类参考基因组；对比对后的测序数据去除结果中的低比对质量序列和PCR重复序列，得到唯一比对、去重后的测序数据；所述的计算单元用于计算基因组每500kb的窗口中比对的测序数据数目；所述的归一化单元用于归一化测序深度，计算每个窗口矫正后的测序数据数目；所述的样本筛选单元用于在得到归一化处理的矫正测序数据数目后，计算每个样本的离散程度，以挑选离散程度满足要求的单倍体样本，在此基础上进一步筛选出具有染色体缺失的单倍体样本用于下一步计算；所述的断点位置确定单元用于动态划分窗口并确定相互易位染色体断点位置，在500kb窗口基础上完成样本离散程度筛选和具有染色体缺失样本的筛选后，以高分辨率5kb窗口为基准，确定染色体易位断点的最优精确判定位置，当在此分辨率下无法确定最优位置时，每次将窗口大小增加5kb，直至得到断点的最优判定位置。

优选地，所述的测序数据清洗比对单元使用trim_galore软件默认参数对获得的测序数据去除二代测序接头序列和低质量碱基，保留序列长度大于36 bp的测序数据；使用BWA软件以默认参数将处理后的测序数据比对至UCSC hg38人类参考基因组；对比对后的测序数据使用samtools软件去除结果中的低比对质量序列和PCR重复序列，得到唯一比对、去重后的测序数据。

优选地，所述的计算单元在得到比对结果后，将参考基因组划分为500kb大小的窗口，根据序列在基因组上的位置，使用samtools 软件统计每个样本在基因组上每个500kb窗口的比对测序数据数目。

优选地，所述的归一化单元通过归一化不同窗口所得的测序数据数目来消除由于扩增和建库所带来的偏倚；具体为：i为一批测序中第i个样本，j为基因组中第j个窗口，为i样本的测序深度，为i样本所测碱基总数，为人类基因组序列长度，为i样本中j窗口比对的测序数据数目，为i样本j窗口的归一化测序深度，n为样本总数， 为j窗口的归一化因子，为i样本j窗口矫正后的比对测序数据数目；

①计算每个样本中每个窗口的归一化测序深度：

= ；

= ；

②计算每个窗口的归一化因子：

= ；

③计算每个样本中每个窗口的矫正后的测序数据数目：

= 。

优选地，所述的样本筛选单元计算每个样本的离散程度，具体过程为：i为一批测序中第i个样本，j为基因组中第j个窗口，N为该样本全部染色体划分的窗口总数，i样本的离散程度用表示，为i样本j窗口矫正后的比对测序数据数目，为i样本每个窗口矫正后测序数据数目的平均值；SD的值表示样本的离散程度，与样本质量呈反比；

①计算每个样本矫正后拷贝数的平均值：

= ；

②计算每个样本的SD值，表示每个样本的离散程度：

= ；

③筛选出SD＜0.25的单倍体样本。

优选地，所述的样本筛选单元还用于进一步根据矫正测序数据数目所得拷贝数变异（CNV）结果，筛选出染色体有缺失的单细胞测序样本。

优选地，所述的断点位置确定单元首先将参考基因组划分为5kb大小的窗口，使用samtools统计所选样本每个窗口的比对测序数据数目，使用所述的归一化单元对每个样本每个窗口的测序数据数目进行矫正，得到每个样本每个窗口的矫正后的测序数据数目；由于所用样本是单倍体细胞，因此使用0表示缺失，1表示正常单倍体，2表示重复；将所有样本的每个窗口归为0、1、2三个拷贝数数值，如果矫正后测序数据数目小于0.4，则将该窗口归为拷贝数0；如果矫正后测序数据数目在0.4到1.4之间，则将该窗口归为拷贝数1；如果矫正后测序数据数目大于1.4，则将该窗口归为拷贝数2；当拷贝数由连续的0转变为连续的1或2时，0拷贝转变为1或2拷贝的窗口即为相互易位染色体断点所在的区域，根据拷贝数变化记录断点位置所在窗口。

优选地，所述的断点位置确定单元在划分窗口大小为5kb的基准上，进行样本间断点区域的比较，对样本所得的断点区域进行汇总，计算出现次数最多的区域的频数，当该频数占总频数2/3以上时，确定该区域为相互易位染色体的断点位置。

优选地，所述的断点位置确定单元在划分窗口为5kb大小的分辨率下无法判定断点位置时，逐步增加窗口大小，每次将窗口增加5kb，直至区域最大频数与总频数的比值达到2/3以上时，确定该区域为断点位置。

本发明还提供一种使用单倍体低深度测序鉴定染色体相互易位断点位置的方法，包括以下步骤：

第一步、获得相互易位携带者的单倍体细胞（如精子等），进行单细胞DNA扩增，扩增后对DNA进行纯化，继而构建二代测序文库，获得测序数据（reads）；

第二步、对获得的测序数据去除二代测序接头序列和低质量碱基，保留序列长度大于36 bp的测序数据；将处理后的测序数据比对至UCSC hg38人类参考基因组；对比对后的测序数据去除结果中的低比对质量序列和PCR重复序列，得到唯一比对、去重后的测序数据；

第三步、计算基因组每500kb的窗口中比对的测序数据数目；

第四步、归一化测序深度，计算每个窗口矫正后的测序数据数目；

第五步、在得到归一化处理的矫正测序数据数目后，计算每个样本的离散程度，以挑选离散程度满足要求的单倍体样本，在此基础上进一步筛选出具有染色体部分缺失的单倍体样本用于下一步计算；

第六步、动态划分窗口并确定相互易位染色体断点位置，在500kb窗口基础上完成样本离散程度筛选和具有染色体缺失样本的筛选后，以高分辨率5kb窗口为基准，确定染色体易位断点的最优精确判定位置，当在此分辨率下无法确定最优位置时，每次将窗口大小增加5kb，直至得到断点的最优判定位置。

有益效果

与现有技术相比，本发明的有益效果是：

本发明提供了一种采取低测序深度达到高精确度的定位相互易位染色体断点位置的方法，对单倍体精子细胞测序后进行筛选，去掉数据离散程度大（SD值较大）的样本，选取染色体有部分缺失的样本，继而在划分窗口长度为5kb的基础上进行动态划分窗口处理，选取划分窗口的最优大小，在确保断点位置准确性的同时达到最精确的程度。本发明利用染色体缺失部分的数据特性，精准确定相互易位染色体的断点位置。由于单倍体精子的易位相关染色体在缺失部分数据相对二倍体及重复区域离散程度最小，得到的结果更加准确。同时，通过动态划分窗口大小的方法，确定最优计算结果，因此可在高分辨率下得到易位的精确断点位置。本发明在降低操作难度的同时提高了定位相互易位染色体断点的精确性。

附图说明

附图用来提供对本发明技术方案的进一步理解，并且构成说明书的一部分，与本申请的具体实施方式一起用于解释本发明的技术方案，并不构成对本发明技术方案的限制。

图1是本发明的技术流程图。

图2A至图2K是窗口大小为500kb，6号染色体有缺失的11个精子细胞cnv图。

图3A至图3M是窗口大小为500kb，7号染色体有缺失的13个精子细胞cnv图。

图4是窗口大小为5kb,6号染色体断点所在区域分布情况示意图。

图5是窗口大小为10kb，6号染色体断点所在区域分布情况示意图。

图6是窗口大小为5kb，7号染色体断点所在区域分布情况示意图。

图7A至图7X是窗口大小为500kb，24个精子细胞的异常6号染色体cnv

图8A至图8X是窗口大小为500kb，24个精子细胞的异常7号染色体cnv图。

图9是窗口大小为500kb,6号染色体断点所在区域分布情况示意图。

图10是窗口大小为500kb,7号染色体断点所在区域分布情况示意图。

具体实施方式

在下文中更详细地描述了本发明以有助于对本发明的理解。

本发明的目的是通过对单倍体细胞进行二代测序，通过一系列优化的计算方法，对所得样本进行筛选，使用动态划分窗口长度的高分辨率计算方法，更为精确地定位相互易位染色体的断点位置。

如图1所示，本发明的基本方法是首先对获得的男性单倍体细胞（如精子）进行扩增、测序、比对，筛选出数据离散程度SD值小且有染色体部分缺失的单倍体细胞样本，将基因组划分为长度5kb的窗口，确定每个单倍体样本的易位断点位置，通过多个单倍体样本细胞的相互比较，准确判定该男性的染色体相互易位的断点位置。在5kb标准下所得结果达不到期望的集中程度时，增加窗口长度，每次以5kb大小递增划分的窗口长度，再次对每个单倍体样本的断点位置和多个样本的断点位置集中度进行比较，直至得到最优的断点位置，以此准确判定该男性的染色体相互易位的断点位置。

在一个具体实施例中，本申请选取一个6号染色体（chr6）和7号染色体（chr7）相互易位携带者作为说明案例，其核型为46,XY,t(6;7)(p23;q36)。本技术方案的目的是使用该携带者的单倍体精子细胞，通过筛选样本和动态划分窗口长度，确定相互易位染色体断点位置。该批次一共获得47个单倍体精子细胞。

1. 单细胞扩增及测序

获得相互易位携带者的精子细胞，首先通过MALBAC方法进行单细胞DNA扩增，扩增后对DNA进行纯化，继而构建二代测序文库，并测序获得测序数据（reads）。测序使用Illumina Hiseq X10测序仪，测序方式采用双端150 bp，平均每个样本测序深度为0.3x。

2.测序数据清洗、比对至参考基因组及比对后处理

所得测序数据（reads）使用trim_galore软件默认参数去除二代测序接头序列、低质量碱基，保留序列长度大于36 bp的reads。经过处理后的reads，使用BWA软件以默认参数比对至UCSC hg38人类参考基因组。比对后的reads使用samtools软件去除结果中的低比对质量序列、PCR重复序列，得到唯一比对、去重后的reads。

3.计算基因组每500kb的窗口中比对的reads数

得到比对结果后，将参考基因组划分为500kb大小的窗口，根据序列在基因组上的位置，使用samtools 软件统计每个样本在基因组上每个500kb窗口的比对reads数。

4.归一化测序深度，计算每个窗口矫正后的reads数

由于单细胞二代测序数据比对至每个窗口的reads数受到扩增和建库过程的影响，因此通过归一化不同窗口所得的reads数来消除由于扩增和建库所带来的偏倚。具体为：i为一批测序中第i个样本，j为基因组中第j个窗口， 为i样本的测序深度，为i样本所测碱基总数，为人类基因组序列长度，为i样本中j窗口比对的reads数，为i样本j窗口的归一化测序深度，n为样本总数，为j窗口的归一化因子，为i样本j窗口矫正后的比对reads数。

①计算每个样本中每个窗口的归一化测序深度：

= ；

= ；

②计算每个窗口的归一化因子：

= ；

③计算每个样本中每个窗口的矫正后reads数：

= ；

5.样本筛选

（1）得到归一化处理的矫正reads数后，计算每个样本的离散程度，以挑选离散程度小的高质量样本用于下一步计算。i为一批测序中第i个样本，j为基因组中第j个窗口，N为该样本全部染色体划分的窗口总数，i样本的离散程度用表示，为i样本j窗口矫正后的比对reads数，为i样本每个窗口矫正后reads数的平均值。

①计算每个样本矫正后拷贝数的平均值：

= 

②计算每个样本的SD值，表示每个样本的离散程度：

= 

③筛选出SD＜0.25的单倍体样本。

SD值表示样本的离散程度，与样本质量呈反比，故剔除SD值较大的样本以提高后续分析的准确性。本实施例共47个精子细胞，根据多样本SD值的汇总，将SD=0.25设为界值，在此标准下，剔除1个不符合标准的样本WN-38，SD=0.2701，筛选出46个符合标准的精子细胞，结果如[表1]所示。

[表1]：46个SD<0.25的样本及其SD值及1个剔除样本的SD值

（2）进一步根据矫正reads数所得CNV结果，筛选出染色体有缺失的单细胞测序样本，并进行后续分析。

由于单倍体细胞正常拷贝数为1，染色体易位时会发生部分缺失或重复，缺失区域的拷贝数为0，重复区域的拷贝数为2。染色体存在重复区域而没有缺失时，拷贝数发生1-2转变的区域由于测序reads比对及计数导致区域附近离散程度较大，其推断结果没有0-1或0-2转变准确，故筛选出染色体部分缺失的样本进行后续分析，筛选结果如下：

10个精子细胞仅6号染色体有缺失：WN-1、WN-6、WN-13、WN-28、WN-29、WN-42、WN-45、WN-46、WN-47、WN-48；

12个精子细胞仅7号染色体有缺失：WN-2、WN-4、WN-5、WM-9、WN-11、WN-16、WM-23、WN-32、WN-33、WN-34、WN-37、WN-41；

1个精子细胞6号和7号染色体都有缺失：WN-39。

①11个6号染色体有缺失的精子细胞：WN-1、WN-6、WN-13、WN-28、WN-29、WN-39、WN-42、WN-45、WN-46、WN-47、WN-48，其cnv如图2A至图2K所示。

②13个7号染色体有缺失的精子细胞：WN-2、WN-4、WN-5、WM-9、WN-11、WN-16、WM-23、WN-32、WN-33、WN-34、WN-37、WN-39、WN-41，其cnv如图3A至图3M所示。

6.动态划分窗口确定相互易位染色体断点位置

在500kb窗口基础上完成高质量样本筛选和具有染色体缺失样本的筛选后，以高分辨率5kb窗口为基准，确定染色体易位断点的最优精确判定位置，当在此分辨率下无法确定最优位置时，每次将窗口大小增加5kb，直至得到断点的最优判定位置。具体为：

（1）首先将参考基因组划分为5kb大小的窗口，使用samtools统计所选样本每个窗口的比对reads数，使用上述步骤4所述方法对每个样本每个窗口的reads数进行矫正，得到每个样本每个窗口的矫正后reads数。 由于所用样本是单倍体细胞，因此使用0表示缺失，1表示正常单倍体，2表示重复。将所有样本的每个窗口归为0,1,2三个拷贝数数值，如果矫正后reads数小于0.4，则将该窗口归为拷贝数0；如果矫正后reads数在0.4到1.4之间，则将该窗口归为拷贝数1；如果矫正后reads数大于1.4，则将该窗口归为拷贝数2。当拷贝数由连续的0转变为连续的1或2时，0拷贝转变为1或2拷贝的窗口即为相互易位染色体断点所在的区域，根据拷贝数变化记录断点位置所在窗口。

本实施例中，对于筛选出的6号染色体有缺失的11个样本WN-1、WN-6、WN-13、WN-28、WN-29、WN-39、WN-42、WN-45、WN-46、WN-47、WN-48，记录其拷贝数发生0-1/0-2变化时上下游共4个窗口的位置及相应的矫正后拷贝数，如[表2]所示。

[表2]：5kb分辨率下，6号染色体拷贝数发生0-1/0-2转变时上下游4个窗口的位置及拷贝数。

（2）记录筛选出的11个样本WN-1、WN-6、WN-13、WN-28、WN-29、WN-39、WN-42、WN-45、WN-46、WN-47、WN-48中6号染色体矫正后拷贝数发生0-1/0-2转变的区域，如[表3]所示。

[表3]：窗口大小为5kb，6号染色体易位断点所在区域

（3）在划分窗口大小为5kb的基准上，进行样本间断点区域的比较，对样本所得的断点区域进行汇总，计算出现次数最多的区域的频数，当该频数占总频数2/3以上时，可确定该区域为相互易位染色体的断点位置。

本实施例中，对所得的区域进行汇总，可以发现6号染色体断点位置集中在19675001-19685000区域，该区域频数为6，总频数为11，其分布情况如图4所示，频数之比6/11未达到2/3，在5kb基准下尚不能确定断点位置，因此进一步分析。

（4）在划分窗口为5kb大小的分辨率下无法判定断点位置时，逐步增加窗口大小，每次将窗口增加5kb，直至区域最大频数与总频数的比值达到2/3以上时，即可确定该区域为断点位置。本实施例中，在窗口大小为5kb的分辨率下，6号染色体断点位置分布在19675000-19685000区域的频数比值达不到2/3，因此，需要划分10kb窗口再次对6号染色体断点位置进行分析。

①对于筛选出的6号染色体有缺失的11个样本WN-1、WN-6、WN-13、WN-28、WN-29、WN-39、WN-42、WN-45、WN-46、WN-47、WN-48，划分窗口大小为10kb，记录其矫正后拷贝数发生0-1/0-2变化时上下游共4个窗口的位置及相应的矫正后拷贝数，如[表4]所示。

[表4]：10kb分辨率下，6号染色体矫正后拷贝数发生0-1/0-2转变时上下游共4个窗口的位置及矫正后拷贝数。

②记录11个样本WN-1、WN-6、WN-13、WN-28、WN-29、WN-39、WN-42、WN-45、WN-46、WN-47、WN-48中6号染色体矫正后拷贝数发生0-1/0-2转变的区域，如[表5]所示。

[表5]：窗口大小为10kb，6号染色体断点所在区域汇总

③在10kb窗口大小的基准下，对所得的区域进行汇总，可以发现6号染色体断点位置集中在19670001-19690000区域，该区域频数为10，总频数为11，频数之比10/11达到2/3以上，可以确定6号染色体断点位置在19670001-19690000区域，其分布情况如图5所示。

（5）本实施例中，对于同样存在易位的7号染色体，使用和6号染色体相同的判定策略。

①筛选出7号染色体有缺失的13个样本WN-2、WN-4、WN-5、WM-9、WN-11、WN-16、WM-23、WN-32、WN-33、WN-34、WN-37、WN-39、WN-41，记录其矫正后拷贝数发生0-1/0-2变化时上下游共4个窗口的位置及相应的矫正后拷贝数，如[表6]所示。

[表6]：5kb分辨率下，7号染色体矫正后拷贝数发生0-1/0-2转变时上下游共4个窗口的位置及矫正后拷贝数。

②记录13个样本WN-2、WN-4、WN-5、WM-9、WN-11、WN-16、WM-23、WN-32、WN-33、WN-34、WN-37、WN-39、WN-41中7号染色体矫正后拷贝数发生0-1/0-2转变的区域，如[表7]所示。

[表7]：窗口大小为5kb，7号染色体断点所在区域汇总

③在划分窗口大小为5kb的基准下，对所得的区域进行汇总，可以发现7号染色体断点位置全部分布在154670001-154680000区域，该区域频数为13，总频数13，频数之比大于2/3，可以确定7号染色体断点位置在154670001-154680000区域，其分布情况如图6所示。

7.对比实验：使用已有专利文件CN106929595B检测本实施例中相互易位染色体断点位置。

（1）对于本批次精子细胞，同时采用已有专利文件CN106929595B提到的划分500kb窗口的检测方法判断相互易位染色体断点位置。本实施例中有24个精子细胞中的6号和7号染色体存在缺失或重复，24个异常6号染色体cnv如图7所示，24个异常7号染色体cnv如图8所示。

（2）①筛选出6号染色体有异常的24个精子细胞，记录其矫正后拷贝数发生转变时上下游共4个窗口的位置及相应的矫正后拷贝数，如[表8]所示。

[表8]： 6号染色体矫正拷贝数发生转变时上下游共4个窗口的位置及矫正后拷贝数。

②筛选出7号染色体有异常的24个精子细胞，记录其矫正后拷贝数发生转变时上下游共4个窗口的位置及相应的矫正后拷贝数，如[表9]所示。

[表9] ：7号染色体矫正拷贝数发生转变时上下游共4个窗口的位置及矫正后拷贝数。

如[表10]所示,7号染色体断点所在区域如[表11]所示。6号染色体断点所在区域分布情况如[图9]所示，7号染色体断点所在区域分布情况如[图10]所示。通过已有专利文件CN106929595B公开的技术方案，可粗略确定相互易位染色体断点位置在chr6的19000001-20000000区域。7号染色体有14个断点出现在154000001-155000000区域，有10个断点出现在154500001-155500000区域，无法确定7号染色体断点位置。[表10]示出了窗口大小为500kb，6号染色体相互易位断点所在窗口位置。[表11]示出了窗口大小为500kb，7号染色体相互易位断点所在窗口位置。

[表10]：窗口大小为500kb，6号染色体断点所在区域汇总。

[表11]：窗口大小为500kb，7号染色体断点所在区域汇总。

（4）使用已有的专利文件CN106929595B公开的技术方案仅能确定6号染色体断点位置位于19000001-20000000区域，范围较大；另尚不能确定7号染色体断点位置，而本申请提出的方法得到的断点位置更为精确，6号染色体断点位置为19670001-19690000区域，7号染色体断点位置为154670001-154680000区域，得到的断点位置比已有专利所得结果更加精确。综上所述，本申请提出的方法所得的结果优于已有的专利文件CN106929595B公开的技术方案。

基于以上实施例和对照结果，本发明提供一种使用单倍体低深度测序鉴定染色体相互易位断点位置的系统，所述的鉴定染色体相互易位断点位置的系统包括单细胞扩增及测序单元、测序数据清洗比对单元、计算单元、归一化单元、样本筛选单元和断点位置确定单元，所述的单细胞扩增及测序单元用于获得相互易位携带者的精子细胞，进行单细胞DNA扩增，扩增后对DNA进行纯化，继而构建二代测序文库，获得测序数据（reads）；所述的测序数据清洗比对单元用于对获得的测序数据去除二代测序接头序列和低质量碱基，保留序列长度大于36 bp的测序数据；将处理后的测序数据比对至UCSC hg38人类参考基因组；对比对后的测序数据去除结果中的低比对质量序列和PCR重复序列，得到唯一比对、去重后的测序数据；所述的计算单元用于计算基因组每500kb的窗口中比对的测序数据数目；所述的归一化单元用于归一化测序深度，计算每个窗口矫正后的测序数据数目；所述的样本筛选单元用于在得到归一化处理的矫正测序数据数目后，计算每个样本的离散程度，以挑选离散程度满足要求的单倍体样本，在此基础上进一步筛选出具有染色体缺失的单倍体样本用于下一步计算；所述的断点位置确定单元用于动态划分窗口并确定相互易位染色体断点位置，在500kb窗口基础上完成高质量样本筛选和具有染色体缺失样本的筛选后，以高分辨率5kb窗口为基准，确定染色体易位断点的最优精确判定位置，当在此分辨率下无法确定最优位置时，每次将窗口大小增加5kb，直至得到断点的最优判定位置。

优选地，所述的测序数据清洗比对单元使用trim_galore软件默认参数对获得的测序数据去除二代测序接头序列和低质量碱基，保留序列长度大于36 bp的测序数据；使用BWA软件以默认参数将处理后的测序数据比对至UCSC hg38人类参考基因组；对比对后的测序数据使用samtools软件去除结果中的低比对质量序列和PCR重复序列，得到唯一比对、去重后的测序数据。

优选地，所述的计算单元在得到比对结果后，将参考基因组划分为500kb大小的窗口，根据序列在基因组上的位置，使用samtools 软件统计每个样本在基因组上每个500kb窗口的比对测序数据数目。

优选地，所述的归一化单元通过归一化不同窗口所得的测序数据数目来消除由于扩增和建库所带来的偏倚；具体为：i为一批测序中第i个样本，j为基因组中第j个窗口，为i样本的测序深度，为i样本所测碱基总数，为人类基因组序列长度，为i样本中j窗口比对的测序数据数目，为i样本j窗口的归一化测序深度，n为样本总数，为j窗口的归一化因子，为i样本j窗口矫正后的比对测序数据数目；

①计算每个样本中每个窗口的归一化测序深度：

= ；

= ；

②计算每个窗口的归一化因子：

= ；

③计算每个样本中每个窗口的矫正后的测序数据数目：

=。

优选地，所述的样本筛选单元用于计算每个样本的离散程度，i为一批测序中第i个样本，j为基因组中第j个窗口，N为该样本全部染色体划分的窗口总数，i样本的离散程度用表示，为i样本j窗口矫正后的比对测序数据数目，为i样本每个窗口矫正后测序数据数目的平均值；SD的值表示样本的离散程度，与样本质量呈反比。具体过程为：

①计算每个样本矫正后拷贝数的平均值：

= ；

②计算每个样本的SD值，表示每个样本的离散程度：

=；

③筛选出SD＜0.25的单倍体样本；

SD值表示样本的离散程度，与样本质量呈反比，故剔除SD值较大的样本以提高后续分析的准确性。

优选地，所述的样本筛选单元还用于进一步根据矫正测序数据数目所得拷贝数变异（CNV）结果，筛选出23个染色体有缺失的单细胞测序样本。

优选地，所述的断点位置确定单元首先将参考基因组划分为5kb大小的窗口，使用samtools统计所选样本每个窗口的比对测序数据数目，使用所述的归一化单元对每个样本每个窗口的测序数据数目进行矫正，得到每个样本每个窗口的矫正后的测序数据数目；由于所用样本是单倍体细胞，因此使用0表示缺失，1表示正常单倍体，2表示重复；将所有样本的每个窗口归为0、1、2三个拷贝数数值，如果矫正后测序数据数目小于0.4，则将该窗口归为拷贝数0；如果矫正后测序数据数目在0.4到1.4之间，则将该窗口归为拷贝数1；如果矫正后测序数据数目大于1.4，则将该窗口归为拷贝数2；当拷贝数由连续的0转变为连续的1或2时，0拷贝转变为1或2拷贝的窗口即为相互易位染色体断点所在的区域，根据拷贝数变化记录断点位置所在窗口。

优选地，所述的断点位置确定单元在划分窗口大小为5kb的基准上，进行样本间断点区域的比较，对样本所得的断点区域进行汇总，计算出现次数最多的区域的频数，当该频数占总频数2/3以上时，确定该区域为相互易位染色体的断点位置。

优选地，所述的断点位置确定单元在划分窗口为5kb大小的分辨率下无法判定断点位置时，逐步增加窗口大小，每次将窗口增加5kb，直至区域最大频数与总频数的比值达到2/3以上时，确定该区域为断点位置。

本发明还提供一种使用单倍体低深度测序鉴定染色体相互易位断点位置的方法，包括以下步骤：

第一步、获得相互易位携带者的精子细胞，进行单细胞DNA扩增，扩增后对DNA进行纯化，继而构建二代测序文库，获得测序数据（reads）；

第二步、对获得的测序数据去除二代测序接头序列和低质量碱基，保留序列长度大于36 bp的测序数据；将处理后的测序数据比对至UCSC hg38人类参考基因组；对比对后的测序数据去除结果中的低比对质量序列和PCR重复序列，得到唯一比对、去重后的测序数据；

第三步、计算基因组每500kb的窗口中比对的测序数据数目；

第四步、归一化测序深度，计算每个窗口矫正后的测序数据数目；

第五步、在得到归一化处理的矫正测序数据数目后，计算每个样本的离散程度，以挑选离散程度满足要求的单倍体样本，进一步筛选出具有染色体部分缺失的单倍体样本用于下一步计算；

第六步、动态划分窗口并确定相互易位染色体断点位置，在500kb窗口基础上完成高质量样本筛选和具有染色体缺失样本的筛选后，以高分辨率5kb窗口为基准，确定染色体易位断点的最优精确判定位置，当在此分辨率下无法确定最优位置时，每次将窗口大小增加5kb，直至得到断点的最优判定位置。

优选地，第二步中，使用trim_galore软件默认参数对获得的测序数据去除二代测序接头序列和低质量碱基，保留序列长度大于36 bp的测序数据；使用BWA软件以默认参数将处理后的测序数据比对至UCSC hg38人类参考基因组；对比对后的测序数据使用samtools软件去除结果中的低比对质量序列和PCR重复序列，得到唯一比对、去重后的测序数据。

优选地，第三步中，在得到比对结果后，将参考基因组划分为500kb大小的窗口，根据序列在基因组上的位置，使用samtools 软件统计每个样本在基因组上每个500kb窗口的比对测序数据数目。

优选地，第四步中，通过归一化不同窗口所得的测序数据数目来消除由于扩增和建库所带来的偏倚；具体为：i为一批测序中第i个样本，j为基因组中第j个窗口，为i样本的测序深度，为i样本所测碱基总数，为人类基因组序列长度，为i样本中j窗口比对的测序数据数目，为i样本j窗口的归一化测序深度，n为样本总数，为j窗口的归一化因子，为i样本j窗口矫正后的比对测序数据数目；

①计算每个样本中每个窗口的归一化测序深度：

=；

=；

②计算每个窗口的归一化因子：

=

③计算每个样本中每个窗口的矫正后的测序数据数目：

=。

优选地，第五步中，计算每个样本的离散程度，i为一批测序中第i个样本，j为基因组中第j个窗口，N为该样本全部染色体划分的窗口总数，i样本的离散程度用表示，为i样本j窗口矫正后的比对测序数据数目，为i样本每个窗口矫正后测序数据数目的平均值；SD的值表示样本的离散程度，与样本质量呈反比。具体过程为：

①计算每个样本矫正后拷贝数的平均值：

=；

②计算每个样本的SD值，表示每个样本的离散程度：

=；

③筛选出SD＜0.25的单倍体样本；

SD值表示样本的离散程度，与样本质量呈反比，故剔除SD值较大的样本以提高后续分析的准确性。

优选地，第六步中，首先将参考基因组划分为5kb大小的窗口，使用samtools统计所选样本每个窗口的比对测序数据数目，使用所述的归一化单元对每个样本每个窗口的测序数据数目进行矫正，得到每个样本每个窗口的矫正后的测序数据数目；由于所用样本是单倍体细胞，因此使用0表示缺失，1表示正常单倍体，2表示重复；将所有样本的每个窗口归为0、1、2三个拷贝数数值，如果矫正后测序数据数目小于0.4，则将该窗口归为拷贝数0；如果矫正后测序数据数目在0.4到1.4之间，则将该窗口归为拷贝数1；如果矫正后测序数据数目大于1.4，则将该窗口归为拷贝数2；当拷贝数由连续的0转变为连续的1或2时，0拷贝转变为1或2拷贝的窗口即为相互易位染色体断点所在的区域，根据拷贝数变化记录断点位置所在窗口；在划分窗口大小为5kb的基准上，进行样本间断点区域的比较，对样本所得的断点区域进行汇总，计算出现次数最多的区域的频数，当该频数占总频数2/3以上时，确定该区域为相互易位染色体的断点位置；在划分窗口为5kb大小的分辨率下无法判定断点位置时，逐步增加窗口大小，每次将窗口增加5kb，直至区域最大频数与总频数的比值达到2/3以上时，确定该区域为断点位置。

本发明采用单倍体细胞测序、分步筛选样本及动态划分窗口及判定技术，能够精确定位相互易位的染色体断点位置，同时，利用单倍体精子基因组测序技术来研究相互易位染色体断点位置具有显著优势。首先，精子细胞取材方便，来源广泛。其次，由于精子细胞是单倍体细胞，缺失的染色体片段拷贝数为0，与二倍体缺失区域拷贝数为1相比能更明确地显示缺失的区域，从而高效定位染色体断点位置，同时，单精子细胞测序周期更短，测序深度较低，成本更低，能够更迅速定位相互易位染色体断点位置。

本发明的优势包括：

（1）使用单倍体精子细胞，染色体缺失的部分没有比对的reads，通过阈值划分后拷贝数为0，拷贝数变化在数据计算上更加明显，获得相互易位染色体断点的位置更加精确。

（2）通过筛选预处理精子细胞数据，可以消除个别样本的离散偏倚，同时可以筛选出最适合用来分析断点位置的精子细胞。

（3）低深度测序（0.3×）即可达到高分辨率，从而精确定位断点位置，周期短，成本低。

（4）动态划分窗口，划分窗口的大小以5kb水平递增，同时要求频数之比大于2/3，可以在确保断点位置准确性的同时达到最精确的程度。

本发明提出的技术方案的关键技术点包括：

（1）通过对样本的标准差（SD）进行样本筛选，减少高离散度样本对结果造成的偏倚。

（2）通过筛选染色体存在缺失的样本达到更高的准确度，减少后续区域统计的离散程度。

（3）动态划分窗口，划分窗口的大小每次增加5kb，在确保断点位置准确性的同时达到最高的精确度。

以上描述了本发明优选实施方式，然其并非用以限定本发明。本领域技术人员对在此公开的实施方案可进行并不偏离本发明范畴和精神的改进和变化。

Claims (10)

1. 一种鉴定染色体相互易位断点位置的系统，其特征在于，所述的鉴定染色体相互易位断点位置的系统包括单细胞扩增及测序单元、测序数据清洗比对单元、计算单元、归一化单元、样本筛选单元和断点位置确定单元，所述的单细胞扩增及测序单元用于获得相互易位携带者的精子细胞，进行单细胞DNA扩增，扩增后对DNA进行纯化，继而构建二代测序文库，获得测序数据；所述的测序数据清洗比对单元用于对获得的测序数据去除二代测序接头序列和低质量碱基，保留序列长度大于36 bp的测序数据；将处理后的测序数据比对至UCSC hg38人类参考基因组；对比对后的测序数据去除结果中的低比对质量序列和PCR重复序列，得到唯一比对、去重后的测序数据；所述的计算单元用于计算基因组每500kb的窗口中比对的测序数据数目；所述的归一化单元用于归一化测序深度，计算每个窗口矫正后的测序数据数目；所述的样本筛选单元用于在得到归一化处理的矫正测序数据数目后，计算每个样本的离散程度，以挑选离散程度满足要求的单倍体样本，在此基础上进一步筛选出具有染色体部分缺失的单倍体样本用于下一步计算；所述的断点位置确定单元用于动态划分窗口并确定相互易位染色体断点位置，在500kb窗口基础上完成样本离散程度筛选和具有染色体缺失样本的筛选后，以高分辨率5kb窗口为基准，确定染色体易位断点的最优精确判定位置，当在此分辨率下无法确定最优位置时，每次将窗口大小增加5kb，直至得到断点的最优判定位置。

2.根据权利要求1所述的鉴定染色体相互易位断点位置的系统，其特征在于，所述的测序数据清洗比对单元使用trim_galore软件默认参数对获得的测序数据去除二代测序接头序列和低质量碱基，保留序列长度大于36 bp的测序数据；使用BWA软件以默认参数将处理后的测序数据比对至UCSC hg38人类参考基因组；对比对后的测序数据使用samtools软件去除结果中的低比对质量序列和PCR重复序列，得到唯一比对、去重后的测序数据。

3. 根据权利要求1所述的鉴定染色体相互易位断点位置的系统，其特征在于，所述的计算单元在得到比对结果后，将参考基因组划分为500kb大小的窗口，根据序列在基因组上的位置，使用samtools 软件统计每个样本在基因组上每个500kb窗口的比对测序数据数目。

4. 根据权利要求1所述的鉴定染色体相互易位断点位置的系统，其特征在于，所述的归一化单元通过归一化不同窗口所得的测序数据数目来消除由于扩增和建库所带来的偏倚；具体为：i为一批测序中第i个样本，j为基因组中第j个窗口， 为i样本的测序深度，为i样本所测碱基总数，为人类基因组序列长度，为i样本中j窗口比对的测序数据数目，为i样本j窗口的归一化测序深度，n为样本总数，为j窗口的归一化因子，为i样本j窗口矫正后的比对测序数据数目；

①计算每个样本中每个窗口的归一化测序深度：

 =  ；

 =  ；

②计算每个窗口的归一化因子：

 =  ；

③计算每个样本中每个窗口的矫正后的测序数据数目：

 =  。

5.根据权利要求1所述的鉴定染色体相互易位断点位置的系统，其特征在于，i为一批测序中第i个样本，j为基因组中第j个窗口，N为该样本全部染色体划分的窗口总数，i样本的离散程度用表示，为i样本j窗口矫正后的比对测序数据数目，为i样本每个窗口矫正后测序数据数目的平均值；SD的值表示样本的离散程度，与样本质量呈反比；所述的样本筛选单元计算每个样本的离散程度的具体过程为：

①计算每个样本矫正后拷贝数的平均值：

 = ；

②计算每个样本的SD值，表示每个样本的离散程度：

 =  ；

③筛选出SD＜0.25的单倍体样本。

6.根据权利要求1所述的鉴定染色体相互易位断点位置的系统，其特征在于，所述的样本筛选单元还用于进一步根据矫正测序数据数目所得拷贝数变异结果，筛选出易位染色体有缺失的单细胞测序样本。

7.根据权利要求1所述的鉴定染色体相互易位断点位置的系统，其特征在于，所述的断点位置确定单元首先将参考基因组划分为5kb大小的窗口，使用samtools统计所选样本每个窗口的比对测序数据数目，使用所述的归一化单元对每个样本每个窗口的测序数据数目进行矫正，得到每个样本每个窗口的矫正后的测序数据数目；由于所用样本是单倍体细胞，因此使用0表示缺失，1表示正常单倍体，2表示重复；将所有样本的每个窗口归为0、1、2三个拷贝数数值，如果矫正后测序数据数目小于0.4，则将该窗口归为拷贝数0；如果矫正后测序数据数目在0.4到1.4之间，则将该窗口归为拷贝数1；如果矫正后测序数据数目大于1.4，则将该窗口归为拷贝数2；当拷贝数由连续的0转变为连续的1或2时，0拷贝转变为1或2拷贝的窗口即为相互易位染色体断点所在的区域，根据拷贝数变化记录断点位置所在窗口。

8.根据权利要求1所述的鉴定染色体相互易位断点位置的系统，其特征在于，所述的断点位置确定单元在划分窗口大小为5kb的基准上，进行样本间断点区域的比较，对样本所得的断点区域进行汇总，计算出现次数最多的区域的频数，当该频数占总频数2/3以上时，确定该区域为相互易位染色体的断点位置。

9.根据权利要求1所述的鉴定染色体相互易位断点位置的系统，其特征在于，所述的断点位置确定单元在划分窗口为5kb大小的分辨率下无法判定断点位置时，逐步增加窗口大小，每次将窗口增加5kb，直至区域最大频数与总频数的比值达到2/3以上时，确定该区域为断点位置。

10.一种根据权利要求1至9任一项所述的系统鉴定染色体相互易位断点位置的方法，其特征在于，包括以下步骤：

第一步、获得相互易位携带者的单倍体细胞，进行单细胞DNA扩增，扩增后对DNA进行纯化，继而构建二代测序文库，获得测序数据；

第二步、对获得的测序数据去除二代测序接头序列和低质量碱基，保留序列长度大于36 bp的测序数据；将处理后的测序数据比对至UCSC hg38人类参考基因组；对比对后的测序数据去除结果中的低比对质量序列和PCR重复序列，得到唯一比对、去重后的测序数据；

第三步、计算基因组每500kb的窗口中比对的测序数据数目；

第四步、归一化测序深度，计算每个窗口矫正后的测序数据数目；

第五步、在得到归一化处理的矫正测序数据数目后，计算每个样本的离散程度，以挑选离散程度满足要求的单倍体样本，在此基础上进一步筛选出具有染色体缺失的样本用于下一步计算；

第六步、动态划分窗口并确定相互易位染色体断点位置，在500kb窗口基础上完成样本离散程度筛选和具有染色体缺失样本的筛选后，以高分辨率5kb窗口为基准，确定染色体易位断点的最优精确判定位置，当在此分辨率下无法确定最优位置时，每次将窗口大小增加5kb，直至得到断点的最优判定位置。

Images