CN108573125B - 一种基因组拷贝数变异的检测方法及包含该方法的装置 - Google Patents

一种基因组拷贝数变异的检测方法及包含该方法的装置 Download PDF

Info

Publication number
CN108573125B
CN108573125B CN201810353495.4A CN201810353495A CN108573125B CN 108573125 B CN108573125 B CN 108573125B CN 201810353495 A CN201810353495 A CN 201810353495A CN 108573125 B CN108573125 B CN 108573125B
Authority
CN
China
Prior art keywords
window
data
correction
copy number
genome
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
CN201810353495.4A
Other languages
English (en)
Other versions
CN108573125A (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.)
Yikon Genomics Shanghai Co ltd
Original Assignee
Yikon Genomics Shanghai 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 Yikon Genomics Shanghai Co ltd filed Critical Yikon Genomics Shanghai Co ltd
Priority to CN201810353495.4A priority Critical patent/CN108573125B/zh
Publication of CN108573125A publication Critical patent/CN108573125A/zh
Application granted granted Critical
Publication of CN108573125B publication Critical patent/CN108573125B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations

Landscapes

  • Bioinformatics & Cheminformatics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Genetics & Genomics (AREA)
  • Biotechnology (AREA)
  • Biophysics (AREA)
  • Chemical & Material Sciences (AREA)
  • Molecular Biology (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Analytical Chemistry (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明提供了一种基因组拷贝数变异的检测方法及包含该方法的装置，所述方法包括输入原始数据、质控清理、将序列比对到参考基因组、利用大小不同的窗口计算唯一比对序列数、GC矫正、参考矫正、屏蔽不可检测区域、CBS分段、核型报告整合与生成报告的步骤，通过实验摸索优化，建立了一整套完整的检测方法和装置，通过特定顺序步骤的连用，创造性地采用参考校正的步骤，并选用大小不同的窗口进行比对整合，各步骤相互配合，最终提高敏感性和特异性，使检测准确度和结果形式能够符合临床需求，自动化程度高、易于扩展，检测的准确度高，能够降低数据分析的成本，具有极高的应用价值。

Description

一种基因组拷贝数变异的检测方法及包含该方法的装置
技术领域
本发明涉及生物信息学技术领域，尤其涉及一种基因组拷贝数变异的检测方法及包含该方法的装置。
背景技术
拷贝数变异(Copy number variation，CNV)主要指基因组的DNA片段大小从1kb到几个Mb范围内的缺失、插入、重复等，包括数目异常和结构异常。拷贝数嵌合指的是染色体拷贝数变异数量介于整数之间，比如2.5倍的拷贝数变异为50％的三倍体嵌合。基因组拷贝数变异检测的应用领域除了科学研究外，在临床应用领域，拷贝数变异检测可用于单细胞拷贝数变异检测、流产物组织拷贝数变异检测、拷贝数变异的遗传病和肿瘤的检测等；其中单细胞拷贝数变异指的是分析单个细胞的拷贝数变异，比如单精子、单个受精卵细胞、单个卵细胞等。
单细胞拷贝数变异检测的主要应用领域之一是胚胎植入前遗传学筛查(Preimplantation Genetic Screening，PGS)和胚胎植入前遗传学诊断(PreimplantationGenetic diagnostics，PGD)，主要对早期胚胎的单个细胞进行染色体拷贝数异常的检测，通过一次性检测胚胎23对染色体的结构和数目，分析胚胎是否有遗传物质异常的一种早期产前筛查/诊断方法，从而挑选正常的胚胎植入子宫，以期获得正常的妊娠，提高患者的临床妊娠率，降低胎儿患病风险。流产物组织拷贝数变异检测的目的是通过检测流产胎儿或胎盘样本，检测是不是胎儿基因组拷贝数变异导致的导致流产的。拷贝数变异的遗传病和肿瘤的检测主要是对病人进行基因组拷贝数变异检测用于发现导致遗传病/肿瘤的原因，并根据这些信息选择合适的药物进行精准治疗，可见，基因组拷贝数变异检测在整个生物医学领域有着重要的应用价值。
目前全基因组拷贝数变异检测的主要方法有：基因芯片，如比较基因组杂交(comparative Genomic Hybridization，array CGH)和二代测序(Next GenerationSequencing，NGS)，其它方法只是针对单个基因或目标片段分析，无法覆盖整个基因组；基因芯片的方法通量低、分辨率低，不能检测到精确的断点；二代测序的方法具有更高通量、更精细的分辨率、能够更精确地检测到断点以及更低的价格。然而，二代测序技术也有一定的缺陷；虽然二代测序技术可以产生大量的数据，但如何处理和分析这些数据成为制约二代测序技术用于拷贝数变异检测临床应用的主要瓶颈；目前二代测序技术的数据分析方面还存在的问题主要是检测结果的假阳性高、准确度低。
CN104133914A提供了一种消除高通量测序引入的GC偏差及对染色体拷贝数变异的检测方法，通过对人类基因组进行处理，并结合高通量测序得到的基因序列进行比对，对基因序列进行校正后，在染色体间做T-test，从而判断混合样本中染色体是否存在整倍体变异，很好地解决了高通量测序引入的GC偏差的技术问题，从而使得高通量测序在混合样本中染色体拷贝数变异的检测上的应用成为可能。CN106845154A涉及一种FFPE样本拷贝数变异检测装置包括测序数据获取模块、序列比对模块、前期数据处理模块、归一化模块、背景库筛选模块、数据波动消除模块、GC校正模块以及输出模块；CN105574361A涉及一种检测基因组拷贝数变异的方法，具体包括以下步骤：对样本基因组进行测序，以获得基因组序列；将序列比对到参考基因组，得到序列在基因组上的位置；将参考基因组分成一定长度的窗口，统计落在每个窗口的序列及碱基；根据每个窗口的序列及碱基GC含量，对每个窗口做校正；确定拷贝数正常的阈值，扫描每个窗口，确定窗口拷贝数是否变异；精确扫描异常的窗口，以确定精确的断点，来确定拷贝数变异的具体位置；但上述现在技术的检测结果假阳性高、准确度低，检测流程繁琐，步骤冗余，有待进一步提高和优化。
因此，需要高敏感性和特异性的生物信息学方法，以便于将基于二代测序的拷贝数变异检测技术更广泛地应用于临床。
发明内容
针对现有技术的不足及实际的需求，本发明提供一种基因组拷贝数变异的检测方法及包含该方法的装置，经过实验摸索优化，建立了一整套完整的检测方法和装置，通过特定顺序步骤的连用，创造性地采用参考校正的步骤，并选用大小不同的窗口进行比对整合，各步骤相互配合，最终提高敏感性和特异性，使检测准确度和结果形式能够符合临床需求，自动化程度高、易于扩展，检测的准确度高，能够降低数据分析的成本，具有极高的应用价值。
为达上述目的，本发明采用以下技术方案：
第一方面，本发明提供一种基因组拷贝数变异的检测方法，所述检测方法包括如下步骤：
(1)获取样本的原始数据并进行质控和清理；
(2)将步骤(1)得到的数据与参考基因组进行比对，然后排序并去重复；
(3)将参考基因组分成至少两个大小不同的窗口，计算落入窗口的唯一比对的序列数；
(4)统计落入步骤(3)所述窗口的GC含量，并进行GC校正；
(5)将步骤(4)得到的GC矫正后的窗口计数结果的中位数的倒数作为权重进行参考校正；
(6)将步骤(5)得到的数据进行CBS算法进行分段，找到基因组拷贝数相同的区域和不同的区域；
(7)将步骤(6)得到的CBS分段结果针对单个窗口进行全基因组核型分析，得到初步核型结果；
(8)将步骤(7)得到的不同大小窗口的核型结果进行整合，得到最终核型结果；
(9)对步骤(1)-(8)得到的数据进行加工处理整合，生成报告。
发明人在长期生产研究中，总结现有技术的优缺点，通过大量试验进行摸索优化，建立了一整套完整的检测方法和装置，通过特定顺序步骤的连用，创造性地采用参考校正的步骤，并选用大小不同的窗口进行比对整合，各步骤相互配合，最终提高敏感性和特异性，使检测准确度和结果形式能够符合临床需求，自动化程度高、易于扩展，检测的准确度高，能够降低数据分析的成本，具有极高的应用价值。
本发明中，发明人经实验论证后发现，整个检测方法的顺序要按照数据清洗与转换的合理逻辑为准则，不能颠倒，从而才能够保证检测方法的流程完整性，即数据处理的过程与过程的顺序：数据处理的过程和过程的顺序为对原始数据质控，去除掉低质量的数据、将数据比对到基因组上获得唯一比对序列、GC矫正、参考矫正。
不仅如此，发明人通过在检测流程中引入参考矫正的步骤，可以将测序仪、试剂盒等技术噪音屏蔽掉，因此可以减少这些方面的假阳性和假阴性，而具体的参考矫正采用正常的样本的每个GC矫正后的拷贝数值的中位数的倒数作为该窗口参考矫正的权重来进行。
另外，发明人吸收归纳现有技术的优点，采用多窗口扫描整合的方法，同时利用两个或多个不同大小的窗口对基因组扫描，计算拷贝数变异，然后对这些不同大小窗口的结果进行整合，由于每个样本的数据量是一定的，窗口越小则落入每个窗口的数据越少，数据方差越大越不稳定，因此一些大的拷贝数会在小窗口的时候检测不出来，因此需要多种窗口并行运行，最后对这些结果进行整合，以减少假阴性的发生。
与此同时，发明人提供的方法，能够根据CBS分段的结果自动化报告出符合临床需求的核型结果，能够推算性别、染色体数目和结构异常以及染色体嵌合及嵌合比例。
优选地，步骤(1)所述样品包括单细胞、少量混合的细胞、痕量DNA或组织的基因组DNA中的任意一种或至少两种的组合。
所述样本经过DNA提取、全基因组扩增、建库、测序等实验步骤获得测序数据。
优选地，步骤(1)所述获取数据的测序仪包括Illumina平台、Ion Torrent平台或DA8600平台的测序仪。
优选地，步骤(1)所述数据的格式包括FASTQ和/或BAM格式。
优选地，步骤(1)所述质控和清理的软件包括Trimmomatic、cutadapt、FASTQC或fastp中的任意一种或至少两种的组合。
本发明中，质控和清理的目的是去掉接头序列、低质量的序列、切掉序列中低质量的部分和去掉长度太短的序列。
优选地，步骤(2)所述参考基因组包括UCSC hg19、UCSC hg38、GRCh37或GRCh38中的任意一种或至少两种的组合。
优选地，步骤(2)所述比对软件包括TMAP(torrent mapping alignmentprogram)，BWA(Burrows-Wheeler Aligner)，Bowtie/Bowtie2，SOAP/SOAP2(ShortOligonucleotide Analysis Package)中的任意一种或至少两种的组合，所用参数为默认参数。
优选地，步骤(2)所述排序和去重复所用软件包括Samtools、Picard或GATK中的任意一种或至少两种的组合。
优选地，步骤(2)所述排序和去重复后得到的数据格式为BAM。
优选地，步骤(3)所述窗口大小包括1000K、500K、100K、50K、10K或1K中的至少两种的组合，例如可以是1000K和500K，500K和50K，100K和10K，1000K和100K或1000K和50K。
不同大小窗口的选择依赖于要检测的最小拷贝数变异的大小限制和数据量的多少，为了避免误差，会选择大小不同的多个窗口进行组合从而检测。计数时去掉比对质量低的序列、PCR重复序列、未比对到参考基因组上的序列以及比对到参考基因组上面多个位置不唯一比对的序列，只统计能比对到参考基因组上面唯一比对的序列。
优选地，步骤(4)所述GC校正的方法包括局部加权回归散点平滑法或GC梯度变化的倒数权重法。
1)局部加权回归散点平滑法(locally weighted scatterplot smoothing，LOWESS)：LOWESS主要思想是取一定比例的局部数据，在这部分子集中拟合多项式回归曲线，这样我们便可以观察到数据在局部展现出来的规律和趋势，再将局部范围从左往右依次推进，最终一条连续的曲线就被计算出来了。
2)GC梯度变化的倒数权重法：该方法的核心思想是将窗口GC含量从小到大按照一个固定的梯度如0.1％增加，形成一系列的GC含量梯度值，将含有相同GC梯度含量的窗口分到一类，将该类实际计算到的每个窗口的序列数的平均数的倒数作为该类窗口矫正的权重，计算步骤如下：首先，根据参考基因组GC含量，确定参考基因组GC含量的最大值和最小值；第二，从最小值到最大值按照0.1％的递增得到GC含量的浓度梯度值；第三，把有相同GC含量浓度梯度的窗口内的序列数取平均值Mi，i为不同特定GC浓度梯度；第四，计算不同特定GC浓度梯度的权重，w＝avg(M)/Mi，avg(M)为所有窗口里面的序列数的平均值；最后，每个窗口矫正后的值为该窗口的权重乘以每个窗口统计出来的序列值。
本发明中，GC含量定义为一段序列中G+C/序列长度；由于每个窗口的GC含量不同，会引起测序数据相应的窗口中的序列分布不均匀，引起GC偏好。GC偏好对于拷贝数变异检测的准确性影响很大，因此需要做GC矫正。
优选地，步骤(5)所述参考校正的值为该窗口GC校正后的值乘以参考数据每个窗口的权重。
本发明中，除了GC含量影响数据的偏好性外，不同的基因组区域也会影响检测的准确性；正常人基因组里面有一些异染色质区域和着丝粒区域或者重复区域导致该区域无法矫正到正常水平，因此需要进行参考矫正。将一批正常人的单细胞样本数据进行上述步骤分析，得到GC矫正后的窗口计数结果，再计算每个窗口的中位数的倒数作为参考矫正的权重，每个窗口参考矫正后的值为该窗口GC矫正后的值乘以参考数据每个窗口的权重。
优选地，步骤(5)所述参考校正的步骤还包括屏蔽掉异常不可检测区域的步骤，所述异常不可检测区域的权重定义为0，通过将不可检测区域的权重定义为0会降低假阳性，减少这些区域的检测的随机性的影响。
优选地，所述异常不可检测区域包括窗口GC含量为0、参考计数的中位数为0或参考计数的变异系数大于0.2的数据窗口。
优选地，步骤(6)所述CBS分段采用的R语言软件包包括DNAcopy、seqCBS或PSCBS中的任意一种或至少两种的组合，根据CBS分段结果，自动化报告出符合临床需求的核型结果，能够推算性别、染色体数目和结构异常以及染色体嵌合及嵌合比例。
优选地，步骤(6)所述CBS分段之前还包括将步骤(5)得到的数据的基因组每个窗口在染色体上位置按照p、q臂分开的步骤，避免着丝粒对分段的影响，有利于结果的报告和解读。
优选地，步骤(7)所述初步核型结果包括：染色体总数、性别、异常染色体数目、嵌合比例和异常片段的位置。
具体地，步骤(7)所述核型分析的步骤包括：核型报告与画图，根据质控结果、分段结果、染色体区段(cytoband)信息对拷贝数变异(CNV)加工处理，确定拷贝数增加或减少、拷贝数变化的多少、拷贝数变化的基因组位置和染色体区段、嵌合比例以及是否为整条染色体/染色体臂/小片段的变异等信息，报告出核型结果并画出全基因组和每条染色体拷贝数变异图谱。
优选地，步骤(8)所述核型整合的方法包括：通过比较不同大小窗口核型结果的异常片段的染色体位置、起始位置、终止位置、拷贝数变异数目和嵌合比例，从而保留大的拷贝数变异，舍去被包含的小的拷贝数变异，保留小窗口得到的更高分辨率的拷贝数变异。
具体地，步骤(8)的步骤包括：对1000K和100K窗口的核型报告结果进行整合然后得到一个最终的整合后的核型报告结果；首先，程序读取1000K和100K窗口的核型报告结果；第二，对于每个1000K和100K的变异的核型进行两两比较，比较二者是否为同一条染色体，大小范围是否相同/包含关系，最后，根据核型比较结果进行取舍；
优选地，步骤(9)的具体步骤包括：将整合后的核型结果和质控信息等信息写入预先准备好的报告模板里面，自动生成报告。
作为优选技术方案，一种基因组拷贝数变异的检测方法，具体包括如下步骤：
(1)获取样本的原始数据并进行质控和清理；
(2)将步骤(1)得到的数据与参考基因组进行比对，排序并去重复；
(3)将参考基因组分成至少两个大小不同的窗口，计算落入窗口的唯一比对的序列数；
(4)统计落入步骤(3)所述窗口的GC含量，并进行GC校正；
(5)将步骤(4)得到的GC矫正后的窗口计数结果的中位数的倒数作为权重进行参考校正，并屏蔽掉异常不可检测区域；
(6)将步骤(5)得到的数据的基因组每个窗口在染色体上位置按照p、q臂分开，然后进行CBS算法进行分段，找到基因组拷贝数相同的区域和不同的区域；
(7)将步骤(6)得到的CBS分段结果针对单个窗口进行全基因组核型分析，得到初步核型结果，包括染色体总数、性别、异常染色体数目、嵌合比例和异常片段的位置；
(8)将步骤(7)得到的不同大小窗口的核型结果进行整合，通过比较不同大小窗口核型结果的异常片段的染色体位置、起始位置、终止位置、拷贝数变异数目和嵌合比例，从而保留大的拷贝数变异，舍去被包含的小的拷贝数变异，保留小窗口得到的更高分辨率的拷贝数变异，得到最终核型结果；
(9)对步骤(1)-(8)得到的数据进行加工处理整合，生成报告。
第二方面本发明提供一种包含第一方面所述方法的装置，包括如下模块：
(1)测序数据获取模块：用于获取样本的原始数据并进行质控和清理；
(2)序列比对模块：用于将测序数据与参考基因组进行比对，排序并去重复；
(3)数据处理模块：用于将参考基因组分成至少两个大小不同的窗口，计算落入窗口的唯一比对的序列数；
(4)GC校正模块：用于统计落入数据处理模块所述窗口的GC含量，并进行GC校正；
(5)参考校正模块：用于将GC矫正后得到的窗口计数结果的中位数的倒数作为权重进行参考校正，并屏蔽掉异常不可检测区域；
(6)CBS分段模块：用于将参考校正后得到的数据的基因组每个窗口在染色体上位置按照p、q臂分开，然后进行CBS算法进行分段，找到基因组拷贝数相同的区域和不同的区域；
(7)核型分析模块：将步骤(6)得到的CBS分段结果针对单个窗口进行全基因组核型分析，得到初步核型结果；
(8)核型整合模块：将步骤(7)得到的不同大小窗口的核型结果进行整合，得到最终核型结果；
(9)报告输出模块：对步骤(1)-(8)得到的数据进行加工处理整合，生成报告。
第三方面，本发明提供一种计算机可读存储介质，所述存储介质存储有计算机可执行指令，用于执行第一方面所述方法和/或第二方面所述装置的指令。
与现有技术相比，本发明具有如下有益效果：
本发明提供的基于高通量测序背景下的基因组拷贝数变异检测的方法和装置，能够自动化报告出用于临床的核型结果，提高敏感性和特异性，降低假阳性和假阴性，使检测准确度和结果形式能够符合临床需求，装置的自动化程度高、易于扩展，检测的准确度高，能够降低数据分析的成本，具有极高的应用价值。
附图说明
图1为本发明的数据流程图；
图2为本发明的1000K被屏蔽窗口分布图；
图3为本发明的100K被屏蔽窗口分布图；
图4为本发明的1000K全基因组CNV图；
图5为本发明的100K全基因组CNV图；
图6为本发明的1000K异常染色体chr21拷贝数变异图；
图7为本发明的100K异常染色体chr21拷贝数变异图；
图8为本发明对比例的有参考矫正和无参考矫正的拷贝数变异系数(CV)的对比图。
具体实施方式
为更进一步阐述本发明所采取的技术手段及其效果，以下结合附图并通过具体实施方式来进一步说明本发明的技术方案，但本发明并非局限在实施例范围内。
实施例1检测基因拷贝数异常的装置的组装
将以下模块组装为检测基因组拷贝数变异的装置：
(1)测序数据获取模块：用于获取样本的原始数据并进行质控和清理；
(2)序列比对模块：用于将测序数据与参考基因组进行比对，排序并去重复；
(3)数据处理模块：用于将参考基因组分成至少两个大小不同的窗口，计算落入窗口的唯一比对的序列数；
(4)GC校正模块：用于统计落入数据处理模块所述窗口的GC含量，并进行GC校正；
(5)参考校正模块：用于将GC矫正后得到的窗口计数结果的中位数的倒数作为权重进行参考校正，并屏蔽掉异常不可检测区域；
(6)CBS分段模块：用于将参考校正后得到的数据的基因组每个窗口在染色体上位置按照p、q臂分开，然后进行CBS算法进行分段，找到基因组拷贝数相同的区域和不同的区域；
(7)核型分析模块；将步骤(6)得到的CBS分段结果针对单个窗口进行全基因组核型分析，得到初步核型结果；
(8)核型整合模块；将步骤(7)得到的不同大小窗口的核型结果进行整合，得到最终核型结果；
(9)报告输出模块：对步骤(1)-(8)得到的数据进行加工处理整合，生成报告。
实施例2
本发明中，采用实施例1中的装置进行拷贝数变异的检测，完整数据流程图见图1所示，具体步骤如下；
1.对样本进行全基因组扩增、建库、测序
在本实施例中，检测样本为胚胎植入前染色体非整倍体国家参考品，该参考品用于高通量测序法胚胎植入前染色体非整倍体检测试剂盒的性能评价，评价高通量测序法用于囊胚筛选中染色体不同大小CNV的检测能力。
全基因组扩增方法选用MALBAC-LAB早期胚胎植入前染色体非整倍体检测文库制备试剂盒，扩增建库方法按照上海亿康医学检验所有限公司提供的产品说明书操作。
上机测序采用达安基因公司的DA8600高通量测序平台(Ion Torrent平台)，按照达安基因公司提供的说明书操作，测序类型为单端(Single End)测序，测序数据量为2M左右，最终下机数据的格式为BAM。
2.原始数据格式转换和质控
将原始BAM格式数据转换成FASTAQ格式，对FASTAQ数据进行质控和清理，去掉接头序列、低质量的序列、切掉序列中低质量的部分和去掉长度太短的序列，所采用的软件为Trimmomatic，所用的参数为：ILLUMINACLIP:ADAPTERS:2:20:6SLIDINGWINDOW:4:15LEADING:3TRAILING:3MINLEN:25HEADCROP:12。
3.将序列比对到参考基因组，排序和去除重复
将上一步输出的高质量纯净数据和人参考基因组(UCSC hg19)进行比对，所用软件为BWA(Burrows-Wheeler Alignment tool)，再用Picard软件将比对生成的BAM文件排序、去除重复最后生成最终的BAM文件，所用参数均为默认参数。
4.将参考基因组分成1000K和100K大小的窗口，同时进行计算落入窗口中的唯一比对的序列数
将人的基因组按1000K和100K两个不同大小的窗口分割，分别统计落入每个窗口的序列数，两种窗口并行计算，计数时去掉比对质量低的序列、PCR重复序列、未比对到参考基因组上的序列以及对对到参考基因组上面多个位置不唯一比对的序列，只统计能比对到参考基因组上面唯一比对的序列。
5.GC矫正
GC矫正采用GC梯度变化的倒数权重法，梯度值定位0.1％，步骤如下：首先，根据参考基因组GC含量，确定参考基因组GC含量的最大值和最小值；第二，从最小值到最大值按照0.1％的递增得到GC含量的浓度梯度值；第三，把有相同GC含量浓度梯度的窗口内的序列数取平均值Mi，i为不同特定GC浓度梯度；第四，计算不同特定GC浓度梯度的权重，w＝avg(M)/Mi，avg(M)为所有窗口里面的序列数的平均值；最后，每个窗口矫正后的值为该窗口的权重乘以每个窗口统计出来的序列值。
6.参考矫正
将国家参考品种的正常人的单细胞样本数据进行上述步骤分析，得到GC矫正后的窗口计数结果，再计算每个窗口的中位数的倒数作为参考矫正的权重，每个窗口参考矫正后的值为该窗口GC矫正后的值乘以参考数据每个窗口的权重，实施例1的参考矫正的参考文件及权重的部分数据如表1所示，由于数据量过于庞大，表1展示前20行的部分数据，能够代表完整结果；
表1：100K窗口的参考矫正的参考文件及权重(部分)
bin_ID chr start end bin_mean bin_sd bin_median bin_CV mask weight
1 chr1 1 100000 0.275915 0.072238 0.279693 0.261813 1 0
2 chr1 100001 200000 0.06241 0.032531 0.057854 0.521252 1 0
3 chr1 200001 300000 0.054662 0.030796 0.050616 0.563387 1 0
4 chr1 300001 400000 0.000629 0.003219 0 5.116599 1 0
5 chr1 400001 500000 0.000903 0.004118 0 4.560376 1 0
6 chr1 500001 600000 0.10713 0.048043 0.105956 0.448455 1 0
7 chr1 600001 700000 0.063918 0.031424 0.06496 0.49163 1 0
8 chr1 700001 800000 0.581337 0.111411 0.573815 0.191645 0 1.742721
9 chr1 800001 900000 0.842661 0.13989 0.809979 0.16601 0 1.234599
10 chr1 900001 1000000 1.001536 0.001653 1.001149 0.00165 0 0.998852
11 chr1 1000001 1100000 1.099251 0.13463 1.092497 0.122475 0 0.915335
12 chr1 1100001 1200000 0.903832 0.128415 0.898397 0.142078 0 1.113093
13 chr1 1200001 1300000 1.001536 0.001653 1.001149 0.00165 0 0.998852
14 chr1 1300001 1400000 0.677938 0.113173 0.675723 0.166937 0 1.479897
15 chr1 1400001 1500000 0.844607 0.139653 0.836856 0.165347 0 1.194949
16 chr1 1500001 1600000 0.562343 0.096328 0.560479 0.171298 0 1.784189
17 chr1 1600001 1700000 0.681637 0.102328 0.681244 0.150121 0 1.467902
18 chr1 1700001 1800000 0.822244 0.12059 0.804463 0.146659 0 1.243065
19 chr1 1800001 1900000 0.76811 0.109957 0.767425 0.143152 0 1.303059
20 chr1 1900001 2000000 0.917051 0.126308 0.909784 0.137733 0 1.099162
由表1可知，基因组不同窗口的权重不同，mask＝1的区域的权重为0，是不可检测区域，该区域由于重复序列多、包含N的未知碱基多或者位于着丝粒区域，同时该区域的变异系数CV(bin_cv)要比正常区域的大。
7.屏蔽掉异常不可检测区域。
本实施例中定义窗口GC含量为0，参考计数的中位数为0，参考计数的变异系数(CV)大于0.2的数据均为被屏蔽的窗口，将这些区域的权重定义为0。本实施例中的1000K窗口的不可检测区域在全基因组上的位置和分布如图2所示，100K窗口的不可检测区域分布在全基因组上的位置和分布如图3所示。
由图2可知，位于着丝粒区域和染色体末端的端粒区域通常是不可检测区域，此外还有其它区域。
由图3可知，由于窗口变小，可检测的拷贝数变异的分辨率提高，可检测的拷贝数变异更小，位于着丝粒区域和染色体末端的端粒区域通常是不可检测区域，此外还有其它区域。
8.CBS分段
对上述GC矫正和参考矫正好的数据进行分段(segment)，运用R语言软件包DNAcopy对基因组拷贝数数据进行分段，找到基因组拷贝数相同的区域和不同的区域，所用的参数为：alpha＝0.05，nperm＝10000，p.method＝"hybrid"，undo.splits＝"sdundo"，undo.SD＝2.5，verbose＝1，min.width＝2
9.核型报告与画图
根据质控结果、分段结果、染色体区段(cytoband)信息对拷贝数变异(CNV)加工处理，确定拷贝数增加或减少、拷贝数变化的多少、拷贝数变化的基因组位置和染色体区段、嵌合比例以及是否为整条染色体/染色体臂/小片段的变异等信息，报告出核型结果并画出全基因组和每条染色体拷贝数变异图谱。
10.核型报告整合
对1000K和100K窗口的核型报告结果进行整合然后得到一个最终的整合后的核型报告结果；
首先，程序读取1000K和100K窗口的核型报告结果；然后，对于每个1000K和100K的变异的核型进行两两比较，比较二者是否为同一条染色体，大小范围是否相同/包含关系，最后，根据核型比较结果进行取舍。
11.报告生成
最终，将整合后的核型结果和质控信息等信息写入预先准备好的报告模板里面，自动生成报告。
12.检测结果
对107例已知基因组拷贝数变异核型的胚胎植入前染色体非整倍体国家参考品细胞系样本进行上述方法的分析，得到核型结果见表2。
表2胚胎植入前染色体非整倍体国家参考品细胞系样本检测结果
Figure BDA0001634015550000161
Figure BDA0001634015550000171
Figure BDA0001634015550000181
Figure BDA0001634015550000191
Figure BDA0001634015550000201
由表2可知，本实施例的所有样本的检测结果均与参考品的标准答案一致。
挑选其中一例21三体的阳性样本作为例子来展示1000K窗口全基因组拷贝数变异图谱(图4)、100K窗口拷贝数变异图谱(图5)、1000K窗口特定染色体拷贝数变异图谱(图6)以及100K窗口特定染色体拷贝数变异图谱(图7)。
由图4可知，本方法及装置可以将1000K窗口的全基因组拷贝数信息输出到图示中。
由图5可知，本方法及装置可以将100K窗口的全基因组拷贝数信息输出到图示中。
由图6可知，本方法和装置可以将每个染色体的1000K窗口的拷贝数信息输出到图示中，并且将染色体位置和条带信息输出作为横坐标。用于更加清晰的展示该染色体的详细情况。
由图7可知，本方法和装置可以将每个染色体的100K窗口的拷贝数信息输出到图示中，并且将染色体位置和条带信息输出作为横坐标。用于更加清晰的展示该染色体的详细情况。
对比例
与实施例2相比，除不采用参考校正外，其他步骤与实施例2相同，检测结果见图8；
如图8可知，该图反映了在采用参考矫正和不采用参考矫正的情况下，这些样本的基因组的拷贝数变异系数(CV)的分布，由于CV越小数据越稳定，可以看出有参考矫正要比没参考矫正使得拷贝数的数据更加稳定、波动更小，从而使得拷贝数变异分析更加准确。
综上所述，本实施例的所有样本的检测结果均与参考品的标准答案一致，在2M数据量下可以检测1M以上非整倍体和10M以上嵌合，表明本方法的准确度高，假阳性和假阴性为0，报告核型形式和核型图谱能直接应用于临床。
申请人声明，本发明通过上述实施例来说明本发明的详细方法，但本发明并不局限于上述详细方法，即不意味着本发明必须依赖上述详细方法才能实施。所属技术领域的技术人员应该明了，对本发明的任何改进，对本发明产品各原料的等效替换及辅助成分的添加、具体方式的选择等，均落在本发明的保护范围和公开范围之内。

Claims (18)

1.一种基因组拷贝数变异的检测方法，其特征在于，所述检测方法包括如下步骤：
（1）获取样本的高通量测序原始数据并进行质控和清理；
（2）将步骤（1）得到的数据与参考基因组进行比对，排序并去重复；
对获得的数据并行地应用至少两种不同的窗口大小来运行以下步骤(3)-(7)；
（3）针对所述至少两种不同大小的窗口中的每一种，分别进行参考基因组的窗口分割，和计算落入窗口的唯一比对的序列数；
（4）统计落入步骤（3）所述窗口的GC含量，并进行GC校正，得到GC校正后的窗口计数结果；
（5）将步骤（4）得到的样本GC校正后的窗口计数结果，以参考数据的GC校正后的窗口计数结果的中位数的倒数作为权重，进行参考校正，其中通过所述样本GC校正后的窗口计数结果乘以所述权重，得到样本该窗口的参考校正的值，其中所述参考数据通过采用正常样本数据进行上述步骤(1)-(4)分析获得；
所述参考校正还包括：将异常不可检测区域的权重定义为0，屏蔽掉所述异常不可检测区域的步骤；其中，所述异常不可检测区域包括窗口GC含量为0、参考计数的中位数为0或参考计数的变异系数大于0.2的数据窗口；
（6）将步骤（5）得到的数据进行CBS算法进行分段，找到基因组拷贝数相同的区域和不同的区域；
（7）将步骤（6）得到的CBS分段结果针对每种窗口大小分别进行全基因组核型分析，得到初步核型结果；
（8）将步骤（7）得到的不同大小窗口的核型结果进行整合，得到最终核型结果；
（9）对步骤（1）-（8）得到的数据进行加工处理整合，生成报告。
2.根据权利要求1所述的检测方法，其特征在于，步骤（1）所述样本包括单细胞DNA、少量混合的细胞DNA、痕量DNA或组织的基因组DNA中的任意一种或至少两种的组合。
3.根据权利要求1所述的检测方法，其特征在于，步骤（1）所述获取样本的高通量测序原始数据的测序仪包括Illumina平台、Ion Torrent平台或DA8600平台的测序仪。
4.根据权利要求1所述的检测方法，其特征在于，步骤（1）所述数据的格式包括FASTQ和/或BAM格式。
5.根据权利要求1所述的检测方法，其特征在于，步骤（1）所述质控和清理的软件包括Trimmomatic、cutadapt、FASTQC或fastp中的任意一种或至少两种的组合。
6.根据权利要求1所述的检测方法，其特征在于，步骤（2）所述参考基因组包括UCSChg19、UCSC hg38、GRCh37或GRCh38中的任意一种或至少两种的组合。
7.根据权利要求1所述的检测方法，其特征在于，步骤（2）所述比对使用的软件包括TMAP、BWA、Bowtie/Bowtie2、SOAP或SOAP2中的任意一种或至少两种的组合。
8.根据权利要求1所述的检测方法，其特征在于，步骤（2）所述排序和去重复所用软件包括Samtools、Picard或GATK中的任意一种或至少两种的组合。
9.根据权利要求1所述的检测方法，其特征在于，步骤（2）所述排序和去重复后得到的数据格式为BAM。
10.根据权利要求1所述的检测方法，其特征在于，步骤（2）所述窗口大小包括1000K、500K、100K、50K、10K或1K中的至少两种的组合。
11.根据权利要求1所述的检测方法，其特征在于，步骤（4）所述GC校正的方法包括局部加权回归散点平滑法或GC梯度变化的倒数权重法。
12.根据权利要求1所述的检测方法，其特征在于，步骤（6）所述CBS分段采用的R语言软件包包括DNAcopy、seqCBS或PSCBS中的任意一种或至少两种的组合。
13.根据权利要求1所述的检测方法，其特征在于，步骤（6）所述CBS分段之前还包括将步骤（5）得到的数据根据基因组每个窗口在染色体上的位置按照p、q臂分开的步骤。
14.根据权利要求1所述的检测方法，其特征在于，步骤（7）所述初步核型结果包括：染色体总数、性别、异常染色体数目、嵌合比例和异常片段的位置。
15.根据权利要求1所述的检测方法，其特征在于，步骤（8）所述核型整合的方法包括：通过比较不同大小窗口核型结果的异常片段的染色体位置、起始位置、终止位置、拷贝数变异数目和嵌合比例，从而保留大的拷贝数变异，舍去被包含的小的拷贝数变异，保留小窗口得到的更高分辨率的拷贝数变异。
16.一种基因组拷贝数变异的检测方法，其特征在于，具体包括如下步骤：
（1）获取样本的高通量测序原始数据并进行质控和清理；
（2）将步骤（1）得到的数据与参考基因组进行比对，排序并去重复；
对获得的数据并行地应用至少两种不同的窗口大小来运行以下步骤(3)-(7)；
（3）针对所述至少两种不同大小的窗口中的每一种，分别进行参考基因组的窗口分割，和计算落入窗口的唯一比对的序列数；
（4）统计落入步骤（3）所述窗口的GC含量，并进行GC校正，得到GC校正后的窗口计数结果；
（5）将步骤（4）得到的样本GC校正后的窗口计数结果，以参考数据的GC校正后的窗口计数结果的中位数的倒数作为权重，进行参考校正，其中通过所述样本GC校正后的窗口计数结果乘以所述权重，得到样本该窗口的参考校正的值，其中所述参考数据通过采用正常样本数据进行上述步骤(1)-(4)分析获得，
所述参考校正还包括：将异常不可检测区域的权重定义为0，屏蔽掉所述异常不可检测区域；其中，所述异常不可检测区域包括窗口GC含量为0、参考计数的中位数为0或参考计数的变异系数大于0.2的数据窗口；
（6）将步骤（5）得到的数据根据基因组每个窗口在染色体上位置按照p、q臂分开，然后进行CBS算法进行分段，找到基因组拷贝数相同的区域和不同的区域；
（7）将步骤（6）得到的CBS分段结果针对每种窗口大小分别进行全基因组核型分析，得到初步核型结果，包括染色体总数、性别、异常染色体数目、嵌合比例和异常片段的位置；
（8）将步骤（7）得到的不同大小窗口的核型结果进行整合，通过比较不同大小窗口核型结果的异常片段的染色体位置、起始位置、终止位置、拷贝数变异数目和嵌合比例，从而保留大的拷贝数变异，舍去被包含的小的拷贝数变异，保留小窗口得到的更高分辨率的拷贝数变异，得到最终核型结果；
（9）对步骤（1）-（8）得到的数据进行加工处理整合，生成报告。
17.一种用于执行权利要求1-16中任一项所述方法的装置，包括如下模块：
（1）测序数据获取模块：用于获取样本的高通量测序原始数据并进行质控和清理；
（2）序列比对模块：用于将测序数据与参考基因组进行比对，排序并去重复；
（3）数据处理模块：用于将参考基因组按照至少两种大小不同的窗口并行地进行窗口划分，和计算落入窗口的唯一比对的序列数；
（4）GC校正模块：用于统计落入数据处理模块所述窗口的GC含量，并进行GC校正，得到GC校正后的窗口计数结果；
（5）参考校正模块：用于将参考数据GC校正后得到的窗口计数结果的中位数的倒数作为权重，对GC校正模块得到的样本GC校正后的窗口计数结果，进行参考校正，并屏蔽掉异常不可检测区域；
（6）CBS分段模块：用于将参考校正后得到的数据根据基因组的每个窗口在染色体上位置按照p、q臂分开，进行CBS算法进行分段，找到基因组拷贝数相同的区域和不同的区域；
（7）核型分析模块：将步骤（6）得到的CBS分段结果针对每种窗口大小分别进行全基因组核型分析，得到初步核型结果；
（8）核型整合模块：将步骤（7）得到的不同大小窗口的核型结果进行整合，得到最终核型结果；
（9）报告输出模块：对步骤（1）-（8）得到的数据进行加工处理整合，生成报告。
18.一种计算机可读存储介质，其特征在于，所述存储介质存储有计算机可执行指令，用于执行权利要求1-16中任一项所述方法。
CN201810353495.4A 2018-04-19 2018-04-19 一种基因组拷贝数变异的检测方法及包含该方法的装置 Active CN108573125B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810353495.4A CN108573125B (zh) 2018-04-19 2018-04-19 一种基因组拷贝数变异的检测方法及包含该方法的装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810353495.4A CN108573125B (zh) 2018-04-19 2018-04-19 一种基因组拷贝数变异的检测方法及包含该方法的装置

Publications (2)

Publication Number Publication Date
CN108573125A CN108573125A (zh) 2018-09-25
CN108573125B true CN108573125B (zh) 2022-05-13

Family

ID=63575233

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810353495.4A Active CN108573125B (zh) 2018-04-19 2018-04-19 一种基因组拷贝数变异的检测方法及包含该方法的装置

Country Status (1)

Country Link
CN (1) CN108573125B (zh)

Families Citing this family (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109545279B (zh) * 2018-11-29 2023-12-29 深圳市第二人民医院 染色体微阵列数据的分析方法、装置、设备及存储介质
CN110129419B (zh) * 2018-12-18 2023-03-31 华联生物科技股份有限公司 拷贝数变异的检测方法
CN109935275B (zh) * 2018-12-29 2021-09-07 北京安诺优达医学检验实验室有限公司 序列变异校验方法和装置、生产变异序列的方法和装置及电子设备
CN109801677B (zh) * 2018-12-29 2023-05-23 浙江安诺优达生物科技有限公司 测序数据自动化分析方法、装置和电子设备
CN110797088B (zh) * 2019-10-17 2020-09-15 南京医基云医疗数据研究院有限公司 全基因组重测序分析及用于全基因组重测序分析的方法
CN110797081B (zh) * 2019-10-17 2020-11-10 南京医基云医疗数据研究院有限公司 激活区域识别方法及装置、存储介质及电子设备
CN110628890B (zh) * 2019-11-07 2020-11-13 中国人民解放军军事科学院军事医学研究院 测序质控标准品及其应用与产品
CN111243666B (zh) * 2020-01-08 2023-04-07 华南理工大学 一种基于Nextflow的环状核糖核酸自动化分析方法及系统
CN111755069A (zh) * 2020-07-10 2020-10-09 苏州科贝生物技术有限公司 基于高通量测序分析拷贝数变异的方法
CN113113085B (zh) * 2021-03-15 2022-08-19 杭州杰毅生物技术有限公司 基于智能宏基因组测序数据肿瘤检测的分析系统及方法
CN112967756B (zh) * 2021-03-30 2022-07-26 上海欧易生物医学科技有限公司 基于snakemake语言快速批量可自动邮件反馈结果的高通量测序质控分析方法
CN113270138B (zh) * 2021-04-13 2023-09-22 杭州博圣医学检验实验室有限公司 基于生物信息学富集胎儿游离dna用于拷贝数变异的分析方法
CN113823353B (zh) * 2021-08-12 2024-02-09 上海厦维医学检验实验室有限公司 基因拷贝数扩增检测方法、装置及可读介质
CN114420208B (zh) * 2022-02-28 2023-04-18 上海亿康医学检验所有限公司 一种用于鉴定核酸样本中cnv的方法和装置
CN114758720B (zh) * 2022-06-14 2022-09-02 北京贝瑞和康生物技术有限公司 用于检测拷贝数变异的方法、设备和介质
CN114792548B (zh) * 2022-06-14 2022-09-09 北京贝瑞和康生物技术有限公司 校正测序数据、检测拷贝数变异的方法、设备和介质
CN114864000B (zh) * 2022-07-05 2022-09-09 北京大学第三医院（北京大学第三临床医学院） 一种动态鉴定人类单细胞染色体拷贝数的方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105574361A (zh) * 2015-11-05 2016-05-11 上海序康医疗科技有限公司 一种检测基因组拷贝数变异的方法
CN106520940A (zh) * 2016-11-04 2017-03-22 深圳华大基因研究院 一种染色体非整倍体和拷贝数变异检测方法及其应用

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2593708C2 (ru) * 2012-01-20 2016-08-10 БГИ Диагносис Ко., Лтд. Способ и система выявления вариации числа копий в геноме

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105574361A (zh) * 2015-11-05 2016-05-11 上海序康医疗科技有限公司 一种检测基因组拷贝数变异的方法
CN106520940A (zh) * 2016-11-04 2017-03-22 深圳华大基因研究院 一种染色体非整倍体和拷贝数变异检测方法及其应用

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于多元方差分析的成对肿瘤SNP array数据分段算法;习佳宁等;《科学通报》;20141120(第32期);全文 *
基于烟草基因组重测序数据的SNP提取软件组合比较;余世洲等;《烟草科技》;20171015(第10期);全文 *

Also Published As

Publication number Publication date
CN108573125A (zh) 2018-09-25

Similar Documents

Publication Publication Date Title
CN108573125B (zh) 一种基因组拷贝数变异的检测方法及包含该方法的装置
CN108920899B (zh) 一种基于目标区域测序的单个外显子拷贝数变异预测方法
CN112669901A (zh) 基于低深度高通量基因组测序的染色体拷贝数变异检测装置
CN108256289B (zh) 一种基于目标区域捕获测序基因组拷贝数变异的方法
CN111341383B (zh) 一种检测拷贝数变异的方法、装置和存储介质
CN111755068B (zh) 基于测序数据识别肿瘤纯度和绝对拷贝数的方法及装置
CN110268044B (zh) 一种染色体变异的检测方法及装置
JP2015506684A (ja) ゲノムのコピー数変異の有無を判断する方法、システム及びコンピューター読み取り可能な記憶媒体
WO2021232388A1 (zh) 确定胚胎细胞染色体中预定位点碱基类型的方法及其应用
CN106096330B (zh) 一种无创产前生物信息检测分析方法
CN110993029B (zh) 一种检测染色体异常的方法及系统
EP4152334A1 (en) Gene sequencing analysis method and apparatus, and storage medium and computer device
CN110016497B (zh) 一种检测肿瘤单细胞基因组拷贝数变异的方法
CN116030892B (zh) 一种鉴定染色体相互易位断点位置的系统和方法
Lun et al. From reads to regions: a Bioconductor workflow to detect differential binding in ChIP-seq data
CN106795551B (zh) 单细胞染色体的cnv分析方法和检测装置
US20160154931A1 (en) Method and device for detecting chromosomal aneuploidy
CN106591451B (zh) 测定胎儿游离dna含量的方法及其用于实施该方法的装置
CN107239676B (zh) 一种针对胚胎染色体的序列数据处理装置
CN114067908B (zh) 一种评估单样本同源重组缺陷的方法、装置和存储介质
TW202300656A (zh) 基因組序列上之拷貝數變異之候選斷點之機械性檢測
US20210366569A1 (en) Limit of detection based quality control metric
WO2024140881A1 (zh) 胎儿dna浓度的确定方法及装置
Hernandez-Lopez et al. Lossy compression of quality scores in differential gene expression: A first assessment and impact analysis
RU2712175C1 (ru) Способ неинвазивного пренатального скрининга анеуплоидий плода

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant