CN115083529B - 一种检测样本污染率的方法及装置 - Google Patents

一种检测样本污染率的方法及装置 Download PDF

Info

Publication number
CN115083529B
CN115083529B CN202210811098.3A CN202210811098A CN115083529B CN 115083529 B CN115083529 B CN 115083529B CN 202210811098 A CN202210811098 A CN 202210811098A CN 115083529 B CN115083529 B CN 115083529B
Authority
CN
China
Prior art keywords
site
rate
sample
pollution
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
CN202210811098.3A
Other languages
English (en)
Other versions
CN115083529A (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.)
Shenzhen Guiinga Medical Laboratory
Beijing Jiyinjia Medical Laboratory Co ltd
Original Assignee
Shenzhen Guiinga Medical Laboratory
Beijing Jiyinjia 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 Shenzhen Guiinga Medical Laboratory, Beijing Jiyinjia Medical Laboratory Co ltd filed Critical Shenzhen Guiinga Medical Laboratory
Priority to CN202210811098.3A priority Critical patent/CN115083529B/zh
Publication of CN115083529A publication Critical patent/CN115083529A/zh
Application granted granted Critical
Publication of CN115083529B publication Critical patent/CN115083529B/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
    • G16B50/00—ICT programming tools or database systems specially adapted for bioinformatics
    • G16B50/30—Data warehousing; Computing architectures
    • 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/30—Detection of binding sites or motifs

Landscapes

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

Abstract

一种检测样本污染率的方法及装置，该方法包括：位点MAF提取步骤，包括提取待测样本的测序数据中的位点在数据库中的MAF；过滤步骤，包括过滤去除不符合条件的SNP位点；错误率计算步骤，包括计算不同碱基替换的错误率；似然值计算步骤，包括计算待测样本在不同污染率下的似然值；候选污染率计算步骤，包括根据每个SNP位点计算的似然值对数与位点深度计算加权平均值，选择加权平均值最大的似然值对应的污染率为候选污染率；优化步骤，包括根据优化函数优化候选污染率，获得最终的样本污染率。该方法的分析结果可信度高。

Description

一种检测样本污染率的方法及装置
技术领域
本发明涉及生物信息学领域，具体涉及一种检测样本污染率的方法及装置。
背景技术
DNA甲基化是最早发现的基因表观修饰方式之一，可能存在于所有高等生物中，能够在不改变基因序列的前提下，改变遗传现象。它是基因调控的手段之一，即通过对位于启动子及第一外显子区的CpG岛的甲基化而抑制基因的表达，对生命活动非常重要。相比于普通的测序技术，甲基化测序会导致DNA上面的碱基信息改变。主流的甲基化测序方法是使用重亚硫酸盐处理，会导致非甲基化的C碱基变换成T。而新的甲基化测序方法，TET酶和吡啶硼烷结合处理的方法(TAPS)会导致甲基化的C碱基变换成T。
目前，甲基化测序在肿瘤基因组中的应用越来越多，因此对于甲基化测序污染率的探索极其重要。甲基化测序样本污染主要存在三种情况，即个体内、不同个体间以及跨物种间的污染。跨物种间的污染比较容易解决，因为可以通过评估样本比对到物种参考基因组的情况，推测污染率的占比。但是同一物种不同个体间的污染则难以发现，因为污染可能发生在许多意想不到的场合，比如：样品存储过程中，收集样本的容器被污染；样品运输过程中，容器密封不严导致样本外溢；以及实验室人为操作制备过程中，不同样本移液时忘记更换枪尖或未使用带滤芯枪尖。
在现有的甲基化测序过程中，检测和计算样本间污染率对于甲基化测序的下游分析非常重要，即使只是少量的污染，也会导致分析结果出现许多假阳性或者假阴性，特别是在肿瘤与正常甲基化测序样本的比对研究中。因此，甲基化测序数据的样本间污染需要严格控制，但是遗憾的是，目前并没有相关软件或者流程能够实现对甲基化测序数据的污染率评估。如何实现对甲基化测序数据的污染率评估是目前亟待解决的问题。
发明内容
根据第一方面，在一实施例中，提供一种检测样本污染率的方法，包括：
位点MAF提取步骤，包括提取待测样本的测序数据中的SNP位点在数据库中的最小等位基因频率(MAF)；
过滤步骤，包括过滤去除不符合条件的SNP位点；
错误率计算步骤，包括计算不同碱基替换的错误率；
似然值计算步骤，包括计算待测样本在不同污染率下的似然值；
候选污染率计算步骤，包括根据每个SNP位点计算的似然值对数与位点深度计算加权平均值，选择加权平均值最大的似然值对应的污染率为候选污染率；
优化步骤，包括根据优化函数优化所述候选污染率，获得最终的样本污染率。
根据第二方面，在一实施例中，提供一种检测样本污染率的装置，包括：
位点MAF提取模块，用于提取待测样本的测序数据中的位点在数据库中的最小等位基因频率(MAF)；
过滤模块，用于过滤去除不符合条件的位点；
错误率计算模块，用于计算不同碱基替换的错误率；
似然值计算模块，用于计算待测样本在不同污染率下的似然值；
候选污染率计算模块，用于根据每个SNP位点计算的似然值对数与位点深度计算加权平均值，选择加权平均值最大的似然值对应的污染率为候选污染率；
优化模块，用于根据优化函数优化所述候选污染率，获得最终的样本污染率。
根据第三方面，在一实施例中，提供一种检测样本污染率的装置，包括：
存储器，用于存储程序；
处理器，用于通过执行所述存储器存储的程序以实现第一方面的方法。
根据第四方面，在一实施例中，提供一种计算机可读存储介质，所述介质上存储有程序，所述程序能够被处理器执行以实现第一方面的方法。
依据上述实施例的一种检测样本污染率的方法及装置，该方法的分析结果可信度高。
在一实施例中，该方法在实际应用过程中，可通过总体评估甲基化测序样本在不同污染率下的评估污染率集合，在总体情况下，确定污染率阈值，能进一步增加下游分析结果的可靠性。
具体实施方式
下面通过具体实施方式对本发明作进一步详细说明。其中不同实施方式中类似元件采用了相关联的类似的元件标号。在以下的实施方式中，很多细节描述是为了使得本申请能被更好的理解。然而，本领域技术人员可以毫不费力的认识到，其中部分特征在不同情况下是可以省略的，或者可以由其他元件、材料、方法所替代。在某些情况下，本申请相关的一些操作并没有在说明书中显示或者描述，这是为了避免本申请的核心部分被过多的描述所淹没，而对于本领域技术人员而言，详细描述这些相关操作并不是必要的，他们根据说明书中的描述以及本领域的一般技术知识即可完整了解相关操作。
另外，说明书中所描述的特点、操作或者特征可以以任意适当的方式结合形成各种实施方式。同时，方法描述中的各步骤或者动作也可以按照本领域技术人员所能显而易见的方式进行顺序调换或调整。因此，说明书中的各种顺序只是为了清楚描述某一个实施例，并不意味着是必须的顺序，除非另有说明其中某个顺序是必须遵循的。
本文中为部件所编序号本身，例如“第一”、“第二”等，仅用于区分所描述的对象，不具有任何顺序或技术含义。
如本文所用，dbSNP是指NCBI于1998年建立的主要存储单核苷酸多态性(SNP)的免费公共数据库。该数据库包含多种模式生物。虽然其名称为dbSNP，但该数据库实际上包括多种分子变异，具体如下：
单核苷酸多态性SNP；
短缺失和插入多态性short deletion and insertion polymorphisms(indels/DIPs)；
微卫星标记或短串联重复microsatellite markers or short tandem repeats(STRs)；
多核苷酸多态性multinucleotide polymorphisms(MNPs)；
杂合序列heterozygous sequences；
命名变体named variants。
如本文所用，“MAF”是指最小等位基因频率，通常是指在给定人群中的不常见的等位基因发生频率，例如TT、TC、CC三个基因型，在人群中C的频率＝0.36，T的频率＝0.64，则等位基因C就为最小等位基因频率，MAF＝0.36。
如本文所用，“杂合子”是指同一位点上的两个等位基因不相同的基因型。
如本文所用，“纯合子”是指同一位点上的两个等位基因相同的基因型。
根据第一方面，在一实施例中，提供一种检测样本污染率的方法，包括：
位点MAF提取步骤，包括提取待测样本的测序数据中的SNP位点在数据库中的最小等位基因频率(MAF)；
过滤步骤，包括过滤去除不符合条件的SNP位点；
错误率计算步骤，包括计算不同碱基替换的错误率；
似然值计算步骤，包括计算待测样本在不同污染率下的似然值；
候选污染率计算步骤，包括根据每个SNP位点计算的似然值对数与位点深度计算加权平均值，选择加权平均值最大的似然值对应的污染率为候选污染率；
优化步骤，包括根据优化函数优化所述候选污染率，获得最终的样本污染率。
在一实施例中，所述位点MAF提取步骤中，将SNP位点映射到数据库的文件中，如果该位点存在，则保留该位点，并提取数据库中该位点的等位基因信息；如果该位点不存在，则删除该位点。
在一实施例中，所述位点MAF提取步骤中，如果该位点存在，则提取其在数据库中的等位基因信息以及人群频率。
在一实施例中，所述位点MAF提取步骤中，如果该位点存在，则提取其在数据库中的最小等位基因频率(MAF)。
在一实施例中，所述位点MAF提取步骤中，所述数据库包括但不限于dbSNP数据库、HapMap数据库中的至少一种。
在一实施例中，所述过滤步骤中，包括过滤并确定SNP位点深度、基因型、位点先验污染概率和背景噪音读段(read)数。
在一实施例中，所述过滤步骤中，如果SNP位点的等位基因不是T(胸腺嘧啶)、C(胞嘧啶)、G(鸟嘌呤)、A(腺嘌呤)或者N(未知碱基)，则过滤去除该位点。
在一实施例中，所述过滤步骤中，计算该位点主要等位基因的读段数占该位点主要等位基因的读段数与该位点次要等位基因的读段数的加和的比例，记为主要占比；计算该位点的次要等位基因的读段数占该位点主要等位基因的读段数与该位点次要等位基因的读段数的加和的比例，记为次要占比；根据主/次要占比判断样本基因型。
在一实施例中，所述过滤步骤中，如果主要占比小于杂合子位点的限制阈值，则将相应基因型判定为“1/1”；如果次要占比小于杂合子位点的限制阈值，则将相应基因型判定为“0/0”；其他则为“0/1”。
在一实施例中，所述限制阈值可以为0.25。
在一实施例中，所述过滤步骤中，如果基因型为“0/0”，背景噪音read计数为次要占比乘以位点深度，该位点的先验污染率为1-[1-(MAF)]2；
如果基因型为“1/1”，背景噪音read计数为主要占比乘以位点深度，该位点的先验污染率为1-(MAF)2；
如果基因型为“0/1”，背景噪音read计数为主要占比、次要占比中额最小值乘以位点深度，该位点的先验污染率为1-[2×(MAF)×(1-MAF)]。
在一实施例中，所述错误率计算步骤中，对待测样本的测序数据中纯合子基因型进行分析，计算所述纯合子基因型中主等位基因的碱基为X、次等位基因不为Y的读段数占参考碱基为X的读段数的2/3的比值，即为X>Y的错误率，X与Y为互补配对的碱基。例如，X为A时，Y为T；再例如，X为C时，Y为G；再例如，X为T时，Y为A。
在一实施例中，所述似然值计算步骤中，根据给定的初始污染率集合，计算不同初始污染率下每个SNP位点的对数似然值的比值，即用该位点的污染的似然值除以非污染的似然值(即为SNP位点的污染似然值与非污染似然值的比值)，再取对数。
在一实施例中，所述似然值计算步骤中，给定样本初始污染率集合，根据所述样本初始污染率集合中的每个值分别计算待测样本的污染率的概率值。
在一实施例中，所述似然值计算步骤中，首先移除杂合子SNP位点；使用以下公式(1)分别计算每一个SNP位点的污染似然值：
Figure BDA0003739064110000041
式(1)中，μ为位点先验污染率，ε为该位点的碱基替换错误率，α为样本污染率，n为位点深度，k为位点背景噪音read计数；
再使用公式(2)计算每一个SNP位点非污染似然值，最后将SNP位点的污染似然值除以非污染似然值，再取对数，即为SNP位点的对数似然比值；
公式(2)如下：
Figure BDA0003739064110000042
公式(2)中各符号的含义同公式(1)。
在一实施例中，所述对数以自然数e为底。
在一实施例中，所述候选污染率计算步骤中，取所有SNP位点的对数似然比值与位点深度，计算加权平均值，则每个给定的初始污染率均会计算出对应的污染率概率值，选取其中最大的污染率概率值对应的污染率作为候选污染率。
在一实施例中，所述优化步骤中，选取所述候选污染率，将所述候选污染率前后各一位的初始污染率作为优化函数的取值范围，代入公式(1)，并根据公式(2)，计算得到对数似然比值，计算全局最优解，即为待测样本最终的污染率。
在一实施例中，所述待测样本的测序数据包括但不限于全基因组测序数据、全外显子组测序数据、靶向捕获测序数据。
在一实施例中，所述待测样本的测序数据包括甲基化测序数据。
在一实施例中，所述待测样本的测序数据为甲基化测序数据时，则从所述测序数据中提取未转化的SNP位点的碱基信息，并提取所述SNP位点的最小等位基因频率。
在一实施例中，“未转化”是指未进行甲基化转化，没有提供甲基化信息的碱基。甲基化转化方法包括TET辅助吡啶硼烷的甲基化测序法(TET-assisted pyridine boranesequencing)，该方法简称TAPS，无需亚硫酸氢盐，利用TET(Ten-Eleven Translocation，简称TET)酶将5mC(5-甲基胞嘧啶)和5hmC(5-羟甲基胞嘧啶)氧化为5caC(5-羧甲基胞嘧啶)，随后使用有机硼烷(包括但不限于吡啶硼烷、2-甲基吡啶硼烷等等)将5caC还原为二氢尿嘧啶(dihydrouracil，DHU)，而后的PCR再将DHU转化为胸腺嘧啶(T)，对目标序列直接进行DNA甲基化测序。本发明中，在进行污染率检测时，先通过中国专利《一种基于甲基化测序数据进行变异检测的方法及装置》(申请号：202110960293.8，公开号：CN113674802A)中的方法，将转化后的碱基修正为转化之前的碱基，具体是将T修正为C。
在一实施例中，甲基化测序方法不受限制，现有的亚硫酸氢盐测序法等其他甲基化测序方法获得的测序数据也适用于本发明。
根据第二方面，在一实施例中，提供一种检测样本污染率的装置，包括：
位点MAF提取模块，用于提取待测样本的测序数据中的SNP位点在dbSNP中的最小等位基因频率(MAF)；
过滤模块，用于过滤去除不符合条件的SNP位点；
错误率计算模块，用于计算不同碱基替换的错误率；
似然值计算模块，用于计算待测样本在不同污染率下的似然值；
候选污染率计算模块，用于根据每个SNP位点计算的似然值对数与位点深度计算加权平均值，选择加权平均值最大的似然值对应的污染率为候选污染率；
优化模块，用于根据优化函数优化所述候选污染率，获得最终的样本污染率。
根据第三方面，在一实施例中，提供一种检测样本污染率的装置，包括：
存储器，用于存储程序；
处理器，用于通过执行所述存储器存储的程序以实现第一方面任意一项的方法。
根据第四方面，在一实施例中，提供一种计算机可读存储介质，所述介质上存储有程序，所述程序能够被处理器执行以实现第一方面任意一项的方法。
在一实施例中，提供一种基于甲基化测序数据进行样本污染率检测的方法，包括：快速提取甲基化测序样本中的未转化位点的碱基信息；确定未转化位点在dbSNP中的最小等位基因频率；过滤并确定SNP位点深度、基因型、位点污染先验概率和背景噪音read数；计算不同碱基替换的错误率；计算样本不同污染率下的似然值；选取最大似然值；根据优化函数优化最大似然对应的污染率。该方法可以应用于全基因组甲基化测序数据。
在一实施例中，提供一种基于甲基化测序数据进行样本污染率检测的方法，包括如下步骤：
(1)快速提取甲基化测序样本中的未转化位点的碱基信息。使用软件处理甲基化测序的压缩比BAM文件，该软件能够检测甲基化测序样本的SNP位点，并且能够将SNP位点的碱基区分为经过转化或者未经过转化的碱基(具体方法参见申请号为202110960293.8的中国专利《一种基于甲基化测序数据进行变异检测的方法及装置》说明书第105～109段)，对于经过甲基化转化的碱基，修正为转化之前的碱基，即为未转化的碱基，通过提取未转化的碱基信息作为污染率检测软件的输入文件，评估样本污染率。
(2)确定未转化位点在dbSNP中的最小等位基因频率(MAF)。将未转化的位点映射到dbSNP数据库的vcf文件中，如果该位点存在，则保留该位点并提取dbSNP中该位点的MAF；如果该位点不存在，则删除该位点。
(3)过滤并确定SNP位点深度、基因型、位点先验污染概率和背景噪音read数。如果SNP位点的等位基因不是T、C、G、A或者N，则过滤去除该位点。
(4)计算不同碱基替换的错误率。每种碱基均有3种碱基替换形式，比如A碱基，可能的碱基替换形式有：A>T，A>G，A>C。例如要计算碱基替换A>T的错误率，仅需对样本中纯合子基因型进行分析，计算这部分基因型中主等位基因为A、次等位基因不为T的reads占参考碱基为A的reads的比值，即为A>T的错误率。
(5)计算样本在不同污染率下的似然值。根据给定的初始污染率集合，计算不同初始污染率下每个SNP位点的对数似然值的比值，即用该位点的污染的似然值比去非污染的似然值，再取对数。
(6)选取最大似然值。将每个SNP位点计算的似然值对数与位点深度计算加权平均值，选择加权平均值最大的似然值对应的污染率为候选污染率。
(7)根据优化函数优化最大似然值对应的污染率。取候选污染率前后各一位污染率为优化函数的取值参数，再求全局最优解作为最终的样本对应的污染率。
实施例1
下述实施例中涉及的肿瘤样本由北京吉因加医学检验实验室有限公司提供。样本进行TAPS测序，下机数据预处理和基因组比对(使用BWA软件)得到压缩比对文件BAM。BAM文件为本方法的输入文件。
检测方法如下：
(1)读取BAM文件，使用变异检测软件(具体方法参见申请号为202110960293.8的中国专利《一种基于甲基化测序数据进行变异检测的方法及装置》说明书第105～109段)处理，该软件能够检测甲基化测序样本的SNP位点，并且能够将SNP位点的碱基区分为经过转化或者未经过转化的碱基，通过提取未转化的碱基信息作为污染率软件的输入文件。该文件每一行为一个SNP位点上未经过转化的碱基信息，具体包括：染色体、位点位置、位点深度、参考碱基、A、T、G、C、N碱基个数、INS个数、DEL个数。
(2)将上述文件SNP位点映射到dbSNP数据库中。dbSNP数据库是单核苷酸多态性数据库，收录了SNP、短插入缺失多态性和短重复序列等数据。将映射到的SNP位点保留，并提取其在dbSNP数据库中的等位基因信息以及人群频率；在dbSNP数据库中未收录的位点，则删除该位点。
(3)过滤SNP位点，如果SNP位点的等位基因都不是单碱基或者N，则将这种类型的SNP位点过滤掉。计算该位点主要等位基因的读段数(即reads个数)占该位点主要等位基因的读段数与该位点次要等位基因的读段数的加和的比例，记为主要占比；计算该位点的次要等位基因的读段数占该位点主要等位基因的读段数与该位点次要等位基因的读段数的加和的比例，记为次要占比。根据主/次要占比判断样本基因型。如果主要占比小于杂合子位点的限制阈值(0.25)，则相应基因型判定为“1/1”；如果次要占比小于杂合子位点的限制阈值(0.25)，则相应基因型判定为“0/0”；其他则为“0/1”。如果基因型为“0/0”，背景噪音read计数为次要占比乘以位点深度，该位点的先验污染率为1-[1-(MAF)]2；如果基因型为“1/1”，背景噪音read计数为主要占比乘以位点深度，该位点的先验污染率为1-(MAF)2；如果基因型为“0/1”，背景噪音read计数为主要占比、次要占比中的最小值乘以位点深度，该位点的先验污染率为1-[2×(MAF)×(1-MAF)]。例如，等位基因A/T，A的计数是3，T的计数是2，那么A为主要等位基因，B为次要等位基因。主要占比＝A count/(A count+B count)＝3/(3+2)＝0.6，次要占比＝0.4。因此与0.25比较，两者都大于0.25，因此基因型为0/1。
(4)计算不同碱基替换的错误率。因为要计算碱基替换的错误率，所以仅需考虑基因型是纯合子的位点。假设计算A>T的错误率，错误率＝(样本中主等位基因是A、次等位基因不是T的read数和)÷(样本中参考碱基为A的全部read数×2/3)。这里之所以分母要乘以2/3，是因为计算A>T的错误率会和A>C或者A>G重复。其他的碱基替换方法相同。
(5)给定样本初始污染率集合，具体为0.001、0.003、0.005、0.007、0.01、0.03、0.05、0.07、0.1、0.2、0.3。根据污染率集合中的每个值分别计算样本的污染率的概率值。首先移除杂合子SNP位点；然后使用以下公式(1)分别计算每一个SNP位点的污染似然值：
Figure BDA0003739064110000071
其中μ为位点先验污染率，ε为该位点的碱基替换错误率，α为样本污染率，n为位点深度，k为位点背景噪音read计数。
再使用公式(2)计算每一个SNP位点非污染似然值。最后将SNP位点的污染似然值比去非污染似然值，再取对数(以自然数e为底)，即为SNP位点的对数似然比值。
公式(2)如下：
Figure BDA0003739064110000072
公式(2)中各符号的含义同公式(1)。
(6)将上述所有SNP位点的对数似然比值与位点深度计算加权平均值，则相应的11个初始污染率均会计算出对应的污染率概率值，选取其中最大的污染率概率值对应的污染率作为候选污染率。
(7)挑选上述候选污染率。因为该污染率还不够精确，真实的污染率可能在候选污染率的前后区间内，因此将其前后各一位的初始污染率作为优化函数的取值范围，代入上述公式(1)，并根据公式(2)继续计算对数似然比值，计算全局最优解，即作为样本最终的污染率。
本实施例对6例样本进行干实验模拟污染率情况，具体是将其他甲基化测序样本按一定比例混到当前样本中，其中每一例样本都分别模拟了0.001、0.005、0.01、0.05、0.1、0.2、0.3污染率的情况。按上述所示方法，对模拟样本进行污染率检测，具体检测结果如表1所示。
表1六例模拟污染率样本的污染率检测结果
污染率 Mix1 Mix2 Mix3 Mix4 Mix5 Mix6 误差
0.001 0.0014 0.0015 0.0016 0.0014 0.0015 0.0014 0.047％
0.005 0.0051 0.0052 0.0052 0.0050 0.0053 0.0048 0.017％
0.010 0.0100 0.0099 0.0099 0.0095 0.0100 0.0092 0.029％
0.050 0.0565 0.0587 0.0528 0.0532 0.0551 0.0524 0.477％
0.100 0.1291 0.1289 0.1237 0.1243 0.1302 0.1189 2.585％
0.200 0.2310 0.2273 0.2311 0.2275 0.2339 0.2235 2.904％
0.300 0.2802 0.2760 0.2828 0.2778 0.2845 0.2758 2.048％
通过上述结果可知，当模拟污染率为0.001时，6个样本计算的污染率与模拟污染率的平均误差仅为0.047％；
当模拟污染率为0.005时，6个样本计算的污染率与模拟污染率的平均误差仅为0.017％；
当模拟污染率为0.01时，6个样本计算的污染率与模拟污染率的平均误差仅为0.029％；
当模拟污染率为0.05时，6个样本计算的污染率与模拟污染率的平均误差仅为0.477％；
当模拟污染率为0.1时，6个样本计算的污染率与模拟污染率的平均误差仅为2.585％；
当模拟污染率为0.2时，6个样本计算的污染率与模拟污染率的平均误差仅为2.904％；
当模拟污染率为0.3时，6个样本计算的污染率与模拟污染率的平均误差仅为2.048％。
7种模拟污染率的平均误差仅为1.158％。显然，本方法的分析结果可信度高，并且在实际应用过程中，可通过总体评估甲基化测序样本在不同污染率下的评估污染率集合，在总体情况下，确定污染率阈值，能进一步增加下游分析结果的可靠性。
在一实施例中，本发明可应用于甲基化测序样本中。
本领域技术人员可以理解，上述实施方式中各种方法的全部或部分功能可以通过硬件的方式实现，也可以通过计算机程序的方式实现。当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘、光盘、硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述实施方式中全部或部分功能。
以上应用了具体个例对本发明进行阐述，只是用于帮助理解本发明，并不用以限制本发明。对于本发明所属技术领域的技术人员，依据本发明的思想，还可以做出若干简单推演、变形或替换。

Claims (13)

1.一种检测样本污染率的方法，其特征在于，包括：
位点MAF提取步骤，包括提取待测样本的测序数据中的SNP位点在数据库中的最小等位基因频率；
过滤步骤，包括过滤去除不符合条件的SNP位点；
错误率计算步骤，包括计算不同碱基替换的错误率；
似然值计算步骤，包括计算待测样本在不同污染率下的似然值；
候选污染率计算步骤，包括根据每个SNP位点计算的似然值对数与位点深度计算加权平均值，选择加权平均值最大的似然值对应的污染率为候选污染率；
优化步骤，包括根据优化函数优化所述候选污染率，获得最终的样本污染率。
2.如权利要求1所述的方法，其特征在于，所述位点MAF提取步骤中，将SNP位点映射到数据库的文件中，如果该位点存在，则保留该位点，并提取数据库中该位点的等位基因信息；如果该位点不存在，则删除该位点；
或，所述位点MAF提取步骤中，如果该位点存在，则提取其在数据库中的等位基因信息以及人群频率；
或，所述位点MAF提取步骤中，如果该位点存在，则提取其在数据库中的最小等位基因频率；
或，所述位点MAF提取步骤中，所述数据库包括dbSNP数据库、HapMap数据库中的至少一种。
3.如权利要求1所述的方法，其特征在于，所述过滤步骤中，包括过滤并确定SNP位点深度、基因型、位点先验污染概率和背景噪音读段数；
或，所述过滤步骤中，如果SNP位点的等位基因不是T、C、G、A或者N，则过滤去除该位点；
或，所述过滤步骤中，计算该位点主要等位基因的读段数占该位点主要等位基因的读段数与该位点次要等位基因的读段数的加和的比例，记为主要占比；计算该位点的次要等位基因的读段数占该位点主要等位基因的读段数与该位点次要等位基因的读段数的加和的比例，记为次要占比；根据所述主/次要占比判断样本基因型；
或，如果主要占比小于杂合子位点的限制阈值，则将相应基因型判定为“1/1”；如果次要占比小于杂合子位点的限制阈值，则将相应基因型判定为“0/0”；其他则为“0/1”；
或，所述限制阈值为0.25。
4.如权利要求3所述的方法，其特征在于，所述过滤步骤中，如果基因型为“0/0”，背景噪音读段计数为次要占比乘以位点深度，该位点的先验污染率为1-[1-(MAF)]2；
如果基因型为“1/1”，背景噪音读段计数为主要占比乘以位点深度，该位点的先验污染率为1-(MAF)2；
如果基因型为“0/1”，背景噪音读段计数为主要占比、次要占比中的最小值乘以位点深度，该位点的先验污染率为1-[2×(MAF)×(1-MAF)]。
5.如权利要求1所述的方法，其特征在于，所述错误率计算步骤中，对待测样本的测序数据中纯合子基因型进行分析，计算所述纯合子基因型中主等位基因的碱基为X、次等位基因的对应位置的碱基不为Y的读段数占待测样本中参考碱基为X的读段数的2/3的比值，即为X>Y的错误率，X与Y为互补配对的碱基。
6.如权利要求1所述的方法，其特征在于，所述似然值计算步骤中，根据给定的初始污染率集合，计算不同初始污染率下每个SNP位点的对数似然值的比值，即用该位点的污染的似然值除以非污染的似然值，再取对数；
或，所述似然值计算步骤中，给定样本初始污染率集合，根据所述样本初始污染率集合中的每个值分别计算待测样本的污染率的概率值。
7.如权利要求1所述的方法，其特征在于，所述似然值计算步骤中，首先移除杂合子SNP位点；然后使用以下公式(1)分别计算每一个SNP位点的污染似然值：
Figure QLYQS_1
(1)
公式（1）中，μ为位点先验污染率，ε为该位点的碱基替换错误率，α为样本污染率，n为位点深度，k为位点背景噪音读段计数；
再使用公式（2）计算每一个SNP位点非污染似然值，最后将SNP位点的污染似然值除以非污染似然值，再取对数，即为SNP位点的对数似然比值；
公式（2）如下：
Figure QLYQS_2
（2）
公式（2）中各符号的含义同公式（1）。
8.如权利要求7所述的方法，其特征在于，所述对数以自然数e为底。
9.如权利要求7所述的方法，其特征在于，
所述优化步骤中，选取所述候选污染率，将所述候选污染率前后各一位的初始污染率作为优化函数的取值范围，代入公式（1），并根据公式（2），计算得到对数似然比值，计算全局最优解，即为待测样本最终的污染率。
10.如权利要求1所述的方法，其特征在于，所述待测样本的测序数据包括全基因组测序数据、全外显子组测序数据；
或，所述待测样本的测序数据包括甲基化测序数据；
或，所述甲基化测序包括TET辅助吡啶硼烷的甲基化测序；
或，所述待测样本的测序数据为甲基化测序数据时，则从所述测序数据中提取未转化的SNP位点的碱基信息，并提取所述SNP位点的最小等位基因频率。
11.一种检测样本污染率的装置，其特征在于，包括：
位点MAF提取模块，用于提取待测样本的测序数据中的SNP位点在数据库中的最小等位基因频率；
过滤模块，用于过滤去除不符合条件的SNP位点；
错误率计算模块，用于计算不同碱基替换的错误率；
似然值计算模块，用于计算待测样本在不同污染率下的似然值；
候选污染率计算模块，用于根据每个SNP位点计算的似然值对数与位点深度计算加权平均值，选择加权平均值最大的似然值对应的污染率为候选污染率；
优化模块，用于根据优化函数优化所述候选污染率，获得最终的样本污染率。
12.一种检测样本污染率的装置，其特征在于，包括：
存储器，用于存储程序；
处理器，用于通过执行所述存储器存储的程序以实现如权利要求1~10任意一项所述的方法。
13.一种计算机可读存储介质，其特征在于，所述介质上存储有程序，所述程序能够被处理器执行以实现如权利要求1~10任意一项所述的方法。
CN202210811098.3A 2022-07-11 2022-07-11 一种检测样本污染率的方法及装置 Active CN115083529B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210811098.3A CN115083529B (zh) 2022-07-11 2022-07-11 一种检测样本污染率的方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210811098.3A CN115083529B (zh) 2022-07-11 2022-07-11 一种检测样本污染率的方法及装置

Publications (2)

Publication Number Publication Date
CN115083529A CN115083529A (zh) 2022-09-20
CN115083529B true CN115083529B (zh) 2023-03-14

Family

ID=83259782

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210811098.3A Active CN115083529B (zh) 2022-07-11 2022-07-11 一种检测样本污染率的方法及装置

Country Status (1)

Country Link
CN (1) CN115083529B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115410649B (zh) * 2022-04-01 2023-03-28 北京吉因加医学检验实验室有限公司 一种同时检测甲基化和突变信息的方法及装置
CN116153400B (zh) * 2022-12-20 2023-11-21 深圳吉因加信息科技有限公司 一种用于检测同源污染的模型构建方法与装置

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20180373832A1 (en) * 2017-06-27 2018-12-27 Grail, Inc. Detecting cross-contamination in sequencing data
CN113136422A (zh) * 2020-01-19 2021-07-20 北京圣谷同创科技发展有限公司 通过成组snp位点检测高通量测序样本污染的方法
CN111370065B (zh) * 2020-03-26 2022-10-04 北京吉因加医学检验实验室有限公司 一种检测rna跨样本交叉污染率的方法和装置
CN114530198A (zh) * 2020-11-23 2022-05-24 福建和瑞基因科技有限公司 一种用于检测样本污染水平的snp位点的筛选方法及样本污染水平的检测方法
CN112746097A (zh) * 2021-01-29 2021-05-04 深圳裕康医学检验实验室 一种检测样本交叉污染的方法以及预测交叉污染源的方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
James E.Barrett etc.."The DNA methylome of cervical cells can predict the presence of ovarian cancer".2022,全文. *
Yidong Chen etc.."DNA methylome reveals cellular origin of cell-free DNA in spent medium of human preimplantation embryos".2021,全文. *

Also Published As

Publication number Publication date
CN115083529A (zh) 2022-09-20

Similar Documents

Publication Publication Date Title
CN115083529B (zh) 一种检测样本污染率的方法及装置
US12002544B2 (en) Determining progress of chromosomal aberrations over time
AU2022203114A1 (en) Detecting mutations for cancer screening and fetal analysis
JP7113838B2 (ja) 配列バリアントコールのための有効化方法およびシステム
CN102682224B (zh) 检测拷贝数变异的方法和装置
US20220130488A1 (en) Methods for detecting copy-number variations in next-generation sequencing
Chen et al. Using Mendelian inheritance to improve high-throughput SNP discovery
Katzman et al. GC-biased evolution near human accelerated regions
US20220336051A1 (en) Method for Determining Relatedness of Genomic Samples Using Partial Sequence Information
WO2018150378A1 (en) Detecting cross-contamination in sequencing data using regression techniques
AU2017279575B2 (en) Detection of genetic or molecular aberrations associated with cancer
CN112639129A (zh) 确定新发突变在胚胎中的遗传状态的方法和装置
WO2022061189A1 (en) Detecting cross-contamination in sequencing data
CN110462063B (zh) 一种基于测序数据的变异检测方法、装置和存储介质
Swaegers et al. Restricted X chromosome introgression and support for Haldane's rule in hybridizing damselflies
CN112735519B (zh) 一种定位偏分离性状的方法、装置及存储介质
이선호 New Methods for SNV/InDel Calling and Haplotyping from Next Generation Sequencing Data
Hu Methods and Analyses in the Study of Human DNA Methylation
CN114708905A (zh) 基于ngs的染色体非整倍体检测方法、装置、介质和设备
Wagner et al. Computational Analysis of Whole-Genome Differential Allelic Expression Data in
Borucki et al. Viral shift and drift, report part 1

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