肿瘤TMB的计算原理和数学模型
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了肿瘤TMB的计算原理和数学模型相关的知识，希望对你有一定的参考价值。
参考技术A美国FDA目前批准了3个检测肿瘤TMB（ tumor mutation burden）的方法：Omisc Core WES、MSK-IMPACT、FMI。如果按实验流程及操作复杂度来讲，这可能是有史以来最复杂的病理诊断方法。流程一般如下：手术取样、提取DNA、建测序文库、靶向捕获、NGS测序、数据比对、变异检测、变异注释、结果解读。复杂方法下的简单目的：数一数肿瘤特有的变异有几个，多数情况下平均每100万个位点还不超过10个变异。
TMB定义是：每Mb区域里非同义的体细胞变异个数。
有效覆盖区域L容易知道，如果覆盖度足够好，可以用bed文件区间长度代替，那么计算TMB主要的工作就是怎么计算somatic。在肿瘤基因检测中计算somatic处于最核心的位置，靶向治疗、免疫治疗等都围绕着somatic展开。
Somatic定义为个体生长过程中产生的变异。目标就是要找在测序中肿瘤组织测到而正常组织没有测到的变异。因为测序有错误，如果测序足够深，肿瘤及正常组织大多数位点都会有变异。那么通过判断正常组织有无变异将不适合。我们可以模型化为肿瘤组织与正常组织某一变异有相同的频率为种系变异germline，否则为somatic。方法是对肿瘤组织及正常组织同时外显子靶向高通量测序，对每个位点进行比较。通过测序深度、频率等参数比较判断是否somatic。业内有很多计算软件可以用来计算somatic。下表是业内使用最多的开源软件[1]，现在选取业界使用最多的VarScan2 的算法来介绍somatic及其TMB的计算原理。
表1 热门软件采用的数学模型
以下是VarScan2输出VCF格式文件位点的结果
SPV=8.8337E-2是用Fisher’s Exact Test计算的P-value值，用来判断是否somatic。那这值是怎么计算得到的，源文献没有具体说明，本人在这里给出计算过程。Fisher’s Exact Test接受2×2行列表作为计算对象，从上述结果提取数据可以得到下表。
表2. 肿瘤样本与正常样本覆盖深度2×2行列表。RD: 支持参考序列的碱基数；AD: 支持变异的碱基数。
Fisher’s Exact Test接受2×2行列表作为输入，并按如下方式计算P-value值：
1、计算行列表中每一行，每一列的总和以及观察总数。
2、给定行和列总和，如果原假设为真，则使用超几何概率函数计算条件概率，以观察行列表中的准确结果。条件概率为
其中R1和R2是行总和，C1和C2是列总和，N是行列表中观测的总数，nij是表中第i行和第j列的值。
3、查找与行和列之和一致的所有可能的非负整数矩阵。对于每个矩阵，使用P公式计算相关的条件概率。
4、根据感兴趣的替代假设，使用这些值来计算检验的p值。
(1)对于双边测试，对于观察到的行列表，将所有小于或等于Pcutoff的条件概率求和。
(2)对于左侧测试，将（1,1）象限频率小于或等于n11的所有矩阵的条件概率求和。
(3)对于右侧测试，在观察到的行列表中求和（1,1）单元频率大于或等于n11的所有矩阵的条件概率。
VarScan2 somaitc计算的是单侧测试，如果按下表排列数据是左侧测试。
表2. 肿瘤样本与正常样本覆盖深度2×2行列表。
那怎么计算P-value呢？根据P-value的定义：在原假设为真的前提下，出现该样本或比该样本更极端的结果的概率之和。上表就是所有tumor Ref depth小于等于47的概率之和。
(P根据公式1计算)
P-value=0.0883，这P的意思是假设tumor Alt与normal Alt有同样的频率，但这的假设可能性为0.0883。按通常0.05或0.01阈值来cutoff的话，应该接受这一假设。肿瘤样本中测到的位点是种系变异而非somatic。
再看一个测序位点
表3.肿瘤样本与正常样本覆盖深度2×2行列表。
P-value=1.3226e-15，这P值就足够小，认为是种系变异就不合适了，可以判断为肿瘤特有的变异somatic。
按照上述方法对肿瘤样本与正常样本所有同时覆盖的位点都计算一个P-value值，当P-value值小于设定值时判定位肿瘤体细胞变异。因影响NGS测序结果的因素多，reads覆盖的概率分布并不跟理想模型一致。如果按照常见0.05值来判定somatic，结果会比预期多。可以通过更低的P-value值或其他条件过滤掉可能的假阳性位点。按本人计算经验，P-value设0.001跟mutect2设的默认参数差不多。如果somatic是10个，覆盖到的区域2M，那么TMB=10/2=5。
参考资料
肿瘤增长数学模型
曲面偏微分方程的一个实际应用——肿瘤增长
文章目录
简介
基本想法是利用表面有限元数值方法求解演化曲面上的反应扩散方程。不断增长的生物表面上的模式形成，是通过求解曲面上的反应扩散方程得到的，可以参考图灵的经典文章 1。这种方程表现出空间均匀结构的扩散驱动不稳定性，导致空间不均匀的图案。这项工作的进一步的发展，有了更多的数学模型，其中的一些应用包括九头蛇的图案形成、脊椎动物的连续图案形成、动物的皮肤表面的颜色标记、蝴蝶翅膀上的色素沉着图案等等。可以参考综述文章 2。
图灵的原始工作没有考虑几何的增长和变化。然而，这些因素在自然界中所观察到的模式的发展中很重要。增长的区域上的模式形成已被广泛研究。例如 Crampin 等人 3 考虑一维区域的区域增长，并表明它可能是一种机制，用于增加模式形成时，相对于初始条件的稳健性。在了解在一种名叫 Pomacanthus 的鱼中观察到的条纹进化的背景下4，二维增长平面区域被考虑5。也可以看看6，其中图灵扩散驱动的不稳定性条件被推广到具有缓慢、各向同性增长的反应扩散系统。本文中提出的数值方法是一种在演化表面上探索此类反应动力学的工具。
许多有趣的工作都是关于增长平面区域的。然而，在生物的应用中，很自然地考虑在不断发展的空间三维区域的表面上形成图案，例如我们参考7，和8等，其中研究了曲率
对固定球体和生长锥体的图案形成的影响。另一项初步研究是在这篇文章 9 中，它提出了一种数值方法，用于解决球形表面图案形成问题，并应用于实体肿瘤生长。他们的数值方法基于参数化球体上的表面，然后使用线方法。
下面提到的方法基于表面有限元方法，它是有限元方法6的自然延伸，并且能够处理复杂的几何形状和形状101112131415。这个想法是对表面进行三角剖分，并使用基于三角剖分的分段线性表面有限元空间来逼近偏微分方程组。在这种情况下，三角剖分的顶点以规定的或受某些进化规律支配的速度移动。为了做到这一点，我们需要在表面给出一个适当的守恒定律。当以适当的变分形式给出时，演化表面有限元方法利用了该守恒定律的特殊特征。
接下来，我们要介绍的内容陈述如下。
首先，我们给出演化表面上的反应扩散方程。在这个框架内，我们给出了表面梯度的符号表达，这是演化表面上的第一原理反应扩散系统推导出来的关键。我们考虑经过充分研究的活化剂耗尽衬底模型161718，它也被称为 Brusselator 模型。同样，在此框架中可以轻松处理任何其他合理的反应动力学。我们还讨论了表面生长模型的建模。
其次，我们看看一个演化表面有限元方法，该方法应用于改变形状的演化表面上的反应扩散系统。变分公式的一个特点是不会出现某些几何量，例如曲面的平均曲率和法线。然后通过有限元方法以自然的方式利用这一点。在这种方法中，在导出反应扩散系统或利用表面有限元方法时不需要转换或参数化。在空间中离散会产生一个非自治常微分方程组，该方程组在时间上是离散的。特别是因为这个应用涉及半线性抛物方程，我们线性化非线性动力系统[^马兹瓦缪斯]，以获得线性代数系统。然后使用 GMRes 算法19求解。
最后，在固定和进化的表面上呈现出图灵模式图案。我们展示了表面有限元方法的普遍适用性，该方法具有处理生物模型可能导致的各向异性生长的能力。结果证实，结合区域增长增强了模式选择过程。该方法也适用于将表面演化与表面上的反应扩散系统耦合的模型。这在再下一个章节中有说明，它应用于模拟实体瘤的生长。
部分结论性意见在最后一个章节。其他也考虑表面偏微分方程数值解的文章的例子包括2021222324。这篇文章25 对计算演化表面的不同方法的回顾。
演化表面上反应扩散方程的推导
曲面梯度
我们假设
Γ
\\Gamma
Γ 可以全局表示为函数
d
d
d 的零水平集，即存在开集
U
⊂
R
3
\\mathcal{U} \\subset \\mathbb{R}^{3}
U⊂R3和函数
d
∈
C
2
(
U
)
d \\in C^{2}(\\mathcal{U})
d∈C2(U)，使得
U
∩
Γ
=
{
x
∈
U
:
d
(
x
)
=
0
}
,
and
∇
d
(
x
)
≠
0
,
∀
x
∈
U
∩
Γ
\\mathcal{U} \\cap \\Gamma=\\{x \\in \\mathcal{U}: d(x)=0\\}, \\quad \\text { and } \\quad \\nabla d(x) \\neq 0, \\quad \\forall x \\in \\mathcal{U} \\cap \\Gamma
U∩Γ={x∈U:d(x)=0}, and ∇d(x)=0,∀x∈U∩Γ
我们定义单位法向量场为，
v
(
x
)
=
∇
d
(
x
)
∣
∇
d
(
x
)
∣
\\boldsymbol{v}(x)=\\frac{\\nabla d(x)}{|\\nabla d(x)|}
v(x)=∣∇d(x)∣∇d(x)
我们定义切向梯度为，
∇
Γ
η
(
x
)
=
∇
η
(
x
)
−
∇
η
(
x
)
⋅
v
(
x
)
v
(
x
)
,
x
∈
Γ
\\nabla_{\\Gamma} \\eta(x)=\\nabla \\eta(x)-\\nabla \\eta(x) \\cdot \\boldsymbol{v}(x) \\boldsymbol{v}(x), \\quad x \\in \\Gamma
∇Γη(x)=∇η(x)−∇η(x)⋅v(x)v(x),x∈Γ
它只依赖于
η
\\eta
η 在曲面上的值。它具有三个分量，表示为：
∇
Γ
η
(
x
)
=
(
D
‾
1
η
(
x
)
,
D
‾
2
η
(
x
)
,
D
‾
3
η
(
x
)
)
\\nabla_{\\Gamma} \\eta(x)=\\left(\\underline{D}_{1} \\eta(x), \\underline{D}_{2} \\eta(x), \\underline{D}_{3} \\eta(x)\\right)
∇Γη(x)=(D1η(x),D2η(x),D3η(x))
Laplace-Beltrami 算子就定义为切向梯度的切向散度：
Δ
Γ
η
(
x
)
=
∇
Γ
⋅
∇
Γ
η
(
x
)
=
∑
i
=
1
3
D
‾
i
D
‾
i
η
(
x
)
\\Delta_{\\Gamma} \\eta(x)=\\nabla_{\\Gamma} \\cdot \\nabla_{\\Gamma} \\eta(x)=\\sum_{i=1}^{3} \\underline{D}_{i} \\underline{D}_{i} \\eta(x)
ΔΓη(x)=∇Γ⋅∇Γη(x)=i=1∑3DiDiη(x)
演化表面上的反应扩散系统
让
Γ
(
t
)
\\Gamma(t)
Γ(t) 是嵌入到三维空间
Ω
(
t
)
\\Omega(t)
Ω(t) 的二维曲面。曲面上每一点的速度表示为，
v
:
=
V
v
+
v
T
\\boldsymbol{v}:=V \\boldsymbol{v}+\\boldsymbol{v}_{T}
v:=Vv+vT
u = { u i } i = 1 m \\boldsymbol{u}=\\left\\{u_{i}\\right\\}_{i=1}^{m} u={ui}i=1m 是个矢量场， R ( t ) \\mathcal{R}(t) R(t) 是曲面的某一部分，由物质守恒定律，
d d t ∫ R ( t ) u i = − ∫ ∂ R ( t ) q i ⋅ μ + ∫ R ( t ) f i ( u ) \\frac{d}{d t} \\int_{\\mathcal{R}(t)} u_{i}=-\\int_{\\partial \\mathcal{R}(t)} \\boldsymbol{q}_{i} \\cdot \\boldsymbol{\\mu}+\\int_{\\mathcal{R}(t)} f_{i}(\\boldsymbol{u}) dtd∫R(t)ui=−∫∂R(t)qi⋅μ+∫R(t)fi(u)
这里边的 q i \\boldsymbol{q}_{i} qi 表示曲面通量， f i ( u ) f_{i}(\\boldsymbol{u}) fi(u) 表示曲面内的净生产率。
d d t ∫ R ( t ) u i = − ∫ ∂ R ( t ) q i ⋅ μ + ∫ R ( t ) f i ( u ) \\frac{d}{d t} \\int_{\\mathcal{R}(t)} u_{i}=-\\int_{\\partial \\mathcal{R}(t)} \\boldsymbol{q}_{i} \\cdot \\boldsymbol{\\mu}+\\int_{\\mathcal{R}(t)} f_{i}(\\boldsymbol{u}) dtd∫R(t)ui=−∫∂R(t)qi⋅μ+∫R(t)f以上是关于肿瘤TMB的计算原理和数学模型的主要内容，如果未能解决你的问题，请参考以下文章