CN108804876A - 用于计算癌症样本纯度和染色体倍性的方法和装置 - Google Patents
用于计算癌症样本纯度和染色体倍性的方法和装置 Download PDFInfo
- Publication number
- CN108804876A CN108804876A CN201710312237.7A CN201710312237A CN108804876A CN 108804876 A CN108804876 A CN 108804876A CN 201710312237 A CN201710312237 A CN 201710312237A CN 108804876 A CN108804876 A CN 108804876A
- Authority
- CN
- China
- Prior art keywords
- peak
- tre
- formula
- indicate
- window
- 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.)
- Granted
Links
- 206010028980 Neoplasm Diseases 0.000 title claims abstract description 124
- 201000011510 cancer Diseases 0.000 title claims abstract description 111
- 238000000034 method Methods 0.000 title claims abstract description 63
- 238000004364 calculation method Methods 0.000 claims abstract description 17
- 238000009826 distribution Methods 0.000 claims description 61
- 239000012634 fragment Substances 0.000 claims description 42
- 108090000623 proteins and genes Proteins 0.000 claims description 39
- 108700028369 Alleles Proteins 0.000 claims description 35
- 210000004027 cell Anatomy 0.000 claims description 27
- 238000012163 sequencing technique Methods 0.000 claims description 27
- 238000004458 analytical method Methods 0.000 claims description 10
- 238000012937 correction Methods 0.000 claims description 10
- 210000001519 tissue Anatomy 0.000 claims description 10
- 238000000605 extraction Methods 0.000 claims description 8
- 238000007476 Maximum Likelihood Methods 0.000 claims description 7
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 6
- 230000014509 gene expression Effects 0.000 claims description 6
- 238000012545 processing Methods 0.000 claims description 6
- 210000000349 chromosome Anatomy 0.000 claims description 5
- 239000000284 extract Substances 0.000 claims description 5
- 238000013467 fragmentation Methods 0.000 claims description 5
- 238000006062 fragmentation reaction Methods 0.000 claims description 5
- UYTPUPDQBNUYGX-UHFFFAOYSA-N guanine Chemical compound O=C1NC(N)=NC2=C1N=CN2 UYTPUPDQBNUYGX-UHFFFAOYSA-N 0.000 claims description 4
- 239000011159 matrix material Substances 0.000 claims description 4
- 238000002864 sequence alignment Methods 0.000 claims description 4
- 238000004043 dyeing Methods 0.000 claims description 3
- 238000001914 filtration Methods 0.000 claims description 3
- 230000001186 cumulative effect Effects 0.000 claims description 2
- 230000001419 dependent effect Effects 0.000 claims description 2
- 238000009499 grossing Methods 0.000 claims 3
- JUJWROOIHBZHMG-UHFFFAOYSA-N Pyridine Chemical compound C1=CC=NC=C1 JUJWROOIHBZHMG-UHFFFAOYSA-N 0.000 claims 2
- 239000000203 mixture Substances 0.000 claims 2
- 235000013399 edible fruits Nutrition 0.000 claims 1
- UMJSCPRVCHMLSP-UHFFFAOYSA-N pyridine Natural products COC1=CC=CN=C1 UMJSCPRVCHMLSP-UHFFFAOYSA-N 0.000 claims 1
- 238000011160 research Methods 0.000 description 8
- 230000004075 alteration Effects 0.000 description 4
- 201000010099 disease Diseases 0.000 description 4
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 238000012070 whole genome sequencing analysis Methods 0.000 description 4
- 238000004422 calculation algorithm Methods 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 238000012268 genome sequencing Methods 0.000 description 3
- 238000007481 next generation sequencing Methods 0.000 description 3
- 238000001712 DNA sequencing Methods 0.000 description 2
- 108700019961 Neoplasm Genes Proteins 0.000 description 2
- 102000048850 Neoplasm Genes Human genes 0.000 description 2
- 230000002776 aggregation Effects 0.000 description 2
- 238000004220 aggregation Methods 0.000 description 2
- 238000000205 computational method Methods 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 239000003814 drug Substances 0.000 description 2
- 238000010191 image analysis Methods 0.000 description 2
- 239000002773 nucleotide Substances 0.000 description 2
- 125000003729 nucleotide group Chemical group 0.000 description 2
- 238000011002 quantification Methods 0.000 description 2
- 230000011218 segmentation Effects 0.000 description 2
- 230000000392 somatic effect Effects 0.000 description 2
- 241000208340 Araliaceae Species 0.000 description 1
- 208000026310 Breast neoplasm Diseases 0.000 description 1
- 241001269238 Data Species 0.000 description 1
- 206010068052 Mosaicism Diseases 0.000 description 1
- 101001024425 Mus musculus Ig gamma-2A chain C region secreted form Proteins 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 229910002056 binary alloy Inorganic materials 0.000 description 1
- 230000004663 cell proliferation Effects 0.000 description 1
- 230000001413 cellular effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000002759 chromosomal effect Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 229940079593 drug Drugs 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 244000144992 flock Species 0.000 description 1
- 210000004602 germ cell Anatomy 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 230000003834 intracellular effect Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 230000000869 mutational effect Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000007170 pathology Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 210000003765 sex chromosome Anatomy 0.000 description 1
- 235000015170 shellfish Nutrition 0.000 description 1
- 238000010186 staining Methods 0.000 description 1
- 230000004797 therapeutic response Effects 0.000 description 1
- 210000004881 tumor cell Anatomy 0.000 description 1
- 239000002023 wood Substances 0.000 description 1
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
- G16B15/00—ICT specially adapted for analysing two-dimensional or three-dimensional molecular structures, e.g. structural or functional relations or structure alignment
-
- 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
-
- 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
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Engineering & Computer Science (AREA)
- General Health & Medical Sciences (AREA)
- Biophysics (AREA)
- Organic Chemistry (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Biotechnology (AREA)
- Theoretical Computer Science (AREA)
- Zoology (AREA)
- Evolutionary Biology (AREA)
- Analytical Chemistry (AREA)
- Bioinformatics & Computational Biology (AREA)
- Medical Informatics (AREA)
- Wood Science & Technology (AREA)
- Molecular Biology (AREA)
- Crystallography & Structural Chemistry (AREA)
- Microbiology (AREA)
- Immunology (AREA)
- Biochemistry (AREA)
- General Engineering & Computer Science (AREA)
- Genetics & Genomics (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供了一种全自动、高效率、高准确性的计算癌症样本纯度和染色体倍性的方法和装置。通过本发明提供的层次混合高斯模型，实现了对癌症样本纯度和染色体倍性的快速和准确计算，节约了纯度估算的时间和经济成本，同时提高了计算结果的准确性。本发明在癌症样本纯度和染色体倍性计算上具有广阔的运用前景。
Description
技术领域
本发明属于癌症研究领域，具体涉及用于计算癌症样本中癌症细胞纯度和细胞内染色体倍性的方法和装置。
背景技术
癌症的研究是生命医学中的重要研究领域，并对人类健康生活有重大影响。癌症是一类细胞恶性增殖的疾病，因其病理十分复杂，人类尚无法攻克这类疾病。二代测序(next generation sequencing)为快速检测病人遗传信息提供了可能。然而测序需要从病人组织中提取样本，但通常癌症组织并不是只单纯地包含癌症细胞，它还有非常丰富的微环境。癌症细胞微环境指包围或伴随癌症细胞的正常细胞(non-cancerous cells)环境。癌细胞样本提取时，这些微环境会和癌细胞一起被提取，并会伴随癌细胞一起被测序[1]。癌症细胞在癌症样本中的比例被定义为癌症样本的纯度。癌症基因组通常包含着大量体细胞序列拷贝数变异，这些变异主要由基因组片段扩增或删除造成。识别特定肿瘤基因组的基因组片段拷贝数变化，是癌症基因组研究的一个重要课题。要准确鉴定基因组片段拷贝数具有一定挑战，因为癌症片段拷贝数主要由两个因素混合决定，一是癌症样本纯度，即癌症细胞在癌症样本中所占比例，二是染色体倍性[2,3]。传统中鉴定癌症样本纯度和染色体倍性的方法是使用实验技术，如定量图像分析[4]或单细胞测序[5]。但是在大型项目中，这样的方法会耗费大量人力、资金和时间。随着测序技术地发展，测序数据地快速增长，以及测序数据分析技术的积累，各种各样的癌症样本纯度算法被提出，并开发出了相对应的软件。
基于基因组片段拷贝数变异和基于等位基因频率(突变位点的B-等位基因(B-allele))的计算方法被相继提出。基于等位基因频率的方法有PurityEst[6]和PurBayes[7]，主要是依赖于随着肿瘤样本纯度和肿瘤基因组倍性的不同，等位基因的频率会有所不同。基于拷贝数变异的方法有CNAnorm[8]、THetA[9]、和ABSOLUTE[10]等。然而这两种方法都有不同程度的问题，使用等位基因频率的方法由于数据量的问题会有较大的误差，而运用拷贝数变异的方法虽然较稳定，但却无法区分样本纯度和染色体倍性的补偿效应，即存在识别问题。以上基于片段拷贝数的软件都没有解决这一问题，CNAnorm倾向于选择染色体倍性离二倍体最近的解决方案，ABSOLUTE结合了其他的经验数据，THetA则直接将所有可能结果都列出来了。
更优的方案应该是结合等位基因频率信息和片段拷贝数信息共同计算肿瘤样本纯度。PyLOH[11]，patchwork[12]使用了基因组上杂合SNV(单核苷酸变异)位点的频率信息，和基因组片段的拷贝数。PyLOH一定程度上解决了“识别困难”的问题，可以更合理给出唯一解决方案。但是其准确性较差，特别是遇到基因组中存在亚克隆(subclone)的情况下。Patchwork同时使用了两种信息，但是在计算基因型的中间步骤中，需要人工识别，人工判断的结果缺乏准确性，并且这种半自动化软件给应用带来很多不便。
如何充分利用现有的二代测序数据准确计算癌症样本纯度和癌症细胞基因组倍性问题仍然是一项具有挑战性的工作。
参考文献
1、Junttila M R,de Sauvage F J.Influence of tumour micro-environmentheterogeneity on therapeutic response[J].Nature,2013,501(7467):346-354.
2、Carter S L,Cibulskis K,Helman E,et al.Absolute quantification ofsomatic DNA alterations in human cancer[J].Nature Biotechnology,2012,30(5):413-21.
3、Oesper L,Mahmoody A,Raphael B J.Inferring intra-tumor heterogeneityfrom high-throughput DNA sequencing data[J].Genome Biology,2013,14(7):R80.
4、Yuan Y,Failmezger H,Rueda O M,et al.Quantitative Image Analysis ofCellular Heterogeneity in Breast Tumors Complements Genomic Profiling[J].Science Translational Medicine,2012,4(157):157ra143.
5、Navin N,Kendall J,Troge J,et al.Tumour evolution inferred bysingle-cell sequencing[J].Nature,2011,472(7341):90-4.
6、Su X,Zhang L,Zhang J,et al.PurityEst:estimating purity of humantumor samples using next-generation sequencing data.[J].Bioinformatics,2012,28(17):2265-2266.
7、Larson N B.PurBayes:estimating tumor cellularity and subclonalityin next-generation sequencing data[J].Bioinformatics,2013,29(15):1888-9.
8、Gusnanto A,Wood H M,Pawitan Y,et al.Correcting for cancer genomesize and tumour cell content enables better estimation of copy numberalterations from next-generation sequence data[J].Bioinformatics,2012,28(1):40-47.
9、Oesper L,Mahmoody A,Raphael B J.Inferring intra-tumor heterogeneityfrom high-throughput DNA sequencing data[J].Genome Biology,2013,14(7):R80.
10、Carter S L,Cibulskis K,Helman E,et al.Absolute quantification ofsomatic DNA alterations in human cancer[J].Nature Biotechnology,2012,30(5):413-21.
11、Li Y,Xie X.Deconvolving tumor purity and ploidy by integratingcopy number alterations and loss of heterozygosity[J].Bioinformatics,2015,30(4):2121.
12、Mayrhofer M,Dilorenzo S,Isaksson A.Patchwork:allele-specific copynumber analysis of whole-genome sequenced tumor tissue[J].Genome Biology,2013,14(3):R24.
发明内容
术语定义
为了更好地理解本发明，下面提供相关的解释和说明：
Whole Genome Sequencing(WGS)：使用二代测序技术的全基因组测序。
read：高通量测序平台产生的测序序列。
测序深度：测序得到的碱基(bp)总量与基因组(Genome)大小的比值，它是评价测序量的指标之一。
window(窗口)：按照一定长度划分的基因组片段，该长度代表window大小。本方法中window大小可由使用者自由设置，通常设置为几百碱基。一个大基因组片段S可以包含大量window。
Tumor Read Enrichment(TRE)：癌症片段读长富集程度es，指癌症样本中某一片段S内read数量与相应正常样本中对应片段read数量的比值，定义公式如下：
公式(1)中，和分别表示在癌症样本中覆盖片段s的read数量和相匹配的正常样本中覆盖片段s的read数量，Nt表示癌症样本全基因组测序获得read总数量，Nn表示相应正常样本全基因组测序获得read总数量。
Heterozygous Germline Single Nucleotide Variants(HGSNV)：杂合生殖系细胞单碱基变异，由于人类染色体属于二倍体，体细胞均由胚胎细胞发育而来，而生殖细胞中HGSNV位点只有两种碱基类型A和B，其中一种来源于父本，另一种来源于母本。
Major Allele Fraction(MAF)：主要等位基因(allele)分数，本发明中使用的HGSNV只有两种等位基因，一种等位基因与参考基因组相同，另一种与参考基因组不同。这两种等位基因型分数的计算方法为覆盖某一等位基因的read数量除以覆盖该位点总read数量的比值，MAF就是两种等位基因分数中的较大值。计算公式如(1.1)所示，nr为包含与参考基因组相同等位基因的read数量，na为包含另一种等位基因的read的数量，nt表示覆盖该HGSNV位点的总read数量，C为该HGSNV的MAF值。MAF是相对于HGSNV的概念，本发明中“片段的MAF”指片段内所有HGSNV的MAF均值，“peak的MAF”指peak中所有片段包含的HGSNV的MAF均值。
major allele copy number：主要等位基因拷贝数，指在拷贝数为i的片段中，主要等位基因拷贝数的取值，它的取值范围为大于等于的整数。
peak：指基因组所有window的TRE分布中，聚集在一起的TRE簇。如附图1所示，图A表示基因组上所有window的TRE分布，纵轴表示对应某TRE位点的window总数量，该图为基因组GC含量校正之前的TRE分布，图B表示GC含量校正之后的TRE分布，图B中可以看到window明显以簇聚集，本方法将通过类自回归模型鉴定出来的TRE簇定义为peak，本质上是具有相同拷贝的基因组片段内window的聚集。
癌症样本：指从某患癌症的个体身上取下的癌症组织，它包含了一部分癌症细胞和一部分正常细胞。
P：指两个相邻peak之间的间距，由于peak是一个簇，这里peak的TRE由peak的TRE均值来表示，所以实际上是两个相邻peak的TRE均值的差。由于peak是具有相同拷贝数基因组片段内window的聚集，所以这里也表述为相邻拷贝数片段TRE的差值。
发明目的
本发明的目的在于克服现有技术的缺陷，提供一种全自动、高效率、高准确性的计算癌症样本纯度和染色体倍性的方法和装置。本发明在癌症样本纯度和染色体倍性计算上具有广阔的运用前景。
技术方案
为实现上述发明目的，本发明采取的技术方案为：通过癌症样本和匹配的正常样本的全基因组测序数据，对不同拷贝数片段的TRE和HGSNV的MAF分布构建混合高斯模型，计算癌症样本纯度和染色体倍性。
本发明主要运用了全基因组测序数据的TRE信息和HGSNV的MAF信息。TRE基本反映了癌症样本拷贝数变异情况，HGSNV的MAF信息基本反映了癌症样本的基因型。
TRE的差别主要来源于基因组片段的拷贝数差异，高拷贝数基因组片段内测序获得的read数量一定大于低拷贝数基因组片段测序获得的read数量，通过片段内read数量差异计算片段拷贝数差异是基因组拷贝数检测中的常用方法。但是大多数研究中，直接使用癌症样本片段内read数量除以正常样本该片段内read数量的比值(ratio)来进行read数量差异评估。本发明使用公式(1)所示的TRE来评估片不同片段的read数量差异。传统方法计算所得ratio不仅受癌症样本纯度与染色体倍性的影响，还受到癌症样本和正常样本测序深度的影响，而TRE不会受到样本测序深度的影响。
单独依赖read数量差异，无法确定各拷贝数片段的基因型，更重要的是无法区分样本纯度与样本倍性的补偿效应。而结合拷贝数差异片段内的HGSNV可以提供基因型信息，并帮助解决纯度与倍性的补偿效应，然而在前人的研究中，并没有高效利用HGSNV的方法，大部分方法采用枚举的方式将不同拷贝数片段可能对应的基因型一一列出，然后对排列组合的结果进行计算，挑选最可信的结果。这些方法共同的特点是，方法计算时间长，准确性差，对拷贝数很高或基因组变化较大的样本效果很差。本发明依据HGSNV的MAF与TRE的混合高斯模型计算癌症样本纯度和染色体倍性，能显著减少计算时间，并提高计算结果准确率。
假设某癌症样本的纯度为γ，那么癌症样本中的正常细胞比例为1–γ。癌症样本中正常细胞的染色体倍性为2，癌症细胞的染色体倍性为κ。那么癌症样本的染色体倍性ω如下公式(2)所示。
ω＝(1-γ)×2+γ×κ (2)
假设在癌症细胞中某一片段S的拷贝数为CS。那么癌症样本的片段S的拷贝数Ct应该为如下公式(3)所示。
Ct＝(1-γ)×2+γ×Cs (3)
对于基因组片段S，TRE的计算方式如公式(1)所示。而TRE的期望值(expectation)E(es)的推导公式如下公式(4)所示，式中的Nn和Nt与公式(1)中含义相同。
为了更进一步引出es，本方法定义了一些帮助理解的参数。片段S的长度LS，人类参考基因组的长度Lgw，癌症样本的测序深度正常样本的测序深度那么片段S在癌症样本中测序深度为片段S在正常样本中测序深度为λS是指与片段S特性(如GC含量等引起测序深度偏好的特性)有关的参数，所以在癌症和正常样本中是一样的。进一步通过γ，κ，Cs来表示es，如公式(5)所示。
公式(5)中Cs表示在癌症细胞中片段S的拷贝数，那么当片段S的拷贝数为i时和i+1时对应的TRE均值Si和Si+1分别如公式(6)和公式(7)所示：
通过公式(6)和公式(7)，对于相邻的拷贝数对应的片段，它们的TRE的差值P如公式(8)所示，可见P值的大小与片段具体拷贝数没有关系，它只决定于癌症样本纯度和染色体倍性。通过附图2可以直观的看到在TRE分布图中，peak之间的距离是恒定的。
此外，对于i＝2即拷贝数为2的片段，它们的TRE值Q如公式9所示。附图2中Q对应的peak的TRE值略大于1。
通过上述公式(8)和(9)，可以解得癌症样本的纯度(γ)和染色体倍性(κ)分别为：
通过以上分析可以得知，通过确定P和Q可以计算出癌症样本纯度γ和染色体倍性κ。
如附图2所示，计算出全基因组所有片段的TRE分布后，可以计算peak间的间距得到P。在前人的研究方法中，patchwork[12]使用相邻拷贝数片段的对应的read数量的比值的间距来辅助计算癌症样本纯度，但是该研究无法自动识别read数量比值之间的间距，需要人工识别图像确定read数量比值间距来进行下一步计算，效率和准确性都比较低。本发明开创性的使用类自回归模型鉴定TRE之间的间距，如公式(12)、(13)所示。公式(12)和(13)中，Xt表示0到Mt之间的TRE值；t表示扩大了1000倍的TRE值；Mt表示TRE的最大值；P表示两个TRE位点的间隔；C(Xt)表示在TRE为Xt的位点，对应的window数量；C(Xt+1000×P)表示在TRE为Xt+1000×P的位点，对应的window数量；Y(P)表示在P下，类自回归模型的函数值；显而易见当P＝0时，Y(P)的取值最大，但这时候的P并不是实际peak之间的间距。
以0.001为分辨率，遍历0到1之间的所有P值，然后求Y(P)。Y(P)的值分布如附图3所示。根据公式(13)的特点，我们可以知道，当P等于0时，Y(P)的值会是最大的，但此时的P并不是peak之间的间距。我们选择图中第二高峰中Y(P)的最大值对应的x轴坐标值P作为peak之间的间距P的计算结果。
如图1中B图所示的peak之所以簇状分布，是因为具有相同拷贝数的基因组片段的TRE值(指片段内所有window的TRE的均值)并不完全相等，同拷贝数片段TRE相互之间存在误差。该误差服从高斯分布，所以图B中的簇状分布被认为是高斯分布。
如附图2所示，P确定以后，peak会被识别出来，但是有部分基因组片段没有落在识别出的peak上，这些片段被称作亚克隆片段(subclone segmentation)。在考虑亚克隆片段的情况下，会对后面公式(17)和(18)所示的高斯模型取值有影响，进而会影响最终混合高斯模型的取值。由于在后续的分析中，本发明只需考虑落在peak位点的片段，由此排除了亚克隆片段的干扰。
在如图2所示的TRE分布图中，Q位点表示拷贝数为2的片段对应的TRE值。首先我们可以推测，如果癌症细胞基因组中，存在部分片段的拷贝数为1，部分片段的拷贝数为0，那么在Q位点之前应该存在两个peak分别对应拷贝数1和0。如果不存在拷贝数为1的片段，只存在拷贝数为0的片段，那么在Q位点之前距离2P的位点存在一个peak，而在Q位点之前距离P的位点peak的window数为0，这也就是图2所示的情形。另一种情况假若拷贝数为1和0的片段都没有，那么在Q位点之前距离P和距离2P的位点对应的window数都为0。那么对于TRE的分布图，对于Xf，即第一个出现的peak，它可能对应的拷贝数为2(拷贝数1和0的peak对应的window就是0)，也可能对应的拷贝数是1(拷贝数为0的片段对应peak的window为0)，也有可能对应的拷贝数是0。
通过上述分析，我们可以知道，图2中的第一个peak即Xf对应的片段的拷贝数有几种不同的可能，而每一种可能都会使Q对应不同的peak。本发明通过混合高斯模型计算出了最可能的Xf对应的片段的拷贝数，从而确定了Q的取值，最终得到了癌症样本纯度和染色体倍性。首先我们需要鉴定出Xf对应片段的拷贝数的所有可能值。本发明通过如下公式(13.1)来确定Xf的值。(13.1)中，C(Xf+P)表示在TRE为Xf+P的位点，对应的window数量，n表示Mt以内peak的最大数量。当f(Xf)取最大值时Xf为第一个peak的TRE均值。
然后使用公式(13.2)求Xf之前最多可能有几个peak。其中Xf表示第一个peak的TRE均值，P表示相邻拷贝数片段对应的peak之间的间距，floor表示向下取整数，当N＝0时，表示Xf之前没有peak，Xf对应的片段拷贝数为0；当N＝1时，表示Xf之前最多可能有1个peak，也可能没有peak；当N＝2时，表示Xf之前最多可能有2个peak，也可能只有一个peak或者没有peak；
对于Xf之前可能有(1，2，3…N)个peak的情形，每一种情形下都可以通过如下公式(13.3)计算出一个对应的Q值。根据Q的定义，我们知道Q为拷贝数为2的片段的peak对应的TRE值。首先可以推断拷贝数为0的片段对应的TRE值为Xf-n×P，其中n表示Xf之前peak的个数，取值范围是0到N之间的整数，P表示相邻拷贝数片段对应的peak之间的间距，Xf含义与公式(13.1)中相同，公式(13.3)如下所示，其中Qn表示在Xf之前存在n个peak时，Q的取值。
Qn＝Xf-n×P+2×P＝Xf+(2-n)×P，n∈[0，N] (13.3)
根据以上分析，可以得到对于Xf之前可能有(0，1，2，3…N)个peak的情形时，Qn取值可能为(Q0，Q1，Q2，Q3…QN)。而前面的类自回归模型已经计算出了P，那么对于每一个可能的Q，我们可以通过公式(10)和(11)计算出相应的γ和κ。本发明通过混合高斯模型计算出Xf之前最可能的peak数量n，从而确定了Xf对应的片段的拷贝数，进而确定了Q的取值，最终得到了癌症样本纯度和染色体倍性。具体方法如下说明。
对于Qn的每一个可能的取值，结合P，我们可以计算得到相应的γ，然后计算各个拷贝数片段内HGSNV的MAF的理论值。其理论计算方式如公式(14)所示。Cmcp表示主要等位基因的拷贝数(major allele copy number)，Ccp表示peak的整体拷贝数。f表示该peak内MAF的理论值。
但是在实际情况下当测序深度比较低时，f与真实值(期望值)会有较大的误差。这里需要进一步校正f，校正方法如公式(15)、(16)所示。式(15)中m是某peak内所有window中read数量的均值，v是peak内所有window中read数量的方差，所求得的p是用于负二项分布的随机变量成功的概率，r为随机变量失败的次数，这里的随机变量d为测序获取到的测序深度(read coverage)。
在某一个测序深度即d下，MAF实际上是服从以f为概率、d为实验次数的二项分布。本发明以如下公式(16)对f进行校正，得到MAF的期望值fb。公式中k表示在某个HGSNV位点，某一种等位基因(A或B)的数量，测得的等位基因总量为d(与测序深度相等)。
公式(13.2)表明每个基因组片段的拷贝数有N种可能。公式(14)表明，某一拷贝数片段可以有多种主要等位基因拷贝数，于是对每个基因组peak，可以算出多个f，也可以算出多个fb，取距离peak实际观测MAF均值最近的fb作为该peak的MAF期望。而全基因组有多个peak，每个peak的MAF期望不同，对应算出多个MAF期望值{fb}。考虑到某一peak内HGSNV的MAF会存在一定误差但也近似服从高斯分布，peak内所有HGSNV的MAF期望值可以由实际数据直接计算获得。假设某peak的基因型一定，那么通过比较peak的MAF观测值与peak的{fb}的值就可以判断该peak的拷贝数和基因型。也就可以计算出拷贝数为2的peak对应的TRE值即Q的位置。同时也为了进一步校正P，本发明提出了一种混合高斯模型对TRE和HGSNV的观测数据进行拟合。
通过前面的分析可以得知因为公式(12)中εt的存在，Xt并不能十分准确的代表各个peak的TRE均值。TRE的高斯分布模型如公式17所示：
其中，L(es；γ,κ)表示基因组片段TRE的似然函数。N表示基因组上的所有window的数量。I表示基因组中所有片段的最大的拷贝数。σi表示拷贝数为i的所有片段的TRE的标准差。es为第s个window的TRE观测值，Si表示第i个peak的TRE均值。pi表示第s个window的拷贝数为i的权重，本公式中对所有的i，pi均取值为1。该公式表明似然函数的大小与Si的取值相关，当Si与eS越接近，似然函数的取值越大，同时也表示P值越接近真实值。使用L(es；γ,κ)的极大似然估计可以计算出较合理的P值。
然而在部分情况下，P在较小的区间内波动时(如[-0.005,0.005])，对应的似然函数值可能相同。本发明通过结合HGSNV的高斯分布模型如公式(18)所示，来更近一步确定P值，同时确定Q值。
其中，L(fs；γ,κ)表示HGSNV的似然函数。M表示基因组中所有HGSNV。S表示第S个HGSNV。I表示基因组中所有片段的最大的拷贝数。Fi,j为拷贝数为i、主要等位基因的拷贝数为j的片段内HGSNV的MAF期望值，即公式(16)计算出的fb。fs表示该片段内MAF的观测值均值，σi,j表示该片段内HGSNV的MAF观测值的标准差。pi,j表示在主要等位基因的拷贝数为j时，高斯分布的权重，对所有的i和j，pi,j取值均为1。pi表示第S个HGSNV所在片段的拷贝数为i的权重，对所有的i，pi取值均为1。公式(18)表明，似然函数的大小与Fi,j值相关，当Fi,j与fs越接近，似然函数越大，表明fs越准确，同时表明公式(14)中的f越准确，从而可以得到各片段对应的Ccp与Cmcp。于是确定了Q的取值。为了得到最准确的P与Q，本方法将公式(17)和公式(18)相加，得到混合高斯模型。
但对该混合模型统计计算容易发生模型过拟合现象。本发明进一步使用了贝叶斯信息准则(Bayesian Information Criterion，BIC)方法，给混合高斯模型一个罚分函数，用于控制模型的过拟合，最终混合高斯模型如公式(19)所示：
BIC(es,fs；γ,κ)＝-2×logL(fs；γ,κ)-2×logL(es；γ,κ)+I×log(N)+J×log(M) (19)
其中，BIC(Ss,fs；γ,κ)表示混合模型的似然函数，I是公式(17)中高斯分布的个数，J是公式(18)中高斯分布的个数。N是基因组中窗口的数量，M是基因组中HGSNV的个数。
通过遍历[P-0.02,P+0.02]的区间，对所有的(P，Qn)求极大似然估计，可以得到最合适的P值与Q值然后根据公式(10)和(11)可计算癌症样本的纯度和染色体倍性。
因此，本发明一方面提供了一种用于计算癌症样本中癌症细胞纯度和染色体倍性的方法，所述方法包括以下步骤：
步骤A：
获取配对的(来自同一癌症病人的)癌症组织样本和正常组织样本的全基因组测序(WGS)数据，并将测序数据比对到参考基因组；
步骤B：
从步骤A得到的比对结果文件中，提取read位置和长度信息，HGSNV位点和覆盖该位点的read数量信息，计算所有HGSNV的MAF，其中，计算公式如(1.1)所示：
公式(1.1)中，nr为包含与参考基因组相同等位基因的read数量，na为包含另一种等位基因的read的数量，nt表示覆盖该HGSNV位点的总read数量，C为该HGSNV的MAF值；
步骤C：
根据步骤B得到的read位置和长度信息，以window为单位统计各window内包含的read数量，使用基因组GC含量校正所有window内read数量；
步骤D：
使用步骤C校正后的read数量，使用公式(1)计算每一个window的TRE，然后运用TRE，通过BIC-seq软件对基因组进行片段化，获得以拷贝数划分的基因组片段：
公式(1)中，和分别表示在癌症样本中覆盖片段s(这里指window)的read数量和在正常样本中覆盖片段s的read数量，Nt表示癌症样本总read数量，Nn表示相应正常样本总read数量，es为TRE值；
步骤E：
以步骤D中BIC-seq处理后的基因组片段为单位，统计片段内所有window的TRE的均值、方差和该片段内window数量，根据均值和方差对基因组每个片段的window数量进行平滑化(smooth)处理，使TRE的分布更均匀，然后将平滑化处理后所有片段的window分布汇总，得到基因组上window随TRE变化的分布结果；同时以片段为单位，计算片段中所有HGSNV的MAF的均值和方差；
步骤F：
使用如公式(12)、(13)所示的类自回归模型，计算相邻拷贝数片段内TRE的差值即P，具体方法为遍历一定范围的P，计算Y(P)，在Y(P)的分布中，选择第二高峰内Y(P)的最大值对应的P作为P的计算结果：
公式(12)和(13)中，Xt表示0到Mt之间的TRE值；t表示扩大了1000倍的TRE值；Mt表示TRE的最大值；变量P表示两个TRE位点的间隔；C(Xt)表示在TRE为Xt的位点，对应的window数量；C(Xt+1000×P)表示在TRE为Xt+1000×P的位点，对应的window数量；Y(P)表示在变量P下，类自回归模型的函数值；
步骤G：
根据步骤F得到的P，计算TRE分布中第一个实际观测peak的TRE均值，然后计算在第一个实际peak之前最多可能存在理论peak的数量N，最后当第一个实际peak之前存在n个理论peak时，计算Q的值，以Qn表示，其中步骤G可以包括：
G1：
根据步骤F计算的P，使用公式(13.1)，选取使公式(13.1)取最大值的Xf作为第一个实际观测peak的TRE均值：
公式(13.1)中，i表示第i个peak，C(Xf+P×i)表示在TRE为Xf+P×i的位点，对应的window数量，n表示Mt以内peak的最大数量，Mt表示TRE的最大值；
G2：
使用公式(13.2)，根据步骤F计算的P和步骤G1计算的Xf，计算在Xf之前最多可能存在的peak数量N：
公式(13.2)中，Xf表示第一个peak的均值，P表示相邻拷贝数片段对应的peak之间的间距，floor表示向下取整数；
G3：
利用步骤G2计算的N值，当n取0到N之间的整数时，使用公式(13.3)计算Qn的值：
Qn＝Xf-n×P+2×P＝Xf+(2-n)×P，n∈[0，N] (13.3)
公式(13.3)中，n表示Xf之前peak的数量，取值范围是0到N之间的整数，P表示相邻拷贝数片段对应的peak之间的间距，Xf表示第一个实际观测peak的TRE均值，Qn表示在Xf之前理论上存在n个peak时的Q值；
步骤H：
使用步骤F计算的P与步骤G计算的所有可能的Qn，使用公式(10)、(11)计算癌症样本纯度γ和染色体倍性κ：
公式(10)、(11)中，γ表示样本纯度，κ表示染色体倍性，那么对所有的(P，QN)都可以得到对应的(γ，κ)；
步骤I：
当n取[0,N]之间的某个整数值时，使用公式(13.4)计算第i个peak的TRE均值：
Ti＝Xf-n×P+i×P＝Xf+(i-n)×P，n∈[0，N] (13.4)
公式(13.4)中，n表示Xf之前peak的数量，取值范围是0到N之间的整数，P表示相邻拷贝数片段对应的peak之间的间距，Xf表示第一个实际观测peak的TRE均值，Ti表示第i个peak的TRE均值，
对于落在Ti附近的片段，认为该片段具有拷贝数i；对于没有落在Ti附近的片段，将其归类为亚克隆片段，在后续分析中剔除所有亚克隆片段；然后根据步骤H计算的癌症样本纯度γ和peak对应的拷贝数，可计算peak的MAF的期望fb，不同peak的MAF期望不同，对基因组上的所有peak，最终得到MAF期望的集合{fb}；同时计算各个peak的TRE均值和方差(或标准差)；
步骤J：
根据步骤F计算的P和步骤I计算的{fb}构建如公式(19)所示的用“贝叶斯信息准则”校正后的混合高斯分布模型，然后对模型极大似然估计；其中，步骤J可以包括如下几步：
J1：
以步骤F计算的P构建如公式(17)所示的高斯分布模型：
公式(17)中，L(es；γ,κ)表示基因组片段TRE的似然函数，N表示基因组上的所有window的数量，I表示基因组中所有片段的最大的拷贝数，σi表示拷贝数为i的所有片段的TRE的标准差由步骤I得到，es为第s个window的TRE观测值，Si表示第i个peak的TRE均值即步骤I中的Ti，pi表示第s个window的拷贝数为i的权重，对所有的i，pi均取值为1；
J2：
以步骤I计算的fb构建如公式(18)所示的高斯分布模型：
公式(18)中，L(fs；γ,κ)表示HGSNV的似然函数，M表示基因组中所有HGSNV数量，S表示第S个HGSNV，I表示基因组中所有片段的最大的拷贝数；Fi,j表示拷贝数为i，主要等位基因的拷贝数为j的片段内HGSNV的MAF期望值，由步骤I得到；fs表示该片段内所有HGSNV的MAF的观测值均值，由步骤E得到，σi,j表示该片段内所有HGSNV的MAF观测值的标准差，由步骤E得到；pi,j表示在主要等位基因的拷贝数为j时，高斯分布的权重，对所有的i和j，pi,j取值均为1，pi表示第S个HGSNV所在片段的拷贝数为i的权重，对所有的i，pi取值均为1；
J3：
将(17)与(18)相加得到混合高斯模型，然后对混合模型进行BIC(BayesianInformation Criterion)校正得到最终混合模型如公式(19)：
BIC(es,fs；γ,κ)＝-2×logL(fs；γ,κ)-2×logL(es；γ,κ)+I×log(N)+J×log(M) (19)
公式(19)中，BIC(es,fs；γ,κ)表示混合模型的似然函数，I表示基因组中所有片段的最大的拷贝数，J是公式(18)中j的取值个数，N是基因组中window的数量，M是基因组中HGSNV的个数，
对[0,N]范围内的每一个整数值n，可以通过步骤G得到Qn，也可以通过步骤I得到所有peak的MAF期望的集合{fb}，而一对(P，{fb})可以构建一个公式(19)所示的模型，实质上是对每一对(P，Qn)，可以构建一个公式(19)所示的模型；
步骤K：
以0.001为分辨率，对[P-m,P+m]区间的所有P值，重复步骤G～J，可以得到一系列不同的(P，Qn)与对应的似然函数值，取最大的似然函数值对应的(P，Qn)作为最合适的P和Q值，m是0到0.5之间的一个值；
步骤L：
查询步骤H的结果，可以找到在步骤K得到的(P，Q)下，对应的癌症样本纯度和染色体倍性。
另一方面，本发明提供了一种用于计算癌症样本中癌症细胞纯度和染色体倍性的装置，其包括处理器，所述处理器用于运行程序，所述程序运行时执行以下步骤：
步骤A：
获取配对的(来自同一癌症病人的)癌症组织样本和正常组织样本的全基因组测序(WGS)数据，并将测序数据比对到参考基因组；
步骤B：
从步骤A得到的比对结果文件中，提取read位置和长度信息，HGSNV位点和覆盖该位点的read数量信息，计算所有HGSNV的MAF，其中，计算公式如(1.1)所示：
公式(1.1)中，nr为包含与参考基因组相同等位基因的read数量，na为包含另一种等位基因的read的数量，nt表示覆盖该HGSNV位点的总read数量，C为该HGSNV的MAF值；
步骤C：
根据步骤B得到的read位置和长度信息，以window为单位统计各window内包含的read数量，使用基因组GC含量校正所有window内read数量；
步骤D：
使用步骤C校正后的read数量，使用公式(1)计算每一个window的TRE，然后运用TRE，通过BIC-seq软件对基因组进行片段化，获得以拷贝数划分的基因组片段：
公式(1)中，和分别表示在癌症样本中覆盖片段s(这里指window)的read数量和在正常样本中覆盖片段s的read数量，Nt表示癌症样本总read数量，Nn表示相应正常样本总read数量，es为TRE值；
步骤E：
以步骤D中BIC-seq处理后的基因组片段为单位，统计片段内所有window的TRE的均值、方差和该片段内window数量，根据均值和方差对基因组每个片段的window数量进行平滑化(smooth)处理，使TRE的分布更均匀，然后将平滑化处理后所有片段的window分布汇总，得到基因组上window随TRE变化的分布结果；同时以片段为单位，计算片段中所有HGSNV的MAF的均值和方差；
步骤F：
使用如公式(12)、(13)所示的类自回归模型，计算相邻拷贝数片段内TRE的差值即P，具体方法为遍历一定范围的P，计算Y(P)，在Y(P)的分布中，选择第二高峰内Y(P)的最大值对应的P作为P的计算结果：
公式(12)和(13)中，Xt表示0到Mt之间的TRE值；t表示扩大了1000倍的TRE值；Mt表示TRE的最大值；变量P表示两个TRE位点的间隔；C(Xt)表示在TRE为Xt的位点，对应的window数量；C(Xt+1000×P)表示在TRE为Xt+1000×P的位点，对应的window数量；Y(P)表示在变量P下，类自回归模型的函数值；
步骤G：
根据步骤F得到的P，计算TRE分布中第一个实际观测peak的TRE均值，然后计算在第一个实际peak之前最多可能存在理论peak的数量N，最后当第一个实际peak之前存在n个理论peak时，计算Q的值，以Qn表示，其中步骤G可以包括：
G1：
根据步骤F计算的P，使用公式(13.1)，选取使公式(13.1)取最大值的Xf作为第一个实际观测peak的TRE均值：
公式(13.1)中，i表示第i个peak，C(Xf+P×i)表示在TRE为Xf+P×i的位点，对应的window数量，n表示Mt以内peak的最大数量，Mt表示TRE的最大值；
G2：
使用公式(13.2)，根据步骤F计算的P和步骤G1计算的Xf，计算在Xf之前最多可能存在的peak数量N：
公式(13.2)中，Xf表示第一个peak的均值，P表示相邻拷贝数片段对应的peak之间的间距，floor表示向下取整数；
G3：
利用步骤G2计算的N值，当n取0到N之间的整数时，使用公式(13.3)计算Qn的值：
Qn＝Xf-n×P+2×P＝Xf+(2-n)×P，n∈[0，N] (13.3)公式(13.3)中，n表示Xf之前peak的数量，取值范围是0到N之间的整数，P表示相邻拷贝数片段对应的peak之间的间距，Xf表示第一个实际观测peak的TRE均值，Qn表示在Xf之前理论上存在n个peak时的Q值；
步骤H：
使用步骤F计算的P与步骤G计算的所有可能的Qn，使用公式(10)、(11)计算癌症样本纯度γ和染色体倍性κ：
公式(10)、(11)中，γ表示样本纯度，κ表示染色体倍性，那么对所有的(P，QN)都可以得到对应的(γ，κ)；
步骤I：
当n取[0,N]之间的某个整数值时，使用公式(13.4)计算第i个peak的TRE均值：
Ti＝Xf-n×P+i×P＝Xf+(i-n)×P，n∈[0，N] (13.4)
公式(13.4)中，n表示Xf之前peak的数量，取值范围是0到N之间的整数，P表示相邻拷贝数片段对应的peak之间的间距，Xf表示第一个实际观测peak的TRE均值，Ti表示第i个peak的TRE均值。
对于落在Ti附近的片段，认为该片段具有拷贝数i；对于没有落在Ti附近的片段，将其归类为亚克隆片段，在后续分析中剔除所有亚克隆片段；然后根据步骤H计算的癌症样本纯度γ和peak对应的拷贝数，可计算peak的MAF的期望fb，不同peak的MAF期望不同，对基因组上的所有peak，最终得到MAF期望的集合{fb}；。同时计算各个peak的TRE均值和方差(或标准差)；
步骤J：
根据步骤F计算的P和步骤I计算的{fb}构建如公式(19)所示的用“贝叶斯信息准则”校正后的混合高斯分布模型，然后对模型极大似然估计；其中，步骤J可以包括如下几步：
J1：
以步骤F计算的P构建如公式(17)所示的高斯分布模型：
公式(17)中，L(es；γ,κ)表示基因组片段TRE的似然函数，N表示基因组上的所有window的数量，I表示基因组中所有片段的最大的拷贝数，σi表示拷贝数为i的所有片段的TRE的标准差由步骤I得到，es为第s个window的TRE观测值，Si表示第i个peak的TRE均值即步骤I中的Ti，pi表示第s个window的拷贝数为i的权重，对所有的i，pi均取值为1；
J2：
以步骤I计算的fb构建如公式(18)所示的高斯分布模型：
公式(18)中，L(fss；γ,κ)表示HGSNV的似然函数，M表示基因组中所有HGSNV数量，S表示第S个HGSNV，I表示基因组中所有片段的最大的拷贝数；Fi,j表示拷贝数为i，主要等位基因的拷贝数为j的片段内HGSNV的MAF期望值，由步骤I得到；fs表示该片段内所有HGSNV的MAF的观测值均值，由步骤E得到，σi,j表示该片段内所有HGSNV的MAF观测值的标准差，由步骤E得到；pi,j表示在主要等位基因的拷贝数为j时，高斯分布的权重，对所有的i和j，pi,j取值均为1，pi表示第S个HGSNV所在片段的拷贝数为i的权重，对所有的i，pi取值均为1；
J3：
将(17)与(18)相加得到混合高斯模型，然后对混合模型进行BIC(BayesianInformation Criterion)校正得到最终混合模型如公式(19)：
BIC(es,fs；γ,κ)＝-2×logL(fs；γ,κ)-2×ogL(es；γ,κ)+I×log(N)+J×log(M) (19)
公式(19)中，BIC(es,fs；γ,κ)表示混合模型的似然函数，I表示基因组中所有片段的最大的拷贝数，J是公式(18)中j的取值个数，N是基因组中window的数量，M是基因组中HGSNV的个数。
对[0,N]范围内的每一个整数值n，可以通过步骤G得到Qn，也可以通过步骤I得到所有peak的MAF期望的集合{fb}，而一对(P，{fb})可以构建一个公式(19)所示的模型，实质上是对每一对(P，Qn)，可以构建一个公式(19)所示的模型；
步骤K：
以0.001为分辨率，对[P-m,P+m]区间的所有P值，重复步骤G～J，可以得到一系列不同的(P，Qn)与对应的似然函数值，取最大的似然函数值对应的(P，Qn)作为最合适的P和Q值，m是0到0.5之间的一个值；
步骤L：
查询步骤H的结果，可以找到在步骤K得到的(P，Q)下，对应的癌症样本纯度和染色体倍性。
作为一种优选方案，上述计算癌症样本纯度和染色体倍性的方法和装置中，所述步骤A中，采用1000基因组计划第三期(phase 3)项目使用的参考基因组hs37d5(ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequ ence/hs37d5.fa.gz)作为本发明的参考基因组,它包含了GRCh37中的所有染色体和零散序列(decoy sequences)。比对软件使用Burrows-Wheeler Aligner(BWA)，比对方法使用其中的bwa mem，最终获得癌症和正常样本的比对结果bam格式文件。
作为一种优选方案，上述计算癌症样本纯度和染色体倍性的方法和装置中，所述步骤B中，采用了samtools软件提取read的位置和长度信息，HGSNV的位点和覆盖该位点的read数量信息。使用samtools view命令提取read信息时，过滤掉序列比对质量(MAPQ)低于31的序列(参数-q 31,q表示过滤掉测序质量差的序列)，同时过滤掉未能正确匹配的read(参数-f 0x2-F 0x18，f表示提取符合一定要求的序列，F表示过滤符合一定要求的序列)。使用samtools mpileup命令提取HGSNV信息时，过滤掉序列比对质量(MAPQ)低于20的序列(参数-q 20)，并过滤掉碱基质量小于20的序列(参数-Q 20，Q表示过滤掉碱基质量差的序列)。选取等位基因频率(allele frequence)时，本发明使用samtools mpileup的-l参数。使用该参数需要提前准备一个包含SNP位点信息的bed格式文件。本发明方法提前收集了1000基因组(genome)计划(http://www.internationalgenome.org/)中，根据大量样本统计出来的杂合等位基因位点，并且过滤掉B-等位基因频率(B-allele frequence)小于0.05的位点，然后做成bed文件。使用“-l”参数在确保能提供充足的HGSNV位点基础上，大大加快了HGSNV位点的提取速度，提高了装置运行效率。
作为一种优选方案，上述计算癌症样本纯度和染色体倍性的方法和装置中，所述步骤C中，步骤C可以包括4步：
C1、将全基因组按照一定碱基长度的window为单位进行划分，对每个window统计覆盖该window的read数量，统计时以每条read的中点代表该read的位置；
C2、对参考基因组创建索引文件，提高GC含量的统计速度；
C3、以每个window的GC含量为自变量，以每个window的read数量为因变量，拟合read数量随GC含量变化的函数；
C4、使用拟合出的模型对全基因组read数量进行调整。
作为一种优选方案，上述计算癌症样本纯度和染色体倍性的方法和装置中，所述步骤C2中，本发明为参考基因组创建GC含量索引文件。对每一条染色体分别统计1、5、25、125个碱基间隔的区域内，鸟嘌呤(G)和胞嘧啶(C)的累积数量。那么在统计某一个window中的GC含量时，可以用a*125+b*25+c*5+d*1(其中a,b,c,d表示系数变量)的快速算法提取。例如想要统计某380bp区域内的GC含量时，可以分解为3*125+1*5形式，那么只需要读取特定索引文件的某个区域5碱基中的GC含量和某个区域125碱基中的GC含量即可。同时本发明将索引文件存储为二进制的格式，极大的加快了对特定区域的GC含量的提取。
作为一种优选方案，上述计算癌症样本纯度和染色体倍性的方法和装置中，所述步骤C3中，本发明使用步骤C1和步骤C2提取的各window GC含量，通过如下弹性网络模型拟合read数量随GC含量变化。本发明使用window的GC含量为变量x，使用x,x2,x3,x4,x5,x6作为弹性网络模型的输入变量，以read数量为输出变量，构建弹性网络模型如公式(20)所示。式中，y表示window内观测到的read数量，X表示输入变量矩阵，β表示变量系数矩阵，j表示变量系数下标，P表示系数总数，λ1和λ2表示罚分系数。
作为一种优选方案，上述计算癌症样本纯度和染色体倍性的方法和装置中，所述步骤C4中，使用步骤C3中的模型预测每一个window理论上的read数量μgc，基因组的平均GC含量定义为μ，window内观测到的read数量定义为y，window内校正后的read数量为Y。那么校正公式如下(21)所示：
作为一种优选方案，上述计算癌症样本纯度和染色体倍性的方法和装置中，所述步骤D中，本发明使用公式(1)计算每个window的TRE取值。然后运用TRE的值，使用BIC-seq软件对全基因组进行片段化(segmentation)。BIC-seq的思路是使用BayesianInformation Criterion(BIC)算法，统计相邻窗口的BIC值，值越小说明两个窗口越相似，然后将BIC值小于0的window合并，最终BIC-seq会按照片段拷贝数的差异，将全基因组分割为不同片段。每一个片段与相邻片段有不同的TRE均值，即拷贝数存在差异。
作为一种优选方案，上述计算癌症样本纯度和染色体倍性的方法和装置中，所述步骤E中，使用步骤D中BIC-seq处理后的基因组片段为单位，计算片段所包含的window数量，TRE的平均值以及方差。然后对片段的TRE进行smooth处理。处理方式如公式(22)所示。针对每一个基因组片段，以TRE的均值作为正态分布的均值μ，以TRE的方差作为正态分布的方差σ，计算出TRE在[μ-2σ,μ+2σ]范围内window数量的分布，定义v为TRE坐标，取值范围为[μ-2σ,μ+2σ]，分辨率为0.000，Cwin为该片段分配到v位点的window数量，CT表示该片段内window的总数。将所有片段的window根据TRE值smooth后，可使片段内的window数量呈现正态分布，对所有片段各TRE位点对应的window数求和汇总，可以得到基因组范围的window随TRE变化的分布。
作为一种优选方案，上述计算癌症样本纯度和染色体倍性的方法和装置中，所述步骤F中，以0.001为分辨率，遍历(0,1]范围内的所有P，使用类自回归模型，计算Y(P)的值。Y(P)表现为多峰分布，类似图3所示，图中横轴为P，纵轴表示Y(P)，本发明使用第二高峰内Y(P)的最大值对应的P作为P的计算结果，Mt是TRE的最大取值，这里将Mt设置为3。
作为一种优选方案，上述计算癌症样本纯度和染色体倍性的方法和装置中，所述步骤G中，步骤G包括3个步骤，步骤G1中，[0,1]的TRE区间作为变量Xf的取值范围，过滤掉C(Xf)小于1000的TRE位点，计算使公式(13.1)取最大值时的Xf作为第一个peak的均值点。
作为一种优选方案，上述计算癌症样本纯度和染色体倍性的方法和装置中，所述步骤I中，根据步骤H计算的癌症样本纯度γ和peak对应的拷贝数，可计算peak的MAF的期望fb。其中步骤I可以包含I1、I2、I3三个步骤。
I1，使用公式(14)计算peak内HGSNV的MAF理论值，公式(14)中，Cmcp表示主要等位基因的拷贝数(major allele copy number)，Ccp表示peak的整体拷贝数，由步骤I得到，f表示该peak内MAF的理论值,可见当Ccp较大时，f有多种不同的可能值。
I2，利用负二项分布估计覆盖每个HGSNV位点的read总数的概率，使用公式(15)计算负二项分布的概率p和失败次数r。公式(15)中，m是片段内所有window中read数量的均值，v是片段内所有window中read数量的方差，所求得的p是用于负二项分布的随机变量成功的概率，r为随机变量失败的次数，随机变量为覆盖某个HGSNV中的read数量。
I3，利用二项分布求得的覆盖某个HGSNV的read数的概率。结合在一定read数量下，HGSNV只有两种基因型，服从二项分布规律，利用公式(16)计算f的校正值fb(即f的期望)。同一个peak中，不同的Cmcp可以计算得到不同的fb，选择与该peak的MAF观测均值最接近的fb作为该peak的fb。公式(16)中，k表示在某个HGSNV位点，某一种等位基因(A或B)的数量，d为覆盖该HGSNV的read数量，r为随机变量失败的次数，p是用于负二项分布的随机变量成功的概率；
对每一个Qn，可推断获得基因组所有peak对应的拷贝数和癌症样本纯度，从而对每一个peak可求fb，进而可以得到所有peak的MAF的期望值得集合{fb}。
作为一种优选方案，上述计算癌症样本纯度和染色体倍性的方法和装置中，所述步骤K中，m取0.02，P值的遍历区间为[P-0.02,P+0.02]。
通过本发明提供的层次混合高斯模型，实现了对癌症样本纯度的快速和准确计算，节约了纯度估算的时间和经济成本，同时提高了计算结果的准确性。
附图说明
图1表示全基因组中window数量在TRE上的分布。其中，图A表示的是未进过GC校正的TRE分布，图B表示经过GC含量校正后的TRE分布图。
图2表示一种癌症细胞中的TRE分布的模型，图为smooth处理以后，图中的peak满足以P为周期的分布，少量不满足周期性分布的小peak被认为是亚克隆片段。Q表示拷贝数为2的peak，不存在拷贝数为1的片段，所以在大约0.6的位置的peak的window数量为0。
图3表示横轴为P纵轴为类自回归模型计算值得分布。
图4表示本发明方法和装置的流程图。
具体实施方式
为更好地说明本发明的目的、技术方案和优点，下面将结合附图和具体实施例对本发明作进一步说明。但是提供实施例仅用于说明目的，而本发明的范围不限于实施例。
使用本发明所述的装置计算癌症样本纯度和染色体倍性的流程图如图4所示。
实施例中，所使用的实验材料为TCGA(https://cancergenome.nih.gov/)数据库下载的样本(TCGA-AD-A5EJ)的正常组织TCGA-AD-A5EJ-10A和癌症组织TCGA-AD-A5EJ-01A的全基因组测序数据。计算平台为ubuntu 16.04，方法的具体实现为C++，python，R程序。
实施例：根据样本TCGA-CM-4746的癌症组织和正常组织的全基因组测序数据，使用层次性混合高斯模型计算癌症样本的纯度和染色体倍性。
一、收集样本数据，在TCGA中下载TCGA-CM-4746-01A的肿瘤样本和正常样本的全基因组测序数据。癌症样本bam文件大小为12.6G，正常样本bam文件大小为10.1G。将bam文件用PICARD软件处理为fastq文件。将fastq使用bwa mem比对到参考基因组hs37d5得到新的癌症样本和正常样本bam文件，文件大小分别为12.4G和9.9G。
二、下载1000genome项目提供的1到22号染色体的vcf文件(ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/)，使用GATK的SelectVariants方法，提取参考基因组hs37d5.fa中的等位基因频率大于5％的BIALLELIC位点作为潜在的HGSNV位点，最终得到5633774个biallele位点。
三、提取正常样本和癌症样本的read信息，同时提取癌症样本的HGSNV信息。使用samtools提取癌症样本的序列覆盖度和HGSNV，得到HGSNV 67732个。提取HGSNV时，采用上述步骤一获得的biallele位点作为备选位点表单。使用samtools view方法，直接从备选表单中提取HGSNV，加快提取速度。
四、以500bp为window对参考基因组建立GC含量的索引文件，对上述步骤一下载得到的参考基因组hs37d5.fa文件，建立1,5,25,125区段中的GC含量索引文件。存储为二进制格式。
五、以500bp为一个窗口，统计全基因组范围内每个窗口中的read数量。同时使用步骤四中产生的索引文件计算每个窗口中的GC含量。通过弹性网络模型对read数量进行GC含量校正。
六、对每一个窗口，使用矫正后的read数量计算TRE。并依据TRE，通过BIC-seq对基因组进行片段化。片段化的结果如表一所示，每一列数据表示了一个基因组片段的位置信息和TRE的均值，方差和片段内window的数量。
表1 BIC-seq对基因组片段化后的结果
七、通过步骤六得到了每个片段的TRE的均值和方差，以及该片段内包含的窗口数量。使用正态分布的方法，以每个片段的TRE均值和方差作为正态分布的均值和方差，将片段中的窗口按照正态分布进行平滑化。汇总所有片段smooth后的TRE以及对应窗口数的信息。
八、对smooth后的TRE的窗口数进行自回归分析，得到P的取值为0.386。
九、P等于0.386时，第一个实际观测peak的TRE均值为0.562，第一个实际观测peak前最多可以存在1个理论peak，即N＝1。可能的Q为：Q0＝1.334，Q1＝0.948，这两种Q的混合高斯模型的似然函数值分别为1.77E+07，1.78E+07。
十、计算P在取值范围[P-0.02,P+0.02]内，BIC校正后的混合高斯模型的极大似然值，计算结果如表2所示。
表2在P的取值范围内混合高斯模型的结果
十一、步骤十中的结果显示，P为0.382时，混合模型取极大值，此时的Q为0.948，据此可计算获得癌症样本纯度为0.80，癌症细胞染色体倍性为2.14。
Claims (13)
1.一种用于计算癌症样本中癌症细胞纯度和染色体倍性的方法，所述方法包括以下步骤：
步骤A：
获取配对的癌症组织样本和正常组织样本的全基因组测序数据，并将测序数据比对到参考基因组；
步骤B：
从步骤A得到的比对结果文件中，提取read位置和长度信息，HGSNV位点和覆盖该位点的read数量信息，计算所有HGSNV的MAF，其中，计算公式如(1.1)所示：
公式(1.1)中，nr为包含与参考基因组相同等位基因的read数量，na为包含另一种等位基因的read的数量，nt表示覆盖该HGSNV位点的总read数量，C为该HGSNV的MAF值；
步骤C：
根据步骤B得到的read位置和长度信息，以window为单位统计各window内包含的read数量，使用基因组GC含量校正所有window内read数量；
步骤D：
使用步骤C校正后的read数量，使用公式(1)计算每一个window的TRE，然后运用TRE，通过BIC-seq软件对基因组进行片段化，获得以拷贝数划分的基因组片段：
公式(1)中，和分别表示在癌症样本中覆盖片段s的read数量和在正常样本中覆盖片段s的read数量，Nt表示癌症样本总read数量，Nn表示相应正常样本总read数量，es为TRE值；
步骤E：
以步骤D中BIC-seq处理后的基因组片段为单位，统计片段内所有window的TRE的均值、方差和该片段内window数量，根据均值和方差对基因组每个片段的window数量进行平滑化处理，使TRE的分布更均匀，然后将平滑化处理后所有片段的window分布汇总，得到基因组上window随TRE变化的分布结果；同时以片段为单位，计算片段中所有HGSNV的MAF的均值和方差；
步骤F：
使用如公式(12)、(13)所示的类自回归模型，计算相邻拷贝数片段内TRE的差值即P，其中，遍历一定范围的P，计算Y(P)，在Y(P)的分布中，选择第二高峰内Y(P)的最大值对应的P作为P的计算结果：
公式(12)和(13)中，Xt表示0到Mt之间的TRE值；t表示扩大了1000倍的TRE值；Mt表示TRE的最大值；变量P表示两个TRE位点的间隔；C(Xt)表示在TRE为Xt的位点，对应的window数量；C(Xt+1000×P)表示在TRE为Xt+1000×P的位点，对应的window数量；Y(P)表示在变量P下，类自回归模型的函数值；
步骤G：
根据步骤F得到的P，计算TRE分布中第一个实际观测peak的TRE均值，然后计算在第一个实际peak之前最多可能存在理论peak的数量N，最后当第一个实际peak之前存在n个理论peak时，计算Q的值，以Qn表示，其中步骤G包括：
G1：
根据步骤F计算的P，使用公式(13.1)，选取使公式(13.1)取最大值的Xf作为第一个实际观测peak的TRE均值：
公式(13.1)中，i表示第i个peak，C(Xf+P×i)表示在TRE为Xf+P×i的位点，对应的window数量，n表示Mt以内peak的最大数量，Mt表示TRE的最大值；
G2：
使用公式(13.2)，根据步骤F计算的P和步骤G1计算的Xf，计算在Xf之前最多可能存在的peak数量N：
公式(13.2)中，Xf表示第一个peak的均值，P表示相邻拷贝数片段对应的peak之间的间距，floor表示向下取整数；
G3：
利用步骤G2计算的N值，当n取0到N之间的整数时，使用公式(13.3)计算Qn的值：
Qn＝Xf-n×P+2×P＝Xf+(2-n)×P，n∈[0，N] (13.3)
公式(13.3)中，n表示Xf之前peak的数量，取值范围是0到N之间的整数，P表示相邻拷贝数片段对应的peak之间的间距，Xf表示第一个实际观测peak的TRE均值，Qn表示在Xf之前理论上存在n个peak时的Q值；
步骤H：
使用步骤F计算的P与步骤G计算的Qn，使用公式(10)、(11)计算癌症样本纯度γ和染色体倍性κ：
公式(10)、(11)中，γ表示样本纯度，κ表示染色体倍性，由此对(P，QN)得到对应的(γ，κ)；
步骤I：
当n取[0,N]之间的某个整数值时，使用公式(13.4)计算第i个peak的TRE均值：
Ti＝Xf-n×P+i×P＝Xf+(i-n)×P，n∈[0，N] (13.4)
公式(13.4)中，n表示Xf之前peak的数量，取值范围是0到N之间的整数，P表示相邻拷贝数片段对应的peak之间的间距，Xf表示第一个实际观测peak的TRE均值，Ti表示第i个peak的TRE均值，
对于落在Ti附近的片段，认为该片段具有拷贝数i；对于没有落在Ti附近的片段，将其归类为亚克隆片段，在后续分析中剔除所有亚克隆片段；然后根据步骤H计算的癌症样本纯度γ和peak对应的拷贝数，计算peak的MAF的期望fb，不同peak的MAF期望不同，对基因组上的所有peak，最终得到MAF期望的集合{fb}；同时计算各个peak的TRE均值和方差或标准差；
步骤J：
根据步骤F计算的P和步骤I计算的{fb}构建如公式(19)所示的用“贝叶斯信息准则”校正后的混合高斯分布模型，然后对模型极大似然估计；其中，步骤J包括如下几步：
J1：
以步骤F计算的P构建如公式(17)所示的高斯分布模型：
公式(17)中，L(es；γ,κ)表示基因组片段TRE的似然函数，N表示基因组上的所有window的数量，I表示基因组中所有片段的最大的拷贝数，σi表示拷贝数为i的所有片段的TRE的标准差由步骤I得到，es为第s个window的TRE观测值，Si表示第i个peak的TRE均值即步骤I中的Ti，pi表示第s个window的拷贝数为i的权重，对所有的i，pi均取值为1；
J2：
以步骤I计算的fb构建如公式(18)所示的高斯分布模型：
公式(18)中，L(fs；γ,κ)表示HGSNV的似然函数，M表示基因组中所有HGSNV数量，S表示第S个HGSNV，I表示基因组中所有片段的最大的拷贝数；Fi,j表示拷贝数为i，主要等位基因的拷贝数为j的片段内HGSNV的MAF期望值，由步骤I得到；fs表示该片段内所有HGSNV的MAF的观测值均值，由步骤E得到；σi,j表示该片段内所有HGSNV的MAF观测值的标准差，由步骤E得到；pi,j表示在主要等位基因的拷贝数为j时，高斯分布的权重，对所有的i和j，pi,j取值均为1，pi表示第S个HGSNV所在片段的拷贝数为i的权重，对所有的i，pi取值均为1；
J3：
将(17)与(18)相加得到混合高斯模型，然后对混合模型进行BIC(BayesianInformation Criterion)校正得到最终混合模型如公式(19)：
BIC(es,fs；γ,κ)＝-2×log L(fs；γ,κ)-2×log L(es；γ,κ)+I×log(N)+J×log(M) (19)
公式(19)中，BIC(es,fs；γ,κ)表示混合模型的似然函数，I表示基因组中所有片段的最大的拷贝数，J是公式(18)中j的取值个数，N是基因组中window的数量，M是基因组中HGSNV的个数，
对[0,N]范围内的每一个整数值n，通过步骤G得到Qn，或者通过步骤I得到所有peak的MAF期望的集合{fb}，由一对(P，{fb})构建一个公式(19)所示的模型；
步骤K：
以0.001为分辨率，对[P-m,P+m]区间的所有P值，重复步骤G～J，得到一系列不同的(P，Qn)与对应的似然函数值，取最大的似然函数值对应的(P，Qn)作为最合适的P和Q值，m是0到0.5之间的一个值；
步骤L：
查询步骤H的结果，找到在步骤K得到的(P，Q)下，对应的癌症样本纯度和染色体倍性。
2.一种用于计算癌症样本中癌症细胞纯度和染色体倍性的装置，其包括处理器，所述处理器用于运行程序，所述程序运行时执行以下步骤：
步骤A：
获取配对的癌症组织样本和正常组织样本的全基因组测序数据，并将测序数据比对到参考基因组；
步骤B：
从步骤A得到的比对结果文件中，提取read位置和长度信息，HGSNV位点和覆盖该位点的read数量信息，计算所有HGSNV的MAF，其中，计算公式如(1.1)所示：
公式(1.1)中，nr为包含与参考基因组相同等位基因的read数量，na为包含另一种等位基因的read的数量，nt表示覆盖该HGSNV位点的总read数量，C为该HGSNV的MAF值；
步骤C：
根据步骤B得到的read位置和长度信息，以window为单位统计各window内包含的read数量，使用基因组GC含量校正所有window内read数量；
步骤D：
使用步骤C校正后的read数量，使用公式(1)计算每一个window的TRE，然后运用TRE，通过BIC-seq软件对基因组进行片段化，获得以拷贝数划分的基因组片段：
公式(1)中，和分别表示在癌症样本中覆盖片段s的read数量和在正常样本中覆盖片段s的read数量，Nt表示癌症样本总read数量，Nn表示相应正常样本总read数量，es为TRE值；
步骤E：
以步骤D中BIC-seq处理后的基因组片段为单位，统计片段内所有window的TRE的均值、方差和该片段内window数量，根据均值和方差对基因组每个片段的window数量进行平滑化处理，使TRE的分布更均匀，然后将平滑化处理后所有片段的window分布汇总，得到基因组上window随TRE变化的分布结果；同时以片段为单位，计算片段中所有HGSNV的MAF的均值和方差；
步骤F：
使用如公式(12)、(13)所示的类自回归模型，计算相邻拷贝数片段内TRE的差值即P，其中，遍历一定范围的P，计算Y(P)，在Y(P)的分布中，选择第二高峰内Y(P)的最大值对应的P作为P的计算结果：
公式(12)和(13)中，Xt表示0到Mt之间的TRE值；t表示扩大了1000倍的TRE值；Mt表示TRE的最大值；变量P表示两个TRE位点的间隔；C(Xt)表示在TRE为Xt的位点，对应的window数量；C(Xt+1000×P)表示在TRE为Xt+1000×P的位点，对应的window数量；Y(P)表示在变量P下，类自回归模型的函数值；
步骤G：
根据步骤F得到的P，计算TRE分布中第一个实际观测peak的TRE均值，然后计算在第一个实际peak之前最多可能存在理论peak的数量N，最后当第一个实际peak之前存在n个理论peak时，计算Q的值，以Qn表示，其中步骤G包括：
G1：
根据步骤F计算的P，使用公式(13.1)，选取使公式(13.1)取最大值的Xf作为第一个实际观测peak的TRE均值：
公式(13.1)中，i表示第i个peak，C(Xf+P×i)表示在TRE为Xf+P×i的位点，对应的window数量，n表示Mt以内peak的最大数量，Mt表示TRE的最大值；
G2：
使用公式(13.2)，根据步骤F计算的P和步骤G1计算的Xf，计算在Xf之前最多可能存在的peak数量N：
公式(13.2)中，Xf表示第一个peak的均值，P表示相邻拷贝数片段对应的peak之间的间距，floor表示向下取整数；
G3：
利用步骤G2计算的N值，当n取0到N之间的整数时，使用公式(13.3)计算Qn的值：
Qn＝Xf-n×P+2×P＝Xf+(2-n)×P，n∈[0，N](13.3) 公式(13.3)
中，n表示Xf之前peak的数量，取值范围是0到N之间的整数，P表示相邻拷贝数片段对应的peak之间的间距，Xf表示第一个实际观测peak的TRE均值，Qn表示在Xf之前理论上存在n个peak时的Q值；
步骤H：
使用步骤F计算的P与步骤G计算的Qn，使用公式(10)、(11)计算癌症样本纯度γ和染色体倍性κ：
公式(10)、(11)中，γ表示样本纯度，κ表示染色体倍性，由此对(P，QN)得到对应的(γ，κ)；
步骤I：
当n取[0,N]之间的某个整数值时，使用公式(13.4)计算第i个peak的TRE均值：
Ti＝Xf-n×P+i×P＝Xf+(i-n)×P，n∈[0，N] (13.4)
公式(13.4)中，n表示Xf之前peak的数量，取值范围是0到N之间的整数，P表示相邻拷贝数片段对应的peak之间的间距，Xf表示第一个实际观测peak的TRE均值，Ti表示第i个peak的TRE均值，
对于落在Ti附近的片段，认为该片段具有拷贝数i；对于没有落在Ti附近的片段，将其归类为亚克隆片段，在后续分析中剔除所有亚克隆片段；然后根据步骤H计算的癌症样本纯度γ和peak对应的拷贝数，计算peak的MAF的期望fb，不同peak的MAF期望不同，对基因组上的所有peak，最终得到MAF期望的集合{fb}；同时计算各个peak的TRE均值和方差或标准差；
步骤J：
根据步骤F计算的P和步骤I计算的{fb}构建如公式(19)所示的用“贝叶斯信息准则”校正后的混合高斯分布模型，然后对模型极大似然估计；其中，步骤J包括如下几步：
J1：
以步骤F计算的P构建如公式(17)所示的高斯分布模型：
公式(17)中，L(es；γ,κ)表示基因组片段TRE的似然函数，N表示基因组上的所有window的数量，I表示基因组中所有片段的最大的拷贝数，σi表示拷贝数为i的所有片段的TRE的标准差由步骤I得到，es为第s个window的TRE观测值，Si表示第i个peak的TRE均值即步骤I中的Ti，pi表示第s个window的拷贝数为i的权重，对所有的i，pi均取值为1；
J2：
以步骤I计算的fb构建如公式(18)所示的高斯分布模型：
公式(18)中，L(fs；γ,κ)表示HGSNV的似然函数，M表示基因组中所有HGSNV数量，S表示第S个HGSNV，I表示基因组中所有片段的最大的拷贝数；Fi,j表示拷贝数为i，主要等位基因的拷贝数为j的片段内HGSNV的MAF期望值，由步骤I得到；fs表示该片段内所有HGSNV的MAF的观测值均值，由步骤E得到，σi,j表示该片段内所有HGSNV的MAF观测值的标准差，由步骤E得到；pi,j表示在主要等位基因的拷贝数为j时，高斯分布的权重，对所有的i和j，pi,j取值均为1，pi表示第S个HGSNV所在片段的拷贝数为i的权重，对所有的i，pi取值均为1；
J3：
将(17)与(18)相加得到混合高斯模型，然后对混合模型进行BIC校正得到最终混合模型如公式(19)：
BIC(es,fs；γ,κ)＝-2×logL(fs；γ,κ)-2×logL(es；γ,κ)+I×log(N)+J×log(M) (19)
公式(19)中，BIC(es,fs；γ,κ)表示混合模型的似然函数，I表示基因组中所有片段的最大的拷贝数，J是公式(18)中j的取值个数，N是基因组中window的数量，M是基因组中HGSNV的个数，
对[0,N]范围内的每一个整数值n，通过步骤G得到Qn，或者通过步骤I得到所有peak的MAF期望的集合{fb}，由一对(P，{fb})构建一个公式(19)所示的模型；
步骤K：
以0.001为分辨率，对[P-m,P+m]区间的所有P值，重复步骤G～J，得到一系列不同的(P，Qn)与对应的似然函数值，取最大的似然函数值对应的(P，Qn)作为最合适的P和Q值，m是0到0.5之间的一个值；
步骤L：
查询步骤H的结果，找到在步骤K得到的(P，Q)下，对应的癌症样本纯度和染色体倍性。
3.根据权利要求1所述的方法或者权利要求2所述的装置，其中，所述步骤A中，采用1000基因组计划第三期(phase 3)项目使用的参考基因组hs37d5(ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz)作为所述参考基因组；和/或，比对软件使用Burrows-Wheeler Aligner(BWA)，比对方法使用其中的bwa mem，最终获得癌症和正常样本的比对结果bam格式文件。
4.根据权利要求1所述的方法或者权利要求2所述的装置，其中，所述步骤B中，采用samtools软件提取read的位置和长度信息，HGSNV的位点和覆盖该位点的read数量信息，其中，使用samtools view命令提取read信息时，使用参数-q 31过滤掉序列比对质量(MAPQ)低于31的序列，其中q表示过滤掉测序质量差的序列，同时使用参数-f0x2-F 0x18过滤掉未能正确匹配的read，其中f表示提取符合一定要求的序列，F表示过滤符合一定要求的序列，使用samtools mpileup命令提取HGSNV信息时，使用参数-q20过滤掉序列比对质量低于20的序列，并使用参数-Q 20过滤掉碱基质量小于20的序列，其中Q表示过滤掉碱基质量差的序列；选取等位基因频率时，使用samtools mpileup的-l参数；使用该参数需要提前准备一个包含SNP位点信息的bed格式文件。
5.根据权利要求1所述的方法或者权利要求2所述的装置，其中，
所述步骤C包括4步：
C1、将全基因组按照一定碱基长度的window为单位进行划分，对每个window统计覆盖该window的read数量，统计时以每条read的中点代表该read的位置；
C2、对参考基因组创建索引文件，提高GC含量的统计速度；
C3、以每个window的GC含量为自变量，以每个window的read数量为因变量，拟合read数量随GC含量变化的函数；
C4、使用拟合出的模型对全基因组read数量进行调整。
6.根据权利要求5所述的方法或者装置，其中，所述步骤C2中，为参考基因组创建GC含量索引文件，对每一条染色体分别统计1、5、25、125个碱基间隔的区域内，鸟嘌呤(G)和胞嘧啶(C)的累积数量，其中，在统计某一个window中的GC含量时，用a*125+b*25+c*5+d*1的快速算法提取，其中a,b,c,d表示系数变量。
7.根据权利要求5所述的方法或者装置，其中，所述步骤C3中，使用步骤C1和步骤C2提取的各window的GC含量，通过如下弹性网络模型拟合read数量随GC含量变化，其中，使用window的GC含量为变量x，使用x,x2,x3,x4,x5,x6作为弹性网络模型的输入变量，以read数量为输出变量，构建弹性网络模型如公式(20)所示：
公式(20)中，y表示window内观测到的read数量，X表示输入变量矩阵，β表示变量系数矩阵，j表示变量系数下标，P表示系数总数，λ1和λ2表示罚分系数。
8.根据权利要求5所述的方法或者装置，其中，所述步骤C4中，使用步骤C3中的模型预测每一个window理论上的read数量μgc，基因组的平均GC含量定义为μ，window内观测到的read数量定义为y，window内校正后的read数量为Y，那么校正公式如下(21)所示：
9.根据权利要求1所述的方法或者权利要求2所述的装置，其中，所述步骤E中，使用步骤D中BIC-seq处理后的基因组片段为单位，计算片段所包含的window数量，TRE的平均值以及方差，然后对片段的TRE进行平滑化处理，处理方式如公式(22)所示，针对每一个基因组片段，以TRE的均值作为正态分布的均值μ，以TRE的方差作为正态分布的方差σ，计算出TRE在[μ-2σ,μ+2σ]范围内window数量的分布，定义v为TRE坐标，取值范围为[μ-2σ,μ+2σ]，分辨率为0.001，Cwin为该片段分配到v位点的window数量，CT表示该片段内window的总数，将所有片段的window根据TRE值平滑化后，可使片段内的window数量呈现正态分布，对所有片段各TRE位点对应的window数求和汇总，得到基因组范围的window随TRE变化的分布：
10.根据权利要求1所述的方法或者权利要求2所述的装置，其中，所述步骤F中，以0.001为分辨率，遍历[0,1]范围内的所有P，使用类自回归模型，计算Y(P)的值，Y(P)表现为多峰分布，使用第二高峰内Y(P)的最大值对应的P作为P的计算结果，Mt是TRE的最大取值，这里将Mt设置为3。
11.根据权利要求1所述的方法或者权利要求2所述的装置，其中，所述步骤G中，步骤G包括3个步骤，步骤G1中，遍历[0,1]的TRE区间作为Xf，过滤掉C(Xf)小于1000的TRE位点，计算使公式(13.1)取最大值时的Xf作为第一个实际观测peak的均值。
12.根据权利要求1所述的方法或者权利要求2所述的装置，其中，所述步骤I中，然后根据步骤H计算的癌症样本纯度γ和peak对应的拷贝数，计算peak的MAF的期望fb，其中步骤I包括：
I1，使用公式(14)计算peak内HGSNV的MAF理论值：
公式(14)中，Cmcp表示主要等位基因的拷贝数，Ccp表示peak的整体拷贝数，由步骤I得到，f表示该peak内MAF的理论值，可见当Ccp较大时，f有多种不同的可能值；
I2，利用负二项分布估计覆盖每个HGSNV位点的read总数的概率，使用公式(15)计算负二项分布的概率p和失败次数r：
公式(15)中，m是peak内所有window中read数量的均值，v是peak内所有window中read数量的方差，所求得的p是用于负二项分布的随机变量成功的概率，r为随机变量失败的次数，随机变量为覆盖某个HGSNV中的read数量；
I3，利用二项分布求得的覆盖某个HGSNV的read数的概率，结合在一定read数量下，HGSNV只有两种基因型，服从二项分布规律，利用公式(16)计算f的校正值fb，同一个peak中，不同的Cmcp计算得到不同的fb，选择与该peak的MAF观测均值最接近的fb作为该peak的fb：
公式(16)中，k表示在某个HGSNV位点，某一种等位基因A或B的数量，d为覆盖该HGSNV的read数量，r为随机变量失败的次数，p是用于负二项分布的随机变量成功的概率；
对每一个Qn，可推断获得基因组所有peak对应的拷贝数和癌症样本纯度，从而对每一个peak求fb，进而得到所有peak的MAF的期望值得集合{fb}。
13.根据权利要求1所述的方法或者权利要求2所述的装置，其中，所述步骤K中，m取0.02作为P值得遍历区间为[P-0.02,P+0.02]。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710312237.7A CN108804876B (zh) | 2017-05-05 | 2017-05-05 | 用于计算癌症样本纯度和染色体倍性的方法和装置 |
| PCT/CN2018/078908 WO2018201805A1 (zh) | 2017-05-05 | 2018-03-14 | 用于计算癌症样本纯度和染色体倍性的方法和装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710312237.7A CN108804876B (zh) | 2017-05-05 | 2017-05-05 | 用于计算癌症样本纯度和染色体倍性的方法和装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108804876A true CN108804876A (zh) | 2018-11-13 |
| CN108804876B CN108804876B (zh) | 2022-03-15 |
Family
ID=64016930
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710312237.7A Expired - Fee Related CN108804876B (zh) | 2017-05-05 | 2017-05-05 | 用于计算癌症样本纯度和染色体倍性的方法和装置 |
Country Status (2)
| Country | Link |
|---|---|
| CN (1) | CN108804876B (zh) |
| WO (1) | WO2018201805A1 (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111627498A (zh) * | 2020-05-21 | 2020-09-04 | 北京吉因加医学检验实验室有限公司 | 一种测序数据gc偏向性校正的方法及其装置 |
| CN112216344A (zh) * | 2020-09-05 | 2021-01-12 | 西安翻译学院 | 肿瘤纯度和平均倍体信息的预测方法、系统、存储介质 |
| CN112767999A (zh) * | 2021-01-05 | 2021-05-07 | 中国科学院上海药物研究所 | 一种全基因组测序数据的分析方法及装置 |
| CN113808009A (zh) * | 2021-09-24 | 2021-12-17 | 熵智科技(深圳)有限公司 | 一种峰值初相位估计方法、装置、计算机设备及存储介质 |
| CN115948521A (zh) * | 2022-12-29 | 2023-04-11 | 东北林业大学 | 一种检测非整倍体缺失染色体信息的方法 |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110808084B (zh) * | 2019-09-19 | 2023-02-28 | 西安电子科技大学 | 一种基于单样本二代测序数据的拷贝数变异检测方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106460070A (zh) * | 2014-04-21 | 2017-02-22 | 纳特拉公司 | 检测染色体片段中的突变和倍性 |
| CN106520940A (zh) * | 2016-11-04 | 2017-03-22 | 深圳华大基因研究院 | 一种染色体非整倍体和拷贝数变异检测方法及其应用 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107075509B (zh) * | 2014-05-23 | 2021-03-09 | 数字基因公司 | 通过数字化转座子的单倍体组测定 |
| CN104560697A (zh) * | 2015-01-26 | 2015-04-29 | 上海美吉生物医药科技有限公司 | 一种基因组拷贝数不稳定性的检测装置 |
-
2017
- 2017-05-05 CN CN201710312237.7A patent/CN108804876B/zh not_active Expired - Fee Related
-
2018
- 2018-03-14 WO PCT/CN2018/078908 patent/WO2018201805A1/zh active Application Filing
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106460070A (zh) * | 2014-04-21 | 2017-02-22 | 纳特拉公司 | 检测染色体片段中的突变和倍性 |
| US20170107576A1 (en) * | 2014-04-21 | 2017-04-20 | Natera, Inc. | Detecting mutations and ploidy in chromosomal segments |
| CN106520940A (zh) * | 2016-11-04 | 2017-03-22 | 深圳华大基因研究院 | 一种染色体非整倍体和拷贝数变异检测方法及其应用 |
Non-Patent Citations (1)
| Title |
|---|
| YI LI等: ""Deconvolving tumor purity and ploidy by integrating copy number alterations and loss of heterozygosity"", 《BIOINFORMATICS》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111627498A (zh) * | 2020-05-21 | 2020-09-04 | 北京吉因加医学检验实验室有限公司 | 一种测序数据gc偏向性校正的方法及其装置 |
| CN111627498B (zh) * | 2020-05-21 | 2022-10-04 | 北京吉因加医学检验实验室有限公司 | 一种测序数据gc偏向性校正的方法及其装置 |
| CN112216344A (zh) * | 2020-09-05 | 2021-01-12 | 西安翻译学院 | 肿瘤纯度和平均倍体信息的预测方法、系统、存储介质 |
| CN112767999A (zh) * | 2021-01-05 | 2021-05-07 | 中国科学院上海药物研究所 | 一种全基因组测序数据的分析方法及装置 |
| CN113808009A (zh) * | 2021-09-24 | 2021-12-17 | 熵智科技(深圳)有限公司 | 一种峰值初相位估计方法、装置、计算机设备及存储介质 |
| CN115948521A (zh) * | 2022-12-29 | 2023-04-11 | 东北林业大学 | 一种检测非整倍体缺失染色体信息的方法 |
| CN115948521B (zh) * | 2022-12-29 | 2024-06-25 | 东北林业大学 | 一种检测非整倍体缺失染色体信息的方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108804876B (zh) | 2022-03-15 |
| WO2018201805A1 (zh) | 2018-11-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108804876A (zh) | 用于计算癌症样本纯度和染色体倍性的方法和装置 | |
| Tarabichi et al. | A practical guide to cancer subclonal reconstruction from DNA sequencing | |
| Gawad et al. | Dissecting the clonal origins of childhood acute lymphoblastic leukemia by single-cell genomics | |
| Cheema et al. | Computational approaches and software tools for genetic linkage map estimation in plants | |
| CN109689891A (zh) | 用于无细胞核酸的片段组谱分析的方法 | |
| Van de Wiel et al. | Preprocessing and downstream analysis of microarray DNA copy number profiles | |
| CN109994200A (zh) | 一种基于相似度融合的多组学癌症数据整合分析方法 | |
| Riester et al. | A differentiation-based phylogeny of cancer subtypes | |
| RU2013140708A (ru) | Способ для оценки потока информации в биологических сетях | |
| CN110016497B (zh) | 一种检测肿瘤单细胞基因组拷贝数变异的方法 | |
| CN110010195A (zh) | 一种探测单核苷酸突变的方法及装置 | |
| Yang et al. | Detecting recent positive selection with a single locus test bipartitioning the coalescent tree | |
| CN111223525A (zh) | 一种肿瘤外显子测序数据分析方法 | |
| CN106021986B (zh) | 超低频突变分子一致性序列简并算法 | |
| Roman et al. | A simplicial complex-based approach to unmixing tumor progression data | |
| CN110223732A (zh) | 多类生物序列注释的整合方法 | |
| WO2019236842A1 (en) | Difference-based genomic identity scores | |
| Moravec et al. | Testing for phylogenetic signal in single-cell RNA-seq data | |
| CN114974432A (zh) | 一种生物标志物的筛选方法及其相关应用 | |
| CN107563152A (zh) | 基于生物云平台的甲基化数据分析应用系统 | |
| Liu et al. | CRSCNV: A cross-model-based statistical approach to detect copy number variations in sequence data | |
| Wang et al. | A graph-based algorithm for estimating clonal haplotypes of tumor sample from sequencing data | |
| CN111850124A (zh) | 一种特征lincRNA表达谱组合及肺鳞癌早期预测方法 | |
| CN109754843A (zh) | 一种探测基因组小片段插入缺失的方法及装置 | |
| Dong et al. | AeQTL: eQTL analysis using region-based aggregation of rare genomic variants |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20220315 |