CN104221022B - 一种拷贝数变异检测方法和系统 - Google Patents

一种拷贝数变异检测方法和系统 Download PDF

Info

Publication number
CN104221022B
CN104221022B CN201280066929.3A CN201280066929A CN104221022B CN 104221022 B CN104221022 B CN 104221022B CN 201280066929 A CN201280066929 A CN 201280066929A CN 104221022 B CN104221022 B CN 104221022B
Authority
CN
China
Prior art keywords
sequence
window
sequence label
candidate
sample
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.)
Active
Application number
CN201280066929.3A
Other languages
English (en)
Other versions
CN104221022A (zh
Inventor
李旭超
陈盛培
陈芳
谢伟伟
汪建
王俊
杨焕明
张秀清
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
BGI Shenzhen Co Ltd
Original Assignee
BGI Shenzhen Co Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by BGI Shenzhen Co Ltd filed Critical BGI Shenzhen Co Ltd
Publication of CN104221022A publication Critical patent/CN104221022A/zh
Application granted granted Critical
Publication of CN104221022B publication Critical patent/CN104221022B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • 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
    • 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
    • G16B15/00—ICT specially adapted for analysing two-dimensional or three-dimensional molecular structures, e.g. structural or functional relations or structure alignment
    • 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
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • 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
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/10—Sequence alignment; Homology search
    • 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
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C99/00—Subject matter not provided for in other groups of this subclass
    • 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
    • C12Q2535/00—Reactions characterised by the assay type for determining the identity of a nucleotide base or a sequence of oligonucleotides
    • C12Q2535/122—Massive parallel sequencing
    • 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
    • C12Q2537/00—Reactions characterised by the reaction format or use of a specific feature
    • C12Q2537/10—Reactions characterised by the reaction format or use of a specific feature the purpose or use of
    • C12Q2537/165—Mathematical modelling, e.g. logarithm, ratio

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Organic Chemistry (AREA)
  • Biophysics (AREA)
  • Biotechnology (AREA)
  • General Health & Medical Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Theoretical Computer Science (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Zoology (AREA)
  • Wood Science & Technology (AREA)
  • Medical Informatics (AREA)
  • Evolutionary Biology (AREA)
  • Molecular Biology (AREA)
  • Microbiology (AREA)
  • Immunology (AREA)
  • Biochemistry (AREA)
  • General Engineering & Computer Science (AREA)
  • Genetics & Genomics (AREA)
  • Computing Systems (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioethics (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Epidemiology (AREA)
  • Evolutionary Computation (AREA)
  • Public Health (AREA)
  • Software Systems (AREA)
  • Artificial Intelligence (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
  • Apparatus Associated With Microorganisms And Enzymes (AREA)

Abstract

本发明公开一种基因组拷贝数变异检测方法和系统，涉及生物信息学技术领域。该方法包括：获得读序；根据读序确定序列标签；统计落入各个窗口的序列标签数目；对各个窗口的序列标签数目进行GC校正并根据以对照样本集修正的期望序列标签数目进行修正获得调整后的序列标签数目；选取显著性值较小的分界点为候选的CNV断点；每次剔除最不显著的候选CNV断点并更新被剔除的候选CNV断点左右两个候选CNV断点的差异显著性值，循环迭代，直至所有候选CNV断点的差异显著性值都小于终止阈值，从而确定CNV断点。本发明的方法和系统，具有临床可行性，在使用50M左右的数据情况下，可精确检测到0.5M的微缺失/微重复区域。

Description

一种拷贝数变异检测方法和系统
技术领域
本发明涉及生物信息学技术领域，尤其涉及一种拷贝数变异(Copy NumberVariation，CNV)检测方法和系统。
背景技术
CNV是基因组结构变异的一种，狭义的CNV通常指染色体上DNA片段的拷贝数发生了变化。基因组结构变异的类型及原因包括：1.缺失(末端缺失、中间缺失)；2.易位(相互易位、罗伯逊易位)；3.倒位；4.环状染色体；5.双着丝粒染色体；6.插入等。广义的CNV也包括例如染色体非整倍性和部分非整倍性的结构变异。
目前拷贝数变异的检测方法，主要有高分辨率染色体核型分析、FISH(荧光原位杂交)、Array CGH(比较基因组杂交)、MLPA(多重连接探针扩增技术)和PCR(PolymeraseChain Reaction，多聚酶链式反应)的方法等，其中遗传学诊断以FISH检测为金标准，可以有效地检测出大部分已知染色体缺失或重复。但这些方法普遍存在效率较低，尤其在全基因组水平全面扫描的情况下，资源消耗较大或无法检测未知CNV等缺点。
因此，开发一种新的拷贝数变异检测方法，对已知的位点进行鉴定，并对未知的位点进行发现性探索，已经迫在眉睫。
发明内容
本发明要解决的一个技术问题是提供一种拷贝数变异检测方法和系统，能够精确检测包括微缺失/微重复的拷贝数变异。
根据本发明的一个方面，提供一种拷贝数变异检测方法，其特征在于，包括：获取受试样品中至少一部分核酸分子的读序信息；根据读序信息确定唯一比对至(基因组)参考序列的序列标签；将基因组参考序列划分窗口,统计落入各个窗口的序列标签数目；对各个窗口的序列标签数目进行GC校正并根据以对照样本集修正的期望序列标签数目进行修正获得调整后的序列标签数目；以窗口的起点或终点作为分界点，计算两侧由调整后的序列标签数目组成数值群体的差异显著性值，选取显著性值较小的分界点为候选的CNV断点；对于每个CNV断点至在先CNV断点和在后CNV断点的两段序列，计算两段序列包含的窗口的调整后的序列标签数目组成的两个数值群体的差异显著性值，每次剔除最不显著的候选CNV断点并更新被剔除的候选CNV断点左右两个候选CNV断点的差异显著性值，循环迭代，直至所有候选CNV断点的差异显著性值都小于终止阈值，从而确定CNV断点。
可选地，还包括对受试样本中的至少一部分核酸分子进行测序以获得读序信息的步骤。
可选地，每个窗口具有相同的参考序列标签数目(reference unique reads)，或者，每个窗口具有相同的长度。
可选地，终止阈值根据由正常样本组成对照样本集获得。
可选地，对各个窗口的序列标签数目进行GC校正包括以下的步骤：将窗口按GC含量分组，基于组内序列标签的平均数与所有窗口的序列标签平均数获得修正系数，对该窗口序列标签数进行修正获得GC校正的序列标签数。
可选地，对照样本集修正的期望序列标签数目由以下方法获得：计算对照集中每个窗口经过GC校正的序列标签数与标签总数的比值；基于该比值，获得所有对照样本对应窗口比值的平均数；基于上述平均数和待测样本的序列标签总数，计算待测样本各个窗口序列标签数的期望值。
可选地，在确定CNV断点后还包括:对根据CNV断点之间的片段进行置信选择的步骤；所述置信选择包括以下步骤：根据调整后的序列标签数目的分布规律，利用对照集确定调整后的序列标签数目正常的置信区间；当片段内调整后的序列标签数目的均值处于置信区间之外时，认为该CNV断点之间的片段确实存在异常。
可选地，序列标签数目符合正态分布，所述置信区间为95％置信区间。
可选地，在选择候选的CNV断点时进行单条染色体环化或者全基因组水平环化处理。
可选地，还包括：所述受试样品是来自人类样本的样品，所述受试样品包括是羊膜腔穿刺获得的羊水、绒毛活检得到的绒毛、经腹脐静脉穿刺得到的脐带血、自然流产的胎儿组织或人体的外周血；和/或所述受试样品的基因组DNA分子采用盐析法、柱层析法、磁珠法、或SDS法的DNA提取方法获得；和/或对所述受试样品的基因组DNA分子采用酶切、雾化、超声、或者HydroShear法随机打断得到DNA片段；和/或对所述受试样品的基因组DNA分子片段进行单向测序或双向测序获得所述DNA片段读序信息。
可选地，该方法还包括：在每个受试样品的DNA片段上加上不同的标签序列以区别不同的受试样品。
根据本发明的另一方面，提供一种拷贝数变异检测系统，包括：读序获取单元，用于获取据受试样品中至少一部分核酸分子的读序信息；序列标签确定单元，用于根据读序信息确定唯一比对至(基因组)参考序列的序列标签；标签数目统计单元，用于将基因组参考序列划分窗口,统计落入各个窗口的序列标签数目；标签数目调整单元，用于对各个窗口的序列标签数目进行GC校正并根据以对照样本集修正的期望序列标签数目进行修正获得调整后的序列标签数目；候选断点选择单元，用于以窗口的起点或终点作为分界点，计算两侧由调整后的序列标签数目组成数值群体的差异显著性值，选取显著性值较小的分界点为候选的CNV断点；断点确定单元，用于对于每个CNV断点至在先CNV断点和在后CNV断点的两段序列，计算所述两段序列包含的窗口的调整后的序列标签数目组成的两个数值群体的差异显著性值，每次剔除最不显著的候选CNV断点并更新被剔除的候选CNV断点左右两个候选CNV断点的差异显著性值，循环迭代，直至所有候选CNV断点的差异显著性值都小于终止阈值，从而确定CNV断点。
可选地，每个窗口具有相同的参考序列标签数目，或者，每个窗口具有相同的长度。
可选地，终止阈值根据由正常样本组成对照样本集获得。
可选地，标签数目调整单元包括：GC校正模块，用于将窗口按GC含量分组，基于组内序列标签的平均数与所有窗口的序列标签平均数获得修正系数，对该窗口序列标签数进行修正获得GC校正的序列标签数；窗口修正模块，用于计算对照集中每个窗口经过GC校正的序列标签数与标签总数的比值；基于该比值，获得所有对照样本对应窗口比值的平均数；基于上述平均数和待测样本的序列标签总数，计算待测样本各个窗口序列标签数的期望值，对GC校正后的序列标签数以对照样本集修正的期望序列标签数目进行修正获得调整后的序列标签数目。
可选地，该系统还包括断点过滤单元，用于在所述断点确定单元确定CNV断点后根据序列标签数目的分布规律，利用对照集确定序列标签数目正常的置信区间；当片段内序列标签数目的均值处于置信区间之外时，认为该CNV断点之间的片段确实存在异常。
可选地，调整后的序列标签数目符合正态分布，所述置信区间为95％置信区间。
可选地，候选断点选择单元在选择候选的CNV断点时进行单条染色体环化或者全基因组水平环化处理。
可选地，受试样品是来自人类样本的样品，所述受试样品包括是羊膜腔穿刺获得的羊水、绒毛活检得到的绒毛、经腹脐静脉穿刺得到的脐带血、自然流产的胎儿组织或人体的外周血；和/或所述受试样品的基因组DNA分子采用盐析法、柱层析法、磁珠法、或SDS法的DNA提取方法获得；和/或对所述受试样品的基因组DNA分子采用酶切、雾化、超声、或者HydroShear法随机打断得到DNA片段；
和/或
对所述受试样品的基因组DNA分子片段进行单向测序或双向测序获得所述DNA片段读序信息。
可选地，通过受试样品DNA片段上所加的不同标签序列以区别不同的受试样品。
本发明的拷贝数变异检测方法和系统的一个优点在于，具有临床可行性，可精确检测到包括较小微缺失/微重复区域的拷贝数变异。
附图说明
图1示出本发明的拷贝数变异检测方法的一个实施例的流程图；
图2示出本发明的拷贝数变异检测方法的另一个实施例的流程图；
图3示出本发明的拷贝数变异检测方法的又一个实施例的流程图；
图4示出本发明一个实现中对染色体CNV分析的简要流程图；
图5示出本发明的拷贝数变异检测系统的一个实施例的结构图；
图6示出本发明的拷贝数变异检测系统的另一个实施例的结构图；
图7A-7H示出本发明的一个应用例中8个样品的检测结果示意图。
具体实施方式
对本文中用到的术语进行说明:
拷贝数变异(copy number variation,CNV):是指待测样品核酸序列与正常样品核酸序列相比，大于1kb的核酸分子的拷贝数发生了变化。拷贝数变异的情形及原因包括：缺失，例如微缺失；插入，例如微插入，微重复，重复，倒置，转座和复杂的多位点变异。
非整倍体：是指与正常样本相比，遗传物质中存在染色体数目的增加或减少，进一步包括整条或部分染色体的增加或减少。本发明中拷贝数变异也包括非整倍体的情形。
测序：获得样品核酸序列信息的过程。测序可通过各种测序方法进行，包括但不限于双脱氧链终止法；优选高通量的测序方法，包括但不限于第二代测序技术或者是单分子测序技术。
第二代测序平台(Metzker ML.Sequencing technologies-the nextgeneration.Nat Rev Genet.2010Jan；11(1):31-46)包括但不限于Illumina-Solexa(GATM,HiSeq2000TM等)、ABI-Solid和Roche-454(焦磷酸测序)测序平台；单分子测序平台(技术)包括但不限于Helicos公司的真实单分子测序技术(True Single Molecule DNAsequencing),Pacific Biosciences公司单分子实时测序(single molecule real-time(SMRTTM))，以及Oxford Nanopore Technologies公司的纳米孔测序技术等(Rusk,Nicole(2009-04-01).Cheap Third-Generation Sequencing.Nature Methods6(4):2446(4)。
测序类型可以为单向(single-end)测序和双向(Pair-end)测序，测序长度可以为50bp、90bp、或100bp。在本发明的一个实施方案中，测序平台为Illumina/Solexa，测序类型为Pair-end测序，得到具有双向位置关系的100bp大小的DNA序列分子。
本发明的一个实施方案中，测序的测序深度可以依据检测的受试样品染色体变异片段大小确定，测序深度越高，检测的灵敏度越高，即可检出的缺失和重复的片段越小。测序深度可以是0.1-30×，即总数据量为人类基因组长度的0.1-30倍，例如在本发明的一个实施方案中，测序深度为0.1×，(2.5×108bp)。
读序(reads)：一定长度的核酸序列(一般大于20bp),例如由测序仪产生的测序序列的结果；通过序列比对的方法，可定位至参考序列的特定区域或位置。
序列比对(比对)：是指一个或多个核酸序列与参考序列进行比较的过程。具体常见为将一段较短的核酸序列(如读序)与参考基因组序列相比较，以确定较短核酸序列在参考基因组上的位置。在利用计算机进行序列比对时，序列比对可以通过任何一种序列比对程序，ELAND(efficient local alignment of nucleotide data),SOAP(ShortOligonucleotide Analysis Package)和BWA(Burrows-Wheeler Aligner)等程序进行。比对成功的认定标准又分为不容错比对(100％匹配)和部分容错比对(小于100％的匹配)。
序列标签：是指可定位参考序列(例如参考基因组序列)唯一位置的读序(reads)。
参考序列标签(reference unique reads)：是指具有固定长度，且在参考序列(通常是参考基因组)上有唯一比对位置的序列。获得参考序列标签的过程，例如包括，将参考基因组切割为固定长度的序列，将这些序列比对回参考基因组，选择唯一比对到参考基因组的序列为参考唯一比对序列。固定长度依测序仪的测序结果序列长度而定，具体可参考平均长度。不同测序仪得到的测序结果长度是不同的，具体每一次测序，测序结果的长度也可能不同，该长度的选取存在一定主观和经验因素。
标签序列(index)：是一段特定长度的，起标识性作用的核酸序列。当待测的DNA分子来自多个受试样本时，每个样本可以被加上不同的标签序列，以用于在测序过程中进行样品的区分(Micah Hamady,Jeffrey J Walker,J Kirk Harris et al.Error-correctingbarcoded primers forpyrosequencing hundreds of samples in multiplex.NatureMethods,2008,March,Vol.5 No.3)，从而实现同时对多个样品进行测序。标签序列为了区分不同序列，但不影响添加标签序列的DNA分子的其他功能。
GC校正：因为测序批次间/内存在一定的GC偏向性，会使基因组中高GC或低GC区域出现拷贝数偏差，对测序数据基于对照样本集进行GC校正得到每个窗口中校正后的相对测序序列数，可以去除此偏向性，提高拷贝数变异检测的精度。
平均数:本文所指平均数一般为算术平均数或中位数。
序列标签数目：序列标签数目可以基于最初的数量统计得到个数统计量，也可以是序列标签个数经其他参数修正得到的修正值，例如可以是一个比值，一些情况下与“拷贝率”可以互换。
受试样品：在一些情形下可被称为待测样品，是指含有核酸分子，而且核酸分子被怀疑发生变异的样品。核酸的类型并不受特别限制，可以是脱氧核糖核酸(DNA)，也可以是核糖核酸(RNA)，优选DNA。对于RNA，可以通过常规手段将其转换为具有相应序列的DNA，进行后续检测和分析。
对照样本：是相对受试样本而言的，是被认为正常的样本，通常情况下，正常是指表型正常。
对照样本集(对照集)：是指对照样本组成的集合，本发明一个实施中，集合中对照样本的数目要求大于30。
下面参照附图对本发明进行更全面的描述，其中说明本发明的示例性实施例。
随着高通量测序技术的不断发展与测序成本的不断降低，测序技术在染色体变异检测方面得到了越来越广泛的应用。
为了提高目前临床上对拷贝数变异检测的技术，本公开设计了一种基于高通量测序技术的进行全基因组水平拷贝数变异筛查的技术方案，具有通量高、特异性高、定位准确的优点。通过获取受试者样品，提取DNA，进行高通量测序，对获得的数据进行分析，得出检测结果。
图1示出本发明的拷贝数变异检测方法的一个实施例的流程图。
如图1所示，在步骤102，获取受试样品中至少一部分核酸分子的读序(reads)信息。可以对受试样本中的至少一部分核酸分子或者全部核酸分子进行测序以获得读序信息。可以获取受试样品的部分核酸分子的读序信息，或者获取全部核酸分子的读序信息。例如，将来自受试样品的基因组DNA分子随机打断得到DNA片段，对DNA片段进行测序，然后获得一定长度的读序。获得的读序的长度可能在一定范围内，可以通过截取操作获得固定长度的读序。DNA片段的长度可以在50bp～1500bp，例如50bp～150bp、150bp～350bp、350bp～500bp、500bp～700bp、700bp～1000bp、或1000bp～1500bp。例如，可以选择50bp、90bp、100bp、150bp、300bp、350bp、500bp、700bp、1000bp、1500bp。在一个实施例中，优选300bp～700bp，更优选350bp～500bp。读序的长度鉴于测序仪的不同可能会有较大的差别，例如illumina-solexa和life technologies-solid仪器目前普遍读序长度在300bp以内，而roche-454、传统的sanger测序、最新的单分子的测序系统得到的读序长度可接近或超过1000bp。为了适应唯一比对的需要，在基于读序选择序列标签时，一般会选择读序中20bp以上的序列去比对，优选地，在26bp以上。
步骤104，根据读序信息确定唯一比对至(基因组)参考序列的序列标签。例如，将读序的全部或部分序列与基因组参考序列进行比对以获得读序在基因组上的位点信息，得到读序在具体染色体上的位点信息。对于人类受试样本，人类基因组参考序列可以是NCBI数据库中的人类基因组参考序列。在本发明的一个实施例中，人类基因组序列是NCBI数据库中版本36(hg18；NCBI Build36)的人类基因组参考序列，所采用的比对软件是SOAPaligner/soap2。选择与基因组参考序列唯一比对上的DNA片段读序，即只在人类基因组参考序列上比对上一次的读序，也就是唯一比对至(基因组)参考序列的序列标签。
步骤106，将基因组参考序列划分窗口，统计落入各个窗口的序列标签数目。对于窗口采用以下方法确定：将参考基因组打断为与待检测样品测序长度相同的片段，采用相同的比对软件的相同参数进行比对，筛选出染色体上可被唯一比对的位置；然后每间隔一定长度的可比对位置划定一次窗口。窗口之间可以选择是否有交叉滑动。窗口内应包含的可被唯一比对的位点数与待测样品的测序数据量有关。一般要使待测样品落在每个窗口中的期望reads数在300条以上，从而保证落在窗口内的reads数分布符合泊松分布。例如，假设基因组可被唯一比对的点数为N，待检测样品的有效reads数为n，每个窗口的期望reads数为E，则参考基因组的每个窗口中应包含个可被唯一比对位点。
步骤108，对各个窗口的序列标签数目进行GC校正并根据以对照样本集修正的期望序列标签数目进行修正获得调整后的序列标签数目。在一个实施例中，对各个窗口的序列标签数目进行GC校正包括以下的步骤:将窗口按GC含量分组，基于组内序列标签的平均数与所有窗口的序列标签平均数获得修正系数，对该窗口序列标签数进行修正获得GC校正的序列标签数；对照样本集修正的期望序列标签数目通过以下方法获得：计算对照集中每个窗口经过GC校正的序列标签数与标签总数的比值；基于该比值，获得所有对照样本对应窗口比值的平均数；基于上述平均数和待测样本的序列标签总数，计算待测样本各个窗口序列标签数的期望值
步骤110，以窗口的起点或终点作为分界点，计算两侧由调整后的序列标签数目组成数值群体的差异显著性值，选取显著性值较小(即差异显著性较大)的分界点为候选的CNV断点。例如，在全基因组范围内根据表征每个窗口两侧拷贝数变异的显著性水平的p值选取预定数量的窗口作为候选CNV断点，获得每个候选CNV断点的差异显著性值，即p值。
步骤112，对于每个CNV断点至在先CNV断点和在后CNV断点的两段序列，计算所述两段序列包含的窗口的调整后的序列标签数目组成的两个数值群体的差异显著性值，每次剔除最不显著的候选CNV断点并更新被剔除的候选CNV断点左右两个候选CNV断点的差异显著性值，循环迭代，直至所有候选CNV断点的差异显著性值都小于终止阈值，从而确定CNV断点。终止阈值通常预先设定。例如，通过对由正常样本组成对照样本集进行分析处理获得该终止阈值。
上述实施例中，通过将获取的读序与基因组参考序列比对，选择唯一比对的读序进行窗口统计，对每个窗口的序列标签数目进行GC校正和对照集修正后进行基因显著性差异的循环迭代进行CNV断点选择，从而完成CNV的检测，能够精确检测到较小的包括微缺失/微重复区域的拷贝数变异。
对于人类样本，受试样品可以是羊膜腔穿刺获得的羊水、绒毛活检得到的绒毛、经腹脐静脉穿刺得到的脐带血、自然流产的胎儿组织或人体外周血中获取的基因组DNA，可以采用盐析法、柱层析法、磁珠法、SDS法等常规DNA提取方法。在一个实施例中，优选采用柱层析法，该柱层析法是指血液、组织或细胞经过细胞裂解液和蛋白酶K的作用后得到裸露的DNA分子，DNA分子在高盐条件下与硅胶膜结合，在低盐、高pH值时DNA分子从硅胶膜上洗脱下来。具体原理和方法参见Tiangen TIANamp Micro DNA Kit(DP316)产品说明书))。
当待测的DNA片段来自多个受试样本的DNA分子时，每个受试样本的DNA片段可以被加上不同的标签序列(index)，长度可选4bp～12bp，以用于在测序过程中进行受试样品的区分(Micah Hamady,Jeffrey J Walker,J Kirk Harris et al.Error-correctingbarcoded primers forpyrosequencing hundreds of samples in multiplex.NatureMethods,2008,March,Vol.5No.3)。通过这样的方式，能够同时处理多个受试样品的检测，提高了效率，减小了检测成本。
图2示出本发明的拷贝数变异检测方法的另一个实施例的流程图。
步骤202，对受试样品的基因组DNA分子随机打断得到DNA片段。DNA分子随机打断处理可以采用酶切、雾化、超声、或者HydroShear法。优选地，采用超声法，例如Covaris公司的S-series(基于AFA技术，当由传感器释放的声能/机械能通过DNA样品时，溶解气体形成气泡。当能量移除后，气泡破裂并产生断裂DNA分子的能力。通过设置一定的能量强度和时间间隔等条件，可将DNA分子打断至一定范围的大小。具体原理和方法请参见Covaris公司的S-series说明书)，将DNA分子打断为比较集中的一定大小的片段。
步骤204，对DNA片段进行测序获得DNA片段测序序列，即读序。测序获得的读序可以在一定长度范围内，可以通过截取操作从DNA片段测序序列获得固定长度的读序，在本实施例下文中DNA片段测序序列指固定长度的读序。所采用的测序方法可以为高通量测序方法Illumina/Hiseq2000、ABI/SOLiD、Roche/454。测序类型可以为Single-end(单向)测序或Pair-end(双向)测序，测序长度可以为50bp～1500bp。在本发明的一个实施例中，测序平台为Illumina/Hiseq2000，测序类型为Pair-end测序，得到具有双向位置关系的100bp大小的DNA序列分子。对测序深度有一定要求，测序深度可以依据检测的染色体变异片段大小确定，测序深度越高，检测的灵敏度越高，即可检出的缺失和重复的片段越小。在本发明的一个实施例中，人类受试样本的测序量在2～900×108条测序片段之间。
步骤206，将读序与基因组参考序列进行比对以获得读序在基因组上的位点信息。
步骤206，选择与基因组参考序列唯一比对上的读序作为序列标签。
步骤208，统计落入基因组上每个窗口的DNA片段唯一比对的序列标签的个数。对于每个受试样品，统计每个窗口所落的序列标签个数(记为ni,j，下标i和j分别代表窗口编号和样本编号，以作区分，下不复述)。
步骤210，确定基因组上每个窗口的平均GC含量以确定该窗口修正系数，根据修正系数获得每个窗口的序列标签的修正个数。该步骤主要根据每个窗口的GC含量对窗口的序列标签个数进行修正，可称为批次修正或GC校正。
统计每个窗口所落的序列标签的平均GC含量(记为GCi,j)。平均GC含量(GCi,j)的计算，即计算落在统计窗口内的所有序列标签的平均GC含量水平。在计算时，将所有序列标签中的G和C碱基的个数总和即为Ngc，所有序列标签的总长度记为L，则
为了修正样品数据量的差异，将各个统计窗口按其GCi,j进行分组，即GCi,j相同的窗口划为一组。取各分组中序列标签数ng,j的中位数或算术平均数mg,j，除全基因组水平序列标签数的中位数或算术平均数mj，得到修正系数cg,j(cg,j＝mg,j/mj)，下标g代表不同分组的GC含量。将原始的每个窗口的序列标签数ni,j乘以修正系数cg,j得到每个窗口内序列标签数的修正值(记为ni,j′)
步骤212，将每个窗口的序列标签的修正个数除以该窗口的个数期望值获得每个窗口的调整后的序列标签数，即拷贝率。其中每个窗口的个数期望值通过由正常样本组成对照集(control set)获得。该步骤主要根据正常样本数据进行每个窗口的序列标签的个数的校正，可以称为窗口校正。
在对照集样本中，定义相对序列标签数百分比(Pi,j)为窗口内序列标签数(ni,j')和全基因组总序列标签数(Nj)之比，即然后计算得到对照集中每个窗口的平均百分比 在受试样本中，每个窗口的拷贝率ri,j等于修正序列标签数ni,j′除以窗口内的期望值(全基因组总的序列标签数乘以窗口序列标签数百分比)，即
在对照集的选择中，建库方法、测序试剂及测序类型等应尽量与待测样品一致，从而提高对照样品对待测样品的校正效果。对照集中的样品应为正常样本，且样本量大于30。
步骤214，在全基因组范围内根据表征每个窗口两侧拷贝数变异的显著性差异的p值选取预定数量的窗口作为候选CNV断点，获得每个候选CNV断点的p值。
1)候选CNV断点选择：针对全基因组的所有窗口，计算每个窗口左右两侧一定数量窗口(窗口数量通常大于30或满足检验模型的最低样品量限制，以使得检验模型具有统计意义)的拷贝率变化差异，根据全基因组范围的显著性差异水平(p值从小到大)，选取一定数量(例如，总窗口数量的1％)的(与窗口对应)点作为候选CNV断点(Breakpoint，即每个CNV片段的边界点)。
2)初始化：将排序后的所有断点的集合记为Bc＝{b1,b2,....bk..,bs}，则在每个断点k左右两侧，分别存在相邻的断点k+1、k-1。通过计算k-1与k之间的拷贝数集合与k与k+1之间的拷贝率集合之间的显著性差异，得到代表每个断点两侧显著性差异的p值。例如，在一个实施例中，选用非参数检验中的游程检验，利用两个群体元素混合后的分布均匀状态评价此两个群体的差异显著性，详见：Wald,A.&Wolfowitz,J.On a test whether twosamples are from the same population.The Annals of Mathematical Statistics11,147-162(1940)。
步骤216，每次剔除最不显著的候选CNV断点并更新被剔除的候选CNV断点左右两个候选CNV断点的p值，循环迭代，直至所有候选CNV断点的p值都小于终止p值(即，终止阈值)，终止p值根据对照集获得。
迭代合并：通过不断的循环迭代，每次将最不显著的候选突破点剔除，同时更新左右两个断点的p值，直到所有p值都小于终止p值为止。
终止p值的获得：例如，可以在将对照样品进行上述迭代合并的操作，记录每次迭代合并的最大p值，迭代合并至一条片段时终止。此时，根据最大p值的变化趋势，将最大p值变化最剧烈的点(即取p值的变化曲线中斜率变化最明显的点(曲率最大的点))对应或之前的那次合并过程的最大p值作为终止阈值。在一个实施例中，亦可将迭代合并的终止设置为迭代合并后的片段数等于预知的片段数，例如，当对全基因组进行分析时，对照样品迭代合并后的片段数为24个时终止，通过计算此时终止p值的平均值，可以有效地获得终止p值。
上述步骤214和216也可以称为片段化。注意：在步骤214中1)和2)中的窗口和断点选择时，可以考虑进行单条染色体环化或者全基因组水平环化。单条染色体环化是指：在计算染色体起点附近的窗口时，如果左侧的有效窗口数量不足以进行统计检验，则从本条染色体的尾部反向获得足够的窗口进行计算；同理，对于终点附近右侧无法获得足够有效窗口的位置则从染色体前端获得。这一操作使得位于染色体前部和尾部的窗口依然可以计算。全基因组水平环化则是在每条染色体前端有效窗口数不足时索引前一条染色体的尾端，而尾端不足时索引下一条染色体的前端，1号染色体则和Y染色体相连。
在步骤216后还可以包括对根据CNV断点之间的片段进行置信选择的步骤：根据拷贝率的分布规律，利用对照集确定拷贝率正常的置信区间；当片段内拷贝率的均值处于置信区间之外时，认为该CNV断点之间的片段确实存在异常。在一个实施例中，拷贝率符合正态分布，置信区间为95％置信区间。通过该步骤对片段化结果进行过滤，获得可靠的结果。如果片段的平均ri,j′小于阈值下限或者大于阈值上限时，将作为阳性结果输出。
阈值选择：统计每个对照样品中窗口的拷贝率分布，根据中心极限定理，读序在窗口内是随机的，所以拷贝率r是服从正态分布的，取其显著性水平为0.05时的左右分位点。计算其在对照集中的均值，分别作为筛选拷贝率变异时的上下限阈值。
上述实施例中，通过批次修正和窗口校正，提高了检测结果的准确率。通过引入参照集，可以通过扩大参照集来提高精度，以减轻对起始DNA量的压力。
图3示出本发明的基因组拷贝数变异检测方法的又一个实施例的流程图。在图3中，包括了由正常样本组成的对照集的处理流程(3A)和受试样品的处理流程(3B)。对照集主要用于获取对受试样品进行校正的数据以及受试样品处理的迭代终止条件的终止阈值。
如图3所示，流程3A包括：
步骤310A，提取对照样品的DNA分子。
步骤311A，对对照样品的DNA分子随机打断为DNA片段进行测序，获得对照样品的DNA片段测试序列数据，即读序。
步骤312A，将对照样品的读序与参考基因组进行对比。
步骤313A，统计窗口唯一比对读序数，即序列标签数目。
步骤314A，对对照样品进行批次修正。
步骤315A，通过对照样品获得窗口个数期望值，用于受试样品的窗口校正。
步骤316A，断点选择与片段化。选择候选CNV断点，每次剔除最不显著的候选CNV断点并更新被剔除的候选CNV断点左右两个候选CNV断点的p值，循环迭代，使最终的片段剩余预定个数(例如24个)时停止。
步骤317A，确定终止阈值。通过计算此时终止p值的平均值，可以有效地获得终止p值，作为受试样品迭代终止条件的终止阈值。
流程3B包括：
步骤310B，提取受试样品的DNA分子。
步骤311B，对受试样品的DNA分子随机打断为DNA片段进行测序，获得受试样品的DNA片段的读序。
步骤312B，将受试样品的DNA片段的读序与参考基因组进行对比。
步骤313B，去统计窗口唯一比对读序数，即序列标签数。
步骤314B，对受试样品进行批次修正。
步骤315B，根据对照样品获得的窗口个数期望值对受试样品进行窗口校正。
步骤316B，断点选择与片段化。
步骤317B，进行结果过滤。
在对照集的选择中，建库方法、测序试剂及测序类型等应尽量与受试样品一致，从而提高对照样品对受试样品的校正效果。对照集中的样品应为正常样本，且样本量大于30。
图4示出本发明一个实现中染色体CNV分析的简要流程图。
如图4所示，步骤401，DNA提取及测序：根据Tiangen DP327-02Kit操作手册提取基因组DNA后，基本按照Illumina/Hiseq2000标准建库流程进行建库。在这个过程中，本身集中于500bp的DNA分子两端被加上测序所用接头，每个样本被加上不同的标签序列(index)，从而在一次测序得到的数据中可以使多个样本的数据区分开。
步骤402，序列比对：利用测序方法Illumina/Hiseq2000测序(用其它测序方法如ABI/SOLiD能达到相同或相近的效果)，每个样本得到一定长度片段的DNA读序，将其与NCBI数据库中的标准人类基因组参考序列进行SOAP2比对，得到读序定位于基因组相应位置的信息。为避免重复序列对CNV分析的干扰，只选取与人类基因组参考序列唯一比对的读序，即只在人类基因组参考序列上比对上一次的读序，也即序列标签数，作为后续CNV分析的有效数据。
步骤403，PSCC分析。运用本申请人自主开发的针对全基因组拷贝数变异检测的一系列生物信息学方法(下称PSCC)，在对受试样本进行批次性差异修正，根据对照集(control set)，对受试样品进行窗口校正(correction)、标准化(Normalization)和拷贝数片段化(segmentation)。
步骤404，基于步骤403中得到的拷贝数片段进行CNV分析，将受试样品拷贝率≤0.7和≥1.3分别作为片段缺失和重复的检测阈值，分析得到全基因组水平拷贝数变异片段，然后进行结果的可视化。
上述实施例中，实现的软件算法是一种由深圳华大基因研究院开发针对全基因组拷贝数变异检测的一系列程序，统称为PSCC。它能够通过新一代测序技术产生的数据，将受试样本进行批次修正，然后和对照集合进行数据校正、标准化和片段化，估算出受试样本的拷贝数变异程度和大小。在较低测序深度(50M条测序短序列)的时候可检测出0.5Mb左右的单个拷贝数变异(CNV)片段。
图5示出本发明的基因组拷贝数变异检测系统的一个实施例的结构图。如图5所示，该系统包括：读序获取单元51，获取受试样品中至少一部分核酸分子的读序信息；序列标签确定单元52，确定唯一比对至(基因组)参考序列的读序作为序列标签；标签数目统计单元53，将基因组参考序列划分窗口，统计落入各个窗口的序列标签数目；标签数目调整单元54，对各个窗口的序列标签数目进行GC校正并根据以对照样本集修正的期望序列标签数目进行修正获得调整后的序列标签数目；候选断点选择单元55，以窗口的起点或终点作为分界点，计算两侧由调整后的序列标签数目组成数值群体的差异显著性值，选取显著性值较小(即差异显著性较大)的分界点为候选的CNV断点；断点确定单元56，对于每个CNV断点至在先CNV断点和在后CNV断点的两段序列，计算所述两段序列包含的窗口的调整后的序列标签数目组成的两个数值群体的差异显著性值，每次剔除最不显著的候选CNV断点并更新被剔除的候选CNV断点左右两个候选CNV断点的差异显著性值，循环迭代，直至所有候选CNV断点的差异显著性值都小于终止阈值，从而确定CNV断点；该终止阈值可以根据由正常样本组成对照样本集获得。标签数目统计单元53划分窗口时，可以使得每个窗口具有相同的参考序列标签数目(reference unique reads)，或者，使得每个窗口具有相同的长度。在一个实施例中，候选断点选择单元55在选择候选的CNV断点时进行单条染色体环化或者全基因组水平环化处理
上述实施例中，序列标签确定单元根据读序确定唯一比对至(基因组)参考序列的序列标签，标签数目调整单元对各个窗口的序列标签数目进行修正，由候选断点选择单元和断点确定单元进行基因显著性差异的循环迭代进行CNV断点选择，从而完成CNV的检测，能够精确检测到包括较小的微缺失/微重复区域的拷贝数变异。
图6示出本发明的拷贝数变异检测系统的另一个实施例的结构图。如图6所示，该系统包括：读序获取单元51、序列标签确定单元52、标签数目统计单元53、标签数目调整单元64、候选断点选择单元55和断点确定单元56。读序获取单元51、序列标签确定单元52、标签数目统计单元53、候选断点选择单元55、和断点确定单元56可以参见图5的具体描述，为简洁起见在此不再详细描述。标签数目调整单元64包括GC校正模块641和窗口修正模块642。其中，GC校正模块641，将窗口按GC含量分组，基于组内序列标签的平均数与所有窗口的序列标签平均数获得修正系数，对该窗口序列标签数进行修正获得GC校正的序列标签数；窗口修正模块642，计算对照集中每个窗口经过GC校正的序列标签数与标签总数的比值；基于该比值，获得所有对照样本对应窗口比值的平均数；基于上述平均数和待测样本的序列标签总数，计算待测样本各个窗口序列标签数的期望值，对GC校正后的序列标签数以对照样本集修正的期望序列标签数目进行修正获得调整后的序列标签数目，也称拷贝率。
根据本发明的一个实施例，该系统还包括断点过滤单元67，在断点确定单元确定CNV断点后根据拷贝率的分布规律，利用对照集确定拷贝率正常的置信区间；当片段内拷贝率的均值处于置信区间之外时，认为该CNV断点之间的片段确实存在异常。在一个实施例中，序列标签数目符合正态分布，所述置信区间为95％置信区间。
在一个实施例中，受试样品是来自人类样本的样品，所述受试样品的基因组DNA分子是羊膜腔穿刺获得的羊水、绒毛活检得到的绒毛、经腹脐静脉穿刺得到的脐带血、自然流产的胎儿组织或人体外周血中获取的基因组DNA；受试样品的基因组DNA分子采用盐析法、柱层析法、磁珠法、或SDS法的DNA提取方法获得；对受试样品的基因组DNA分子为采用酶切、雾化、超声、或者HydroShear法随机打断得到DNA片段；DNA片段测序序列通过对受试样品的基因组DNA分子片段进行单向测序或双向测序获得。可以通过受试样品DNA片段上所加的不同标签序列以区别不同的受试样品。
对于图5至图6中各个单元的功能，可以参考上文中关于本发明方法的实施例中对应部分的说明，为简洁起见，在此不再详述。
本领域技术人员将意识到硬件、固件和软件配置在这些情况下的可替换性，以及如何最好地实现每个特定应用地该功能。
下面将结合应用例对本发明的实施例进行详细描述，但是本领域技术人员将会理解，下列应用例仅用于说明本发明，而不应视为限定本发明的范围。应用例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市场获得的常规产品。以下括号内为各个试剂或试剂盒的厂家货号。所使用的测序用的接头和标签序列index来源于Illumina公司的Multiplexing Sample PreparationOligonutide Kit。
应用例一、2个染色体数目异常及6个微缺失样品检测
1、DNA提取：
按照Tiangen的TIANamp Micro DNA Kit(DP316)操作流程提取8例样品(以下简称Sample1、Sample2、Sample3…Sample8)的DNA，所提取DNA按照修改后的Illumina/Hiseq2000标准建库流程进行建库，在本身集中于500bp的DNA分子两端被加上测序所用接头，每个样本被加上不同的标签序列(index)，然后与Flowcell表面互补接头杂交，在一定条件下使核酸分子成簇生长，然后在Illumina Hiseq2000上通过双末端测序，得到长度为100bp的DNA片段序列。
具体的，将获自上述羊水样品的约100ng(Quant-IT dsDNA HS Assay kit定量)的DNA，进行修改后的Illumina/Hiseq2000标准流程建库，具体流程参照现有技术(参见http://www.illumina.com/提供的Illumina/Solexa标准建库说明书)。经2100Bioanalyzer(Agilent)确定DNA文库大小及插入片段为约500bp，QPCR精确定量后可上机测序。
2.测序：本实施例中，对于获自上述8例样品的DNA样本按照Illumina/Solexa官方公布的ClusterStation和Hiseq2000(PEsequencing)说明书进行操作，使每个样品得到约5G数据量进行上机测序，每个样本根据所述index标签区进行分开。利用比对软件SOAP2，将测序所得DNA序列与NCBI数据库中版本36(hg18；NCBIBuild36)的人类基因组参考序列进行比对，得到所测DNA序列定位于基因组相应位置的信息。
3.数据分析
a)基本统计：统计每个窗口所落的序列标签个数(即唯一比对的reads数，记为ni,j，下标i和j分别代表窗口编号和样本编号，以作区分，下不复述)和平均GC含量(记为GCi,j)。
b)批次修正：为了修正样品数据量的差异，将各个统计窗口按落入其中的reads平均GC分组，取各分组的中位数或算术平均数除全基因组水平的中位数或算术平均数之间得到修正系数cg，下标g代表不同分组的GC含量。将原始的每个窗口的reads个数ni,j乘以修正系数cg得到每个窗口内序列标签数的修正值(记为ni,j′)。
c)数据校正：选取90例YH群体的测序数据作为对照集，定义相对序列标签数百分比(Pi,j)为窗口内测序序列标签数(ni,j')和全基因组总序列标签数(Nj)之比，即然后计算得到每个窗口的平均百分比 在受试样本中，每个窗口的拷贝率ri,j等于修正修正序列标签数ni,j′除以窗口内的期望值(全基因组总序列标签数乘以窗口序列标签数百分比)，即
d)片段化
①候选断点选择：针对全基因组的所有ri,j，计算其左右两侧100个窗口的拷贝率变化差异，根据全基因组范围的显著水平(P值从小到大)，选取最显著的10000个点作为候选CNV断点。
②初始化：将排序后的所有断点的集合记为Bc＝{b1,b2,....bk..,bs}，则在每个断点k左右两侧，分别存在相邻的断点k+1、k-1。通过计算k-1与k之间的拷贝率集合与k与k+1之间的拷贝率集合之间的显著性差异(本实施例中选用非参数检验中的游程检验)，得到代表每个断点两侧显著性差异的p值。
③迭代合并：通过不断的循环迭代，每次将最不显著的候选断点剔除，同时更新左右两个断点的p值，直到所有p值都小于1E-50时为止。
注意：在步骤①和②中的窗口和断点选择时，进行了全基因组水平环化，从而使得位于染色体前部和尾部的窗口依然可以计算。
e)阈值与过滤：对片段化过后的结果进行过滤。如果片段的平均ri,j′小于0.7或者大于1.3时，作为阳性结果输出。
f)结果可视化。
4.结果统计
对于8个样品，详细的检测结果和验证结果如下表1所示。
其中验证结果通过CGH芯片(比较基因组杂交芯片)得出。实验使用Human GenomeCGH Microarray Kit，(Agilent Technologies Inc.)，完全按照厂商的使用说明进行操作。
表1.8例样品的CNV结果其中chr代表染色体，T7代表7号染色体三体型。XYY未性染色体三体异常。
图7A-7H示出8个样品的检测结果示意图。
从上表1和图7A-7H可以看到，小至0.4M的微缺失片段，大到整条染色体数目异常，本方法都可精确的进行检测和定位，证明其检测效率和精度都较为优秀。
与目前已报道的拷贝数变异检测的分析方法对比，本公开的优越性主要有一下几点：
(1)分辨率，使用50M左右的数据，可精确检测到0.5M的微缺失区域。
(2)可扩展性。除了通过增加测序量之外，可以通过扩大参照集来提高精度，以减轻对起始DNA量的压力。
(3)更稳定，更加全面。已报道文章中，并无明确指出自身的操作细节，而本发明涉及数据批次修正，群体校正，片段化条件优选等的各个方面。
本发明对适用人群进行全基因组拷贝数变异检测分析，有利于提供遗传咨询和提供临床决策依据，对微缺失综合征患者进行准确的病理判定。本发明适用人群可以是所有微缺失患者或潜在微缺失携带者，适用人群举例仅用于说明本发明，而不应为限定本发明的范围。
本发明的描述是为了示例和描述起见而给出的，而并不是无遗漏的或者将本发明限于所公开的形式。很多修改和变化对于本领域的普通技术人员而言是显然的。选择和描述实施例是为了更好说明本发明的原理和实际应用，并且使本领域的普通技术人员能够理解本发明从而设计适于特定用途的带有各种修改的各种实施例。

Claims (17)

1.一种拷贝数变异检测方法，其特征在于，包括：
获取受试样品中至少一部分核酸分子的读序信息；
根据所述读序确定唯一比对至基因组参考序列的序列标签；
将基因组参考序列划分窗口，统计落入各个窗口的序列标签数目，其中，每个窗口具有相同的参考序列标签数目，其中，每间隔一定长度的可比对位置划定一次窗口，窗口之间可以选择是否有交叉滑动；
对各个窗口的序列标签数目进行GC校正并根据以对照样本集修正的期望序列标签数目进行修正获得调整后的序列标签数目；
以窗口的起点或终点作为分界点，计算两侧由调整后的序列标签数目组成数值群体的差异显著性值，选取显著性值较小的分界点为候选的CNV断点；
对于每个CNV断点至在先CNV断点和在后CNV断点的两段序列，计算所述两段序列包含的窗口的调整后的序列标签数目组成的两个数值群体的差异显著性值，每次剔除最不显著的候选CNV断点并更新被剔除的候选CNV断点左右两个候选CNV断点的差异显著性值，循环迭代，直至所有候选CNV断点的差异显著性值都小于终止阈值，从而确定CNV断点。
2.根据权利要求1所述的方法，其特征在于，还包括对受试样本中的至少一部分核酸分子进行测序以获得读序信息的步骤。
3.根据权利要求1或2所述的方法，其特征在于，所述终止阈值根据由正常样本组成对照样本集获得。
4.根据权利要求1所述的方法，其特征在于，所述对各个窗口的序列标签数目进行GC校正包括以下的步骤：
将窗口按GC含量分组，基于组内序列标签的平均数与所有窗口的序列标签平均数获得修正系数，对该窗口序列标签数进行修正获得GC校正的序列标签数；
和/或
对照样本集修正的期望序列标签数目由以下方法获得：
计算对照集中每个窗口经过GC校正的序列标签数与序列标签总数的比值；基于该比值，获得所有对照样本对应窗口比值的平均数；基于上述平均数和待测样本的序列标签总数，计算待测样本各个窗口序列标签数的期望值。
5.根据权利要求4所述的方法，其特征在于，在选择候选的CNV断点时进行单条染色体环化或者全基因组水平环化处理。
6.根据权利要求1或2所述的方法，其特征在于，在确定CNV断点后还包括:
对根据CNV断点之间的片段进行置信选择的步骤；
所述置信选择包括以下步骤：根据调整后的序列标签数目的分布规律，利用对照集确定调整后的序列标签数目正常的置信区间；当片段内调整后的序列标签数目的均值处于置信区间之外时，认为该CNV断点之间的片段确实存在异常。
7.根据权利要求6所述的方法，其特征在于，所述调整后的序列标签数目符合正态分布，所述置信区间为95％置信区间。
8.根据权利要求1所述的方法，其特征在于，还包括：
所述受试样品是来自人类样本的样品，所述受试样品包括是羊膜腔穿刺获得的羊水、绒毛活检得到的绒毛、经腹脐静脉穿刺得到的脐带血、自然流产的胎儿组织或人体外周血；
和/或
所述受试样品的基因组DNA分子采用盐析法、柱层析法、磁珠法、或SDS法的DNA提取方法获得；
和/或
对所述受试样品的基因组DNA分子采用酶切、雾化、超声、或者HydroShear法随机打断得到DNA片段；
和/或
对所述受试样品的基因组DNA分子片段进行单向测序或双向测序获得所述DNA片段读序。
9.根据权利要求1所述的方法，其特征在于，还包括：
在每个受试样品的DNA片段上加上不同的标签序列以区别不同的受试样品。
10.一种拷贝数变异检测系统，其特征在于，包括：
读序获取单元，用于获取受试样品中至少一部分核酸分子的读序信息；
序列标签确定单元，用于根据所述读序信息确定唯一比对至基因组参考序列的序列标签；
标签数目统计单元，用于将基因组参考序列划分窗口，统计落入各个窗口的序列标签数目，其中，每个窗口具有相同的参考序列标签数目，其中，每间隔一定长度的可比对位置划定一次窗口，窗口之间可以选择是否有交叉滑动；
标签数目调整单元，用于对各个窗口的序列标签数目进行GC校正并根据以对照样本集修正的期望序列标签数目进行修正获得调整后的序列标签数目；
候选断点选择单元，用于以窗口的起点或终点作为分界点，计算两侧由调整后的序列标签数目组成数值群体的差异显著性值，选取显著性值较小的分界点为候选的CNV断点；
断点确定单元，用于对于每个CNV断点至在先CNV断点和在后CNV断点的两段序列，计算所述两段序列包含的窗口的调整后的序列标签数目组成的两个数值群体的差异显著性值，每次剔除最不显著的候选CNV断点并更新被剔除的候选CNV断点左右两个候选CNV断点的差异显著性值，循环迭代，直至所有候选CNV断点的差异显著性值都小于终止阈值，从而确定CNV断点。
11.根据权利要求10所述的系统，其特征在于，所述终止阈值根据由正常样本组成对照样本集获得。
12.根据权利要求10所述的系统，其特征在于，所述标签数目调整单元包括：
GC校正模块，用于将窗口按GC含量分组，基于组内序列标签的平均数与所有窗口的序列标签平均数获得修正系数，对该窗口序列标签数进行修正获得GC校正的序列标签数；
窗口修正模块，用于计算对照集中每个窗口经过GC校正的序列标签数与序列标签总数的比值；基于该比值，获得所有对照样本对应窗口比值的平均数；基于上述平均数和待测样本的序列标签总数，计算待测样本各个窗口序列标签数的期望值，对GC校正后的序列标签数以对照样本集修正的期望序列标签数目进行修正获得调整后的序列标签数目。
13.根据权利要求10所述的系统，其特征在于，还包括断点过滤单元，用于在所述断点确定单元确定CNV断点后根据序列标签数目的分布规律，利用对照集确定序列标签数目正常的置信区间；当片段内序列标签数目的均值处于置信区间之外时，认为该CNV断点之间的片段确实存在异常。
14.根据权利要求13所述的系统，其特征在于，所述序列标签数目符合正态分布，所述置信区间为95％置信区间。
15.根据权利要求10至14中任意一项所述的系统，其特征在于，所述受试样品是来自人类样本的样品，所述受试样品包括是羊膜腔穿刺获得的羊水、绒毛活检得到的绒毛、经腹脐静脉穿刺得到的脐带血、自然流产的胎儿组织或人体的外周血；
和/或
所述受试样品的基因组DNA分子采用盐析法、柱层析法、磁珠法、或SDS法的DNA提取方法获得；
和/或
对所述受试样品的基因组DNA分子采用酶切、雾化、超声、或者HydroShear法随机打断得到DNA片段；
和/或
对所述受试样品的基因组DNA分子片段进行单向测序或双向测序获得所述DNA片段的读序。
16.根据权利要求10所述的系统，其特征在于，通过受试样品DNA片段上所加的不同标签序列以区别不同的受试样品。
17.根据权利要求10所述的系统，其特征在于，候选断点选择单元在选择候选的CNV断点时进行单条染色体环化或者全基因组水平环化处理。
CN201280066929.3A 2012-04-05 2012-04-05 一种拷贝数变异检测方法和系统 Active CN104221022B (zh)

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2012/073545 WO2013149385A1 (zh) 2012-04-05 2012-04-05 一种拷贝数变异检测方法和系统

Publications (2)

Publication Number Publication Date
CN104221022A CN104221022A (zh) 2014-12-17
CN104221022B true CN104221022B (zh) 2017-11-21

Family

ID=49299922

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201280066929.3A Active CN104221022B (zh) 2012-04-05 2012-04-05 一种拷贝数变异检测方法和系统

Country Status (10)

Country Link
US (2) US20150056619A1 (zh)
EP (1) EP2835752B8 (zh)
JP (1) JP5972448B2 (zh)
KR (1) KR101795124B1 (zh)
CN (1) CN104221022B (zh)
AU (1) AU2012376134B2 (zh)
IL (1) IL234875B (zh)
RU (1) RU2014144349A (zh)
SG (1) SG11201406250SA (zh)
WO (1) WO2013149385A1 (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109979535A (zh) * 2017-12-28 2019-07-05 安诺优达基因科技(北京)有限公司 一种胚胎植入前遗传学筛查装置

Families Citing this family (36)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105224543A (zh) * 2014-05-30 2016-01-06 国际商业机器公司 用于处理时间序列的方法和装置
EP4092680A1 (en) * 2014-09-12 2022-11-23 Illumina Cambridge Limited Detecting repeat expansions with short read sequencing data
WO2016045106A1 (zh) * 2014-09-26 2016-03-31 深圳华大基因股份有限公司 单细胞染色体的cnv分析方法和检测装置
US11242559B2 (en) * 2015-01-13 2022-02-08 The Chinese University Of Hong Kong Method of nuclear DNA and mitochondrial DNA analysis
CN104560697A (zh) * 2015-01-26 2015-04-29 上海美吉生物医药科技有限公司 一种基因组拷贝数不稳定性的检测装置
CN104694384B (zh) * 2015-03-20 2017-02-08 上海美吉生物医药科技有限公司 线粒体dna拷贝数变异性的检测装置
CN104745718B (zh) * 2015-04-23 2018-02-16 北京中仪康卫医疗器械有限公司 一种检测人类胚胎染色体微缺失和微重复的方法
US10395759B2 (en) 2015-05-18 2019-08-27 Regeneron Pharmaceuticals, Inc. Methods and systems for copy number variant detection
CN105243299B (zh) * 2015-09-30 2018-03-06 深圳华大基因科技服务有限公司 一种检测cnv的精确断点及断点周围特征的方法及装置
KR101848438B1 (ko) 2015-10-29 2018-04-13 바이오코아 주식회사 디지털 pcr을 이용한 산전진단 방법
CA3005791A1 (en) * 2015-11-18 2017-05-26 Sophia Genetics S.A. Methods for detecting copy-number variations in next-generation sequencing
CN105760712B (zh) * 2016-03-01 2019-03-26 西安电子科技大学 一种基于新一代测序的拷贝数变异检测方法
RS62390B1 (sr) 2016-07-20 2021-10-29 BioNTech SE Izbor neoepitopa kao ciljeva specifičnih za bolest za terapiju sa povećanom efikasnošću
CN106520940A (zh) * 2016-11-04 2017-03-22 深圳华大基因研究院 一种染色体非整倍体和拷贝数变异检测方法及其应用
TWI607332B (zh) * 2016-12-21 2017-12-01 國立臺灣師範大學 Correlation between persistent organic pollutants and microRNAs station
WO2018144449A1 (en) * 2017-01-31 2018-08-09 Counsyl, Inc. Systems and methods for identifying and quantifying gene copy number variations
WO2018161245A1 (zh) * 2017-03-07 2018-09-13 深圳华大基因研究院 一种染色体变异的检测方法及装置
CN109097457A (zh) * 2017-06-20 2018-12-28 深圳华大智造科技有限公司 确定核酸样本中预定位点突变类型的方法
US20200327957A1 (en) * 2017-12-14 2020-10-15 Ancestry.Com Dna, Llc Detection of deletions and copy number variations in dna sequences
CN109979529B (zh) * 2017-12-28 2021-01-08 北京安诺优达医学检验实验室有限公司 Cnv检测装置
CN108256289B (zh) * 2018-01-17 2020-10-16 湖南大地同年生物科技有限公司 一种基于目标区域捕获测序基因组拷贝数变异的方法
KR102036609B1 (ko) * 2018-02-12 2019-10-28 바이오코아 주식회사 디지털 pcr을 이용한 산전진단 방법
CN108427864B (zh) * 2018-02-14 2019-01-29 南京世和基因生物技术有限公司 一种拷贝数变异的检测方法、装置以及计算机可读介质
CN108415886B (zh) * 2018-03-07 2019-04-05 清华大学 一种基于生产工序的数据标签纠错方法及装置
CN108664766B (zh) * 2018-05-18 2020-01-31 广州金域医学检验中心有限公司 拷贝数变异的分析方法、分析装置、设备及存储介质
CN114502744B (zh) * 2019-12-11 2023-06-23 深圳华大基因股份有限公司 一种基于血液循环肿瘤dna的拷贝数变异检测方法和装置
CN111261225B (zh) * 2020-02-06 2022-08-16 西安交通大学 一种基于二代测序数据的反转相关复杂变异检测方法
CN113496761B (zh) * 2020-04-03 2023-09-19 深圳华大生命科学研究院 确定核酸样本中cnv的方法、装置及应用
DE102020116178A1 (de) * 2020-06-18 2021-12-23 Analytik Jena Gmbh Verfahren zum Erkennen einer Amplifikationsphase in einer Amplifikation
CN111968701B (zh) * 2020-08-27 2022-10-04 北京吉因加科技有限公司 检测指定基因组区域体细胞拷贝数变异的方法和装置
CN114220481B (zh) * 2021-11-25 2023-09-08 深圳思勤医疗科技有限公司 基于全基因组测序完成待测样本的核型分析的方法、系统和计算机可读介质
CN114999573B (zh) * 2022-04-14 2023-07-07 哈尔滨因极科技有限公司 一种基因组变异检测方法及检测系统
CN114758720B (zh) * 2022-06-14 2022-09-02 北京贝瑞和康生物技术有限公司 用于检测拷贝数变异的方法、设备和介质
CN114864000B (zh) * 2022-07-05 2022-09-09 北京大学第三医院（北京大学第三临床医学院） 一种动态鉴定人类单细胞染色体拷贝数的方法
CN115132271B (zh) * 2022-09-01 2023-07-04 北京中仪康卫医疗器械有限公司 一种基于批次内校正的cnv检测方法
CN117334249A (zh) * 2023-05-30 2024-01-02 上海品峰医疗科技有限公司 基于扩增子测序数据检测拷贝数变异的方法、设备和介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2004044225A2 (en) * 2002-11-11 2004-05-27 Affymetrix, Inc. Methods for identifying dna copy number changes
CN101449161A (zh) * 2006-05-03 2009-06-03 人口诊断公司 评估遗传病缺陷的方法
WO2011017596A2 (en) * 2009-08-06 2011-02-10 University Of Virginia Patent Foundation Compositions and methods for identifying and detecting sites of translocation and dna fusion junctions
WO2011030838A1 (ja) * 2009-09-10 2011-03-17 富士フイルム株式会社 アレイ比較ゲノムハイブリダイゼーション法による核酸変異解析法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7979215B2 (en) * 2007-07-30 2011-07-12 Agilent Technologies, Inc. Methods and systems for evaluating CGH candidate probe nucleic acid sequences

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2004044225A2 (en) * 2002-11-11 2004-05-27 Affymetrix, Inc. Methods for identifying dna copy number changes
CN101449161A (zh) * 2006-05-03 2009-06-03 人口诊断公司 评估遗传病缺陷的方法
WO2011017596A2 (en) * 2009-08-06 2011-02-10 University Of Virginia Patent Foundation Compositions and methods for identifying and detecting sites of translocation and dna fusion junctions
WO2011030838A1 (ja) * 2009-09-10 2011-03-17 富士フイルム株式会社 アレイ比較ゲノムハイブリダイゼーション法による核酸変異解析法

Non-Patent Citations (10)

* Cited by examiner, † Cited by third party
Title
Accurate and exact CNV identification from targeted high-throughput sequence data;Nord A S,et al.,;《 BMC genomics》;20111231;第12卷(第1期);1～10 *
Adjustment of genomic waves in signal intensities from whole-genome SNP genotyping platforms;Diskin S J, et al.,;《Nucleic acids research》;20080910;第36卷(第19期);e126:1-6 *
Array CGH analsis of copy number variation identifies 1284 new genes variant in healthy white males: implications for association studies of complex diseases;Adam J.de Smith,et al.,;《Human Molecular Genetics》;20070731;第27卷(第23期);2783-2794 *
cn.MOPS:mixture of Poissons for discovering copy number variations in next-generation sequencing data with a low false discoverage rate;G.Klambauer,et al.,;《Nucleic Acids Research》;20120201;第40卷(第9期);1-14 *
Desiging a simple multiplex ligation-dependent probe amplification(MLPA) assay for rapid detection of copy number variants in the genome;Yiping Shen,et al.,;《JOURNAL OF GENETICS AND GENOMICS》;20090401;第36卷(第4期);257-265 *
Experimental strategies for studying transcription factor–DNA binding specificities;Geertz M, et al.,;《Briefings in functional genomics》;20100923;第9卷(第5－6期);1－12 *
High-resolution mapping of copy-number alterations with massively parallel sequencing;Derek，Y.Chiang,et al.,;《Nature Methods》;20090131;第6卷(第1期);99～103 *
Sensitive and accurate detection of copy number variants using read depth of coverage;S.YOON ET AL;《Genome Research》;20090805;第19卷(第9期);1586～1592 *
基因组拷贝数变异研究进展;梁燕，等;《现代生物医学进展》;20120228;第12卷(第5期);964－967 *
拷贝数变异及其研究进展;刘欣，等;《畜牧兽医学报》;20101231;第41卷(第8期);927－931 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109979535A (zh) * 2017-12-28 2019-07-05 安诺优达基因科技(北京)有限公司 一种胚胎植入前遗传学筛查装置
CN109979535B (zh) * 2017-12-28 2021-03-02 浙江安诺优达生物科技有限公司 一种胚胎植入前遗传学筛查装置

Also Published As

Publication number Publication date
CN104221022A (zh) 2014-12-17
SG11201406250SA (en) 2014-11-27
KR101795124B1 (ko) 2017-12-01
WO2013149385A1 (zh) 2013-10-10
US20180148765A1 (en) 2018-05-31
AU2012376134B2 (en) 2016-03-03
US20150056619A1 (en) 2015-02-26
EP2835752A1 (en) 2015-02-11
IL234875B (en) 2019-03-31
AU2012376134A1 (en) 2014-11-06
US11371074B2 (en) 2022-06-28
RU2014144349A (ru) 2016-05-27
EP2835752B8 (en) 2018-12-26
EP2835752B1 (en) 2018-09-19
JP2015512264A (ja) 2015-04-27
KR20140140122A (ko) 2014-12-08
EP2835752A4 (en) 2015-11-18
JP5972448B2 (ja) 2016-08-17

Similar Documents

Publication Publication Date Title
CN104221022B (zh) 一种拷贝数变异检测方法和系统
CN104204220B (zh) 一种遗传变异检测方法
JP7051900B2 (ja) 不均一分子長を有するユニーク分子インデックスセットの生成およびエラー補正のための方法およびシステム
CN106834474B (zh) 利用基因组测序诊断胎儿染色体非整倍性
CN105392894B (zh) 确定样本基因组中是否存在拷贝数变异的方法、系统和计算机可读介质
CN104232777B (zh) 同时确定胎儿核酸含量和染色体非整倍性的方法及装置
AU2015314114B2 (en) Detecting repeat expansions with short read sequencing data
CN110800063B (zh) 使用无细胞dna片段大小检测肿瘤相关变体
EP2935625B1 (en) Noninvasive detection of fetal aneuploidy in egg donor pregnancies
CN113366122B (zh) 游离dna末端特征
CN110291212A (zh) 使用核酸片段的诊断应用
KR20200080272A (ko) 비침습적 산전 검사 및 암 검출을 위한 핵산 크기 범위의 용도
CN104169929A (zh) 用于确定胎儿是否存在性染色体数目异常的方法、系统和计算机可读介质
CN104968802B (zh) 作为诊断标志物的新miRNA
EP4004238A1 (en) Systems and methods for determining tumor fraction
TR201815541T4 (tr) Fetüslü gebe bir dişi denekten alınan bir biyolojik numunenin analiz yöntemi.
CN105555970A (zh) 同时进行单体型分析和染色体非整倍性检测的方法和系统
CN108604258B (zh) 染色体异常判断方法
CN104951671A (zh) 基于单样本外周血检测胎儿染色体非整倍性的装置
CN115989544A (zh) 用于在基因组的重复区域中可视化短读段的方法和系统
CN105765076A (zh) 一种染色体非整倍性检测方法及装置
CN110373458B (zh) 一种地中海贫血检测的试剂盒及分析系统
CN114566224B (zh) 一种用于鉴别或区分不同海拔人群的模型及其应用
CN109280697A (zh) 利用孕妇血浆游离dna进行胎儿基因型鉴定的方法
CN105678110B (zh) 一种样本组合分析核酸序列的方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
CB02 Change of applicant information
CB02 Change of applicant information

Address after: Yantian District of Shenzhen City, Guangdong province 518083 Hongan street No. 21 China Comprehensive Park 7 Building 7 layer -14 layer

Applicant after: BGI SHENZHEN CO LTD

Address before: Yantian District of Shenzhen City, Guangdong province 518083 North Road No. 146 North Industrial Zone 11, floor 3, 2

Applicant before: Shenzhen BGI Medicine Co., Ltd.

GR01 Patent grant
GR01 Patent grant