[发明专利]一种基于高通量测序的基因变异检测方法

权利要求书

1.一种基于高通量测序的基因变异检测方法，包括如下步骤：

(A)对高通量测序数据进行预处理，主要包括碱基质量校正、read去重；

(B)计算测序产生的偏好性，通过偏好性减少支持突变的有效read的数量，获得校正后的支持突变型的read数量；

(C)利用校正后的支持突变型的read数量，计算所有单个样本相关质量，然后得出单个样本变异质量，通过变异质量进行判定；

步骤(C)按照如下进行：

(C1)针对每个需要计算quality的位点，采用如下公式计算family quality，single-strand quality，double-strand quality；

Quality(p，B，b，r)＝logr(BIAS(B，p×B，B，b)×(r-1)+1)×10×log₁₀(OR(B，p×B，B，b))；

其中，Bias(B，p×B，B，b)＝max(0，OR(LS(B，p×B，B，b))-1)*(b÷(b+B))+1；

OR代表比值比；OR的计算公式如下：

OR(B，p×B，B，b)＝(B×b)÷(p×B×B)；

LS代表Laplace平滑；LS公式如下：

LS(B，p×B，B，b)＝Laplace₂(Laplace₁(B，p×B，B，b))；

Laplace1公式为：Laplace₁(B，p×B，B，b)＝(B+p′，p×B+p′，B+p′，b+p′)；

Laplace2公式为：Laplace₂(B，p×B，B，b)＝(B+p′×α，p×B+p′×α，B+p′×β，b+p′×β)；

在Laplace₁公式和Laplace₂公式中，p′的默认值是0.5；

在Laplace₂公式中，α和β公式如下：

α＝max(0，p′×((B+p×B)/(B+b)-1))；

β＝max(0，p′×(B+b)/(B+p×B)-1)；

其中：

(1)针对一个family里的reads，让r＝2并且使用以下定义计算family quality：

p：步骤(A)reads去重中计算出的碱基质量的阈值所对应的错误概率；

B：从数据实际观测到支持任何等位基因的read数量；

b：从数据实际观测到支持突变等位基因的read数量；

(2)如果高通量测序实验采用PCR扩增方法建库，则r＝1.5，如果高通量测序实验采用捕获或者WGS方法建库，则r＝1.25；比对到一条染色体链上的families都在一条single-strand上；针对一条single-strand，使用以下定义计算single-strand quality：

p：每个可能的family quality所对应的错误概率；

B：实际观测到支持任何等位基因的family数量；

b：实际观测到支持突变等位基因的family数量；

如果single-strand quality大于最大family quality数值，则让single-strandquality等于最大family quality数值；所述最大family quality如下：对于C：G＞T：A突变形式是44，对于T：A＞C：G突变形式是48，对于其他SNV突变形式是52，对于所有InDel突变形式是60；

在尝试所有的p之后，取最大single-strand quality作为最终single-strandquality；

(3)每一条single-strand在一条double-strand上的正向或反向；针对一条single-strand，使用以下定义计算double-strandquality：

v₁′＝v₁+v₂×min(1，(min(w₁+w₂，w₀)-w₁)/w₁)；

其中v₁′是某一条链上的最终变异质量，v₁是这一条链的变异质量，即最终single-strand quality，v₂是另外一条链上的变异质量，w₁是这一条链的family quality的平均值，w₂是另外一条链的family quality的平均值，w₀等于最大family quality加上10，所述最大family quality如下：对于C：G＞T：A突变形式是44，对于T：A＞C：G突变形式是48，对于其他SNV突变形式是52，对于所有InDel突变形式是60；

交换v₁和v₂并且交换w₁和w₂，根据计算v₁′公式的对称性计算v₂′；

v₂′＝v₂+v₁×min(1，(min(w₁+w₂，w₀)-w₂)/w₂)；

(C2)把double-strand quality定义成单个样本变异质量；

(C3)按照上述(C1)和(C2)，分别计算得到tumor样本和normal对照样本的单样本突变质量；然后计算TLOD；TLOD代表变异既不来源于躯体也不来源于胚胎的概率，如下：

首先，计算以下函数的输出值；

TNreward(a₁，a₂，b₁，b₂)＝maxH(TNrewardBinomial(a₁，a₂，b₁，b₂)，PLQTN(b₂/b₁，a₂/a₁))；

其中，TNrewardBinomial(a₁，a₂，b₁，b₂)＝10/log(10)×(a₁×KLBernoulli(b₂/b₁，a₂/a₁))；其中，KLBernoulli(b₂/b₁，a₂/a₁)＝(a₂/a₁)×log((a₂/a₁)/(b₂/b₁))+(1-(a₂/a₁))×log((1-(a₂/a₁))/(1-(b₂/b₁)))；如果b₂/b₁＝0、a₂/a₁＝0、b₂/b₁是不定式或者a₂/a₁是不定式，那么KLBernoulli(b₂/b₁，a₂/a₁)＝0；如果max(a，-a)＜max(b，-b)，那么maxH(a，b)＝a，否则maxH(a，b)＝b；

PLQTN(t，n)＝3×10×log₁₀(min(t/n，10^2.5/3))；

TNpenal(a₁，a₂，b₁，b₂)＝max(0，min(DSVQnormal，PLQ(a₂/a₁))-12.5×(max(0，OR(a₁，a₂，b₁，b₂)-1))²)；其中DSVQnormal指的是normal样本的单样本突变质量，也就是normal样本的double-strand quality；OR代表比值比，OR的计算公式如下：OR(a₁，a₂，b₁，b₂)＝(a₁×b₂)÷(a₂×b₁)；a₁、a₂、b₁和b₂含义如下：

a₁：配对normal对照样本中支持任何等位基因的reads的family quality总和；

a₂：配对normal对照样本中支持突变等位基因的reads的family quality总和；

b₁：tumor样本中支持任何等位基因的reads的family quality总和；

b₂：tumor样本中支持突变等位基因的reads的family quality总和；

然后，计算TLOD，计算公式如下：

TLOD(a₁，a₂，b₁，b₂)＝min(DSVQtumor，PLQ(b₂/b₁))+TNreward(a₁，a₂，b₁，b₂)-TNpenal(a₁，a₂，b₁，b₂)；

其中，PLQ(f)＝90+3×10×log₁₀(f)；DSVQtumor指的是tumor样本的单样本突变质量，也就是tumor样本的double-strand quality；

(C4)计算NLOD；NLOD代表来变异来源于胚胎的概率，如下：

首先，针对每个样本，计算每个genotype的genotype likelihood；genotype简称GT，genotype likelihood简称GL；纯合野生型、杂合型和纯合突变型这三种GT的GL的计算方法如下：

纯合野生型的GL计算公式：GL(f,homref)＝max(GLPowerLaw(f,homref),GLBinomial(f,homref))；

杂合型的GL计算公式：GL(f,hetero)＝max(GLPowerLaw(f,hetero),GLBinomial(f,hetero))；

纯合突变型的GL计算公式：GL(f,homalt)＝max(GLPowerLaw(f,homalt),GLBinomial(f,homalt))；

其中，GLPowerLaw(f,homref)＝10×3×min(0,log₁₀(0.02×(1-f)/f))；

GLPowerLaw(f,hetero)＝10×3×min(0,log₁₀(min(f/(1-f),(1-f)/f)))；

GLPowerLaw(f,homalt)＝10×3×min(0,log₁₀(0.02×f/(1-f)))；

GLBinomial(f,homref)＝10/log(10)×d×KLBernoulli(min(0.02,f),f)；

GLBinomial(f,hetero)＝10/log(10)×d×KLBernoulli(0.5,f)；

GLBinomial(f,homalt)＝10/log(10)×d×KLBernoulli(max(1–0.02,f),f)；

KLBernoulli(x,y)＝y×log(y/x)+(1-y)×log((1-y)/(1-x))；如果x＝0、y＝0、x是不定式或者y是不定式，那么KLBernoulli(x,y)＝0；

并且a指的是支持突变型基因的read数量，d指的是支持任何基因型的read数量，f＝a/d是突变频率；

然后，计算NLOD；NLOD＝max(Qtumorhomref,Qnormalhomref)；

其中，在tumor样本上，Qtumorhomref＝-3-max(GL(f,homalt),GL(f,hetero))；

其中在配对的normal样本上，Qnormalhomref＝GL(f,homref)-max(GL(f,homalt),GL(f,hetero))；

(C5)根据如下公式计算最终变异质量：min(TLOD,NLOD+31)；

如果最终变异质量至少是60，则判定为阳性并且输出变异形式；否则判定为阴性。