CN108763859B - 一种基于未知cnv样本建立提供cnv检测所需的模拟数据集的方法 - Google Patents

Abstract

本发明公开了一种基于未知CNV样本建立提供CNV检测所需的模拟数据集的方法。根据本发明的方法模拟生成了正常样本的序列数据，减少了样本测序的成本，解决了缺失配对样本和未知原因造成的测序偏好性的问题，提高了检测的准确度。

Description

一种基于未知CNV样本建立提供CNV检测所需的模拟数据集的 方法

技术领域

本发明属于生物信息学领域，公开了一种基于未知CNV样本建立提供CNV检测所需的模拟数据集的方法。

背景技术

拷贝数变异(CopyNumberVariations，CNV)是指与基因组参考序列相比，样本基因组染色体或染色体片段拷贝数异常，包括但不限于染色体非整倍体、缺失、重复，大于1000bp碱基的微缺失、微重复。在生物医学的科学研究及临床应用领域，经常遇到由于基因组拷贝数变异而引起的疾病，如染色体非整倍体、微缺失、微重复造成的流产，胚胎植入失败，各种遗传病以及癌症等等。基因组拷贝数变异的检测可应用到组织检测如肿瘤组织、羊水、流产物组织，液体活检如血液、尿液的细胞、游离核酸，单细胞领域如胚胎植入前遗传学诊断(PreimplantationGeneticDiagnosis，PGD)、胚胎植入前遗传学筛查(PreimplantationGeneticScreening，PGS)、癌症患者血液中游离的循环肿瘤细胞(CirculatingTumorCells，CTC)、孕妇外周血中游离的胎儿细胞、干细胞、单细胞或几个细胞的微生物。

目前基因组拷贝数变异检测的主要方法有：比较基因组杂交(ComparativeGenomicHybridization，CGH)，荧光定量PCR(realtimefluorescencequantitativePCR，RTFQPCR)，荧光原位杂交(FluorescenceInSituHybridization，FISH)，多重连接探针扩增技术(MultiplexLigation-dependentProbeAmplification，MLPA)，高通量测序。其中，比较基因组杂交分辨率比较低，Mb级，通量低，成本高；荧光定量PCR同样通量低，成本高，一次只能测一个拷贝数变异；荧光原位杂交，只针对特定位置，分辨率低，探针杂交效率不稳定；多重连接探针扩增技术，操作复杂，通量低，成本高，覆盖度小，易造成PCR污染。

高通量测序技术(High-throughput sequencing)又称“下一代”测序技术("Next-generation"sequencing technology)，以能一次并行对几十万到几百万条DNA分子进行序列测定和一般读长较短等为标志。高通量测序技术是对传统测序一次革命性的改变，一次对几十万到几百万条DNA分子进行序列测定，因此在有些文献中称其为下一代测序技术(next generation sequencing)足见其划时代的改变，同时高通量测序使得对一个物种的转录组和基因组进行细致全貌的分析成为可能，所以又被称为深度测序(deepsequencing)。高通量测序技术可以检测样本全基因组范围内的CNV。按照单个样本产出的测序数据的多少，可以大致分为高深度测序和低深度测序。考虑到商业运营成本和CNV分辨率的需求，通常使用较低测序的深度进行CNV测序。低深度测序适合检测流产组织、植入前胚胎筛查细胞等。这些组织或细胞的CNV通常较大(1000000个碱基级别及以上)，不要求CNV的精确断点等。其基本原理是，通过对待测样本和无CNV样本(参照样本)的测序，比较，找出测序reads数显著升高或者降低的区域。这段区域就是发生CNV的区域。然而，对于流产物而言，没有十分理想的参照样本，并且如果每批检测时，都加入参照样本，会提高检测成本。同时，由于一些未知的原因(可能是样本来源，样本状态，以及测序建库的批次不同)，会造成不同样本出现不同的测序偏好性，这一点在19号染色体上，体现得最为明显。为了解决缺失配对样本和未知原因造成的偏好性的问题，本发明开发了一种基于未知CNV的样本，按照样本特征进行样本分类，并在各个分类中，建立模拟数据以提供CNV检测所需的模拟数据集的方法。

发明内容

为了解决现有技术中存在的技术问题，本发明的目的在于提供一种基于未知CNV样本，建立模拟数据以提供CNV检测所需的模拟数据集的方法以及根据此方法建立的模拟数据集。本发明的方法不仅大大降低了检测成本，而且解决了未知原因造成的测序偏好性的问题。

为了实现上述目的，本发明采用了如下技术方案：

根据本发明的一个方面，本发明提供了一种用于对未知CNV人类DNA样本进行分类的随机森林模型的建立方法，所述建立方法包括：

(1)收集至少400例未知CNV人类DNA样本作为待检测样本，至少一例参照样本，该参照样本符合以下要求：常染色体无异常，至少有一条X染色体和Y染色体，至多有2条X染色体，2条Y染色体；

(2)在染色体上划分窗口，统计19号染色体所有窗口内reads数量，提取19号染色体所有窗口及其对应的reads，对这些reads进行均一化，形成文件TAB19；

(3)使用聚类方法，将TAB19文件作为输入，依照结果，对样本进行分类，设置分类数为2、3、4、5、6、7；

(4)通过随机森林算法，将TAB19文件和步骤(2)获得的分类进行训练，建立随机森林模型。并使用该模型，对原来样本进行预测，分别统计错误率。随着分类数的上升，验证之后的错误率也随之上升。在分类4后，错误率发生显著上升。同时，对各个分类进行统计画图，其也符合过去的经验。

前面所述的参照样本包括XY、XXY、XYY或XXYY。

本发明中使用的用以根据样本特征，将样本聚类的方法，可采用任何免费或商业的聚类方法，如：kmeans聚类、层次聚类、mean-shift聚类、主成分分析(PrincipalComponent Analysis，PCA)。

在本发明的具体实施方案中，聚类方法采用的是kmeans。

进一步，上面所述步骤(2)的详细过程如下：

(I)提取样本DNA，进行高通量测序；

(II)将序列比对到参考基因组；

(III)比对之后，使用软件samtools对比对后的数据进行排序；

(IV)对步骤(III)处理后得到的文件进行滑动窗口reads数及GC含量统计；按照固定窗口长度，统计样本数据在各个窗口内的reads数量，以及在参考基因组中，相应的窗口中的GC含量百分比，将这一过程得到的文件，成为TAB文件；将19号染色体所有窗口及其对应的reads提取出来，做成一个新的文件TAB19。

进一步，步骤(IV)中选择的窗口长度，依据测序数据量不同、检测精度不同，窗口可选不同。据经验而言，10kb大小窗口对于测序数据量大约是10M reads，检测精度在50kb；20kb窗口对应的数据量大约是6M reads，对应的检测精度是100kb；100kb窗口对应数据量大约是3M reads，对应检测精度是1M。200kb窗口对应测序量是2.5Mreads，精度是2M。

对本发明的样本进行测序采用的是高通量测序平台。测序平台不受特别限制，第二代测序平台：包括但不限于华大基因的BGI-seq测序平台、Illumina公司的GA、GAII、GAIIx、HiSeq1000/2000/2500/3000/4000、XTen、XFive、NextSeq500/550、MiSeq，AppliedBiosystems的SOLiD，Roche的454FLX，ThermoFisherScientific(LifeTechnologies)的IonTorrent、IonPGM、IonProtonI/II；第三代单分子测序平台：包括但不限于HelicosBioSciences公司的HeliScope系统，PacificBioscience的SMRT系统，OxfordNanoporeTechnologies的GridION、MinION。测序类型可为单端(SingleEnd)测序或双端(PairedEnd)测序，测序长度可为30bp、40bp、50bp、100bp、300bp等大于30bp的任意长度，测序深度可为基因组的0.01、0.02、0.1、1、5、10、30倍等大于0.01的任意倍数。

在本发明的具体实施方案中，采用的是BioelectronSeq 4000测序仪进行的低深度全基因组测序。

本发明的可采用的参考基因组可为全基因组、任意染色体、染色体的一部分。参考基因组通常选择已被公认确定的序列，如人的基因组可为NCBI或UCSC的hg18(GRCh18)、hg19(GRCh19)、hg38(GRCh38)，或任意一条染色体及染色体的一部分。

在本发明的具体实施方案中，参考基因组是hg19(GRCh19)。

比对软件可用任何一种免费或商业软件，如BWA(Burrows-WheelerAlignmenttool)、SOAPaligner/soap2(ShortOligonucleotideAnalysisPackage)、Bowtie/Bowtie2。将序列比对到参考基因组，得到序列在基因组上的位置。

在本发明的具体实施方案中，比对软件采用的是BWA。

用以比较待测样本和参照样本，以检测、定位CNV可用任何一种免费或商业软件及算法，如软件CNV-seq、CNAseq、ReadDepth、CNVnator、cnvHMM、seqCNA、XHMM、nbCNV，如算法隐马尔可夫模型(Hidden Markov Model,HMM)、循环二元分割算法(circular binarysegmentation,CBS)、贝叶斯模型(Bayesian model)。

在本发明的具体实施方案中，比较待测样本和参照样本，用以检测、定位CNV采用的是CNAnorm软件包。该软件包使用的算法是循环二元分割算法(CBS)。

根据已分类样本进行训练，将待测样本进行分类的算法，可采用任何一种算法，如支持向量机回归(support vector regression，SVM)、逻辑回归(logistic regression)、线性回归(linear regression)、决策树(decision tree)、随机森林(random forest)。

在本发明的具体实施方案中，根据已分类样本进行训练，将待测样本进行分类的算法，采用的是随机森林。

根据本发明的另一个方面，本发明提供了一种用于对未知CNV人类DNA样本进行分类的随机森林模型，所述随机森林模型是根据前面所述的建立方法建立的。

根据本发明的又一个方面，本发明提供了一种将未知CNV人类DNA样本进行分类的方法，所述方法包括：

(1)统计样本中19号染色体的所有reads数；

(2)根据步骤(1)获得的reads数和前面所述的随机森林模型，对样本进行分类。

进一步，所述步骤(2)的详细过程如下：

(I)提取样本DNA，进行高通量测序；

(II)将序列比对到参考基因组；

(III)比对之后，使用软件samtools对比对后的数据进行排序；

(IV)对步骤(III)处理后得到的文件进行滑动窗口reads数及GC含量统计；按照固定窗口长度，统计样本数据在各个窗口内的reads数量，以及在参考基因组中，相应的窗口中的GC含量百分比，将这一过程得到的文件，成为TAB文件；将19号染色体所有窗口及其对应的reads提取出来，进行均一化后，做成一个新的文件TAB19。

上面所述的步骤(IV)中滑动窗口的长度至少10kb，最大无上限(不超过单条染色体长度)。在现有数据量及精度要求情况下，选择100kb。

对本发明的样本进行测序采用的是高通量测序平台。测序平台不受特别限制，第二代测序平台：包括但不限于华大基因的BGI-seq测序平台、Illumina公司的GA、GAII、GAIIx、HiSeq1000/2000/2500/3000/4000、XTen、XFive、NextSeq500/550、MiSeq，AppliedBiosystems的SOLiD，Roche的454FLX，ThermoFisherScientific(LifeTechnologies)的IonTorrent、IonPGM、IonProtonI/II；第三代单分子测序平台：包括但不限于HelicosBioSciences公司的HeliScope系统，PacificBioscience的SMRT系统，OxfordNanoporeTechnologies的GridION、MinION。测序类型可为单端(SingleEnd)测序或双端(PairedEnd)测序，测序长度可为30bp、40bp、50bp、100bp、300bp等大于30bp的任意长度，测序深度可为基因组的0.01、0.02、0.1、1、5、10、30倍等大于0.01的任意倍数。

在本发明的具体实施方案中，采用的是BioelectronSeq 4000测序仪进行的低深度全基因组测序。

本发明的可采用的参考基因组可为全基因组、任意染色体、染色体的一部分。参考基因组通常选择已被公认确定的序列，如人的基因组可为NCBI或UCSC的hg18(GRCh18)、hg19(GRCh19)、hg38(GRCh38)，或任意一条染色体及染色体的一部分。

在本发明的具体实施方案中，参考基因组是hg19(GRCh19)。

比对软件可用任何一种免费或商业软件，如BWA(Burrows-WheelerAlignmenttool)、SOAPaligner/soap2(ShortOligonucleotideAnalysisPackage)、Bowtie/Bowtie2。将序列比对到参考基因组，得到序列在基因组上的位置。

在本发明的具体实施方案中，比对软件采用的是BWA。

根据本发明的又一个方面，本发明提供了一种用于模拟出染色体窗口reads数的预测模型的建立方法，所述建立方法包括：

(1)收集至少400例未知CNV人类DNA样本作为待检测样本、至少一例参照样本，所述参照样本符合以下要求：常染色体无异常，至少含有1条X染色体和1条Y染色体，至多含有2条X染色体，2条Y染色体；

(2)使用R语言CNAnorm包，逐个将待检测样本和参照样本组成待检测-参照样本对，进行CNV分析，得到每个待检测样本的CNV区域以及性别；

(3)按照CNV的不同类型，前面所述的TAB文件进行如下处理：除了Y染色体外，对于有拷贝数为1的所有窗口，reads数变更为原有reads数的2倍；对于3倍体的窗口，reads数变更为原有reads数的2/3；对于正常二倍体，reads数不做更改；对于Y染色体，将女性Y染色体的所有窗口定为NA，不参与后续计算；对于Y染色体为1倍体的片段，不更改；对于Y染色体2倍体的片段，reads数变为实际的1/2；对于其他异常情况，将对应窗口的reads数变为NA，不参与后续计算；处理后的TAB文件中所有的reads数求和，定义为处理后的该样本的reads总数；

(4)根据前面所述的分类方法得到的分类结果，分别构建每个分类中，所有样本变更后的reads总数和所有样本每个窗口的reads数的线性关系，建立线性模型，并给出检验的p值；

(5)对于p值达到阈值的窗口，使用步骤(4)的线性模型，作为最终的预测模型；对于p值达不到阈值的窗口，使用窗口的中位数作为预测模型；

进一步，所述异常情况包括嵌合、超过3倍体的多倍拷贝，以及0倍体。

优选地，所述阈值是10的-5次方。

根据本发明的又一个方面，本发明提供了根据前面所述的建立方法建立的线性模型或预测模型。

根据本发明的又一个方面，本发明提供了一种未知CNV人类DNA样本的模拟数据集的建立方法，其特征在于，所述建立方法包括：

(1)按照前面所述的方法对未知CNV的新样本进行分类；

(2)针对于不同的分类，使用前面所述的预测模型，生成模拟数据集。

根据本发明的又一个方面，本发明提供了根据前面所述的建立方法建立的模拟数据集。

根据本发明的又一个方面，本发明提供了一种分析人类DNA样本中CNV的方法，所述方法包括：

(1)按照前面所述的建立方法生成模拟数据集；

(2)将待测样本和模拟数据集组成待测样本-模拟样本对，使用CNAnorm进行CNV分析。

本发明对样本的类型不受特别限制，可以是含有大量核酸的样本，如组织、血液、尿液、唾液、羊水，也可以是含有微量核酸的样本，如血液、尿液、唾液中游离的单细胞、游离的核酸、生殖细胞。

对于含有微量核酸的样本，需要首先对单细胞扩增，以获得更多的核酸用于后续测序分析。单细胞扩增的方法不受特别限制，包括但不限于扩增前引物延伸PCR(PrimerextensionpreamplificationPCR，PEP-PCR)、退变寡核苷酸引物PCR(Degenerateoligonucleotideprimer-PCR，DOP-PCR)、多重置换扩增技术(MultipleDisplacementAmplification，MDA)、多次退火环状循环扩增技术(MultipleAnnealingandLoopingBasedAmplificationCycles，MALBAC)。

“Kmeans”算法是硬聚类算法，是典型的基于原型的目标函数聚类方法的代表，它是数据点到原型的某种距离作为优化的目标函数，利用函数求极值的方法得到迭代运算的调整规则。K-means算法以欧式距离作为相似度测度，它是求对应某一初始聚类中心向量V最优分类，使得评价指标J最小。算法采用误差平方和准则函数作为聚类准则函数。

Kmeans算法工作原理：

输入：聚类个数k，以及包含n个数据对象的数据库。

输出：满足方差最小标准的k个聚类。

“随机森林算法”是一种比较新的机器学习模型。经典的机器学习模型是神经网络，有半个多世纪的历史了。神经网络预测精确，但是计算量很大。上世纪八十年代Breiman等人发明分类树的算法(Breiman et al.1984)，通过反复二分数据进行分类或回归，计算量大大降低。2001年Breiman把分类树组合成随机森林(Breiman 2001a)，即在变量(列)的使用和数据(行)的使用上进行随机化，生成很多分类树，再汇总分类树的结果。随机森林在运算量没有显著提高的前提下提高了预测精度。随机森林对多元公线性不敏感，结果对缺失数据和非平衡的数据比较稳健，可以很好地预测多达几千个解释变量的作用(Breiman2001b)，被誉为当前最好的算法之一(Iverson et al.2008)。

随机森林算法工作原理：

随机森林顾名思义，是用随机的方式建立一个森林，森林里面有很多的决策树组成，随机森林的每一棵决策树之间是没有关联的。在得到森林之后，当有一个新的输入样本进入的时候，就让森林中的每一棵决策树分别进行一下判断，看看这个样本应该属于哪一类(对于分类算法)，然后看看哪一类被选择最多，就预测这个样本为那一类。

“CNAnorm”是一个R包，用于在考虑细胞含量和可能的非整倍体基因组大小的前提下，实现下一代序列数据的拷贝数分析。

“reads”指的是测序仪单次测序所得到的碱基序列，也就是一连串的ATCGGGTA……之类的。不同的测序仪器，reads长度不一样。

“samtools”是一个用于操作sam和bam文件(通常是短序列比对工具如bwa，bowtie2，hisat2，tophat2等等产生的)的工具合集，包含有许多命令。

“bwa”是用于将低分叉序列比对到大的参考基因组比如人基因组的软件包。BWA主要是由三种算法组成：BWA-backtrack，BWA-SW和BWA-MEM。第一个算法是针对于illumina测序reads最多100bp的算法。后面两个主要是针对于从70bp到1Mbp的更长序列。BWA-MEM和BWA-SW拥有一些相同的特征例如长reads支持和序列分开模式。但是相对而言，更加推荐更新的BWA-MEM，因为可以更快和更准确在更高质量上比对。BWA-MEM相比较BWA-backtrack在70-100bp illumina reads上有更好的性能。

本发明的优点和有益效果：

本发明通过构建算法，模拟生成了正常样本的序列数据，减少了样本测序的成本，排除了样本特征不同带来的误差，提高了检测的准确度。

附图说明

图1显示错误分类样本数(Y轴)随着分类数(X轴)的变化趋势图；

图2显示染色体19各个窗口(X轴)和经归一化后窗口的reads数(Y轴)的关系图。

具体实施方式

以下结合实施例来进一步说明本发明，本发明的优点和特点将会随着描述而更为清楚。但这些实施例仅是范例性的，并不对本发明的范围构成任何限制。本领域技术人员应该理解的是，在不偏离本发明的精神和范围下可以对本发明技术方案的细节和形式进行修改或替换，但这些修改和替换均落入本发明的保护范围内。下列实施例中未注明具体条件的实验方法，通常按照常规条件，或按照制造厂商所建议的条件。比例和百分比基于重量，除非特别说明。

实施例1模拟数据集的建立

1、样本收集

收集497例流产组织样本，1例不因为染色体异常而流产的男性流产组织样本，实验对象均签署了知情同意书。

2、DNA提取

DNA提取使用的试剂盒为：购自天根公司的血液组织细胞基因组提取试剂盒(DP304)

填写《组织基因组DNA提取操作记录单(天根)》。

2.1实验准备

以下所有实验步骤均在常温下(15–25℃)进行。

室温融化组织样本。

开启恒温振荡器，调至56℃。

Buffer GD及PW使用前，请注意是否按照配置要求加入无水乙醇并颠倒混匀；未开封的Buffer GD、PW有效期为室温保存一年。

观察缓冲液GB中是否有白色絮状沉淀；若有沉淀，可在37℃水浴锅中重新溶解。

2.2实验操作

(1)按照《组织基因组DNA提取操作记录单(天根)》上的样本编号在纯化柱盖和1.5mL离心管盖上写上对应的DNA编号；带编号的1.5mL离心管需要准备2套，带编号的2mL离心管需要准备1套。

(2)打开超净工作台紫外灯照射剪刀、镊子、培养皿等15min后，通风10min。

(3)处理组织

将剪刀镊子浸泡到配置好的生理盐水中待用；

待组织融化后，用清水冲洗，再用生理盐水冲洗至无明显血色；取组织约10mg剪碎后放入对应标记的第一套1.5mL离心管中。并将剩余样本放入2mL离心管中用于留样，注意核对样本编号；并更换样本时更换剪刀和镊子；

瞬时离心，用移液器弃去液体，并用蒸馏水洗两遍，至无明显血色；

10,000rpm(～11,200×g)离心1min，倒尽上清，加200μl缓冲液GA，振荡至彻底悬浮。加入20μl Proteinase K溶液，涡旋1min混匀，瞬离；

在56℃放置，900rpm振荡过夜(此步至少1h，如样本未裂解完全可适当延长裂解时间)，瞬时离心以去除管盖内壁的水珠。

(4)裂解样品

加入200μL缓冲液GB，充分颠倒混匀，瞬时离心后，置于恒温振荡仪上70℃上10min，溶液应变清亮(如溶液未彻底变清亮，请延长裂解时间至溶液清亮为止)；

瞬时离心后，加200μL无水乙醇，充分颠倒混匀，瞬离，此时可能会出现絮状沉淀。

(5)结合DNA

将上述离心管中所得溶液和絮状沉淀全部加入吸附柱中(吸附柱放入收集管中)；

3,400×g离心30秒，将收集管中的废液倒入原1.5mL离心管，暂存。

(6)清洗硅胶膜

向吸附柱中加入500μL缓冲液GD，13,400×g离心30秒，倒掉收集管中的废液，用吸水纸蘸去管口残余废液，将吸附柱放回收集管中；

向吸附柱中加入700μL漂洗液PW，13,400×g离心30秒，倒掉收集管中的废液，用吸水纸蘸去管口残余废液，将吸附柱放回收集管中；

向吸附柱中加入500μL漂洗液PW，13,400×g离心30秒，倒掉收集管中的废液，用吸水纸蘸去管口残余废液；

向吸附柱中加入200μL漂洗液PW，13,400×g离心30秒，倒掉收集管中的废液，用吸水纸蘸去管口残余废液。

(7)干燥硅胶膜

将吸附柱放回收集管中，13,400×g离心2分钟，倒掉废液，用吸水纸蘸去管口残余废液；

将吸附柱放回收集管中，13,400×g离心2分钟，倒掉废液，用吸水纸蘸去管口残余废液；

将吸附柱置于一个新1.5mL收集管中，开盖放置5分钟。

(8)洗脱DNA

向吸附膜中间位置悬空滴加105μL洗脱缓冲液TE，室温放置3分钟，13,400×g离心2分钟，将溶液收集到离心管中。

(9)DNA定量

取2μL DNA样品进行NanoQ定量分析，DNA的A260/280在1.80～2.0，浓度应在30ng/μL以上。

3、打断

3.1准备

将提取好的DNA用Nanodrop检测，A260/280在1.8～2.0，浓度>30ng/ul，体积约100uL；

用TE将样本基因组稀释到40ng/ul，125uL，瞬时离心，备用；

将稀释好的基因组DNA分别转移到标好对应编号的打断管中，并注意核对，打断前样本尽量放置在冰上。

3.2操作步骤

(1)接通电源，按照电脑-仪器-软件的顺序依次打开基因打断仪(Covaris,M220)。

(2)装配打断仪：打开舱门，将蓝色打断架正确安放到打断仪的水槽上，注意缺口在右下角；并将规格为130uL的试管托安放到蓝色底座上，放下可旋转金属块，关好舱门。

(3)加水：打开舱门从右侧进水口加入去离子水，直至电脑软件页面上waterlevel项显示绿色对勾为止；关好舱门，仪器会对去离子水进行预冷，直至电脑界面上watertemperature右侧显示绿色对勾，方可使用。

(4)样本准备：取出相应数量的打断管并在管盖上写好对应编号。将稀释好的基因组DNA 120uL转移到对应编号的打断管中，注意核对样本编号。

(5)打断

将盛有DNA的打断管放在试管托上，放下可旋转金属块，关好舱门，Check左上角三个样本温度、水位线及舱门三个对呈对勾状，选择200s程序，核对参数，参数如表1所示。

表1参数

参数	Peak Incident Power	Duty Factor	Cycles per Burs	打断持续时间(s)
						50	20	200	200

参数无误后，将点击Run，即开始打断。

打断结束后跳出对话框，关闭对话框，开启舱门，拿出打断好的样本，放入待打断样本，关闭舱门，重复上面步骤继续打断。

4、文库构建

4.1实验前准备

将4℃保存的磁珠AgencourtAMPure XP磁珠取出，室温平衡30min，充分震荡混匀备用。

将打断后的DNA从-20℃冰箱取出，将其置于离心管架上室温溶解，振荡混匀，短暂离心备用。

取出末端修复缓冲液置于冰上融化。

4.2操作步骤

(1)末端修复

取出1.5mL的EP管，并做好标记，按照表2配置体系，手持震荡并瞬离。

表2末端修复反应体系

组分	加入体积
		NF水	19.5μL
5X末端修复缓冲液	10μL
		末端修复酶	0.5μL
总量	30μL

将配置好的体系分装到标记好的0.2mL离心管中，根据样本编号加入对应DNA 20μL，手持震荡混合，并瞬离，置于PCR仪上反应，25℃，20min。

(2)纯化末端修复后的文库

将90μL的AgencourtAMPure XP磁珠(1.8X样本体积)分装到写好对应编号的EP管中；

将下机后的DNA全部转入上述离心管，震荡混合，瞬离(转速不得超过3000rpm)，静置5min；

上架5min，直到溶液清澈，弃掉上清，注意不要弃掉磁珠；

保持EP管在磁力架上，向EP管中加入500μL新鲜配置的70％乙醇，在磁力架上转动EP管两次，等溶液澄清后，弃掉上清，注意不要弃掉磁。该步骤重复一次；

从磁力架上取下EP管，瞬离，用20μL的枪将EP管中剩余的液体移除，不要碰到磁珠；

保持EP管在磁力架上，室温晾干磁珠≤5min。待磁珠出现裂痕后(以1-2道裂痕为宜)向管中加入14μL的TE，取下EP管，震荡混匀，瞬时离心，(转速不得超过3000rpm)，静置5min；

上架，待到溶液清澈后，可将其转入下步MIX中。

(3)加接头和缺口修复

取出1.5mL的EP管，并做好标记，按照表3配置体系，手持震荡并瞬离。

表3连接接头缺口修复反应体系

组分	加入量
		NF水	24.5
10X Ligase Buffer	5μL
		dNTP Mix	1μL
Ion P1 Adapter	1μL
		DNA Ligase	1μL
Nick Repair Polymerase	4μL
		Total	36.5μL

将配置好的体系分装到标记好的0.2mL离心管中，根据记录单加入对应Barcode，并将上步所得纯化后DNA 12.5μL转移到对应编号的离心管中，手持震荡混合，并瞬离，置于PCR仪上反应，25℃，15min；72℃，5min。

(4)纯化连接接头和缺口修复后的文库

准备2套EP管分别分装40μL和7.5μL的AgencourtAMPure XP磁珠分装到写好对应编号的EP管中；

将下机后的DNA全部转入装有40μL AgencourtAMPure XP磁珠的离心管中，震荡混合，瞬离(转速不得超过3000rpm)，静置5min；

上架5min，直到溶液清澈，将上清液转到装有7.5μLAgencourtAMPure XP磁珠的离心管中(一次转不完全可将EP管离心后再用20μL移液器再转一遍)，震荡混合，瞬离(转速不得超过3000rpm)，静置5min；

上架5min，直到溶液清澈，弃掉上清，注意不要弃掉磁珠；

保持EP管在磁力架上，向EP管中加入500μL新鲜配置的70％乙醇，在磁力架上转动EP管两次，等溶液澄清后，弃掉上清，注意不要弃掉磁珠。该步骤重复一次；

从磁力架上取下EP管，瞬离，用20μL的枪将EP管中剩余的液体移除，不要碰到磁珠；

保持EP管在磁力架上，室温晾干磁珠≤5min。待磁珠出现裂痕后(以1-2道裂痕为宜)向管中加入14μL的TE，取下EP管，震荡混匀，瞬时离心，(转速不得超过3000rpm)，静置5min；

上架，待到溶液清澈后，可将其转入下步MIX中。

(5)扩增文库

取出1.5mL的EP管，并做好标记，按照表4配置体系，手持震荡并瞬离。

表4文库扩增反应体系

将配置好的体系分装到标记好的0.2mL离心管中，并将上步所得纯化后DNA12.5μL转移到对应编号的离心管中，手持震荡混合，并瞬离，置于PCR仪上反应，95℃，5min；(95℃15s,58℃15s,70℃，1min)x10；4℃Hold。

(6)纯化扩增后文库

将78μL的AgencourtAMPure XP磁珠分装到写好对应编号的EP管中；

将下机后的DNA全部转入上述离心管，震荡混合，瞬离(转速不得超过3000rpm)，静置5min；

上架5min，直到溶液清澈，弃掉上清，注意不要弃掉磁珠；

保持EP管在磁力架上，向EP管中加入500μL新鲜配置的70％乙醇，在磁力架上转动EP管两次，等溶液澄清后，弃掉上清，注意不要弃掉磁珠。该步骤重复一次；

从磁力架上取下EP管，瞬离，用20μL的枪将EP管中剩余的液体移除，不要碰到磁珠；

保持EP管在磁力架上，室温晾干磁珠≤5min。待磁珠出现裂痕后(以1-2道裂痕为宜)向管中加入20μL的TE，取下EP管，震荡混匀，瞬时离心，静置5min；

上架，待到溶液清澈后，将其收集到事先准备好的文库收集管中。

(7)文库收集及储存

将澄清后的文库收集到事先准备好的文库收集管中，注意反复核对样本编号。

使用qPCR对文库进行定量。文库浓度达到≥1nmol/L时，可以上机测序。若文库浓度<1nmol/L时，需要重新构建文库。

5、上机测序

上机测序采用博奥生物集团有限公司基因测序仪器BioelectronSeq 4000。

5.1仪器清洗

开启氮气瓶，调节压力为30psi。开启测序仪，进入主页面单击“Clean”，根据测序仪的使用情况选择水洗或者氯洗程序：机器停用72小时以内，使用水洗；机器停用超过72小时，使用氯洗。

水洗：清空C1和C2清洗瓶，用18MΩ纯水洗两次。向C1和C2清洗瓶个加入100mL 18MΩ纯水，安装在C1和C2位置上。在测序仪屏幕上选择“Clean”选项，装上芯片，根据屏幕提示单击“Next”直至水洗程序开始运行，清洗时间约15分钟。

氯洗：清空两个氯洗瓶，用18MΩ纯水洗两次氯洗瓶和玻璃氯瓶。向装有1L 18MΩ纯水的玻璃氯瓶中加入氯片，溶解10分钟，加入1ml 1M NaOH，颠倒混匀。用0.45μm氯头分别过滤100mL洗液到氯洗瓶中，将两个氯洗瓶安装在C1和C2位置上，在测序仪屏幕上选择“Clean”选项，装上芯片，根据屏幕提示单击“Next”直至氯洗程序开始运行，清洗时间约为15分钟。氯洗完成后，进行一次水洗。

5.2仪器初始化

清空W2试剂瓶，用18MΩ纯水洗三次。加入1920mL的18MΩ纯水，并加入125mL测序溶液II，盖好瓶盖，颠倒混匀4～8次。取出两个新的140mL试剂管，分别标记为W1瓶和W3瓶。向W1瓶中加入32μL 1M NaOH，往W3瓶中加入40～50mL测序溶液III，盖好瓶盖。

在测序仪屏幕上选择“Initialization”和“Ion P1 Hi-Q200 V3 Kit”，更换W1、W2、W3位置的吸管，将试剂瓶安装到对应的位置，拧紧接口盖。装好用于初始化的芯片，根据屏幕提示确定仪器状态，点击“Next”直至程序开始运行，第一阶段约需30分钟。

预先取出dGTP、dCTP、dATP、dTTP，将其置于冰盒上，融化后振荡10秒，离心2秒，取出4个新的140mL试剂管，分别做好“G、C、A、T”标记，并加入70μL对应的dNTP溶液。

第一阶段初始化完成后，更换dNTP位置的吸管，将试剂管安装到对应的位置并拧紧。根据屏幕提示，继续完成初始化，约需30分钟。初始化完成后回到主界面。

5.3设置Plan

登陆测序仪服务器，单击选择“Plan”，在“Plan”下拉菜单栏单击选择“Templates”。

在对应的应用下选择“Plan new run”，或在存在的模板中选择“Plan Run”。

按照提示在“Application”中选择Whole Genome，单击“Next(Kit)”。

按照下表5选择Kits信息，其余留空，然后点击“Next(Monitor)”。

表5 Kits信息

选项	选项类型
		Library Kit Type	Ion Plus Fragment Library Kit
Templating Kit	Ion PI Hi-QOT2 200 Kit
		Sequence Kit	IonPI Hi-Q 200 Sequence Kit
Chip Type	Ion PI Chip V3
		Barcode Set	Ionxpress
Flow	400

Monitor测序实时监控参数为默认值，不做修改。然后点击“Next(Reference)”。

Plugins插件选项留空不选，单击“Next(Project)”。

在Project选择中选择样品所对应的项目，单击“Next(Export)”。

Export报告产生形式选项为默认值，不做修改，单击“Next(Plan)”。

输入Run Plan Name和样本的标签序列使用情况。Run Name命名规则：P30(机器编号)-文库上机浓度-混合文库编号。

标签序列使用情况：在对应的标签序列(Barcode)编号后面输入对应的样本编号，输入后需核对标签序列号和样本编号是否属实。

Reference选项留空不选，单击“Next(Plugins)”。

单击“Plan Run”保存Plan设置。

5.4准备上机文库

预先取出质控微珠溶液、测序引物、测序聚合酶置于冰盒上。

在仪器初始化快结束时进行文库的上机准备。将质控微珠溶液振荡30秒混匀，瞬时离心；取5μL加入样本微珠溶液中，涡旋振荡30秒，15500g离心5分钟；枪头对着非沉淀处，小心吸走上清，在管中保留10μL溶液。

往上步样本微珠溶液加入15μL退火缓冲液，总体积是25μL。

测序引物在冰上融化后，涡旋振荡30秒，混匀离心2秒。往上步ISP样本溶液加入20μL测序引物，确保总体积45μL，将样本微珠溶液涡旋振荡60秒，混匀，瞬时离心。

在PCR仪上按以下程序对样本微珠溶液进行退火：95℃2分钟，37℃2分钟，20℃hold。

退火完毕后，往上步微珠溶液加入10μL上样缓冲液，涡旋振荡10秒，混匀，瞬时离心，室温条件放置。

5.4上样和测序

将芯片取出，置于吊篮中。吸取55μL上一步得到的样本微珠溶液，注入到上步芯片的进样槽中(非进样孔)。

将芯片放置在离心机上，芯片的缺口朝外，与旧芯片配平。离心10分钟。

离心过程中准备以下试剂：

在1.5mL离心管中，分别加入0.5mL的退火缓冲液和0.5mL纯化水混匀，配制成50％退火缓冲液。该溶液配置后，在7天内可以使用。

在1.5mL离心管中，分别加入0.5mL的100％异丙醇溶液和0.5mL退火缓冲液混匀，配制成50％冲洗液。该溶液当天使用当天配置。

在1.5mL离心管中，分别加入60μL 50％的退火缓冲液和60μL的测序聚合酶混匀，配制成酶反应液。该溶液配制后置于冰盒上。

在1.5mL离心管中，分别加入49μL 50％的退火缓冲液和1μL发泡液混匀，配制成发泡剂。

用100μL移液器打入100μL空气进入发泡剂中，快速反复吹打液体5秒钟，使大气泡变成小泡。重复该操作一次，使发泡剂成为密致的小泡。

离心完毕后取出芯片，匀速将100μL细小的泡泡注入加样孔中，并将出样槽溢出的液体吸走；往进样槽加入55μL 50％的退火缓冲液，将芯片放回离心机，离心30秒。重复该操作一次。

往芯片的金阳空垂直缓慢加入100μL冲洗液两次。每一次往芯片加入溶液后，在另一个孔吸走排出来的液体。

往芯片的进样孔垂直缓慢加入100μL 50％的退火缓冲液三次，避免产生气泡。每一次往芯片加入溶液后，在另一孔吸走排出来的液体。

往芯片的进样孔垂直缓慢加入65μL酶反应液，避免产生气泡。在另一孔吸走排出来的液体。

室温酶孵育5分钟后，将芯片安装到测序仪芯片槽上，运行程序，选择提前设置好的Plan，核对无误后，开始测序，运行时间为2.5小时。

程序运行结束后72小时内，返回主菜单，进行水洗。如超过72小时，则对仪器进行氯洗再进行水洗。水洗完毕，在主菜单上关闭仪器，在放dNTP孔上套上管子防尘。如长时间不使用，关闭氮气瓶气阀。

6、测序后原始数据处理

(1)测序的原始下机reads数在3.2M～6M之间不等。使用软件bwa与人类基因组参考序列19版本(hg19)进行比对。比对之后，使用软件samtools对比对后的文件进行排序。

(2)对上述处理后得到的文件进行滑动窗口reads数及GC含量统计。按照固定窗口长度，统计样本数据在各个窗口内的reads数量，以及在参考基因组中，相应的窗口中的GC含量百分比，将这一过程得到的文件，成为TAB文件。将19号染色体所有窗口及其对应的reads提取出来，将所有窗口进行均一化。均一化方法为，该窗口均一化后的数据＝(该窗口原始reads数-19号染色体所有窗口reads的平均数)/(19号染色体所有窗口reads数的标准差)做成一个新的文件TAB19。

(3)使用R语言CNAnorm包，逐个将497例流产样本和1例样本组成497对待检测-参照样本对，进行CNV分析，得到每个流产组织样本的CNV区域以及性别。

(4)按照CNV的不同类型，对第三步得到的TAB文件进行处理。除了Y染色体外，对于有拷贝数为1的所有窗口，reads数变更为原有reads数的2倍。对于3倍体的窗口，reads数变更为原有reads数的2/3。对于其他异常情况，包括嵌合、超过3倍体的多倍拷贝，将对应窗口的reads数变为NA，不参与后续计算。对于正常二倍体，reads数不做更改。对于Y染色体，将女性Y染色体的所有窗口定为NA，不参与后续计算。对于Y染色体为1倍体的片段，不更改。对于Y染色体2倍体的片段，reads数变为实际的1/2。其他异常情况，都将reads数变为NA，不参与后续计算。处理后的TAB文件中所有的reads数求和，定义为处理后的该样本的reads总数。

(5)使用R语言中kmeans函数的聚类方法，将TAB19文件作为输入矩阵，依照结果，对样本进行分类，最佳分类数筛选：首先，设置分类数为2、3、4、5、6、7；通过随机森林算法，将TAB19文件和上面所述的分类进行训练，建立随机森林模型。并使用该模型，对原来样本进行预测，分别统计错误率。图1结果显示，随着分类数的上升，验证之后的错误率也随之上升。在分类4后，错误率发生显著上升。同时，对各个分类进行统计画图(如图2所示，图中4种颜色代表分类1-4，每一条线代表一个样本。每一张小图是所有样本叠加后形成的。4种分类的情况，符合经验上的染色体19分类情况，即1类对应中间稍高，两侧稍低。2类对应中间极高，两侧极低，且存在较大波动。3类对应中间较低，两侧较高。4类对应两侧与中间基本一致的情况)，其也符合过去的经验。因此认为设置分类4是最佳的。之后，人为设定聚为4类，其他参数使用默认参数，依照结果，将样本文件进行分类。

(6)使用R语言randomForest软件包中randomForest方法，将(5)中得到的分类作为训练集，设置参数importance＝T，proximity＝T,ntree＝15000，其他参数为默认，将TAB19文件和上一步得到的分类进行训练，得到随机森林模型。

(7)按照步骤5得到的分类结果，分别构建每个分类中，所有样本变更后的reads总数和所有样本每个窗口的reads数的线性关系，并给出检验的p值。

(8)对于p值达到阈值的窗口，使用上一步骤的线性模型，作为最终的预测模型。对于p值达不到阈值的窗口，使用窗口的中位数作为预测模型。

(9)对于一个待测样本，首先使用步骤(1)和(2)所述方法，进行基本的比对、划分窗口、计算GC含量。而后，19号染色体的所有reads数以及步骤6构建的预测模型，对样本进行分类。而后，针对于不同的分类，使用步骤8建立的预测模型，生成最终的模拟数据集。将待测样本和模拟数据集组成待测样本-模拟样本对，使用CNAnorm进行CNV分析。

实施例2方法的适用性和准确性的验证

1、样本收集

使用了不参与实施例1的预测模型构建的400例流产组织样本；10例有异常生育史的病患的血液样本以及8例国际上通用的标准品DNA样本(来源：coriell公司https://www.coriell.org/生产的，编号分别为NA12662，NA12074，NA20305，NA14943，NA13277，NA06230，NA22569，NA01723)进行验证。

2、检测

400例流产样本中的200例利用本发明的方法进行参照样本数据生成，并用软件CNV-seq进行CNV检测；

另外200例使用不需要参考数据集的CNV检测方法。该方法为Z-score方法。首先，对这200例样本进行窗口划分，reads数及GC含量统计，方法与前述一致。而后，通过均一化，即(该窗口reads数减去该样本所有窗口reads数)除以(该样本所有窗口的标准差)，用以消除不同样本测序量不同而对样本产生的影响。得到所有样本所有窗口均一化数值后，对所有样本同一窗口进行Z检验，以Z值小于-3或者大于3为阈值，确定异常窗口。异常窗口即为包含CNV的窗口；

10例异常生育史的病患样本，使用核型分析方法确定真实的CNV信息。

8例标准品采用核型分析、芯片检测、fish验证等一种或多种方法确定标准CNV。

结果：与本发明的方法相比，对于5M以上的异常区域，400例样本的一致性达到99％以上；10例异常生育史的核型结果，染色体异常区带的一致性达到100％；8例标准品样本，一致性达到100％。

三种来源不同，检测方法不同的样本，其检出的CNV和利用本发明方法检出的CNV，一致性极高，证明了本发明方法的普适性和准确性。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (3)

1.一种用于模拟出染色体窗口reads数的预测模型的建立方法，其特征在于，所述建立方法包括：

（1）收集至少400例未知CNV人类DNA样本作为待检测样本、至少一例参照样本，该参照样本符合以下要求：常染色体无异常，至少有一条X染色体和一条Y染色体，至多有2条X染色体，2条Y染色体；

（2）提取待检测样本和参照样本的DNA，进行高通量测序；将序列与参考基因组进行比对；比对之后，使用软件samtools对比对后的数据进行排序；对排序后得到的文件进行滑动窗口reads数及GC含量统计；按照固定窗口长度，统计样本数据在各个窗口内的reads数量，以及在参考基因组中，相应的窗口中的GC含量百分比，将这一过程得到的文件，称为TAB文件；将19号染色体所有窗口及其对应的reads提取出来，做成一个新的文件TAB19；使用聚类方法，人为设定聚为4类，将TAB19文件作为输入，依照结果，对待检测样本和参照样本进行分类；

（3）通过随机森林算法，将TAB19文件和步骤（2）获得的分类进行训练，建立随机森林模型；

（4）使用R语言CNAnorm包，逐个将待检测样本和参照样本组成待检测-参照样本对，进行CNV分析，得到每个待检测样本的CNV区域以及性别；

（5）按照CNV的不同类型，对步骤（2）得到的TAB文件进行如下处理：除了Y染色体外，对于有拷贝数为1的所有窗口，reads数变更为原有reads数的2倍；对于3倍体的窗口，reads数变更为原有reads数的2/3；对于正常二倍体，reads数不做更改；对于Y染色体，将女性Y染色体的所有窗口定为NA，不参与后续计算；对于Y染色体为1倍体的片段，不更改；对于Y染色体2倍体的片段，reads数变为实际的1/2；对于其他异常情况，将对应窗口的reads数变为NA，不参与后续计算；处理后的TAB文件中所有的reads数求和，定义为处理后的该样本的reads总数；所述异常情况包括嵌合、超过3倍体的多倍拷贝以及0倍体；

（6）根据步骤（2）得到的分类结果，分别构建每个分类中，所有样本变更后的reads总数和所有样本每个窗口的reads数的线性关系，建立线性模型，并给出检验的p值；

（7）对于p值达到阈值的窗口，使用步骤（6）的线性模型，作为最终的预测模型；对于p值达不到阈值的窗口，使用窗口的中位数作为预测模型；所述阈值是10的-5次方。

2.一种未知CNV人类DNA样本的模拟数据集的建立方法，其特征在于，所述建立方法包括：

（1）统计待检测样本和参照样本中的19号染色体的所有reads数；

（2）根据步骤（1）获得的reads数和权利要求1建立的随机森林模型，对待检测样本和参照样本进行分类；

（3）针对于不同的分类，使用权利要求1建立的预测模型，生成模拟数据集。

3.一种分析人类DNA样本中CNV的方法，其特征在于，所述方法包括：

（1）按照权利要求2的建立方法生成模拟数据集；

（2）将待检测样本和模拟数据集组成待测样本-模拟样本对，使用CNAnorm进行CNV分析。

Images