CN117253539B - 基于胚系突变检测高通量测序中样本污染的方法和系统 - Google Patents

基于胚系突变检测高通量测序中样本污染的方法和系统 Download PDF

Info

Publication number
CN117253539B
CN117253539B CN202311540335.8A CN202311540335A CN117253539B CN 117253539 B CN117253539 B CN 117253539B CN 202311540335 A CN202311540335 A CN 202311540335A CN 117253539 B CN117253539 B CN 117253539B
Authority
CN
China
Prior art keywords
snp
genotype
sample
pollution
snp locus
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
CN202311540335.8A
Other languages
English (en)
Other versions
CN117253539A (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.)
Beijing Qiuzhen Medical Laboratory Co ltd
Original Assignee
Beijing Qiuzhen Medical Laboratory 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 Beijing Qiuzhen Medical Laboratory Co ltd filed Critical Beijing Qiuzhen Medical Laboratory Co ltd
Priority to CN202311540335.8A priority Critical patent/CN117253539B/zh
Publication of CN117253539A publication Critical patent/CN117253539A/zh
Application granted granted Critical
Publication of CN117253539B publication Critical patent/CN117253539B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G16B20/20—Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
    • 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
    • G16B20/50—Mutagenesis
    • 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
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • G16B25/20—Polymerase chain reaction [PCR]; Primer or probe design; Probe optimisation
    • 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
    • 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
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B50/00—ICT programming tools or database systems specially adapted for bioinformatics
    • G16B50/30—Data warehousing; Computing architectures
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Theoretical Computer Science (AREA)
  • Medical Informatics (AREA)
  • General Health & Medical Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • Biophysics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Chemical & Material Sciences (AREA)
  • Genetics & Genomics (AREA)
  • Molecular Biology (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Analytical Chemistry (AREA)
  • Databases & Information Systems (AREA)
  • Bioethics (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Chemical Kinetics & Catalysis (AREA)
  • Epidemiology (AREA)
  • Evolutionary Computation (AREA)
  • Public Health (AREA)
  • Software Systems (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明涉及生物信息学技术领域，特别是涉及基于胚系突变检测高通量测序中样本污染的方法和系统。本发明提供的方法具体如下优点：利用基因型一致性和频率相关性共同判断肿瘤和对照样本是否配对，未配对的样本的检测结果不可信；根据纯合突变和杂合突变对最终突变频率影响的差异来确定污染的比例，且可以对污染源进行溯源；仅需肿瘤和对照样本的Vcf文件即可进行后续分析，使用简单，提高了运行效率，节省运行资源；不需要大量数据集来得到训练模型，只需利用panel本身的覆盖范围，并且不需要大量样本来做训练，即可得出准确结果。

Description

基于胚系突变检测高通量测序中样本污染的方法和系统
技术领域
本发明涉及生物信息学技术领域，特别是涉及基于胚系突变检测高通量测序中样本污染的方法和系统。
背景技术
NGS技术的主要步骤是：通过DNA杂交捕获或PCR扩增富集感兴趣的DNA区域，用条形码识别每个患者的DNA，通过NGS测序仪测序，并对原始数据进行生物信息学分析。
NGS技术的进步带来了更高的测序通量和更低的测序成本，可以分析比以前更多的样品，但是这也增加了样品混淆和污染的机会。虽然测序数据的质量在不断提高，但是样本污染仍然是DNA测序研究中的一个常见问题[1]，因为在很多过程中都有可能引入污染，如样本运输过程中（移液错误或干冰不足等），文库制备和扩增过程中（由于凝胶切割片段大小选择或条形码接头之间的意外切换），测序过程中，或者按照条形码拆分样本数据的过程中[1,2]。样本污染影响样本所有下游分析，并可能导致基因分型错误或者增加假阳性结果的产生。
癌症研究通常联合分析匹配的肿瘤和正常（T-N）样本，将测序的DNA与参考基因组进行比较，以检测肿瘤中存在的突变。目前最常用的方法之一是ContEst，这是基因组分析工具GATK软件中的一个模块。ContEst使用贝叶斯方法计算每个污染水平的后验概率，并找到纯合子位点污染水平的最大后验概率（MAP）估计；其次还有VerifyBamID，假设测试DNA样本中不包含超过一种污染物，VerifyBamID实现了基于可能性和基于回归的方法。通过在每个污染级别上使用网格搜索来最大化污染级别的可能性。Conpair专注于肿瘤-正常配对样本中的纯合位点，考虑到纯合标记不受拷贝数变化的影响，结合预先选择的高信息量基因组标记的短列表进行污染检测[3,4]。
目前的方法依赖于人类参考基因组数据，以及至少两个大的内存密集型文件：它们都需要肿瘤和对照样本的BAM文件（Conpair），或者Vcf文件和BAM文件（VerifyBamID和ContEst）[3,4]。另外，现有技术通常先筛选出大量的位点，这些位点需要包含在panel中，就会增加panel的大小和成本[5,6]，还需要大量数据集通过群体训练得到相关参数[6]，现有的大多数方法通常不会判断样本是否可以配对，未配对的样本的检测结果不可信。
参考文献：
[1] Matthew Flickinger, Goo Jun, Goncalo R. Abecasis, MichaelBoehnke, Hyun Min Kan. Correcting for Sample Contamination in GenotypeCalling of DNA Sequence Data [J]. The American Journal of Human Genetics,2016, 97:284-290.
[2] Fan Zhang, Matthew Flickinger, Sarah A. Gagliano Taliun, et, al.Ancestry-agnostic estimation of DNA sample contamination from sequence reads[J]. Genome Research, 2020, 30:185-194.
[3] Tao Jiang, Martin Buchkovich, Alison Motsinger-Reif. Same-SpeciesContamination Detection with Variant Calling Information from Next GenerationSequencing [J]. bioRxiv,2019.
[4] Ewa A. Bergmann, Bo-Juen Chen, Kanika Arora, Vladimir Vacict,Michael C. Zody. Conpair: concordance and contamination estimator for matchedtumor-normal pairs [J]. Bioinformatics, 2016, 32(20):3196-3198.
[5] 用于检测高通量测序中样本污染的SNP位点的筛选方法及检测样本污染的方法。中国发明专利，申请号：CN201810997769.3。
[6] 用于判断样本配对或污染的位点组合及其筛选方法和应用。中国发明专利，申请号：CN202211064680.4。
发明内容
为了解决上述问题，本发明提供了基于胚系突变检测高通量测序中样本污染的方法和系统。本发明提供的方法只需要处理Vcf文件，使用简单，提高了运行效率，节省运行资源，而且利用基因型一致性和频率相关性共同判断肿瘤和对照样本是否配对，未配对的样本的检测结果不可信，另外不需要大量数据集来得到训练模型。
为了实现上述目的，本发明提供如下技术方案：
本发明提供了一种基于胚系突变检测高通量测序中样本污染的方法，包括以下步骤；获取待测样本的Fastp数据和待检测的SNP位点，所述待测样本包括肿瘤样本和对照样本；将所述Fastp数据处理为Vcf文件；根据所述待检测的SNP位点的基因型一致率和频率相关系数判断所述肿瘤样本和对照样本是否配对，若基因型一致率＞60%，则肿瘤样本和对照样本可以配对；若基因型一致率≤60%，且频率相关性系数＞50%，则肿瘤样本和对照样本可以配对；若基因型一致率≤60%，且频率相关性系数≤50%，则肿瘤样本和对照样本不能配对，不进行后续处理；所述基因型一致率根据式Ⅰ计算得到，所述频率相关系数根据式Ⅱ计算得到；
式Ⅰ；
其中，Consistency rate是基因型一致率，S是肿瘤样本中SNP位点的基因型和对照样本中SNP位点的基因型一致的个数，N是筛选的SNP位点总位点数；
式Ⅱ；
其中，r(X,Y)是频率相关系数，X是肿瘤样本中每个SNP位点的突变频率，Y是对照样本中每个SNP位点的突变频率，Cov(X,Y)为X与Y的协方差，Var[X]为X的方差，Var[Y]为Y的方差；
采用式Ⅲ对所述SNP位点在能够配对的肿瘤样本和对照样本中突变频率差的绝对值进行核密度估计：
式Ⅲ；
其中，K为核函数，h为带宽，是概率密度，n是SNP位点个数，x是以xi为中心带宽内的所有点，xi是第i个SNP位点在肿瘤样本和对照样本中突变频率差的绝对值；
将所述核密度结果进行高斯函数进行数据拟合，高斯函数公式如下：
式Ⅳ；
其中，x是式Ⅲ计算得到的概率密度，a为高斯曲线的峰值，b为其对应的横坐标，c为标准差；对拟合曲线中峰值对应的横坐标进行提取，最大的横坐标为污染比例p。
优选的，所述方法还包括：根据对照样本的SNP位点的基因型、污染比例p和肿瘤样本SNP位点的突变频率确定污染源；所述确定污染源包括：推测污染源SNP位点的基因型，同一批样本中其他肿瘤样本记为可能污染样本，根据所述污染源SNP位点的基因型和所述可能污染样本中SNP位点的基因型一致率确定污染源，基因型一致率≥90%的可能污染样本，为待测样本中肿瘤样本的污染源；
所述推测污染源SNP位点的基因型包括：
当对照样本的SNP位点的基因型为纯合时，若肿瘤样本SNP位点的突变频率为100%，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为100%-0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为100%-p，则污染源的基因型为野生型；当对照样本的SNP位点的基因型为杂合时，若肿瘤样本SNP位点的突变频率为0.5×(100%-p)+p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为50%，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为50%-p，则污染源的基因型为野生型；当对照样本的SNP位点的基因型为野生型时，若肿瘤样本SNP位点的突变频率为p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为0，则污染源的基因型为野生型。
优选的，获取所述待检测的SNP位点的方法包括：
1）从千人基因组数据库中筛选包含在待检测的目标区域中且等位频率为0.1~0.9的SNP位点；2）将步骤1）筛选的SNP位点进行第二次筛选；所述第二次筛选的规则包括：a.将每条染色体上的SNP位点按位置顺序进行排序；b.保留每条染色体上的第1个SNP位点，选择距离>第1个SNP位点100kb的SNP位点暂定为第2个SNP位点，若第2个SNP位点之后100kb以内没有SNP位点，则保留第2个SNP位点；若第2个SNP位点之后100kb以内有SNP位点，则保留等位频率最接近0.5的SNP位点为第2个SNP位点；c.根据b的方法从保留的第2个SNP位点起筛选到每条染色体的最后一个SNP位点，保留的所有SNP位点为所述待检测的SNP位点。
在本发明中，所述待检测的SNP位点如下：
表1 169个SNP位点信息
编号 染色体 突变位置 参考碱基 变异碱基 编号 染色体 突变位置 参考碱基 变异碱基
1 chr1 9784423 C T 86 chr7 142460865 T C
2 chr1 11205058 C T 87 chr8 18257795 C T
3 chr1 22214127 A G 88 chr8 29197672 A G
4 chr1 25291010 A T 89 chr8 30999280 G T
5 chr1 36937059 A G 90 chr8 31497672 G A
6 chr1 45293518 A G 91 chr8 32453358 G A
7 chr1 65310489 T C 92 chr8 32621844 G T
8 chr1 91980447 A G 93 chr8 41794934 C T
9 chr1 97981395 T C 94 chr8 41906095 A G
10 chr1 98348885 G A 95 chr8 57078933 G T
11 chr1 120612006 G A 96 chr8 90967711 A G
12 chr1 161479745 A G 97 chr8 95419698 A G
13 chr1 162737116 C G 98 chr8 108970367 A G
14 chr1 181725110 T C 99 chr8 118819578 C T
15 chr2 17962450 A G 100 chr9 5081780 G A
16 chr2 29455267 A G 101 chr9 8389364 C G
17 chr2 29940529 A T 102 chr9 8518143 T C
18 chr2 38298203 C G 103 chr9 21816758 G A
19 chr2 42396722 A G 104 chr9 27202870 A G
20 chr2 42515437 A G 105 chr9 79986057 A G
21 chr2 47601106 T C 106 chr9 93641175 C T
22 chr2 48010488 G A 107 chr9 98209594 G A
23 chr2 74300717 T C 108 chr9 101597549 T G
24 chr2 112765973 A G 109 chr9 139407932 A G
25 chr2 141032088 C T 110 chr10 43613843 G T
26 chr2 141260668 A G 111 chr10 50740876 G C
27 chr2 141457985 T A 112 chr10 90771829 T C
28 chr2 141571329 T C 113 chr10 95279506 A T
29 chr2 141751592 G A 114 chr10 104386934 T C
30 chr2 142567910 T C 115 chr10 123239112 G A
31 chr2 178682603 T C 116 chr11 9607032 A G
32 chr2 178936373 T G 117 chr11 68192690 G A
33 chr2 191874667 A G 118 chr11 69462910 G A
34 chr2 198265526 A G 119 chr11 103058126 C T
35 chr2 204732714 A G 120 chr11 128333503 T C
36 chr2 212251864 T C 121 chr11 128651893 A G
37 chr2 215645464 C G 122 chr12 11992168 G A
38 chr2 225362478 C T 123 chr12 69967862 G A
39 chr2 240003870 G A 124 chr12 77438436 T C
40 chr3 10138069 T G 125 chr12 121416622 C G
41 chr3 12475557 C T 126 chr13 22275394 A G
42 chr3 30713842 C T 127 chr13 28624294 G A
43 chr3 37053568 A G 128 chr13 103504517 T C
44 chr3 124951821 C T 129 chr14 38061742 C T
45 chr3 128614185 A C 130 chr14 62213848 T C
46 chr3 142281612 A G 131 chr14 75483812 T C
47 chr4 55152040 C T 132 chr14 81574959 A G
48 chr4 55602765 G C 133 chr14 92441066 C T
49 chr4 55972974 T A 134 chr15 38643574 T C
50 chr4 169606649 C A 135 chr15 51502844 A C
51 chr4 169799448 A G 136 chr15 88680684 G A
52 chr5 7870973 A G 137 chr15 89858602 T C
53 chr5 37036492 C T 138 chr16 9916204 C G
54 chr5 38955796 G A 139 chr16 14041958 T C
55 chr5 60200665 A G 140 chr16 15811062 C T
56 chr5 67522722 C T 141 chr16 16281007 A G
57 chr5 68531253 C T 142 chr16 68857441 T C
58 chr5 79966029 G A 143 chr16 69143577 A G
59 chr5 80168937 G A 144 chr16 69745145 G A
60 chr5 176520243 G A 145 chr17 29553485 G A
61 chr5 176636882 C T 146 chr17 41244936 G A
62 chr6 29911190 G A 147 chr17 56435885 G T
63 chr6 30864829 T C 148 chr17 59760996 A G
64 chr6 31238009 T C 149 chr17 62007498 A G
65 chr6 31324595 C G 150 chr17 63533768 G A
66 chr6 32190406 A G 151 chr17 64685078 G A
67 chr6 32407709 C A 152 chr18 20577669 G A
68 chr6 32485491 C T 153 chr18 52895531 T C
69 chr6 32552067 C T 154 chr19 10600442 G C
70 chr6 33048694 A G 155 chr19 11105608 T C
71 chr6 43738977 C T 156 chr20 31024274 T C
72 chr6 112382313 G T 157 chr20 40714479 G A
73 chr6 117725578 T A 158 chr20 41306600 A G
74 chr6 138196066 T G 159 chr20 50407162 T C
75 chr6 152129077 T C 160 chr20 52198340 T A
76 chr6 152265522 G C 161 chr20 57478807 C T
77 chr6 157405930 G A 162 chr21 39870310 G A
78 chr6 160113872 A G 163 chr21 42845383 A G
79 chr6 161807855 C G 164 chr22 23631801 T C
80 chr6 162622197 C T 165 chrX 39932907 T C
81 chr7 55214348 C T 166 chrX 44938563 G A
82 chr7 87138645 A G 167 chrX 66765627 G A
83 chr7 106513011 C T 168 chrX 76937963 G C
84 chr7 116436097 G A 169 chrX 100608191 G A
85 chr7 128846328 G C / / / / /
优选的，将所述Fastp数据处理为Vcf文件包括：用Bwa将Fastq数据和人类hg19参考基因组进行比对，得到bam文件；对所述bam文件进行去重和碱基质量校正，再进行变异检测和变异注释，得到肿瘤样本和对照样本的Vcf文件。
优选的，所述峰值为大于前后四个位点的纵坐标。
本发明还提供了一种基于胚系突变检测高通量测序中样本污染的系统，包括数据获取模块、数据分析模块和样本污染检测模块；
所述数据获取模块用于获取待测样本的Fastp数据，并将Fastp数据传输给数据分析模块，所述待测样本包括肿瘤样本和对照样本；
所述数据分析模块用于将Fastp数据处理为Vcf文件，并将Vcf文件传输给样本污染检测模块；
所述样本污染检测模块用于检测污染源，包括位点选择模块、样本配对检测分析模块和污染比例计算模块；
所述位点选择模块用于选择待检测的SNP位点，包括：
1）从千人基因组数据库中筛选包含在待检测的目标区域中且等位频率为0.1~0.9的SNP位点；2）将步骤1）筛选的SNP位点进行第二次筛选；所述第二次筛选的规则包括：a.将每条染色体上的SNP位点按位置顺序进行排序；b.保留每条染色体上的第1个SNP位点，选择距离>第1个SNP位点100kb的SNP位点暂定为第2个SNP位点，若第2个SNP位点之后100kb以内没有SNP位点，则保留第2个SNP位点；若第2个SNP位点之后100kb以内有SNP位点，则保留等位频率最接近0.5的SNP位点为第2个SNP位点；c.根据b的方法从保留的第2个SNP位点起筛选到每条染色体的最后一个SNP位点，保留的所有SNP位点为所述待检测的SNP位点；
所述样本配对检测分析模块用于检测肿瘤样本和对照样本是否配对，包括根据所述位点选择模块选择的SNP位点的基因型一致率和频率相关系数判断，若基因型一致率＞60%，则肿瘤样本和对照样本可以配对；若基因型一致率≤60%，且频率相关性系数＞50%，则肿瘤样本和对照样本可以配对；若基因型一致率≤60%，且频率相关性系数≤50%，则肿瘤样本和对照样本不能配对，不进行后续处理；所述基因型一致率根据式Ⅰ计算得到，所述频率相关系数根据式Ⅱ计算得到；
式Ⅰ；
其中，Consistency rate是基因型一致率，S是肿瘤样本中SNP位点的基因型和对照样本中SNP位点的基因型一致的个数，N是筛选的SNP位点总位点数；
式Ⅱ；
其中，r(X,Y)是频率相关系数，X是肿瘤样本中每个SNP位点的突变频率，Y是对照样本中每个SNP位点的突变频率，Cov(X,Y)为X与Y的协方差，Var[X]为X的方差，Var[Y]为Y的方差；
所述污染比例计算模块用于将能配对的肿瘤样本和对照样本计算污染比例，包括采用式Ⅲ对所述位点选择模块筛选得到的SNP位点在肿瘤和对照样本中突变频率差的绝对值进行核密度估计：
式Ⅲ；
其中，K为核函数，h为带宽，是概率密度，n是SNP位点个数，x是以xi为中心带宽内的所有点，xi是第i个SNP位点在肿瘤样本和对照样本中突变频率差的绝对值；
将所述核密度结果进行高斯函数进行数据拟合，高斯函数公式如下：
式Ⅳ；
其中，x是式Ⅲ计算得到的概率密度，a为高斯曲线的峰值，b为其对应的横坐标，c为标准差；对拟合曲线中峰值对应的横坐标进行提取，最大的横坐标为污染比例p。
优选的，所述样本污染检测模块还包括推测污染源模块；所述推测污染源模块用于根据对照样本的SNP位点的基因型、污染比例p和肿瘤样本SNP位点的突变频率确定污染源，包括：推测污染源SNP位点的基因型，同一批样本中其他肿瘤样本记为可能污染样本，根据所述污染源SNP位点的基因型和所述可能污染样本中SNP位点的基因型一致率确定污染源，基因型一致率≥90%的可能污染样本，为待测样本中肿瘤样本的污染源；所述推测污染源SNP位点的基因型包括：当对照样本的SNP位点的基因型为纯合时，若肿瘤样本SNP位点的突变频率为100%，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为100%-0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为100%-p，则污染源的基因型为野生型；当对照样本的SNP位点的基因型为杂合时，若肿瘤样本SNP位点的突变频率为0.5×(100%-p)+p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为50%，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为50%-p，则污染源的基因型为野生型；当对照样本的SNP位点的基因型为野生型时，若肿瘤样本SNP位点的突变频率为p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为0，则污染源的基因型为野生型。
有益效果：
本发明提供的基于胚系突变检测高通量测序中样本污染的方法具体如下优点：
1.利用基因型一致性和频率相关性共同判断肿瘤和对照样本是否配对，可以提高检测结果的可信度。
2.污染比例的计算方法，根据提取拟合曲线对应峰值的横坐标得出，从对实验数据和模拟数据的分析结果来看，结果准确性高。
3.仅需肿瘤和对照样本的Vcf文件即可进行后续分析，使用简单，提高了运行效率，节省运行资源。
4.胚系突变中杂合突变的频率接近50%，而纯合突变的频率接近100%，如果样本被污染，频率就会出现偏差，具有代表性。
5.不需要大量数据集来得到训练模型，只需利用panel本身的覆盖范围，并且不需要大量样本来做训练，即可得出准确结果。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍。
图1为本发明整体分析流程图；
图2为样本是否配对分析流程图；
图3为污染比例计算流程图；
图4为污染比例计算方法示例图；
图5为污染源基因型推测方法示例图；
图6为分析方案1样本配对判断结果；
图7为分析方案1污染比例计算结果；
图8为分析方案2样本配对判断结果；
图9为分析方案2污染比例计算结果。
具体实施方式
本发明提供了一种基于胚系突变检测高通量测序中样本污染的方法，包括以下步骤；
获取待测样本的Fastp数据和待检测的SNP位点，所述待测样本包括肿瘤样本和对照样本；将所述Fastp数据处理为Vcf文件；
根据所述待检测的SNP位点的基因型一致率和频率相关系数判断所述肿瘤样本和对照样本是否配对，若基因型一致率＞60%，则肿瘤样本和对照样本可以配对；若基因型一致率≤60%，且频率相关性系数＞50%，则肿瘤样本和对照样本可以配对；若基因型一致率≤60%，且频率相关性系数≤50%，则肿瘤样本和对照样本不能配对，不进行后续处理；所述基因型一致率根据式Ⅰ计算得到，所述频率相关系数根据式Ⅱ计算得到；
式Ⅰ；
其中，Consistency rate是基因型一致率，S是肿瘤样本中SNP位点的基因型和对照样本中SNP位点的基因型一致的个数，N是筛选的SNP位点总位点数；
式Ⅱ；
其中，r(X,Y)是频率相关系数，X是肿瘤样本中每个SNP位点的突变频率，Y是对照样本中每个SNP位点的突变频率，Cov(X,Y)为X与Y的协方差，Var[X]为X的方差，Var[Y]为Y的方差；
采用式Ⅲ对所述SNP位点在能够配对的肿瘤样本和对照样本中突变频率差的绝对值进行核密度估计：
式Ⅲ；
其中，K为核函数，h为带宽，是概率密度，n是SNP位点个数，x是以xi为中心带宽内的所有点，xi是第i个SNP位点在肿瘤样本和对照样本中突变频率差的绝对值；所述带宽优选为函数默认带宽：scott” - 1.059 * A * nobs ** (-1/5.), where A is min(std(x),IQR/1.34)；
将所述核密度结果进行高斯函数进行数据拟合，高斯函数公式如下：
式Ⅳ；
其中，x是式Ⅲ计算得到的概率密度，a为高斯曲线的峰值，b为其对应的横坐标，c为标准差；
对拟合曲线中峰值对应的横坐标进行提取，最大的横坐标为污染比例p。
在本发明中，所述峰值优选为为大于前后四个位点的纵坐标。
本发明优选对待测样本进行DNA提取和测序后，获取待测样本的Fastp数据，所述待测样本包括肿瘤（Tumor）样本和正常（Normal）样本。在本发明中，所述肿瘤样本优选包括肿瘤组织切片样品，所述对照样本优选包括血浆样本中的基因组DNA。
得到待测样本的Fastp数据后，本发明将所述Fastp数据处理为Vcf文件。在本发明中，将所述Fastp数据处理为Vcf文件优选包括：用Bwa将Fastq数据和人类hg19参考基因组进行比对，得到bam文件；对所述bam文件进行去重和碱基质量校正，再进行变异检测和变异注释，得到肿瘤样本和对照样本的Vcf文件。本发明提供的方法仅需肿瘤和对照样本的Vcf文件即可进行后续分析，使用简单，提高了运行效率，节省运行资源。
本发明获取所述待检测的SNP位点的方法优选包括：
1）从千人基因组数据库中筛选包含在待检测的目标区域中且等位频率为0.1~0.9的SNP位点；2）将步骤1）筛选的SNP位点进行第二次筛选；所述第二次筛选的规则包括：a.将每条染色体上的SNP位点按位置顺序进行排序；b.保留每条染色体上的第1个SNP位点，选择距离>第1个SNP位点100kb的SNP位点暂定为第2个SNP位点，若第2个SNP位点之后100kb以内没有SNP位点，则保留第2个SNP位点；若第2个SNP位点之后100kb以内有SNP位点，则保留等位频率最接近0.5的SNP位点为第2个SNP位点；c.根据b的方法从保留的第2个SNP位点起筛选到每条染色体的最后一个SNP位点，保留的SNP位点总和为所述待检测的SNP位点。
本发明通过第二次筛选确保每个SNP位点之间的距离>100kb，可以避免连锁不平衡对结果的影响，最终，保留下来169个位点用于后续分析。具体169个位点见表1。
获取待测样本的Fastp数据和待检测的SNP位点后，本发明根据所述待检测的SNP位点的基因型一致率和频率相关系数判断所述肿瘤样本和对照样本是否配对，可以配对的意思为两个样本来自同一个体或除污染部分外来自同一样本，可以得出准确的分析结果；不能配对的意思为两个样本可能不是来自同一个体，无法得出准确的分析结果，需检测人员寻找原因。本发明利用基因型一致性和频率相关性共同判断肿瘤和对照样本是否配对，可以增加后续分析结果的可靠性。
肿瘤样本和对照样本能够配对后，本发明采用式Ⅲ对所述SNP位点在能够配对的肿瘤样本和对照样本中突变频率差的绝对值进行核密度估计，然后将所述核密度结果进行高斯函数进行数据拟合，最后对拟合曲线中峰值对应的横坐标进行提取，最大的横坐标为污染比例p。本发明对肿瘤样本和对照样本的频率差值的绝对值进行核密度估计，对拟合曲线峰值的横坐标进行提取，可以确定污染比例，从而确定样本是否被污染。
得到污染比例p后，本发明优选根据对照样本的SNP位点的基因型、污染比例p和肿瘤样本SNP位点的突变频率确定污染源；
所述确定污染源优选包括：推测污染源SNP位点的基因型，同一批样本中其他肿瘤样本记为可能污染样本，根据所述污染源SNP位点的基因型和所述可能污染样本中SNP位点的基因型一致率确定污染源，基因型一致率≥90%的可能污染样本，为待测样本中肿瘤样本的污染源；
所述推测污染源SNP位点的基因型优选包括：当对照样本的SNP位点的基因型为纯合时，若肿瘤样本SNP位点的突变频率为100%，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为100%-0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为100%-p，则污染源的基因型为野生型；当对照样本的SNP位点的基因型为杂合时，若肿瘤样本SNP位点的突变频率为0.5×(100%-p)+p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为50%，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为50%-p，则污染源的基因型为野生型；当对照样本的SNP位点的基因型为野生型时，若肿瘤样本SNP位点的突变频率为p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为0，则污染源的基因型为野生型。
本发明根据纯合突变和杂合突变对最终突变频率影响的差异结合污染比例可以确定污染源的基因型，从而可以对污染源进行溯源。
在本发明提供的上述技术方案所述的方法基础上，本发明还提供了一种基于胚系突变检测高通量测序中样本污染的系统，包括数据获取模块、数据分析模块和样本污染检测模块；
所述数据获取模块用于获取待测样本的Fastp数据，并将Fastp数据传输给数据分析模块，所述待测样本包括肿瘤样本和对照样本；
所述数据分析模块用于将Fastp数据处理为Vcf文件，并将Vcf文件传输给样本污染检测模块；
所述样本污染检测模块用于检测污染源，包括位点选择模块、样本配对检测分析模块和污染比例计算模块；
所述位点选择模块用于选择待检测的SNP位点，包括：
1）从千人基因组数据库中筛选包含在待检测的目标区域中且等位频率为0.1~0.9的SNP位点；2）将步骤1）筛选的SNP位点进行第二次筛选；所述第二次筛选的规则包括：a.将每条染色体上的SNP位点按位置顺序进行排序；b.保留每条染色体上的第1个SNP位点，选择距离>第1个SNP位点100kb的SNP位点暂定为第2个SNP位点，若第2个SNP位点之后100kb以内没有SNP位点，则保留第2个SNP位点；若第2个SNP位点之后100kb以内有SNP位点，则保留等位频率最接近0.5的SNP位点为第2个SNP位点；c.根据b的方法从保留的第2个SNP位点起筛选到每条染色体的最后一个SNP位点，保留的所有SNP位点为所述待检测的SNP位点；
所述样本配对检测分析模块用于检测肿瘤样本和对照样本是否配对，包括根据所述位点选择模块选择得到的SNP位点的基因型一致率和频率相关系数判断，若基因型一致率＞60%，则肿瘤样本和对照样本可以配对；若基因型一致率≤60%，且频率相关性系数＞50%，则肿瘤样本和对照样本可以配对；若基因型一致率≤60%，且频率相关性系数≤50%，则肿瘤样本和对照样本不能配对，不进行后续处理；所述基因型一致率根据式Ⅰ计算得到，所述频率相关系数根据式Ⅱ计算得到；
式Ⅰ；
其中，Consistency rate是基因型一致率，S是肿瘤样本中SNP位点的基因型和对照样本中SNP位点的基因型一致的个数，N是筛选的SNP位点总位点数；
频率相关系数根据以下公式计算：
式Ⅱ；
其中，r(X,Y)是频率相关系数，X是肿瘤样本中每个SNP位点的突变频率（即每个SNP位点检测到的突变频率），Y是对照样本中每个SNP位点的突变频率（即每个SNP位点检测到的突变频率），Cov(X,Y)为X与Y的协方差，Var[X]为X的方差，Var[Y]为Y的方差；
所述污染比例计算模块用于将能配对的肿瘤样本和对照样本计算污染比例，包括采用式Ⅲ对所述位点选择模块筛选得到的SNP位点在肿瘤和对照样本中突变频率差的绝对值进行核密度估计：
式Ⅲ；
其中，K为核函数，h为带宽，是概率密度，n是SNP位点个数，x是以xi为中心带宽内的所有点，xi是第i个SNP位点在肿瘤样本和对照样本中突变频率差的绝对值；
将所述核密度结果进行高斯函数进行数据拟合，高斯函数公式如下：
式Ⅳ；
其中，x是式Ⅲ计算得到的概率密度，a为高斯曲线的峰值，b为其对应的横坐标，c为标准差；
对拟合曲线中峰值对应的横坐标进行提取，最大的横坐标为污染比例p。
在本发明中，所述数据获取模块优选包括：待测样本测序后将同一个样本所有lane上的数据合并，得到每个样本用于分析的Fastq数据。
在本发明中，所述数据分析模块优选包括：用Bwa将Fastq数据和人类hg19参考基因组进行比对，得到bam文件；对所述bam文件进行去重和碱基质量校正，最后进行变异检测和变异注释，得到肿瘤样本和对照样本的Vcf文件。
在本发明中，所述峰值优选为大于前后四个位点纵坐标的位置；所述横坐标优选为频率差的绝对值。
在本发明中，所述样本污染检测模块优选还包括推测污染源模块；所述推测污染源模块用于根据对照样本的SNP位点的基因型、污染比例p和肿瘤样本SNP位点的突变频率确定污染源，包括：
推测污染源SNP位点的基因型，同一批样本中其他肿瘤样本记为可能污染样本，根据所述污染源SNP位点的基因型和所述可能污染样本中SNP位点的基因型一致率确定污染源，基因型一致率≥90%的可能污染样本，为待测样本中肿瘤样本的污染源；所述推测污染源SNP位点的基因型包括：
当对照样本的SNP位点的基因型为纯合时，若肿瘤样本SNP位点的突变频率为100%，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为100%-0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为100%-p，则污染源的基因型为野生型；
当对照样本的SNP位点的基因型为杂合时，若肿瘤样本SNP位点的突变频率为0.5×(100%-p)+p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为50%，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为50%-p，则污染源的基因型为野生型；
当对照样本的SNP位点的基因型为野生型时，若肿瘤样本SNP位点的突变频率为p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为0，则污染源的基因型为野生型。
本发明提供的系统具有如下优势：
1.利用基因型一致性和频率相关性共同判断肿瘤和对照样本是否配对，增加后续分析结果的可靠性。
2.对频率差值的绝对值进行核密度估计，对拟合曲线的峰值进行提取，将较高的峰值对应的横坐标作为污染比例。
3.根据胚系突变的结果来检测样本是否被污染，根据纯合突变和杂合突变对最终突变频率影响的差异来确定污染的比例，且可以对污染源进行溯源。
为了进一步说明本发明，下面结合实施例对本发明提供的一种基于胚系突变检测高通量测序中样本污染的系统进行详细地描述，但不能将它们理解为对本发明保护范围的限定。
实施例1
一种基于胚系突变检测高通量测序中样本污染的系统，所述系统由数据获取模块、数据分析模块和样本污染检测模块组成；
数据获取模块：待分析样本经过DNA提取、文库构建、文库质检之后通过二代测序仪MGI2000进行上机测序，测序完成后将同一个样本所有lane上的数据合并，得到每个样本用于分析的原始Fastq数据；
数据分析模块：使用Fastp对原始的Fastq数据进行质控，来去除接头和低质量的reads，包括去除ployG序列，ployG长度>10会被去除，去除序列长度<15的过短序列，去除含N碱基>5的reads，过滤低质量碱基（碱基质量<15）占比>40%的reads。然后用Bwa将过滤后的序列和人类hg19参考基因组进行比对，再对得到的bam文件进行去重和碱基质量校正，最后进行变异检测和变异注释即可得到用于分析样本污染的肿瘤（Tumor）样本和其配对的正常（Normal）样本的Vcf文件。
样本污染检测模块：该模块主要包含4个小模块：位点选择模块、样本配对检测分析模块、污染比例计算模块和推测污染源模块。
1）位点选择模块：
a.筛选1000人基因组数据库中等位频率介于0.1~0.9之间（包括0.1和0.9）的SNP位点；
b.筛选出待检测的目标区域中包含的位点；
c.确保每个SNP位点之间的距离>100kb，避免连锁不平衡对结果的影响，保留等位频率接近0.5的位点。
最终，保留下来169个位点用于后续分析。具体169个位点见表1。
2）样本配对检测分析模块：
分析过程如图2所示，根据数据分析模块中Tumor和Normal样本变异检测之后得到的Vcf文件，对上述位点选择模块得到的169个位点，计算基因型一致率和频率相关系数；所述基因型一致率根据式Ⅰ计算得到，所述频率相关系数根据式Ⅱ计算得到。若一致性大于60%，则两个样本可以配对；若一致性≤60%，频率相关性系数＞50%，两个样本可以配对，否则两个样本不能看作配对样本。无法配对的两个样本可能不是来自同一个体，无法得出准确的分析结果，需要查找原因。
3）污染比例计算模块
假设前提：污染源纯合突变或杂合突变对样本频率的影响程度是不同的，理论上来说频率差的绝对值应该分为两个梯度。
分析过程如图3所示，在Tumor和Normal样本中，对位点选择模块得到的169个位点的频率做差值，对差值的绝对值进行核密度估计：采用平滑的峰值函数来拟合观察到的数据点，从而对真实的概率密度分布曲线进行模拟，公式如式Ⅲ所示；本发明选择高斯函数进行数据拟合，高斯函数公式如式Ⅳ所示。
接下来，对拟合曲线的峰值进行提取：获得拟合曲线之后，即可获得拟合曲线各处的坐标，根据拟合曲线的纵坐标选取峰值，为了结果更准确，本发明选择大于前后四个位点纵坐标的位置为峰值，如图4所示，可得到3个峰值；
最后计算污染比例：3个峰值对应的横坐标中，最大的横坐标为污染比例p。在胚系突变中杂合突变的频率接近50%，而纯合突变的频率接近100%，纯合突变比杂合突变对样本的频率影响更大。如图4所示，横坐标是频率差的绝对值，频率差为0的第一个峰是没有被影响的位点，第二个峰为杂合突变影响的位点，第三个峰为纯合突变影响的位点。所以污染比例为横坐标较大的第3个峰所对应的横坐标。
4）推测污染源模块
首先根据对照样本的SNP位点的基因型和污染比例p推测污染源SNP位点的基因型，如图5所示，当对照样本的SNP位点的基因型为纯合时，若肿瘤样本SNP位点的突变频率为100%，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为100%-0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为100%-p，则污染源的基因型为野生型；
当对照样本的SNP位点的基因型为杂合时，若肿瘤样本SNP位点的突变频率为0.5×(100%-p)+p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为50%，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为50%-p，则污染源的基因型为野生型；
当对照样本的SNP位点的基因型为野生型时，若肿瘤样本SNP位点的突变频率为p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为0，则污染源的基因型为野生型。
接着以步骤2）中的公式Ⅰ计算推测的污染源中169个SNP位点基因型和同一批样本中其他Tumor样本中169个所选位点的基因型一致率，基因型一致率可以达到90%以上的Tumor样本，是该样本的污染源，否则无法确定。
实施例2
一、实验流程：
1. DNA提取
按照制造商的方案使用相关试剂盒从患者福尔马林固定石蜡包埋的组织切片样本（FFPE），配对的对照样本的血浆中提取基因组DNA。
2. 样本文库构建
DNA文库采用华大的MGIEasy FS DNA文库准备试剂盒制备。用Covaris S220仪器对DNA进行物理片段化，然后进行末端修复和3’端加A，adapter连接和PCR扩增。将Cot-1DNA阻断试剂、IDT通用阻断寡核苷酸和IDT接头特异性阻断寡核苷酸添加到混合文库中，并在真空离心浓缩仪中干燥。将干燥后的混合物再溶解在杂交缓冲液、杂交增强剂和捕获探针的混合液中。65℃杂交4 h后，用M270链亲和素珠捕获目标区域，65℃孵育45 min，然后用IDT xGen锁定试剂在65℃下洗涤3次，室温下洗涤3次。然后，进行15次捕获后扩增循环以获得捕获的文库。
3. 样本文库质检
文库构建完成后，先使用Qubit2.0进行初步定量，稀释文库至1ng/ul，随后使用Agilent 2100对文库的insert size进行检测，insert size符合预期后，使用Q-PCR方法对文库的有效浓度进行准确定量（文库有效浓度＞2nM），以保证文库质量。
4. 样本上机测序
库检合格后，把不同文库按照有效浓度及目标下机数据量的需求混合（pooling）后通过二代测序仪MGI2000进行上机测序，采用双端测序的测序方式，读长为100bp。
二、数据分析流程：
1、数据合并：对MGI2000测序仪产生的fastq文件按照样本名称将每条lane的同一样本的fastq文件进行合并，得到每个样本的fastq格式的文件。
2、数据质控和比对：用Fastp对上述fastq数据进行质量控制，去除接头、低质量的数据。利用Bwa将过滤后的fastq数据和人类hg19参考基因组序列进行比对，获得bam文件，然后使用GATK中的MarkDuplicates模块对PCR重复进行标记和去除，消除PCR-duplication的影响，得到去重后的bam文件。
3、重比对: 使用GATK的BaseRecalibrator和ApplyBQSR进行碱基质量校正，得到质量校正之后的bam文件。
4、变异检测和注释：利用上述得到的bam文件，首先提取panel的bed文件，随后使用HaplotypeCaller检测出突变位点，再使用GenotypeGVCFs对突变结果进行分型，接着采用SelectVariants分别提取SNP，INDEL和MNP的结果，用VariantFiltration过滤突变结果，最后使用annovar，snpEff和Vep三个软件对结果进行注释。
三、模型构建：
数据集1：以实验的方式模拟10例样本（两例不同个体的白细胞数据以 1:9，2:8，3:7，4:6，5:5，6:4，7:3，8:2，9:1，10:0 的比例，从DNA层面混合）。
数据集2：以生信的方式模拟112例样本（4对样本tumor数据两两混合，按照 1:9，2:8，3:7，4:6，5:5，6:4，7:3，8:2，9:1 的比例，从fastq层面混合）。
对112例样本，比较常用的Conpair，VerifyBamID和本发明实施例1系统（记为New_Method）的运行时间、运行资源、运行时的IO内存，各指标均值结果如表2所示：
表2 3种软件运行时间和运行资源比较
由表2可以看出，在运行时间和运行内存，以及IO读写方面，本发明提出的方法New_Method均优于其他两种方法。虽然其余两种方法花费的时间也不是特别长，但在IO读写方面，本发明提出的方法更好一些，如果同时运行多个样本的分析任务，大的IO内存可能会造成服务器卡顿。
分析方案1：10例混合后样本分别以两例不同个体的白细胞作为Normal进行分析。仅以其中Sample2作为Normal时为例，Tumor和Normal的对应关系如表3所示。
表3 以实验的方式模拟样本的Tumor和Normal的对应关系
样本配对判断结果如图6所示，mix_1_9的意思是Sample2样本为Normal，往Sample2里混入Sample1，使得Sample1:Sample2 = 1:9，以此类推。从图6可以看出，mix_9_1和mix_10_0的基因一致率没有达到60%，但是mix_9_1的频率相关系数达到了50%，因此只有mix_10_0，得出的结果是不能配对。mix_10_0是代表Sample2里全部混合Sample1，且Sample1和Sample2是两个不同的样本，所以mix_10_0的结果应该就是不能配对，因此配对检测的结果正确。
污染比例结果如图7所示，污染的比例是指在Sample2里混合多少Sampel1，由图7可知，对于分析方案1的10例样本，本方法均可以正确判断出是否被污染，以及计算出大致的污染比例。
污染源的溯源结果以mix_3_7为例，用mix_3_7预测的169个位点的基因型分别和Sample1、Sample2的169个位点的基因型作对比，计算基因型一致率，结果如表4所示：
表4 分析方案1污染源溯源结果示例
可能污染源 基因型一致的位点个数 总位点数 基因型一致率
Sample1 161 169 95.27%
Sample2 101 169 59.76%
mix_3_7表示，以Sample2样本为Normal，往Sample2里混入30%的Sample1。从表3可以看出，预测的mix_3_7样本169个位点的基因型，与Sample1的169个位点的基因型的一致率达到了95.27%，大于90%，而与Sample2的一致性小于90%，所以Sample1是mix_3_7样本的污染源，结果正确。
分析方案2：112例样本与各自的Normal进行配对分析。Tumor和Normal的对应关系如表5所示，_T代表是对应样本的Tumor样本，_N代表是对应样本的Normal样本。样本Sample2_T-0.1-Sample1_T-0.9_Sample2_N表示以Sample2_N作为Normal，往Sampel2_T里添加90%的Sample1_T，使得Sample2_T：Sample1_T=9：1。
表5 以生信的方式模拟112例样本的Tumor和Normal的对应关系
混合后样本 Tumor Normal
Sample1_T_Sample1_N Sample1_T Sample1_N
Sample2_T_Sample2_N Sample2_T Sample2_N
Sample3_T_Sample3_N Sample3_T Sample3_N
Sample4_T_Sample4_N Sample4_T Sample4_N
Sample2_T-0.1-Sample1_T-0.9_Sample2_N Sample2_T-0.1-Sample1_T-0.9 Sample2_N
Sample2_T-0.2-Sample1_T-0.8_Sample2_N Sample2_T-0.2-Sample1_T-0.8 Sample2_N
Sample2_T-0.3-Sample1_T-0.7_Sample2_N Sample2_T-0.3-Sample1_T-0.7 Sample2_N
Sample2_T-0.4-Sample1_T-0.6_Sample2_N Sample2_T-0.4-Sample1_T-0.6 Sample2_N
Sample2_T-0.5-Sample1_T-0.5_Sample2_N Sample2_T-0.5-Sample1_T-0.5 Sample2_N
Sample2_T-0.6-Sample1_T-0.4_Sample2_N Sample2_T-0.6-Sample1_T-0.4 Sample2_N
Sample2_T-0.7-Sample1_T-0.3_Sample2_N Sample2_T-0.7-Sample1_T-0.3 Sample2_N
Sample2_T-0.8-Sample1_T-0.2_Sample2_N Sample2_T-0.8-Sample1_T-0.2 Sample2_N
Sample2_T-0.9-Sample1_T-0.1_Sample2_N Sample2_T-0.9-Sample1_T-0.1 Sample2_N
Sample3_T-0.1-Sample1_T-0.9_Sample3_N Sample3_T-0.1-Sample1_T-0.9 Sample3_N
Sample3_T-0.1-Sample2_T-0.9_Sample3_N Sample3_T-0.1-Sample2_T-0.9 Sample3_N
Sample3_T-0.1-Sample4_T-0.9_Sample3_N Sample3_T-0.1-Sample4_T-0.9 Sample3_N
Sample3_T-0.2-Sample1_T-0.8_Sample3_N Sample3_T-0.2-Sample1_T-0.8 Sample3_N
Sample3_T-0.2-Sample2_T-0.8_Sample3_N Sample3_T-0.2-Sample2_T-0.8 Sample3_N
Sample3_T-0.2-Sample4_T-0.8_Sample3_N Sample3_T-0.2-Sample4_T-0.8 Sample3_N
Sample3_T-0.3-Sample1_T-0.7_Sample3_N Sample3_T-0.3-Sample1_T-0.7 Sample3_N
Sample3_T-0.3-Sample2_T-0.7_Sample3_N Sample3_T-0.3-Sample2_T-0.7 Sample3_N
Sample3_T-0.3-Sample4_T-0.7_Sample3_N Sample3_T-0.3-Sample4_T-0.7 Sample3_N
Sample3_T-0.4-Sample1_T-0.6_Sample3_N Sample3_T-0.4-Sample1_T-0.6 Sample3_N
Sample3_T-0.4-Sample2_T-0.6_Sample3_N Sample3_T-0.4-Sample2_T-0.6 Sample3_N
Sample3_T-0.4-Sample4_T-0.6_Sample3_N Sample3_T-0.4-Sample4_T-0.6 Sample3_N
Sample3_T-0.5-Sample1_T-0.5_Sample3_N Sample3_T-0.5-Sample1_T-0.5 Sample3_N
Sample3_T-0.5-Sample2_T-0.5_Sample3_N Sample3_T-0.5-Sample2_T-0.5 Sample3_N
Sample3_T-0.5-Sample4_T-0.5_Sample3_N Sample3_T-0.5-Sample4_T-0.5 Sample3_N
Sample3_T-0.6-Sample1_T-0.4_Sample3_N Sample3_T-0.6-Sample1_T-0.4 Sample3_N
Sample3_T-0.6-Sample2_T-0.4_Sample3_N Sample3_T-0.6-Sample2_T-0.4 Sample3_N
Sample3_T-0.6-Sample4_T-0.4_Sample3_N Sample3_T-0.6-Sample4_T-0.4 Sample3_N
Sample3_T-0.7-Sample1_T-0.3_Sample3_N Sample3_T-0.7-Sample1_T-0.3 Sample3_N
Sample3_T-0.7-Sample2_T-0.3_Sample3_N Sample3_T-0.7-Sample2_T-0.3 Sample3_N
Sample3_T-0.7-Sample4_T-0.3_Sample3_N Sample3_T-0.7-Sample4_T-0.3 Sample3_N
Sample3_T-0.8-Sample1_T-0.2_Sample3_N Sample3_T-0.8-Sample1_T-0.2 Sample3_N
Sample3_T-0.8-Sample2_T-0.2_Sample3_N Sample3_T-0.8-Sample2_T-0.2 Sample3_N
Sample3_T-0.8-Sample4_T-0.2_Sample3_N Sample3_T-0.8-Sample4_T-0.2 Sample3_N
Sample3_T-0.9-Sample1_T-0.1_Sample3_N Sample3_T-0.9-Sample1_T-0.1 Sample3_N
Sample3_T-0.9-Sample2_T-0.1_Sample3_N Sample3_T-0.9-Sample2_T-0.1 Sample3_N
Sample3_T-0.9-Sample4_T-0.1_Sample3_N Sample3_T-0.9-Sample4_T-0.1 Sample3_N
Sample4_T-0.1-Sample1_T-0.9_Sample4_N Sample4_T-0.1-Sample1_T-0.9 Sample4_N
Sample4_T-0.1-Sample2_T-0.9_Sample4_N Sample4_T-0.1-Sample2_T-0.9 Sample4_N
Sample4_T-0.2-Sample1_T-0.8_Sample4_N Sample4_T-0.2-Sample1_T-0.8 Sample4_N
Sample4_T-0.2-Sample2_T-0.8_Sample4_N Sample4_T-0.2-Sample2_T-0.8 Sample4_N
Sample4_T-0.3-Sample1_T-0.7_Sample4_N Sample4_T-0.3-Sample1_T-0.7 Sample4_N
Sample4_T-0.3-Sample2_T-0.7_Sample4_N Sample4_T-0.3-Sample2_T-0.7 Sample4_N
Sample4_T-0.4-Sample1_T-0.6_Sample4_N Sample4_T-0.4-Sample1_T-0.6 Sample4_N
Sample4_T-0.4-Sample2_T-0.6_Sample4_N Sample4_T-0.4-Sample2_T-0.6 Sample4_N
Sample4_T-0.5-Sample1_T-0.5_Sample4_N Sample4_T-0.5-Sample1_T-0.5 Sample4_N
Sample4_T-0.5-Sample2_T-0.5_Sample4_N Sample4_T-0.5-Sample2_T-0.5 Sample4_N
Sample4_T-0.6-Sample1_T-0.4_Sample4_N Sample4_T-0.6-Sample1_T-0.4 Sample4_N
Sample4_T-0.6-Sample2_T-0.4_Sample4_N Sample4_T-0.6-Sample2_T-0.4 Sample4_N
Sample4_T-0.7-Sample1_T-0.3_Sample4_N Sample4_T-0.7-Sample1_T-0.3 Sample4_N
Sample4_T-0.7-Sample2_T-0.3_Sample4_N Sample4_T-0.7-Sample2_T-0.3 Sample4_N
Sample4_T-0.8-Sample1_T-0.2_Sample4_N Sample4_T-0.8-Sample1_T-0.2 Sample4_N
Sample4_T-0.8-Sample2_T-0.2_Sample4_N Sample4_T-0.8-Sample2_T-0.2 Sample4_N
Sample4_T-0.9-Sample1_T-0.1_Sample4_N Sample4_T-0.9-Sample1_T-0.1 Sample4_N
Sample4_T-0.9-Sample2_T-0.1_Sample4_N Sample4_T-0.9-Sample2_T-0.1 Sample4_N
Sample2_T-0.1-Sample1_T-0.9_Sample1_N Sample2_T-0.1-Sample1_T-0.9 Sample1_N
Sample2_T-0.2-Sample1_T-0.8_Sample1_N Sample2_T-0.2-Sample1_T-0.8 Sample1_N
Sample2_T-0.3-Sample1_T-0.7_Sample1_N Sample2_T-0.3-Sample1_T-0.7 Sample1_N
Sample2_T-0.4-Sample1_T-0.6_Sample1_N Sample2_T-0.4-Sample1_T-0.6 Sample1_N
Sample2_T-0.5-Sample1_T-0.5_Sample1_N Sample2_T-0.5-Sample1_T-0.5 Sample1_N
Sample2_T-0.6-Sample1_T-0.4_Sample1_N Sample2_T-0.6-Sample1_T-0.4 Sample1_N
Sample2_T-0.7-Sample1_T-0.3_Sample1_N Sample2_T-0.7-Sample1_T-0.3 Sample1_N
Sample2_T-0.8-Sample1_T-0.2_Sample1_N Sample2_T-0.8-Sample1_T-0.2 Sample1_N
Sample2_T-0.9-Sample1_T-0.1_Sample1_N Sample2_T-0.9-Sample1_T-0.1 Sample1_N
Sample3_T-0.1-Sample1_T-0.9_Sample1_N Sample3_T-0.1-Sample1_T-0.9 Sample1_N
Sample3_T-0.1-Sample2_T-0.9_Sample2_N Sample3_T-0.1-Sample2_T-0.9 Sample2_N
Sample3_T-0.1-Sample4_T-0.9_Sample4_N Sample3_T-0.1-Sample4_T-0.9 Sample4_N
Sample3_T-0.2-Sample1_T-0.8_Sample1_N Sample3_T-0.2-Sample1_T-0.8 Sample1_N
Sample3_T-0.2-Sample2_T-0.8_Sample2_N Sample3_T-0.2-Sample2_T-0.8 Sample2_N
Sample3_T-0.2-Sample4_T-0.8_Sample4_N Sample3_T-0.2-Sample4_T-0.8 Sample4_N
Sample3_T-0.3-Sample1_T-0.7_Sample1_N Sample3_T-0.3-Sample1_T-0.7 Sample1_N
Sample3_T-0.3-Sample2_T-0.7_Sample2_N Sample3_T-0.3-Sample2_T-0.7 Sample2_N
Sample3_T-0.3-Sample4_T-0.7_Sample4_N Sample3_T-0.3-Sample4_T-0.7 Sample4_N
Sample3_T-0.4-Sample1_T-0.6_Sample1_N Sample3_T-0.4-Sample1_T-0.6 Sample1_N
Sample3_T-0.4-Sample2_T-0.6_Sample2_N Sample3_T-0.4-Sample2_T-0.6 Sample2_N
Sample3_T-0.4-Sample4_T-0.6_Sample4_N Sample3_T-0.4-Sample4_T-0.6 Sample4_N
Sample3_T-0.5-Sample1_T-0.5_Sample1_N Sample3_T-0.5-Sample1_T-0.5 Sample1_N
Sample3_T-0.5-Sample2_T-0.5_Sample2_N Sample3_T-0.5-Sample2_T-0.5 Sample2_N
Sample3_T-0.5-Sample4_T-0.5_Sample4_N Sample3_T-0.5-Sample4_T-0.5 Sample4_N
Sample3_T-0.6-Sample1_T-0.4_Sample1_N Sample3_T-0.6-Sample1_T-0.4 Sample1_N
Sample3_T-0.6-Sample2_T-0.4_Sample2_N Sample3_T-0.6-Sample2_T-0.4 Sample2_N
Sample3_T-0.6-Sample4_T-0.4_Sample4_N Sample3_T-0.6-Sample4_T-0.4 Sample4_N
Sample3_T-0.7-Sample1_T-0.3_Sample1_N Sample3_T-0.7-Sample1_T-0.3 Sample1_N
Sample3_T-0.7-Sample2_T-0.3_Sample2_N Sample3_T-0.7-Sample2_T-0.3 Sample2_N
Sample3_T-0.7-Sample4_T-0.3_Sample4_N Sample3_T-0.7-Sample4_T-0.3 Sample4_N
Sample3_T-0.8-Sample1_T-0.2_Sample1_N Sample3_T-0.8-Sample1_T-0.2 Sample1_N
Sample3_T-0.8-Sample2_T-0.2_Sample2_N Sample3_T-0.8-Sample2_T-0.2 Sample2_N
Sample3_T-0.8-Sample4_T-0.2_Sample4_N Sample3_T-0.8-Sample4_T-0.2 Sample4_N
Sample3_T-0.9-Sample1_T-0.1_Sample1_N Sample3_T-0.9-Sample1_T-0.1 Sample1_N
Sample3_T-0.9-Sample2_T-0.1_Sample2_N Sample3_T-0.9-Sample2_T-0.1 Sample2_N
Sample3_T-0.9-Sample4_T-0.1_Sample4_N Sample3_T-0.9-Sample4_T-0.1 Sample4_N
Sample4_T-0.1-Sample1_T-0.9_Sample1_N Sample4_T-0.1-Sample1_T-0.9 Sample1_N
Sample4_T-0.1-Sample2_T-0.9_Sample2_N Sample4_T-0.1-Sample2_T-0.9 Sample2_N
Sample4_T-0.2-Sample1_T-0.8_Sample1_N Sample4_T-0.2-Sample1_T-0.8 Sample1_N
Sample4_T-0.2-Sample2_T-0.8_Sample2_N Sample4_T-0.2-Sample2_T-0.8 Sample2_N
Sample4_T-0.3-Sample1_T-0.7_Sample1_N Sample4_T-0.3-Sample1_T-0.7 Sample1_N
Sample4_T-0.3-Sample2_T-0.7_Sample2_N Sample4_T-0.3-Sample2_T-0.7 Sample2_N
Sample4_T-0.4-Sample1_T-0.6_Sample1_N Sample4_T-0.4-Sample1_T-0.6 Sample1_N
Sample4_T-0.4-Sample2_T-0.6_Sample2_N Sample4_T-0.4-Sample2_T-0.6 Sample2_N
Sample4_T-0.5-Sample1_T-0.5_Sample1_N Sample4_T-0.5-Sample1_T-0.5 Sample1_N
Sample4_T-0.5-Sample2_T-0.5_Sample2_N Sample4_T-0.5-Sample2_T-0.5 Sample2_N
Sample4_T-0.6-Sample1_T-0.4_Sample1_N Sample4_T-0.6-Sample1_T-0.4 Sample1_N
Sample4_T-0.6-Sample2_T-0.4_Sample2_N Sample4_T-0.6-Sample2_T-0.4 Sample2_N
Sample4_T-0.7-Sample1_T-0.3_Sample1_N Sample4_T-0.7-Sample1_T-0.3 Sample1_N
Sample4_T-0.7-Sample2_T-0.3_Sample2_N Sample4_T-0.7-Sample2_T-0.3 Sample2_N
Sample4_T-0.8-Sample1_T-0.2_Sample1_N Sample4_T-0.8-Sample1_T-0.2 Sample1_N
Sample4_T-0.8-Sample2_T-0.2_Sample2_N Sample4_T-0.8-Sample2_T-0.2 Sample2_N
Sample4_T-0.9-Sample1_T-0.1_Sample1_N Sample4_T-0.9-Sample1_T-0.1 Sample1_N
Sample4_T-0.9-Sample2_T-0.1_Sample2_N Sample4_T-0.9-Sample2_T-0.1 Sample2_N
样本配对判断结果如图8所示，虽然有的样本基因型一致性的比例低于60%，但是所有样本的频率相关系数都大于50%，所有样本均可以配对，结果正确。
污染比例计算结果如图9所示。从图9可以看出，最前面4个Tumor样本和自身Normal样本分析时，没有检测出污染，结果正确；除了这4个，其他也均可以正确判断出被污染，且计算出了污染比例。混合比例10%的样本，计算出的污染比例和混合比例相差的相对较大，在5%以内，其余混合比例的样本，计算出的污染比例和混合比例相差在3%以内。
污染源溯源以Sample2_T -0.5-Sample1_T-0.5_Sample1_N为例，表示以Sample1_N（Sample1样本的Normal样本）为Normal样本，Tumor样本是往Sample1_T（Sample1样本的Tumor样本）里混合Sample2_T样本（Sample2样本的Tumor样本），使得两个样本的比例为5：5，结果如表6所示：
表6 分析方案2污染源溯源结果示例
可能污染源 基因型一致的位点个数 总位点数 基因型一致率
Sample1 100 169 59.17%
Sample2 158 169 93.49%
Sample3 96 169 56.80%
Sample4 95 169 56.21%
预测的Sample2_T -0.5-Sample1_T-0.5_Sample1_N样本，169个位点的基因型，与混合前的4个Tumor的169个位点的基因型相比，与Sample2的基因型一致率达到了93.49%，其余都不大于90%，所以Sample2_T -0.5-Sample1_T-0.5_Sample1_N样本的污染源是Sample2，结果正确。
根据分析方案1和分析方案2的结果，不论是实验数据还是模拟数据，本方法都可以正确判断Tumor和Normal两个样本是否正确配对，并且可以正确判断样本是否被污染，也可以计算出污染比例，且污染比例和预期相差不大。污染源的溯源结果也符合预期。综上可知，本方法可以实现预期目的，且结果较可靠。
尽管上述实施例对本发明做出了详尽的描述，但它仅仅是本发明一部分实施例，而不是全部实施例，人们还可以根据本实施例在不经创造性前提下获得其他实施例，这些实施例都属于本发明保护范围。

Claims (8)

1.一种基于胚系突变检测高通量测序中样本污染的方法，其特征在于，包括以下步骤；
获取待测样本的Fastp数据和待检测的SNP位点，所述待测样本包括肿瘤样本和对照样本；将所述Fastp数据处理为Vcf文件；
根据所述待检测的SNP位点的基因型一致率和频率相关系数判断所述肿瘤样本和对照样本是否配对，若基因型一致率＞60%，则肿瘤样本和对照样本可以配对；若基因型一致率≤60%，且频率相关性系数＞50%，则肿瘤样本和对照样本可以配对；若基因型一致率≤60%，且频率相关性系数≤50%，则肿瘤样本和对照样本不能配对，不进行后续处理；所述基因型一致率根据式Ⅰ计算得到，所述频率相关系数根据式Ⅱ计算得到；
式Ⅰ；
其中，Consistency rate是基因型一致率，S是肿瘤样本中SNP位点的基因型和对照样本中SNP位点的基因型一致的个数，N是筛选的SNP位点总位点数；
式Ⅱ；
其中，r(X,Y)是频率相关系数，X是肿瘤样本中每个SNP位点的突变频率，Y是对照样本中每个SNP位点的突变频率，Cov(X,Y)为X与Y的协方差，Var[X]为X的方差，Var[Y]为Y的方差；
采用式Ⅲ对所述SNP位点在能够配对的肿瘤样本和对照样本中突变频率差的绝对值进行核密度估计：
式Ⅲ；
其中，K为核函数，h为带宽，是概率密度，n是SNP位点个数，x是以xi为中心带宽内的所有点，xi是第i个SNP位点在肿瘤样本和对照样本中突变频率差的绝对值；
将所述核密度结果进行高斯函数进行数据拟合，高斯函数公式如下：
式Ⅳ；
其中，y是式Ⅲ计算得到的概率密度，a为高斯曲线的峰值，b为其对应的横坐标，c为标准差；
对拟合曲线中峰值对应的横坐标进行提取，最大的横坐标为污染比例p。
2.根据权利要求1所述的方法，其特征在于，所述方法还包括：根据对照样本的SNP位点的基因型、污染比例p和肿瘤样本SNP位点的突变频率确定污染源；
所述确定污染源包括：推测污染源SNP位点的基因型，同一批样本中其他肿瘤样本记为可能污染样本，根据所述污染源SNP位点的基因型和所述可能污染样本中SNP位点的基因型一致率确定污染源，基因型一致率≥90%的可能污染样本，为待测样本中肿瘤样本的污染源；
所述推测污染源SNP位点的基因型包括：
当对照样本的SNP位点的基因型为纯合时，若肿瘤样本SNP位点的突变频率为100%，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为100%-0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为100%-p，则污染源的基因型为野生型；
当对照样本的SNP位点的基因型为杂合时，若肿瘤样本SNP位点的突变频率为0.5×(100%-p)+p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为50%，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为50%-p，则污染源的基因型为野生型；
当对照样本的SNP位点的基因型为野生型时，若肿瘤样本SNP位点的突变频率为p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为0，则污染源的基因型为野生型。
3.根据权利要求1所述的方法，其特征在于，获取所述待检测的SNP位点的方法包括：
1）从千人基因组数据库中筛选包含在待检测的目标区域中且等位频率为0.1~0.9的SNP位点；
2）将步骤1）筛选的SNP位点进行第二次筛选；所述第二次筛选的规则包括：
a.将每条染色体上的SNP位点按位置顺序进行排序；
b.保留每条染色体上的第1个SNP位点，选择距离>第1个SNP位点100kb的SNP位点暂定为第2个SNP位点，若第2个SNP位点之后100kb以内没有SNP位点，则保留第2个SNP位点；若第2个SNP位点之后100kb以内有SNP位点，则保留等位频率最接近0.5的SNP位点为第2个SNP位点；
c.根据b的方法从保留的第2个SNP位点起筛选到每条染色体的最后一个SNP位点，保留的所有SNP位点为所述待检测的SNP位点。
4.根据权利要求3所述的方法，其特征在于，所述待检测的SNP位点如下：
。
5.根据权利要求1或2所述的方法，其特征在于，将所述Fastp数据处理为Vcf文件包括：用Bwa将Fastq数据和人类hg19参考基因组进行比对，得到bam文件；
对所述bam文件进行去重和碱基质量校正，再进行变异检测和变异注释，得到肿瘤样本和对照样本的Vcf文件。
6.根据权利要求1或2所述的方法，其特征在于，所述峰值为大于前后四个横坐标位点的纵坐标。
7.一种基于胚系突变检测高通量测序中样本污染的系统，其特征在于，包括数据获取模块、数据分析模块和样本污染检测模块；
所述数据获取模块用于获取待测样本的Fastp数据，并将Fastp数据传输给数据分析模块，所述待测样本包括肿瘤样本和对照样本；
所述数据分析模块用于将Fastp数据处理为Vcf文件，并将Vcf文件传输给样本污染检测模块；
所述样本污染检测模块用于检测污染源，包括位点选择模块、样本配对检测分析模块和污染比例计算模块；
所述位点选择模块用于选择待检测的SNP位点，包括：
1）从千人基因组数据库中筛选包含在待检测的目标区域中且等位频率为0.1~0.9的SNP位点；
2）将步骤1）筛选的SNP位点进行第二次筛选；所述第二次筛选的规则包括：
a.将每条染色体上的SNP位点按位置顺序进行排序；
b.保留每条染色体上的第1个SNP位点，选择距离>第1个SNP位点100kb的SNP位点暂定为第2个SNP位点，若第2个SNP位点之后100kb以内没有SNP位点，则保留第2个SNP位点；若第2个SNP位点之后100kb以内有SNP位点，则保留等位频率最接近0.5的SNP位点为第2个SNP位点；
c.根据b的方法从保留的第2个SNP位点起筛选到每条染色体的最后一个SNP位点，保留的所有SNP位点为所述待检测的SNP位点；
所述样本配对检测分析模块用于检测肿瘤样本和对照样本是否配对，包括根据所述位点选择模块选择的SNP位点的基因型一致率和频率相关系数判断，若基因型一致率＞60%，则肿瘤样本和对照样本可以配对；若基因型一致率≤60%，且频率相关性系数＞50%，则肿瘤样本和对照样本可以配对；若基因型一致率≤60%，且频率相关性系数≤50%，则肿瘤样本和对照样本不能配对，不进行后续处理；所述基因型一致率根据式Ⅰ计算得到，所述频率相关系数根据式Ⅱ计算得到；
式Ⅰ；
其中，Consistency rate是基因型一致率，S是肿瘤样本中SNP位点的基因型和对照样本中SNP位点的基因型一致的个数，N是筛选的SNP位点总位点数；
式Ⅱ；
其中，r(X,Y)是频率相关系数，X是肿瘤样本中每个SNP位点的突变频率，Y是对照样本中每个SNP位点的突变频率，Cov(X,Y)为X与Y的协方差，Var[X]为X的方差，Var[Y]为Y的方差；
所述污染比例计算模块用于将能配对的肿瘤样本和对照样本计算污染比例，包括采用式Ⅲ对所述位点选择模块筛选得到的SNP位点在肿瘤和对照样本中突变频率差的绝对值进行核密度估计：
式Ⅲ；
其中，K为核函数，h为带宽，是概率密度，n是SNP位点个数，x是以xi为中心带宽内的所有点，xi是第i个SNP位点在肿瘤样本和对照样本中突变频率差的绝对值；
将所述核密度结果进行高斯函数进行数据拟合，高斯函数公式如下：
式Ⅳ；
其中，y是式Ⅲ计算得到的概率密度，a为高斯曲线的峰值，b为其对应的横坐标，c为标准差；
对拟合曲线中峰值对应的横坐标进行提取，最大的横坐标为污染比例p。
8.根据权利要求7所述的系统，其特征在于，所述样本污染检测模块还包括推测污染源模块；所述推测污染源模块用于根据对照样本的SNP位点的基因型、污染比例p和肿瘤样本SNP位点的突变频率确定污染源，包括：
推测污染源SNP位点的基因型，同一批样本中其他肿瘤样本记为可能污染样本，根据所述污染源SNP位点的基因型和所述可能污染样本中SNP位点的基因型一致率确定污染源，基因型一致率≥90%的可能污染样本，为待测样本中肿瘤样本的污染源；所述推测污染源SNP位点的基因型包括：
当对照样本的SNP位点的基因型为纯合时，若肿瘤样本SNP位点的突变频率为100%，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为100%-0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为100%-p，则污染源的基因型为野生型；
当对照样本的SNP位点的基因型为杂合时，若肿瘤样本SNP位点的突变频率为0.5×(100%-p)+p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为50%，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为50%-p，则污染源的基因型为野生型；
当对照样本的SNP位点的基因型为野生型时，若肿瘤样本SNP位点的突变频率为p，则污染源的基因型为纯合；若肿瘤样本SNP位点的突变频率为0.5p，则污染源的基因型为杂合；若肿瘤样本SNP位点的突变频率为0，则污染源的基因型为野生型。
CN202311540335.8A 2023-11-20 2023-11-20 基于胚系突变检测高通量测序中样本污染的方法和系统 Active CN117253539B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311540335.8A CN117253539B (zh) 2023-11-20 2023-11-20 基于胚系突变检测高通量测序中样本污染的方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311540335.8A CN117253539B (zh) 2023-11-20 2023-11-20 基于胚系突变检测高通量测序中样本污染的方法和系统

Publications (2)

Publication Number Publication Date
CN117253539A CN117253539A (zh) 2023-12-19
CN117253539B true CN117253539B (zh) 2024-02-06

Family

ID=89137250

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311540335.8A Active CN117253539B (zh) 2023-11-20 2023-11-20 基于胚系突变检测高通量测序中样本污染的方法和系统

Country Status (1)

Country Link
CN (1) CN117253539B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117711487B (zh) * 2024-02-05 2024-05-17 广州嘉检医学检测有限公司 胚系SNV、InDel变异的鉴定方法、系统以及可读存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
SG141218A1 (en) * 2003-10-07 2008-04-28 Nanyang Polytechnic Method for prediction of single nucleotide polymorphisms
WO2018150378A1 (en) * 2017-02-17 2018-08-23 Grail, Inc. Detecting cross-contamination in sequencing data using regression techniques
CN110444255A (zh) * 2019-08-30 2019-11-12 深圳裕策生物科技有限公司 基于二代测序的生物信息质控方法、装置和存储介质
WO2022105629A1 (zh) * 2020-11-23 2022-05-27 福建和瑞基因科技有限公司 一种用于检测样本污染水平的snp位点的筛选方法及样本污染水平的检测方法
CN115394357A (zh) * 2022-09-01 2022-11-25 杭州链康医学检验实验室有限公司 用于判断样本配对或污染的位点组合及其筛选方法和应用

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2016065075A1 (en) * 2014-10-21 2016-04-28 uBiome, Inc. Method and system for microbiome-derived diagnostics and therapeutics

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
SG141218A1 (en) * 2003-10-07 2008-04-28 Nanyang Polytechnic Method for prediction of single nucleotide polymorphisms
WO2018150378A1 (en) * 2017-02-17 2018-08-23 Grail, Inc. Detecting cross-contamination in sequencing data using regression techniques
CN110444255A (zh) * 2019-08-30 2019-11-12 深圳裕策生物科技有限公司 基于二代测序的生物信息质控方法、装置和存储介质
WO2022105629A1 (zh) * 2020-11-23 2022-05-27 福建和瑞基因科技有限公司 一种用于检测样本污染水平的snp位点的筛选方法及样本污染水平的检测方法
CN115394357A (zh) * 2022-09-01 2022-11-25 杭州链康医学检验实验室有限公司 用于判断样本配对或污染的位点组合及其筛选方法和应用
CN116805510A (zh) * 2022-09-01 2023-09-26 杭州链康医学检验实验室有限公司 用于判断样本配对或污染的位点组合及其应用

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
EGFR、ALK、KIT及KRAS基因Panel体细胞突变高通量测序检测法室内质量控制的建立;吴小延;李丹;杨鑫华;刘小云;龙亚康;王芳;邓玲;;肿瘤预防与治疗(第10期);10-18 *
SeqSQC: A Bioconductor Package for Evaluating the Sample Quality of Next-generation Sequencing Data;Qian Liu等;Genomics, Proteomics & Bioinformatics;第17卷(第2期);211-218 *

Also Published As

Publication number Publication date
CN117253539A (zh) 2023-12-19

Similar Documents

Publication Publication Date Title
EP3143537B1 (en) Rare variant calls in ultra-deep sequencing
US12018329B2 (en) Diagnosing fetal chromosomal aneuploidy using massively parallel genomic sequencing
CN117253539B (zh) 基于胚系突变检测高通量测序中样本污染的方法和系统
CN108256289B (zh) 一种基于目标区域捕获测序基因组拷贝数变异的方法
US20190338349A1 (en) Methods and systems for high fidelity sequencing
CN113249495A (zh) 一种绵羊液相芯片及其应用
WO2022105629A1 (zh) 一种用于检测样本污染水平的snp位点的筛选方法及样本污染水平的检测方法
CN112746097A (zh) 一种检测样本交叉污染的方法以及预测交叉污染源的方法
CN113584178A (zh) 一种无创亲子检测分析方法和装置
CN115083521B (zh) 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及系统
CN107075565B (zh) 个体单核苷酸多态性位点分型方法及装置
CN115083529A (zh) 一种检测样本污染率的方法及装置
CN116348957A (zh) 检测测序数据中的交叉污染
CN109920480B (zh) 一种校正高通量测序数据的方法和装置
CN109461473B (zh) 胎儿游离dna浓度获取方法和装置
CN115948521B (zh) 一种检测非整倍体缺失染色体信息的方法
WO2019132010A1 (ja) 塩基配列における塩基種を推定する方法、装置及びプログラム
JP2022537442A (ja) ヒト胚におけるコピー数変異を検証するために単一ヌクレオチド変異の密度を使用するシステム、コンピュータプログラム製品及び方法
CN115961054B (zh) 用于华南虎个体识别和/或亲子鉴定的遗传标记及应用
WO2019213810A1 (zh) 检测染色体非整倍性的方法、装置及系统
CN113981070B (zh) 胚胎染色体微缺失的检测方法、装置、设备和存储介质
Jiang et al. DRAMS: A tool to detect and re-align mixed-up samples for integrative studies of multi-omics data
CN110942806A (zh) 一种血型基因分型方法和装置及存储介质
CN116888274A (zh) 一种利用多态性位点和靶位点测序检测胎儿遗传变异的方法
CN108504734B (zh) 一种恶性肿瘤组织特定个体归属的判断方法及其应用

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