Cnv检测方法和装置的制造方法

Cnv检测方法和装置的制造方法
【专利摘要】本发明提供了一种检测CNV的方法，包括：获取目标个体的基因组测序结果；将所述测序结果与参考序列比对，获得比对结果，所述参考序列包括多个窗口；基于所述比对结果，计算每个窗口的初始比对率，窗口的初始比对率＝比对上所述窗口的读段数目/比对上所有窗口的读段数目的平均值，所述比对上所有窗口的读段数目的平均值＝比对上所有窗口的读段总数/窗口个数；合并初始比对率无显著差异的多个相邻窗口，定义合并后的多个相邻窗口为一次区域，剩余的每个独立窗口也分别称为一次区域；基于所述一次区域的比对率与预定比对率不相等，判定所述一次区域存在CNV。
【专利说明】
GNV检测方法和装置
技术领域
[0001] 本发明涉及生物信息领域，具体的，本发明涉及检测CNV的方法和装置。
【背景技术】
[0002] 单细胞测序技术是利用二代测序技术对单个细胞的微量核酸进行测序。该项技术 主要包括单细胞分离，单细胞核酸的提取与扩增和测序三部分。单细胞测序作为一项革命 性技术，近几年来在科研及生物医药领域得以广泛应用。例如，对肿瘤单细胞进行测序，揭 示肿瘤单细胞层面的异质性，推演肿瘤的演化过程；无创产前诊断；组装不能培养的微生 物基因组；痕量细胞（可以运用于法医学等）基因组的获取；单细胞技术也被引入到胚胎 植入前诊断等等。单细胞测序技术解决了痕量细胞基因组获取的难题，为疾病发病机制和 诊断学研究提供新的方法。
[0003] 在单细胞研究中，单细胞拷贝数变异（Copy Number Variants, CNV)扮演了很重要 的角色。CNV-般指染色体上大于1Kb的片段发生丢失或重复的现象。CNV是一种广泛存 在动植物基因组中的遗传多态，它的突变频率远高于SNP，基因组研究已证明CNV与部分人 类疾病相关，比如与肿瘤、肥胖、自闭症、自身免疫性疾病和系统性红斑狼疮等多种复杂疾 病相关。在肿瘤的异质性和进化研究中，肿瘤单细胞的CNV的检测，通过比较单细胞之间， 以及单细胞与对应组织的CNV的差异，揭示了肿瘤在单个细胞层面的异质性，为肿瘤演化 推演提供了依据；无创产前诊断，则需要对微量DNA进行检测是否有染色体非整倍体变异 (CNV的一种）而导致的21-三体综合征（47, +21)、18-三体综合征（47, +18)、13-三体综 合征（47,+13)等；胚胎植入前的诊断和筛查，需要对单个生殖细胞或胚胎细胞进行相关检 测分析；法医取证样本（痕量的血液，精液等），需要进行痕量细胞的分析等等。综合来看， 当前生物医学领域对于痕量细胞，甚至单个细胞的大片段CNV的检测提出了需求和挑战。
[0004] 现有的CNV检测方法大多是针对组织测序数据的，如CNV-seq，PenCNV，CNAseg和 Readdepth等。单细胞测序数据，尤其是低深度测序数据，具有低基因组覆盖度和高扩增偏 向性，在基因组的不同区域短序列比对波动很大，这些CNV检测方法并不太适合对单个细 胞拷贝数变异的检测。

【发明内容】

[0005] 本发明旨在至少解决上述问题至少之一或者提出至少一种商业选择。
[0006] 依据本发明的一方面，本发明提供一种检测CNV的方法，所述方法包括以下步骤： 获取目标个体的基因组测序结果，所述测序结果包括多个读段；将所述测序结果与参考序 列比对，获得比对结果，所述参考序列包括多个窗口，所述比对结果包括比对上每个所述窗 口的读段的数目；基于所述比对结果，计算每个窗口的初始比对率，窗口的初始比对率=比 对上所述窗口的读段数目/比对上所有窗口的读段数目的平均值，所述比对上所有窗口的 读段数目的平均值=比对上所有窗口的读段总数/窗口个数；合并初始比对率无显著差异 的多个相邻窗口，定义合并后的多个相邻窗口为一次区域，剩余的每个独立窗口也称为一 次区域；基于所述一次区域的比对率与预定比对率不相等，判定所述一次区域存在CNV，所 述一次区域的比对率为所述一次区域包含的窗口的初始比对率的均值，所述预定比对率为 所有窗口的比对率中频率最高的窗口的比对率，所述窗口的比对率为其所在的一次区域的 比对率。在本发明的一个实施例中，所述基因组获自所述目标个体的单个细胞。通过构建 单细胞的基因组测序文库，并对所述测序文库进行序列测定获得所述测序结果。任选的， 构建所述测序文库包括对所述基因组进行简并寡核苷酸引物PCR，多重置换扩增和/或多 次退火环状循环扩增，以获得足够建库的核酸量和/或足够上机测序的核酸量。序列测定 可以利用现有测序平台，包括但不限于CG (Complete Genomics)、Illumina/Solexa、Life Technologies ABI SOLiD和Roche 454测序平台，可根据所选用的测序平台进行相应的测 序文库制备，可选择单端或双端测序，由此获得的测序结果由多个短序列组成，将各个短序 列称为读段。所说的比对可以利用已知比对软件进行，比如利用Bowtie、SOAP、BWA和/或 TeraMap等进行。在本发明的一个实施例中，只利用所述比对结果中的比对到所述参考序列 唯一位置的读段进行比对率的计算，以提高数据准确性以利于提高CNV检测的准确性。
[0007] 所称的窗口可以预先确定，也可以在进行目标个体检测时同时确定。在本发明的 一个实施例中，窗口是预先确定的。所述窗口的确定包括：将短序列集与参考序列比对， 确定比对上所述参考序列的短序列的起始位置，所述短序列集包括多个短序列；在所述参 考序列上划定窗口，使每个所述窗口包含相同数目的所述起始位置，任选的，所述窗口之 间没有重叠。在比对过程中，根据比对参数的设置，一条短序列最多允许有m个碱基错配 (mismatch)，m优选为1或2,若一条短序列中有超过m个碱基发生错配，则视为该短序列无 法比对到参考序列。所称的起始位置为比对上参考序列的各个短序列的起始碱基与参考序 列的匹配位置，当有多个短序列的起始碱基比对到参考序列同一位置时，只记录一次，即记 录所说的起始位置为一个。这里，所称的短序列的起始碱基，即涉及的短序列的方向，是以 参考序列的方向为参照的，例如，将一条短序列上的匹配到参考序列最前位置（位置编号 最小）的碱基称为该短序列的起始碱基。使每个所述窗口包含相同的起始位置数目，而不 限制其包含的没有读段匹配的位点数目，所以一般的各个窗口的大小不一样，这样，有利于 减少单细胞基因组扩增带来的偏向性。据此，在这一构思下，也可以使每个窗口包含相同数 目的特定位置来进行窗口划定，所说的特定位置为比对上参考序列的各个短序列的相同位 置碱基与参考序列的匹配位置，例如，使所说的特定位置为比对上参考序列的各个短序列 的末端碱基与参考序列的匹配位置。
[0008] 所称短序列集可来自模拟序列集和/或测序结果，这里所说的测序结果可以是自 己测定的人核酸的测序数据，也可以是他人公开的核酸样本的测序结果，核酸可以是基因 组DNA也可以是游离DNA。较佳的，使所说的模拟序列集中的模拟序列在比对到参考基因组 上能有相对均匀的分布，在本发明的一个实施例中，模拟序列可以这样获得：从所述参考序 列的长度为Q的染色体的一端的碱基开始，拷贝所述染色体的P个碱基，以获得第一条模拟 序列，沿所述染色体的另一端方向移动一个碱基拷贝所述染色体的P个碱基，以获得第二 条模拟序列，沿所述染色体的另一端方向移动两个碱基拷贝所述染色体的P个碱基，以获 得第二条模拟序列，依此获得第Q _P+1条模拟序列，所述第Q_P+1条模拟序列的末端碱基与 所述染色体的另一端的碱基重合，其中，P为模拟序列的长度，较佳的，P多10。在本发明的 一个实施例中，所述窗口之间没有重叠且所述窗口总数不大于100, 〇〇〇。窗口的大小的设置 可以基于CNV检测精度调整，在人参考基因组大小确定的情况下，窗口的大小与窗口数目 成反比。在该实施例中，窗口的总数不少于10, 〇〇〇且不大于100, 〇〇〇,而且之间没有重叠， 利于准确检测出一般定义的不小于1K的CNV。
[0009] 在本发明的一个实施例中，所述目标个体为人类，人类为二倍体生物，其染色体组 数为2,优选对应的所述参考序列为人参考基因组的至少一部分，例如为HG19，HG19可以从 NCBI数据库获取，或者为所有窗口构成的参考序列。在本发明的另一个实施例中，以N替 代所述人参考基因组的Y染色体的拟常染色体区的每个碱基，N表示A、T、C和G中的任一 种，这样，有利于避免性染色体的拟常染色体区域CNV检测的假阳性。
[0010] 在本发明的一个实施例中，在合并初始比对率无显著差异的多个相邻窗口之前， 利用比对率-GC含量的关系对每个所述窗口的初始比对率进行GC校正，获得各个窗口的校 正比对率，以消除或减少GC含量对测序结果、比对率的影响，并且以窗口的校正比对率替 代所述窗口的初始比对率进行后续检测，例如，所述一次区域的比对率变为其包含的所有 窗口的校正比对率的均值，而在确定所称的预定比对率时，将一次区域的比对率赋给其所 包含的窗口，即处于同一一次区域的各个窗口的比对率均相等，为其所在的一次区域的比 对率，这样，对所有窗口的比对率进行计数，确定各个比对率出现的次数，将出现次数最多 即频率最高的窗口比对率定为所称的预定比对率。所述比对率-GC含量的关系可以预先利 用对照样本的测序数据来建立、保存，用以校正各个待测样本测序结果，优选的对照样本为 与目标个体同物种的组织样本，也可以在检测目标样本时同时利用目标样本基因组的测序 结果来建立。在本发明的一个实施例中，直接利用检测所需的目标样本的测序结果来建立 比对率-GC含量的关系，所说的比对率-GC含量的关系的建立如下：获得至少一个样本的核 酸的测序数据，所述测序数据由多个读段组成；将所述测序数据与参考序列进行比对，获得 比对结果，所述参考序列包括多个窗口，所述比对结果包括比对上每个所述窗口的读段的 数目；计算每个所述窗口的初始比对率，窗口的初始比对率=比对上所述窗口的读段数目 /比对上所有窗口的读段数目的平均值，所述比对上所有窗口的读段数目的平均值=比对 上所有窗口的读段总数/窗口个数；基于多组的窗口的初始比对率和该窗口的GC含量的数 值，利用二维回归分析法建立所述比对率-GC含量的关系。在本发明的一个实施例中，利用 的二维回归分析法为局部加权回归散点平滑法（Lowess)。
[0011] 所称的初始比对率无显著差异的多个相邻窗口指初始比对率无实质差异的相邻 窗口，例如，由于初始比对率或者校正比对率是围绕着"1"波动的一组数值，可以以1或 者以1± 1*10%为有无实质差异的界限，如相邻的、校正比对率都在0.9以下，或0.90~ 1. 10,或1. 10以上的窗口为无实质差异的窗口。在本发明的一个实施例中，所述合并初 始比对率无显著差异的多个相邻窗口为合并符合以下描述的相邻窗口，多个相邻窗口 的校正比对率都大于1或者都小于1。进一步的，为确定所检测的CNV的大小和确切发 生位置（断点），该方法还包括：确定所述一次区域中的二次区域，包括，（1)基于公式
算所述一次区域中的子区域Μ与该一次区域中的所有其它窗口的比对 率的差异，获得所有Zi j，取4= max i i j~| |，（2)将Ζ。与第一临界值比较，当Ζ。超过 第一临界值时，相应的子区域Μ为所述二次区域，所述二次区域为CNV区域，所述二次区域 的边界即为CNV的发生位置，（3)去除所述一次区域中的二次区域，更新i、j和n，进行步 骤（1)和（2)，直至无 Z。超过第一临界值；其中，i和j为所述一次区域中的窗口的编号，η 为所述一次区域中的窗口的数目，所述子区域Μ为所述一次区域中的第i+1个窗口到第j 个窗口之间的区域，氏为所述一次区域中第i个窗口的校正比对率，所述第一临界值是Z ^ 分布中的第一预定概率的概率密度，所述第一预定概率多95%，1 < i < j < mSii L+… +民，S,= R片…+R,，Sn= R彳…+Rn。假设子区域Μ为正常非变异区域，Zy分布指Z u服从标 准正态分布，第一预定概率和第一临界值一一对应，一般统计书籍都包含第一预定概率和 第一临界值对应的表格供查阅。在本发明的一个实施例中，当Z。落入拒绝域，即Z。超过第 一预定概率例如为99. 9%对应的第一临界值，可知发生了小概率事件，否定原假设，即子区 域Μ为变异区域。上述过程，依据窗口的校正比对率的同向性，即都大于1或者都小于1， 对窗口进行合并，获得大的一次区域，再在各个一次区域中进行循环判断以确定其中的CNV 的发生边界，即从其中确定二次区域，这样同时在多个一次区域中并行确定二次区域，利于 快速检测CNV。在本发明的一个实施例中，所述二次区域的比对率为所述二次区域包含的所 有窗口的校正比对率的均值。在本发明的一个实施例中，该方法还包括，基于比较所述二次 区域的比对率与所述预定比对率的大小，判定CNV的类型，其中包括，当所述二次区域的比 对率大于所述预定比对率时，判定所述二次区域为拷贝数增加区域，当所述二次区域的比 对率小于所述预定比对率时，判定所述二次区域为拷贝数减少区域。在本发明的另一个实 施例中，利用以下公式计算所述二次区域的拷贝数，二次区域的拷贝数=该二次区域的比 对率/预定比对率*目标个体的染色体组数，所述二次区域的比对率为其包含的所有窗口 的校正比对率的均值。
[0012] 无显著差异，也可以指统计学上的对数据差异性的评价一一差异无显著性，例如 设定预定概率，通常预定概率可以设为不小于95 %，对相邻多个窗口的校正比对率进行统 计检验，例如可以利用ζ检验或t检验，多个校正比对率之间的差异无显著性（ρ > 0. 05)， 即认为达到所说的无显著差异。在本发明的一个实施例中，所述合并初始比对率无显著差 异的多个相邻窗口为合并满足如下描述的相邻窗口一一校正比对率的差异无统计意义，使 合并得的一次区域为CNV区域。合并初始比对率无显著差异的多个相邻窗口具体包括：(a)
计算区域N与其它所有窗口的比对率的差异，获得所有Zxy， 取Zb= max i < x < y J Zxy |，（b)将Zb与临界值比较，当Z b超过所述临界值时，相应的区域N 为所述一次区域，（c)将所述一次区域去除，更新x、y和w，进行步骤（a)和（b)，直至无 Zb 超过所述临界值，其中，X和y为窗口的编号，w为窗口总数，所述区域N为第x+1个窗口到 第y个窗口之间的区域，Rx为第X个窗口的校正比对率，所述临界值是Z xy分布中的预定概 率的概率密度，所述预定概率彡95%，1彡X < y彡w，Sx=心+…+艮，Sy= RJ…+Ry，Sw= ^+???+1。所说的Zxy分布为Z xy服从标准正态分布，预定概率与临界值一一对应。在本发 明的一个实施例中，假设区域N为正常非变异区域，当Z b落入拒绝域，即Z b超过预定概率例 如为99. 9%对应的临界值，可知发生了小概率事件，否定原假设，即区域N为变异区域。上 述过程，基于对所有窗口进行循环判断确定CNV的发生边界，确定出的一次区域即为CNV区 域。在本发明的一个实施例中，该方法还包括：基于比较所述一次区域的比对率与所述预 定比对率的大小，判定所述CNV的类型，其中包括，当所述一次区域的比对率大于所述预定 比对率时，判定所述一次区域为拷贝数增加区域；当所述一次区域的比对率小于所述预定 比对率时，判定所述一次区域为拷贝数减少区域。在本发明的另一个实施例中，该方法还包 括：利用以下公式计算所述一次区域的拷贝数，一次区域的拷贝数=该一次区域的比对率 /预定比对率*目标个体的染色体组数，所述一次区域的比对率为其包含的所有窗口的校 正比对率的均值。
[0013] 利用上述的本发明一方面的或者任一【具体实施方式】中的CNV检测方法，能够解决 上述现有的CNV检测流程中，存在一些不足，例如现有方法中的采用固定长度的窗口，不能 很好地解决单细胞测序中全基因组扩增所带来的偏性问题以及重复序列问题，不能很好的 用于二倍体单细胞的CNV的检测等。上述的本发明一方面的或者任一【具体实施方式】中的 CNV检测方法，非常适用于基于单细胞测序数据的CNV检测，特别是基于单细胞低深度测序 的CNV检测，对于不同测序平台采用不同扩增方法进行单细胞测序或组织测序的数据都是 有效的，适用性广泛。在不同测序平台使用不同全基因组扩增方法进行单细胞测序时，本发 明的方法在检测CNV的敏感性和特异性都很好，尤其是基于成环循环扩增技术（MALBAC)的 Proton平台的测序数据。而且，利用本发明的方法的检测结果具有高重复性，结果可信。与 现有的CNV检测方法比较，本发明的方法采用长度变化的窗口，有利于保持所有窗口比对 上的短序列数的平均值的稳定性，还可以避免重复序列区域带来的影响，使得CNV检测更 准确。
[0014] 依据本发明的另一方面，本发明提供一种检测CNV的装置，所述装置能够用以执 行或完成上述本发明一方面或者任一【具体实施方式】中的CNV检测方法，所述装置包括：数 据输入单元，用以接收数据；数据输出单元，用以输出数据；处理器，用以执行计算机可执 行程序，执行所述计算机可执行程序包括实现上述本发明一方面或者任一【具体实施方式】中 的CNV检测方法；以及，存储单元，用以存储数据，其中包括所述计算机可执行程序。所说的 计算机可执行程序可以保存在存储介质中，所称存储介质可以包括：只读存储器、随机存储 器、磁盘或光盘等。本发明还提供一种计算机可读存储介质，其用于存储供计算机执行的程 序，所述程序的执行包括完成前述本发明一方面的或者其任一【具体实施方式】中的CNV检测 方法。前述对本发明的CNV检测方法的优点和技术特征的描述也适用于该CNV检测装置和 计算机可读存储介质，在此不再赘述。
【附图说明】
[0015] 本发明的上述和/或附加的方面和优点从结合下面附图对实施方式的描述中将 变得明显和容易理解，其中：
[0016] 图1是本发明的一个【具体实施方式】中的窗口合并后的每个窗口的ratio的密度分 布图；
[0017] 图2是本发明的一个【具体实施方式】中的基于MDA扩增的CG平台单细胞测序数据 检测CNV的结果不意图；
[0018] 图3是本发明的一个【具体实施方式】中的基于MDA扩增的Proton平台单细胞测序 数据的CNV检测的结果示意图；
[0019] 图4是本发明的一个【具体实施方式】中的基于MALBAC扩增的Proton平台单细胞测 序数据的CNV检测的结果示意图。
【具体实施方式】
[0020] 以下对本发明方法的一般步骤或者相关信息的获取方式进行介绍。
[0021] 1.首先从 UCSC 的网站（http://hgdownload. cse. ucsc. edu/goldenPath/hgl9/ bigZips/)下载hgl9参考基因组的序列文件chromFa. tar. gz。在此，CNV检测方法将只使 用那些恰好比对到参考基因组一个位置的短序列，我们将Y染色体上拟常染色体区的序列 用N代替。对于Y染色体，后续过程使用的都是这个经修改的版本。拟常染色体区是X染 色体和Y染色体间唯一可发生互换的位置，这也是它的名称由来，由于会发生互换的是常 染色体，X染色体和Y染色体一般是没有互换的现象，只有在拟常染色体区反常地出现了互 换，导致男性和女性都带有两个该区域基因的复本。这使拟常染色体区的基因表达类似常 染色体，而非性染色体的伴性遗传模式，因而得名。
[0022] 2.确定每一个检测窗口（window)的大小。
[0023] a)针对Proton测序平台的下机数据，可采用单端模拟数据来划分。以hgl9参考 基因组为基准，模拟单端测序短序列，从基因组染色体的第一个碱基开始，每50个碱基为 一条读段（reads)，并为其生成ID和质量值生成fastaq格式。然后依次往后移一个碱基， 直到短序列的末端为染色体的最后一个碱基。使用bowtie把模拟数据对到参考基因组上， 结果只保留那些唯一比对的短序列（即去除可重复比对上的短序列），使用samtools把比 对结果转换为BAM格式。
[0024] 比对参数可设置为：bowtie-S-t-n 2_e 7〇-m 1 -best - strata，后续单细胞测序 数据的比对参数也可以一样，参数意义为：-n 2表示高保真区域内错配数不能超过2个，-e 70表示错配位点质量值不能超过70, 一best报告文件中，每个短序列的匹配结果将按匹配 质量从高到低排序，一srtata与一best -起使用报告质量最高的那部分，-m 1表示报告所 有比对的短序列。
[0025] b)对于CG平台的下机数据，通过对正常人的细胞系的一团细胞DNA的测序数据来 划分。CG的下机数据，进行流程信息分析，例如利用Teramap软件比对到参考基因组，然后 将比对结果的格式转换成BAM格式结果。
[0026] 不同平台的数据最后都可使用软件samtools去除重复的短序列。记录参考基因 组上每一个被短序列起始碱基覆盖的位置，并把这些位置划分到10, 000至100, 000个窗 口，每一个窗口内包含的位置数目完全相同，但其区间长度是变化的。然后分别计算每个窗 口所包含参考基因组序列的GC含量。
[0027] 3.提取单个细胞的DNA，进行全基因组扩增，然后建库上机测序，得到下机数据， 并进行相应分析处理得到bam格式比对的结果。
[0028] a)Proton平台，其下机的数据（BAM格式），我们使用BEDTools转换为FASTQ格 式数据，然后使用Trimmomatic软件对长于50bp的短序列从3'端截取有效长度（50bp加 上与全基因组扩增方法的primer等长的短序列序列），如多重退火和成环循环扩增技术 (MALBAC)的引物为35bp，有效长度为85bp，同时过滤掉长度小于有效长度的短序列。使用 bowtie把截取后的短序列比对到参考基因组，并用samtools view转换成bam文件排序后 去除重复短序列。
[0029] b)CG平台，使用其平台研发的分析流程把下机短序列数据与参考基因组hgl9进 行比对。然后把比对结果转成BAM格式并进行排序，samtools的单端模式去除重复短序列。
[0030] 4.对每一个窗口中比对上的短序列进行统计计数并进行标准化处理，即计算每个 窗口的比对率（ratio)=比对上的短序列数目/所有窗口比对上的短序列数的平均值。
[0031] 5.使用L0WESS算法确定的ratio-GC含量关系对每个窗口中标准化处理后得到的 ratio进行GC校正，获得校正ratio。
[0032] 6.每个样本根据所有窗口校正后的ratio值，可使用CBS segment软件对窗口进 行合并形成无重叠的区域（segment)并计算其ratio值，把此ratio值赋给区域（segment) 内的每个窗口。具体的包括，（a)
计算区域N与其它所有窗口 的比对率的差异，获得所有Zxy，其中，区域N为第x+1个窗口到第y个窗口之间的区域，Zxy 呈标准正态分布，取Zb= maXl<x<y<w|Zxy|，（b)将Z b与临界值比较，当Zb超过临界值时， 相应的区域N为预计的窗口合并区域，即区域N为发生CNV的区域（c)将窗口合并区域去 除，更新X、y和《，进行上述两步骤（a)和（b)，直至无 Zb超过临界值，即循环划分窗口，直 到窗口不能再进行合并；其中，X和y为窗口的编号，w为窗口总数，所述R x为第X个窗口的 校正比对率，所述临界值是Zxy分布中的预定概率的概率密度，所述预定概率多95%，1 < X < y彡w，Sx= R片…+RX，Sy= R彳…+Ry，Sw= R片…+RW。所说的Zxy分布为Z xy服从标准正 态分布，预定概率与临界值一一对应。上述过程可理解为，假设区域N为正常非变异区域， 当Z b落入拒绝域，即Z b超过99. 9%对应的临界值，可知发生了小概率事件，否定原假设，即 区域N为变异区域。各合并区域的ratio为其所包括窗口的校正ratio的均值，然后把该 合并区域的ratio值赋值给其包括的所有窗口，即为窗口的比对率。
[0033] 接着，对所有窗口的ratio画密度曲线分布图，如图1所示。对于近二倍体细胞或 者二倍体是所有倍型的众数的细胞，密度分布图中最大峰值对应的ratio值则为该细胞拷 贝数为2的ratio值。
[0034] 7.把每个区域的ratio除以拷贝数为2的ratio,再乘以2,则得到每个区域的拷 贝数。
[0035] 8.计算CNV检测的敏感性和特异性。敏感性=LT/LC，特异性=LT/L，其中，L :指 单细胞测序找到的CNV(彡1Mb)的总长度，LC :表示组织测序中找到的CNV(彡1Mb)的总长 度，LT :表示单细胞测序和组织测序共同找到的CNV(彡1Mb)的总长度。
[0036] 以下结合具体个体样本对依据本发明的检测方法及检测结果进行详细的描述。下 面示例，仅用于解释本发明，而不能理解为对本发明的限制。在本发明的描述中，"一次"、 "二次"等为指代或描述方便，不能理解为有顺序关系或者相对重要性指示，除非另有说明， "多个"的含义是两个或两个以上。
[0037] 除另有交待，以下实施例中涉及的未特别交待的试剂、序列（接头、标签和引物）、 软件及仪器，都是常规市售产品或者公开的，比如购自Illumina公司的hise q2000测序平 台建库相关试剂盒等。
[0038] 实施例一：基于MDA扩增的CG平台低深度测序数据的CNV检测方法测试
[0039] 随着高通量测序技术的蓬勃发展，以Complete Genomics (CG)、IlluminaSolexa 和Roche454为代表的二代测序，以及三代测序技术（即单分子测序技术）所包括的 HelicosGenetic Analysis System、单分子实时测序技术（SMRT)和纳米孔单分子测序技 术等各种测序技术已成为单细胞组学研究的重要工具。CG平台作为一种专注于人类基因 组的二代测序技术，能够完整并精确地测序人类全基因组，其测序通量大，在业内被高度认 可。CG平台其主要包括测序平台，高通量过程自动化技术和完整的数据管理解决方案三个 部分，其测序平台包括DNA纳米阵列（DNANanoball arrays，DNB? arrays)和组合探针锚定 连接测序法（combinatorial probe-anchor ligation，cPAL?)，这两项技术的应用大大减 少试剂的消耗和缩短成像的时间。所以我们首先在CG平台、利用CG平台的下机数据对本 发明的CNV检测方法进行试验验证。
[0040] 从病人的胶质母细胞瘤组织中分离出了 3个单细胞，组织样本来自北京天坛医院 提供，提取每个单细胞的DNA并利用MDA全基因组扩增技术进行扩增，再进行文库构建，然 后在CG平台进行单细胞低深度测序。最后按本发明的方法进行单细胞CNV的检测分析。为 了验证本发明方法的CNV检测效率，我们提取了组织的DNA进行文库构建，然后在CG平台 进行全基因组测序，并使用CG的标准分析流程检测得到组织的CNV结果。3个单细胞样本 (P1-T2-SC#)和组织样本（P1-T2)检测到的CNV如图2所示，平行X轴的粗黑线表示各区域 的拷贝数，其大于2则说明该区域内发生拷贝数增加，小于2则说明此区域内发生拷贝数减 少，等于2则表示拷贝数正常，每个窗口的ratio值用散点表示。
[0041] 进一步对CNV检测方法的敏感性和特异性进行估计，敏感性=LT/LC，特异性= LT/L。估算5个样本的平均敏感性为91. 01%，特异性为74. 47%，结果如表1所示。
[0042] 表 1
[0043]
[0044] 然后，对基于CG测序平台MDA扩增的CNV检测方法的重复性进行统计，发现样本 间的重复性高于0. 7,结果详见表2。
[0045] 表 2
[0046]
[0047] 从敏感性、特异性以及重复性统计计算结果可以得出，本发明的CNV分析检测流 程的检测结果的有效性，其在CG测序平台上是可行的。
[0048] 实施例二：基于MDA扩增的Proton平台的单细胞低深度测序的CNV检测方法测试
[0049] 现有单细胞测序数据多由Illumina测序平台产出。尽管Illumina测序仪的测序 通量大，但是其上机测序时间周期长，测序成本高，这些会限制单细胞CNV检测分析的快速 发展。而一些研究往往对测序通量没有需求，相对的，对测序的时间和成本具有更高的需 求，此时Proton测序平台将是更好的选择。Proton测序平台运行速度快，测序周期只需几 个小时，测序成本低，更适合部署到医院或第三方检测机构，缩短检测时间，降低成本，从而 提尚检测效率。而对于Proton的单细胞测序CNV检测鲜有报道。
[0050] 提取来自北京大学肿瘤医院的人类胃腺癌细胞系（BGC823)5个细胞（MDA-2_ BGC#)，并用多重置换扩增（MDA)技术进行文库构建后在Proton平台进行单细胞低深度测 序。同时提取人类胃腺癌细胞系（BGC823) -团细胞（BGC)的DNA，进行常规文库构建后在 Proton平台进行测序。然后按我们的方案进行CNV的检测分析，在BGC823的5个单细胞样 本和组织样本（BGC)检测到的CNV如图3,每个区域的拷贝数（平行X轴的粗黑线）大于2 则为说明区域发生拷贝数增加，小于2为拷贝数减少，等于2则为拷贝数正常，三种拷贝数 变化区域内窗口的ratio值分别用不同灰度深度的散点表示。
[0051] 在整个基因组上多个染色体上检测到大片段的CNV，与细胞团的CNV检测结果保 持一致，验证了本发明方法检测CNV的有效性。
[0052] 然后，根据五个单细胞和一团细胞CNV检测结果，我们进一步对检测CNV方法的敏 感性和特异性进行估计，敏感性=LT/LC，特异性=LT/L。估算5个单细胞样本的平均敏感 性为85. 86%，特异性为81. 18%，结果如表3所示。
[0053] 表 3
[0054]
[0055] 对基于Proton测序平台MDA扩增的CNV检测方法的重复性进行统计，发现样本间 的重复性高于0. 7,结果详见表4。
[0056] 表 4
[0057]
[0058] 从敏感性、特异性以及重复性统计计算结果可以得出，本发明的CNV检测流程分 析结果的有效性，其对于Proton测序平台MDA扩增方法的测序数据是可行的。
[0059] 实施例三：基于MALBAC扩增的Proton平台的单细胞低深度测序的CNV检测方法 测试
[0060] 提取5个人类胃腺癌细胞系细胞（BGC823)，用MALBAC全基因组扩增方法进行常规 文库构建后在Proton平台进行单细胞低深度测序；同时提取人类胃腺癌细胞系（BGC823) 一团细胞（BGC)的DNA在Proton平台进行测序。得到的下机数据按我们的方案进行CNV 的检测分析，在BGC823五个样本找到CNV，结果如图4所示，其中，横坐标表示染色体；右侧 纵坐标为5个单细胞样本及团细胞样本，左侧纵坐标为拷贝数，图上的黑色粗线代表划分 区域的ratio值，其值大于2的说明此区域拷贝数增加，小于2则说明区域内拷贝数减少， 等于2则说明区域内拷贝数正常。三种拷贝数变化区域分别用不同灰色深度的散点表示窗 口的 ratio 值。
[0061 ] 进一步对本发明的检测CNV方法的敏感性和特异性进行估计，敏感性=LT/LC，特 异性=LT/L。估算5个样本的平均敏感性为84. 72%，特异性为85. 18%，结果如表5。
[0062] 表 5
[0063]
[0064] 对基于Proton测序平台MALBAC扩增的CNV检测方法的重复性进行统计，发现样 本间的重复性高于〇. 92,详见表6。
[0065] 表 6
[0066]
[0067] 从敏感性、特异性以及重复性统计计算结果可以得出，本发明的CNV检测流程分 析结果的有效性，其对于Proton测序平台MALBAC扩增方法的测序数据是可行的。
【主权项】
1. 一种检测CNV的方法，其特征在于，包括W下步骤： 获取目标个体的基因组测序结果，所述测序结果包括多个读段； 将所述测序结果与参考序列比对，获得比对结果，所述参考序列包括多个窗口，所述比 对结果包括比对上每个窗口的读段的数目； 基于所述比对结果，计算每个窗口的初始比对率，窗口的初始比对率=比对上所述窗 口的读段数目/比对上所有窗口的读段数目的平均值，所述比对上所有窗口的读段数目的 平均值=比对上所有窗口的读段总数/窗口个数； 合并初始比对率无显著差异的多个相邻窗口，定义合并后的多个相邻窗口为一次区 域，剩余的每个独立窗口也分别称为一次区域； 基于所述一次区域的比对率与预定比对率不相等，判定所述一次区域存在CNV， 所述一次区域的比对率为所述一次区域包含的窗口的初始比对率的均值， 所述预定比对率为所有窗口的比对率中频率最高的窗口的比对率，所述窗口的比对率 为其所在的一次区域的比对率。2. 权利要求1的方法，其特征在于，所述基因组获自所述目标个体的单个细胞； 任选的，通过构建所述细胞的基因组测序文库，并对所述测序文库进行序列测定获得 所述测序结果； 任选的，构建所述测序文库包括对所述基因组进行简并寡核巧酸引物PCR，多重置换扩 增和/或多次退火环状循环扩增。3. 权利要求1的方法，其特征在于，所述参考序列为人参考基因组； 任选的，W N替代所述人参考基因组的Y染色体的拟常染色体区的每个碱基，N表示A、 T、C和G中的任一种。4. 权利要求1的方法，其特征在于，所述窗口的确定，包括， 将短序列集与参考序列比对，确定比对上所述参考序列的短序列的起始位置，所述短 序列集包括多个短序列，所述短序列集来自模拟序列集和/或测序结果； 在所述参考序列上划定窗口，使每个所述窗口包含相同数目的所述起始位置； 任选的，所述模拟序列集中的模拟序列的获取包括， 从所述参考序列的长度为Q的染色体的一端的碱基开始，拷贝所述染色体的P个碱基， W获得第一条模拟序列， 沿所述染色体的另一端方向移动一个碱基拷贝所述染色体的P个碱基，W获得第二条 模拟序列， 沿所述染色体的另一端方向移动两个碱基拷贝所述染色体的P个碱基，W获得第=条 模拟序列， 依此获得第Q-P+1条模拟序列，所述第Q-P+1条模拟序列的末端碱基与所述染色体的 另一端的碱基重合，其中， P为模拟序列的长度，P > 10 ; 任选的，所述窗口之间没有重叠； 任选的，所述窗口总数不大于100,000。5. 权利要求1的方法，其特征在于，在所述合并初始比对率无显著差异的多个相邻窗 口之前，利用比对率-GC含量的关系对每个所述窗口的初始比对率进行GC校正，获得各个 窗口的校正比对率， W所述窗口的校正比对率替代该窗口的初始比对率。6. 权利要求5的方法，其特征在于，建立所述比对率-GC含量的关系，包括， 获得至少一个样本的核酸的测序数据，所述测序数据由多个读段组成； 将所述测序数据与参考序列进行比对，获得比对结果，所述参考序列包括多个窗口，所 述比对结果包括比对上每个所述窗口的读段的数目； 计算每个所述窗口的初始比对率，窗口的初始比对率=比对上所述窗口的读段数目/ 比对上所有窗口的读段数目的平均值，所述比对上所有窗口的读段数目的平均值=比对上 所有窗口的读段总数/窗口个数； 基于多组的窗口的初始比对率和该窗口的GC含量，利用二维回归分析法建立所述比 对率-GC含量的关系； 任选的，所述二维回归分析法为局部加权回归散点平滑法。7. 权利要求5的方法，其特征在于，所述合并初始比对率无显著差异的多个相邻窗口 是指，合并满足W下的相邻窗口， 校正比对率的差异无统计意义。8. 权利要求7的方法，其特征在于，所述合并初始比对率无显著差异的多个相邻窗口， 包括， (a) 基于公^5 计算区域N与其它所有窗口的比对率的差异，获得 所有 Zxy,取 Zb= max Jz巧 I， (b) 将Zb与临界值比较，当Z b超过所述临界值时，相应的区域N为所述一次区域， (C)将所述一次区域去除，更新x、y和W，进行步骤（a)和化），直至无 Zb超过所述临界 值，其中， X和y为窗口的编号， W为窗口总数， 所述区域N为第X+1个窗口到第y个窗口之间的区域， Rx为第X个窗口的校正比对率， 所述临界值是Zyy分布中的预定概率的概率密度，所述预定概率> 95%， 1《X < y《W， Sx= R 1+. .. +Rx, Sy= R 1+. .. +Ry, S"= R 1+. . . +斬。9. 权利要求8的方法，其特征在于，还包括， 基于比较所述一次区域的比对率与所述预定比对率的大小，判定所述CNV的类型，其 中包括， 当所述一次区域的比对率大于所述预定比对率时，判定所述一次区域为拷贝数增加区 域， 当所述一次区域的比对率小于所述预定比对率时，判定所述一次区域为拷贝数减少区 域。10. 权利要求7-9任一方法，其特征在于，还包括， 利用W下公式计算所述一次区域的拷贝数， 一次区域的拷贝数=该一次区域的比对率/预定比对率*目标个体的染色体组数， 所述一次区域的比对率为其包含的所有窗口的校正比对率的均值。11. 一种检测CNV的装置，其特征在于，包括， 数据输入单元，用W接收数据； 数据输出单元，用W输出数据； 处理器，用W执行可执行程序，执行所述可执行程序包括完成权利要求1-10任一方 法；化S， 存储单元，用W存储数据，其中包括所述可执行程序。
【文档编号】C12Q1/68GK105986008SQ201510039685
【公开日】2016年10月5日
【申请日】2015年1月27日
【发明人】李甫强, 史旭莲, 谢国云, 鲁娜, 赵至坤, 蒋润泽, 梁瀚, 侯勇, 吴逵
【申请人】深圳华大基因科技有限公司