CN113707216A - 一种浸润免疫细胞比例计数方法 - Google Patents

Abstract

本发明提供一种浸润免疫细胞比例计数方法，属于生物信息学与计算生物学技术领域。所述方法包括：将基因表达谱数据转化为免疫细胞的标记基因的线性组合，得到线性组合方程；基于得到的线性组合方程，构建用于反卷积求解免疫细胞比例的支持向量回归模型；输入训练集，利用粒子群算法对构建的支持向量回归模型进行迭代优化，其中，所述训练集中的每个样本包含个体抽检组织中不同基因的表达谱数据；利用优化后的支持向量回归模型，求解免疫细胞比例。采用本发明，能够显著降低未知混合物含量对免疫细胞比例计数，提高免疫细胞比例计数的准确度。

Description

一种浸润免疫细胞比例计数方法

技术领域

本发明涉及生物信息学与计算生物学技术领域，特别是指一种浸润免疫细胞比例计数方法。

背景技术

近年来，细菌、病毒等病原体入侵人体的时候会诱发免疫系统发生防御作用，这就使得生物学家去研究免疫细胞学这一新的方向。免疫组化法、单细胞测序技术和流式细胞术等传统生物方法可以计算人体组织中免疫细胞的比例，然而它们依赖有限的表型标记库。并且使用流式细胞术这种方法的时候，有些细胞粘附性高，骨髓穿刺难以抽出，无法分析，从而导致免疫细胞计数错误。随着计算机技术的高速发展，已经产生了通过基因表达谱数据，利用各种机器学习算法计算免疫细胞数量的方法。

已经有团队开发出一种强大的免疫细胞计数算法CIBERSORT，它基于微阵列数据(RNA-seq数据也可适用)，利用v-SVR(一种改进的支持向量回归模型)线性回归求解线性方程模型，回归模型的系数即代表22种免疫细胞的相对比例。之后在2016年，Severson E等人设计出了TIMER，基于约束最小二乘法去计算包括CD4 T细胞、CD8 T细胞和B细胞在内的6种免疫细胞的丰度，较好的解决了高维特征带来的多重共线性问题，作者计算并分析了TCGA数据库中的所有癌症样本。另一种MCP-counter方法基于最小二乘算法，去计算包括8类免疫细胞在内的10种细胞群的丰度。2017年，De Jonge K等人基于微阵列数据，利用最小二乘法计算包括T细胞、B细胞等在内的6种免疫细胞的数量，然后再将mRNA数量转化为相关免疫细胞的相对比例的EPIC算法，其具有临床适用性。

目前，机器学习算法广泛应用于生物信息学领域，但在免疫细胞比例计算方面，成熟的技术方案较少，主要还是CIBERSORT、TIMER、MCP-counter和EPIC这四种算法。但这几种算法只是在枚举混合物中未知细胞类型时有效性较高，免疫细胞比例估算的准确度受到未知混合物含量的影响较大，并且能准确估算的免疫细胞种类较少。

发明内容

本发明实施例提供了浸润免疫细胞比例计数方法，能够显著降低未知混合物含量对免疫细胞比例计数，提高免疫细胞比例计数的准确度。

本发明实施例提供了一种浸润免疫细胞比例计数方法，包括：

将基因表达谱数据转化为免疫细胞的标记基因的线性组合，得到线性组合方程；

基于得到的线性组合方程，构建用于反卷积求解免疫细胞比例的支持向量回归模型，即：免疫细胞比例计数模型；

输入训练集，利用粒子群算法对构建的支持向量回归模型进行迭代优化，其中，所述训练集中的每个样本包含个体抽检组织中不同基因的表达谱数据；

利用优化后的支持向量回归模型，求解免疫细胞比例。

进一步地，所述线性组合方程表示为：

x_{n x L}＝S_{n x m} xf_{m x L}

其中，x_{n x L}表示L个样本的n个基因的表达谱数据，L表示样本的数目，n表示基因的数目；f_{m x L}为待求的免疫细胞比例，具体表示L个样本的每种免疫细胞在m种免疫细胞中所占比例，m表示免疫细胞的种类；S_{n x m}为签名矩阵，表示m种免疫细胞，每种免疫细胞中的n个基因的表达谱数据。

进一步地，所述支持向量回归模型表示为：

其中，k(x_i，x)＝φ(x_i)φ(x)，并且得到w、b的值分别为：

其中，f(x)表示支持向量回归模型解的形式；w代表线性组合方程中待求的f_{m x L}，表示待求的免疫细胞比例；b表示残差项；ε表示敏感度；a_i、

都表示拉格朗日乘子，a_i∈[0，C]、

C表示惩罚因子；k(x_i，x)为满足Mercer条件的核函数，φ(x)表示非线性变换，x＝{x₁，x₂，...，x_n}，x_i表示不同样本的第i个基因的表达谱数据组成的向量，i＝1，2，...，n。

进一步地，所述输入训练集，利用粒子群算法对构建的支持向量回归模型进行迭代优化包括：

Step 1：将每个基因分别当作一个粒子，每个粒子由三维参数(ε，C，φ)组成，初始化ε、C、φ的值，其中，C表示惩罚因子，ε表示敏感度，φ表示核函数系数；

Step 2：利用训练集中基因的表达谱数据训练支持向量回归模型，将每个粒子的个体最优值P_ibest定为当前位置，计算每个粒子的适应度，取最小适应度的粒子所对应的p_ibest作为全局最优值g_best；

Step 3：更新粒子的位置、速度；

Step 4：再计算适应度值，若粒子的适应度值小于p_ibest，则更新为p_ibest，否则不变；

Step 5：若粒子的p_ibest大于g_best，则更新为g_best，否则不变；

Step 6：若所得的解g_best＝(ε，C，φ)收敛到预设的向量值或者达到最大迭代次数，则终止迭代，当前的解g_best＝(ε，C，φ)为最优参数解，否则返回到Step 2。

进一步地，适应度函数取为均方根误差，其中，适应度计算方法表示为：

其中，

表示式(6)所示的支持向量回归模型计算得出的第i个基因在第j种免疫细胞的表达谱数据，y_ij表示训练集样本中第i个基因在第j种免疫细胞的表达谱数据，RMSE表示均方根误差。

进一步地，粒子的位置、速度更新公式为：

v_i+1＝qv_i+c₁r₁(p_ibest-x_i)+c₂r₂(g_best-x_i)

x_i+1＝x_i+v_i+1

其中，v_i、x_i分别表示第i个粒子的速度向量和位置向量，i＝1，2，...，m，m表示群体的大小，具体指免疫细胞的种数；v_i+1、x_i+1分别表示第i+1个粒子的速度向量和位置向量；q表示为非负的惯性因子；c₁和c₂均为学习因子；r₁和r₂都表示随机系数；P_ibest表示第i个粒子的个体最优值，g_best为全局最优值。

进一步地，非负的惯性因子表示为：

其中，q_max、q_min分别表示q的上、下限，iter_max表示最大迭代次数，iter_min表示当前迭代次数。

进一步地，所述利用优化后的支持向量回归模型，求解免疫细胞比例包括：

将求得的最优参数解作为支持向量回归模型的输入参数，求解免疫细胞比例。

本发明实施例提供的技术方案带来的有益效果至少包括：

本发明实施例中，利用粒子群算法对构建的支持向量回归模型进行迭代优化，利用优化后的支持向量回归模计算肿瘤浸润免疫细胞比例，这样，能够兼顾支持向量回归分类准确率高、泛化能力强和粒子群算法带来的计算结果的准确性等优点，从而显著降低未知混合物含量对免疫细胞比例计数，提高免疫细胞比例计数的准确度。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的浸润免疫细胞比例计数方法的流程示意图；

图2为本发明实施例提供的浸润免疫细胞比例计数方法的详细流程示意图；

图3为本发明实施例提供的PBMC数据免疫细胞的散点图；

图4为本发明实施例提供的CRC数据免疫细胞的散点图；

图5为本发明实施例提供的Melanoma数据免疫细胞的散点图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

为了降低未知混合物含量对免疫细胞比例计数的影响，本发明实施例引进了支持向量回归(Support Vector Regression，SVR)算法。由于SVR的泛化能力和预测精度取决选择的参数，例如核函数系数和惩罚因子，为了增大算法的鲁棒性，本发明再引进了粒子群算法(Particle swarm optimization，PSO)，这是一种强大的参数迭代优化求解算法，能够较大地提高了SVR解的准确性。本发明实施例将二者结合起来形成一种新的计算浸润免疫细胞比例的方法，命名为EPCPS(Estimating the proportion of cells using PSO-SVR)。

为了更好地理解本发明，先对SVR和粒子群算法进行简要说明：

SVR使用SVM来拟合曲线，做回归分析，也是为了找出一个超平面，使得所有数据到这个超平面的距离最小，即残差最小。SVR利用一条固定宽度的条带(宽度由一个参数来控制)覆盖尽可能多的样本点，从而使得总误差尽可能的小。SVR模型当样本点都落在“隔离带”内时分类效果最好，而SVM模型当样本点都在“隔离带”外时回归拟合效果最好，这导致SVR要同时引入两个松弛变量。

粒子群优化算法是一种进化计算技术，被广泛用于函数优化，通过迭代可以找到本文SVR模型参数的最优解，可以为松弛变量的引入提供依据，同时使得分类器可以提高准确度。SVR的复杂度和泛化能力取决于ε、C、φ这三个参数。选择参数的时候对各个参数分布做优化不太现实。针对该问题，本发明实施例提出了将支持向量回归和粒子群算法相结合。

基于上述问题，为了更好的了解免疫系统对癌症的应对机制，达到有效提高免疫细胞比例计算的准确率的目的。本发明实施例基于RNA-Seq的基因表达谱数据综合考虑了SVR模型算法和PSO参数优化方法，对传统SVR算法做了迭代优化，既减少了人体其他细胞的基因表达对免疫细胞估计的影响，又增加了结果的准确性.达到进一步提高癌症预测的准确率的目的，其中，RNA-seq表示转录组测序技术，用高通量测序技术进行测序分析，反映一些核酸分子的表达水平。

接着，对本发明实施例提供的浸润免疫细胞比例计数方法进行说明，该方法包括：

S101，将基因表达谱数据转化为免疫细胞的标记基因的线性组合，得到线性组合方程；

本实施例中，基于大量组织基因表达谱数据结合支持向量回归和纯化白细胞亚群表达谱数据的先验知识(即代表不同免疫细胞中各基因的表达谱数据的“签名矩阵”)，可以准确估计肿瘤活检的免疫细胞比例其关键思路是：首先将基因表达谱数据转化为免疫细胞的标记基因的线性组合，其目的是求解线性组合方程中的f(代表免疫细胞比例)，所述线性组合方程表示为：

x_{n x L}＝S_{n x m}xf_{m x L} (1)

其中，x_{n x L}表示L个样本的n个基因的表达谱数据，用以反卷积，L表示样本的数目，每个样本包含了个体抽检组织中不同基因的表达谱数据，n表示基因的数目；f_{m x L}为待求的免疫细胞比例，具体表示L个样本的每种免疫细胞在m种免疫细胞中所占比例，m表示免疫细胞的种类；S_{n x m}为签名矩阵，表示m种免疫细胞，每种免疫细胞中的n个基因的表达谱数据。

S102，基于得到的线性组合方程，构建用于反卷积求解免疫细胞比例的支持向量回归模型；

本实施例中，使用机器学习技术v支持向量回归(v-SVR)来反卷积求解f。简而言之，SVR定义了一个超平面，该超平面在给定定义的约束的情况下捕获尽可能多的数据点，并通过使用线性“对ε不敏感”的损失函数仅对特定误差半径之外的数据点(称为支持向量)进行惩罚来减少过度拟合，超平面的方向确定f。

本实施例中，对于回归问题，给定训练集D＝{(x₁，y₁)，(x₂，y₂)，...，(x_n，y_n)}，其中，x＝{x₁，x₂，...，x_n}，y＝{y₁，y₂，...，y_n}，x_i表示输入的不同样本的第i个基因的表达谱数据组成的向量，y_i表示签名矩阵S_{n x m}的第i个向量，y_i＝{y_i1，y_i2，...，y_im)，i＝1，2，...，n，y_ij表示训练集样本中第i个基因在第j种免疫细胞的表达谱数据，j＝1，2，...，m。训练的目的是得到一个支持向量回归模型f(x)＝w^Tx+b使得f(x)与y尽可能接近，w和b是支持向量回归模型参数，在本发明实施例中w代表了式(1)中待求的f_{m x L}，b表示残差项。假设允许f(x)与y之间最多有ε的误差，仅当f(x)与y之间的差的绝对值大于ε时才计算SVR的损失。计算免疫细胞相对比例的准确度的衡量标准是损失函数的大小，为了最小化损失函数，引入如下算法原理：

通过一个非线性变换φ(x)，把x_i映射到高维空间，其中x_i是n维空间的第i个元素，在高维空间里，理论存在一个线性函数f(x)，能够很好的描述输入数据和输出数据之间的联系，这样的线性函数成为SVR函数：

f(x)＝W^Tφ(x)+b，φ：Rⁿ→F，w∈F (2)

其中，F是一个数域，Rⁿ代表实数域上的n维空间；

同时根据SVR回归误差定义下述结构风险函数：

其中，R_reg表示回归误差；式(3)中第一项||w||²为描述函数，第二项为复杂度的项；C为惩罚因子，它的作用是折中支持向量回归模型复杂度和经验风险，即式(3)第二项中的

此处|y-f(x)|_ε为ε不敏感损失函数，通过引入它，把偏差控制住从而让估计得到鲁棒性，|y-f(x)|_ε的定义如下：

上述回归问题等价于最小化代价泛函：

其中

和§_i表示两个松弛变量，其目的是让上述式子的解存在。在求解该优化问题时，先转化为求解对偶问题，然后引入了拉格朗日函数，通过对拉格朗日乘子ai，

以及w，b求偏导，并使各偏导等于零，最后利用KKT(Karush-Kuhn-Tucker conditions)条件求解得到用于反卷积求解免疫细胞比例的支持向量回归(SVR)模型(也可以称为免疫细胞比例计数模型)：

其中，k(x_i，x)＝φ(x_i)φ(x)，并且得到w，b的值分别为：

其中，a_i和

都表示拉格朗日(Lagrange)乘子；ε表示敏感度；a_i∈[0，C]、

C表示惩罚因子；k(x_i，x)为满足Mercer条件的核函数，其中，Mercer条件指一种作为核函数的条件，任何半正定对称函数都可以作为核函数；式(6)中f(x)便是SVR模型解的形式，具体指：

计算得出的训练集样本中各个基因在各种免疫细胞的表达谱数据，而w代表了式(1)中待求的f_{m x L}，即：待求的免疫细胞比例。

S103，输入训练集，利用粒子群算法对构建的支持向量回归模型进行迭代优化，其中，所述训练集中的每个样本包含个体抽检组织中不同基因的表达谱数据；

本实施例中，式(6)所示的支持向量回归模型模型计算的免疫细胞比例准确性受到输入参数(惩罚因子C、敏感度ε和核函数系数φ)的影响，为了找到最优输入参数，结合粒子群算法，通过迭代找到最优(ε，C，φ)值，在每一次的迭代中，粒子通过跟踪两个“最优值”(p_ibest，g_best)来更新自己(详见Step1-Step7)，其中，p_ibest表示粒子的个体最优值，g_best为全局最优值。具体原理如下：

粒子群算法主要有两个向量，第(i+1)个粒子的速度向量v_i+1和位置向量x_i+1分别的更新公式为：

v_i+1＝qv_i+c₁r₁(p_ibest-x_i)+c₂r₂(g_best-x_i) (9)

x_i+1＝x_i+v_i+1 (10)

其中，v_i、x_i分别表示第i个粒子的速度向量和位置向量，这里把每个基因分别当作一个粒子，i＝1，2，...，n，n表示群体的大小，具体指基因的数目；q表示为非负的惯性因子，其值越大，全局寻优能力越强，局部寻优能力越弱；c₁和c₂均为学习因子，一般c₁＝c₂＝2；r₁和r₂都表示属于[0，1]的随机系数；p_ibest表示第i个粒子的个体最优值，g_best为全局最优值。

为了加快迭代速度，q值应该随着迭代次数的增加而减小，本实施例中定义为：

其中，q_max、q_min分别表示q的上、下限，iter_max表示最大迭代次数，iter_min表示当前迭代次数。

将基因表达谱数据中的x_i作为式(9)的输入数据(即：式(9)中的x_i)，根据得出的v(x)(指：x_i经过式(9)求出的速度变量v_i+1)代入式(6)所示的支持向量回归模型，拟合样本基因表达谱数据集{x_i，y_i}(i＝1，2，...，n；x_i∈Rⁿ；y_i∈R)得到的新的回归函数为：

其中，w^*、b^*分别为更新后的w、b，a_i，

为拉格朗日乘子。同式(7)(8)的解法一样，粒子群算法优化后的模型参数解为：

得到的新的回归函数用于取代优化前的支持向量回归模型。

本实施例中，将支持向量回归和粒子群算法相结合，每个粒子由三维参数(ε，C，φ)组成，适应度函数取为RMSE，其中，适应度计算方法表示为：

其中，

表示式(6)所示的支持向量回归模型计算得出的第i个基因在第j种免疫细胞的表达谱数据，即：计算得到的f(x)；y_ij表示训练集样本中第i个基因在第j种免疫细胞的表达谱数据，RMSE表示均方根误差。

本实施例中，基于PSO-SVR算法的(ε、C、φ)参数的优化选择过程如下述所示(其中p_ibest和g_best都是指待求向量(ε，C，φ)在当前迭代次数内的值)：

Step 1：每个粒子由三维参数(ε，C，φ)组成，初始化ε、C、φ的值，其中，C表示惩罚因子，ε表示敏感度ε，φ表示核函数系数；

Step 2：利用训练集中基因的表达谱数据训练支持向量回归模型，将每个粒子的个体最优值p_ibest定为当前位置，利用式(12)和(15)计算每个粒子的适应度，取最小适应度的粒子所对应的P_ibest作为全局最优值g_best；

Step 3：利用式(9)到(11)式，更新粒子的位置、速度；

Step 4：再计算适应度值，若粒子的适应度值小于p_ibest，则更新为p_ibest，否则不变；

Step 5：若粒子的p_ibest大于g_best，则更新为g_best，否则不变；

Step6：若所得的解g_best＝(ε，C，φ)收敛到预设的向量值(包含3个值)或者达到最大迭代次数就终止迭代，否则返回到Step 2。

本实施例中，得到的ε、C、φ等参数的大致范围为：ε＝[0，0.2]，C＝[1，108]，φ＝[0.01，2.0]，q_max，q_min一般分别取为0.9，0.4。

S104，利用优化后的支持向量回归模型，求解免疫细胞比例。

本实施例中，将粒子群的迭代优化过程应用到支持向量回归模型的解中，根据适应度大小不断迭代优化出支持向量回归模型的最优输入参数，将求得的最优参数解(ε，C，φ)作为支持向量回归模型的输入参数，得到最优的支持向量回归模型，对得到的最优的支持向量回归模型进行测试，利用测试完成的罪域的支持向量回归模型求解免疫细胞比例，如图2所示。

综上，对于EPCPS的实现，就是把SVR算法嵌入到PSO算法中去计算适应度值的步骤中，然后再将PSO-SVR算法应用到解决计算免疫细胞比例的线性回归模型中去，迭代求解免疫细胞比例f，流程图如图2所示。

与传统的反卷积方法相比，本发明实施例提供的浸润免疫细胞比例计数方法，应用了粒子群算法对支持向量回归做迭代优化求解。另外，EPCPS实行L2范数正则化，从而减轻由于多重共线性带来的问题。

为了验证本发明实施例提供的EPCPS的有效性，将其与其它算法进行比较：

目前计算免疫细胞相对比例有CIBERSORTX、EPIC和MCP-counter等方法，对此，本发明利用三组实际数据的结果验证EPCPS算法，然后再与CIBERSORTX、EPIC、MCP-counter等方法进行比较。第一组数据为20个PBMC样本的基因表达谱数据和用流式细胞术测得的免疫细胞的比例。第二组数据为10个CRC样本的基因表达谱数据和用免疫组化法测得的免疫细胞的比例。第三组数据为19个melanoma样本的基因表达谱数据和用单细胞测序技术得到的免疫细胞的比例。

为了验证本文算法的准确性，本发明使用三组实际数据进行验证，具体信息如下表1。

表1三组验证数据的具体信息

1)DataSet 1：20个患者血液样本的基因表达谱数据，数据集的来源是CIBERSORT，已经通过流式细胞术测量得到了4种免疫细胞的相对比例。

2)DataSet 2：10个结直肠癌症样本的基因表达谱数据。研究人员也已经通过免疫组化法测量得到的6种细胞的相对比例。

3)DataSet 3：19个melanoma样本的基因表达谱数据和用单细胞测序技术得到的免疫细胞的比例，Newman等人通过单细胞测序技术测得的免疫细胞的相对比例。

对于DataSet 1，流式细胞术用于鉴定并确定异质细胞群中的不同细胞类型、评估分离亚群的纯度以及分析细胞大小和容积。对于DataSet 2，免疫组化学是利用抗原抗体反应，标记抗体的显色剂显色的化学反应来对其进行相对定量的研究。对于DataSet 3，scRNA-Seq是能够在单个细胞的水平上，能够分析相同表型细胞的遗传异质性的一项新技术，从而计算分离计算免疫细胞的数量。

首先对于第一组实际数据，流式细胞术已经给出了B细胞、CD8 T细胞、CD4 T细胞的相对比例，本发明以流式细胞术和EPCPS的结果分别为横、纵轴做图3，其中，图3中的PBMC表示外周血单个核细胞，fraction表示各种免疫细胞比例计数方法计算得出的相对得分。对于第二组数据，以免疫组化法和EPCPS计算的结果为分别作为横、纵轴，如图4，其中，图4中的CRC表示结直肠癌。第三组数据集单细胞测序技术给出了Malignant、T cells CD8、Tcells CD4、B cells、NK cells、Macrophages、Endothelial cells、CAF等8种细胞的相对比例，本文以单细胞测序技术得到的免疫细胞比例为x轴，EPCPS测得结果为y轴做散点图5，图5中的melanoma表示黑色素瘤。

如图3第一组PBMC数据，可以看到B细胞、CD4_Tcells和CD8_T cells基本上靠近y＝x标准线，对于CD 8 T细胞偏离程度相对另外两个较大，这可能是因为CD 8T细胞在PBMC中的T淋巴细胞占比较少(5％-20％)。图3中R₁表示三种免疫细胞的平均斯皮尔曼相关系数，R₂表示皮尔逊相关系数。总体来说3种免疫细胞均具有较好的相关性。

第二组实际数据如图4，R₁表示四种免疫细胞的平均斯皮尔曼相关系数，R₂表示皮尔逊相关系数。除了个别样本的B cells结果差距比较大，其它样本的B cells、Monocytes、NK cells、和T cells结果都比较理想。

第三组实际数据如图5，R₁表示五种免疫细胞的平均斯皮尔曼相关系数，R₂表示皮尔逊相关系数。除了B细胞个别样本偏差较大，其余诸如CD4_T cells，CD8_T cells等免疫细胞均呈现较好的相关性。

为了将发明的结果与流式细胞术，免疫组化学以及单细胞测序技术得到的结果做一些具体比较。选取了三个评价指标，分别是RMSE，Pearson相关系数和spearman相关系数来全面的评价EPCPS的算法的准确性。各个指标如表2、表3、表4。

本实施例中，Pearson相关系数是用来衡量两个数据集合是否在一条线上面，它用来衡量定距变量间的线性关系。Spearman相关系数是利用单调方程评价两个统计变量的相关性。

表2实际数据集1

细胞类型	均方根误差	皮尔逊	斯皮尔曼
				B cells	0.1023	0.7471	0.7642
CD4 T cells	0.1254	0.6318	0.6021
				CD8 T cells	0.1523	0.4921	0.5012

表3实际数据集2

细胞类型	均方根误差	皮尔逊	斯皮尔曼
				B cells	0.2431	0.6532	0.6721
T cells	0.1932	0.7232	0.7024
				Monocytes	0.1221	0.8021	0.7623
NK cells	0.0935	0.9076	0.8935

表4实际数据集3

细胞类型	均方根误差	皮尔逊	斯皮尔曼
				B cells	0.1521	0.6842	0.6521
CD4 T cells	0.1325	0.6465	0.6934
				Macrophages	0.1261	0.7296	0.7432
CD8 T cells	0.1432	0.6757	0.6461
				NK cells	0.0917	0.7633	0.9023

从表2、表3、表4中可以看出，B细胞的误差普遍最大，皮尔逊和斯皮尔曼相关系数都相对较低，这可能是因为B淋巴细胞相对于其他淋巴细胞所占的数量少的缘故。对于CD4T细胞和CD8 T细胞，第一组数据和第三组数据呈现的结果差距较小，两者的RMSE和相关系数很接近，间接的证明了本发明的算法的可靠性。对于NK细胞，第二组数据和第三组数据表明在本文计算的免疫细胞中，NK细胞的RMSE均最小，相关系数最大，第二组数据的NK细胞RMSE为0.0935，pearson为0.8935，在第三组数据NK细胞的RMSE为0.0917，pearson为0.9023。

最后本发明将EPCPS计算的结果与其它算法如CIBERSORTX，EPIC，MCP-counter等计算浸润肿瘤免疫细胞比例的方法相比较，比较结果如表5、表6、表7所示。

表5实际数据集1

表6实际数据集2

表7实际数据集3

从表5、表6、表7可以看出，RMSE，皮尔逊相关系数和斯皮尔曼相关系数这三个指标来看EPCPS整体上均优于EPIC，MCP-counter等方法，但是第一组和第二组数据的总体结果不如CIBERSORTX，但是差距小。总体来说，本发明的准确率较高。

本发明对基于基因表达谱数据计算免疫细胞比例的算法进行了改进，最大的关键点在于：利用粒子群算法(PSO)找到SVR模型的最优(ε，C，φ)和回归拟合训练相结合来计算肿瘤浸润免疫细胞比例，这样，能够兼顾支持向量回归分类准确率高、泛化能力强和粒子群算法带来的计算结果的准确性等优点，从而显著降低未知混合物含量对免疫细胞比例计数，提高免疫细胞比例计数的准确度。结果表明EPCPS的性能在RMSE，Pearson和spearman相关系数等三个指标的准确性普遍强于MCP-counter、EPIC等免疫细胞计数法，本发明实现了一种更精确的免疫细胞比例计算方法。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种浸润免疫细胞比例计数方法，其特征在于，包括：

将基因表达谱数据转化为免疫细胞的标记基因的线性组合，得到线性组合方程；

基于得到的线性组合方程，构建用于反卷积求解免疫细胞比例的支持向量回归模型，即：免疫细胞比例计数模型；

输入训练集，利用粒子群算法对构建的支持向量回归模型进行迭代优化，其中，所述训练集中的每个样本包含个体抽检组织中不同基因的表达谱数据；

利用优化后的支持向量回归模型，求解免疫细胞比例。

2.根据权利要求1所述的浸润免疫细胞比例计数方法，其特征在于，所述线性组合方程表示为：

x_nⅹL＝S_nⅹmⅹf_mⅹL

其中，x_nⅹL表示L个样本的n个基因的表达谱数据，L表示样本的数目，n表示基因的数目；f_mⅹL为待求的免疫细胞比例，具体表示L个样本的每种免疫细胞在m种免疫细胞中所占比例，m表示免疫细胞的种类；S_nⅹm为签名矩阵，表示m种免疫细胞，每种免疫细胞中的n个基因的表达谱数据。

3.根据权利要求1所述的浸润免疫细胞比例计数方法，其特征在于，所述支持向量回归模型表示为：

其中，k(x_i,x)＝φ(x_i)φ(x)，并且得到w、b的值分别为：

其中，f(x)表示支持向量回归模型解的形式；w代表线性组合方程中待求的f_mⅹL，表示待求的免疫细胞比例；b表示残差项；ε表示敏感度；a_i、

都表示拉格朗日乘子，a_i∈[0,C]、

C表示惩罚因子；k(x_i,x)为满足Mercer条件的核函数，φ(x)表示非线性变换，x＝{x₁,x₂,…,x_n}，x_i表示不同样本的第i个基因的表达谱数据组成的向量，i＝1,2,…,n。

4.根据权利要求1所述的浸润免疫细胞比例计数方法，其特征在于，所述输入训练集，利用粒子群算法对构建的支持向量回归模型进行迭代优化包括：

Step 1：将每个基因分别当作一个粒子，每个粒子由三维参数(ε，C，φ)组成，初始化ε、C、φ的值，其中，C表示惩罚因子，ε表示敏感度，φ表示核函数系数；

Step 2：利用训练集中基因的表达谱数据训练支持向量回归模型，将每个粒子的个体最优值p_ibest定为当前位置，计算每个粒子的适应度，取最小适应度的粒子所对应的p_ibest作为全局最优值g_best；

Step 3：更新粒子的位置、速度；

Step 4：再计算适应度值，若粒子的适应度值小于p_ibest，则更新为p_ibest，否则不变；

Step 5：若粒子的p_ibest大于g_best，则更新为g_best，否则不变；

Step 6：若所得的解g_best＝(ε，C，φ)收敛到预设的向量值或者达到最大迭代次数，则终止迭代，当前的解g_best＝(ε，C，φ)为最优参数解，否则返回到Step 2。

5.根据权利要求1所述的浸润免疫细胞比例计数方法，其特征在于，适应度函数取为均方根误差，其中，适应度计算方法表示为：

其中，

表示式(6)所示的支持向量回归模型计算得出的第i个基因在第j种免疫细胞的表达谱数据，y_ij表示训练集样本中第i个基因在第j种免疫细胞的表达谱数据，RMSE表示均方根误差。

6.根据权利要求5所述的浸润免疫细胞比例计数方法，其特征在于，粒子的位置、速度更新公式为：

v_i+1＝qv_i+c₁r₁(p_ibest-x_i)+c₂r₂(g_best-x_i)

x_i+1＝x_i+v_i+1

其中，v_i、x_i分别表示第i个粒子的速度向量和位置向量，i＝1,2,…,n，n表示群体的大小，具体指免疫细胞的种数；v_i+1、x_i+1分别表示第i+1个粒子的速度向量和位置向量；q表示为非负的惯性因子；c₁和c₂均为学习因子；r₁和r₂都表示随机系数；p_ibest表示第i个粒子的个体最优值，g_best为全局最优值。

7.根据权利要求6所述的浸润免疫细胞比例计数方法，其特征在于，非负的惯性因子表示为：

其中，q_max、q_min分别表示q的上、下限，iter_max表示最大迭代次数，iter_min表示当前迭代次数。

8.根据权利要求1所述的浸润免疫细胞比例计数方法，其特征在于，所述利用优化后的支持向量回归模型，求解免疫细胞比例包括：

将求得的最优参数解作为支持向量回归模型的输入参数，求解免疫细胞比例。

Images