阅读说明：本技术 一种基于多连体特征的t细胞受体对应表位预测方法 (T cell receptor corresponding epitope prediction method based on multiconnector characteristics ) 是由 王嘉寅 童瑶 杨玲 郑田 刘涛 李敏 张选平 于 2020-03-19 设计创作，主要内容包括：本发明公开了一种基于多连体特征的T细胞受体对应表位预测方法,将CDR3β链以及对应的表位解析为长度3的碱基,统计每种三联体的频次作为初始特征；根据得到的初始特征建立初始特征矩阵,使用主成分分析法对初始特征矩阵进行降维,进行特征提取；设有n个训练样本,输入预测数据x后,训练得到梯度提升决策树模型,通过梯度提升决策树模型将各个决策树的决策结果线性组合起来做出预测；将特征数据输入训练好的模型中进行预测,根据不同的预测目的选择不同的预测指标。本发明仅使用三联体的统计值作为初始特征,结合梯度提升决策树模型能够在极短的时间内完成模型的训练,且预测的准确度更高。(The invention discloses a T cell receptor corresponding epitope prediction method based on multiconnector characteristics, which comprises the steps of resolving a CDR3 β chain and a corresponding epitope into bases with the length of 3, counting the frequency of each triplet as initial characteristics, establishing an initial characteristic matrix according to the obtained initial characteristics, reducing the dimension of the initial characteristic matrix by using a principal component analysis method, extracting the characteristics, setting n training samples, inputting prediction data x, training to obtain a gradient lifting decision tree model, linearly combining decision results of each decision tree through the gradient lifting decision tree model to predict, inputting the characteristic data into the trained model to predict, and selecting different prediction indexes according to different prediction purposes.)

一种基于多连体特征的T细胞受体对应表位预测方法

技术领域

本发明属于以精准医学为应用背景的数据科学技术领域，具体涉及一种基于多连体特征的 T细胞受体对应表位预测方法。

背景技术

T细胞受体(英文名称：T cell receptor，英文缩写：TCR)和表位(英文名称：epitope) 与主要组织相容性复合体(英文名称：major histocompatibility complex，英文缩写：MHC)的 特异性结合能够激活免疫系统，从而引发一系列特异性免疫反应。免疫治疗正是利用特异性免 疫系统这一特性，通过研制对应药剂，人为激活免疫系统，使机体的免疫系统重新工作以消灭 体内的入侵者或者癌细胞。因此TCR对应表位的预测可为探索疾病机制、癌症的免疫治疗、 药物开发和疫苗制作等领域提供重要的理论依据。

虽然第二代测序技术(英文名称：Next Generation Sequencing，英文缩写：NGS)提供了 海量的核苷酸序列与氨基酸序列，但是由于标注成本高耗时长，目前有标注的数据仍然很少。 如果能够从当前的少量的有标注数据中训练一个相对可靠的预测模型，就能够应用于TCR表 位的标注问题，节省大量时间与经济成本。另外，由于TCR的基因片段是通过一系列的非同 源重组得到的，其涉及组合来自变量区(V)，多样性(D)和连接(J)基因区段的TCR基因 座以及核苷酸随机插入和/或删除，因此能够产生大量不同的TCR，规模可达10¹⁵～10⁶¹。此外， 由于交叉反应(英文名称：cross-reactivity)的存在，一个TCR能够识别多种表位，同时一个 表位也能识别多种TCR。通过人工以及统计学方法很难从这样的数据中发现TCR和pMHC的 匹配模式，如果能借助机器学习算法研究TCR和pMHC的特异性结合机制，在免疫治疗的过 程中具有重要意义。

TCR可以分为四个CDR(Complementarity determining region)区域：CDR1，CDR2，CDR2.5 和CDR3，抗原的特异性识别主要取决于CDR区域。其中CDR3区域多样性最高，主要和表 位的肽链结合，CDR1，CDR2和CDR2.5主要和MHC分子结合，不过也可以和肽链结合。已 有研究发现CDR3β链在预测表位的过程中发挥了主要的作用，但尚不清楚是CDR3β链中的物理化学特性或者结构特性或者其他因素占据主导作用。

目前，国内外均已经有研究人员尝试研究CDR3和表位数据之间的关系，可以大致分为两 类：第一类使用方法定义TCR或CDR3序列相似性度量方法，在求得序列之间的相似性之后 再利用简单的分类器如K近邻(英文名称：K-nearest neighbor，英文缩写：K-nn)算法进行分 类。第二类方法基于TCR或CDR3序列提取氨基酸的物理化学特征或者基于BLOSUM矩阵对 氨基酸序列进行编码，再利用机器学习模型训练得到一个预测模型。

但是，这两类方法的预测性能都不是很好，而且主要存在以下问题：其一，第一类方法需 要计算任意两条TCR序列之间的相似性，因此其计算相似性的时间复杂度是O(n²)，训练过程 比较耗时。其二，第二类方法基本都是基于氨基酸进行编码，由于不同的CDR3序列长度不一 定相等，因此需要进行对齐操作以确保每条TCR序列的特征向量具有相同的维度，已有的方 法的对齐操作都比较暴力，并没有统一的标准。其三、第一类方法主要考虑两条TCR序列整 体的相似性，第二类方法主要考虑序列中每个氨基酸的信息，没有方法考虑到TCR序列中相 邻氨基酸提供的信息在TCR和表位的特异性识别过程中发挥的作用。

发明内容

本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种基于多连体特征的 T细胞受体对应表位预测方法，解决特征提取繁琐复杂的问题；解决模型训练耗时的问题，能 够在较短时间内完成模型训练；能够直接进行多分类预测。

本发明采用以下技术方案：

一种基于多连体特征的T细胞受体对应表位预测方法，包括以下步骤：

S1、将CDR3β链以及对应的表位解析为长度3的碱基，统计每种三联体的频次作为初始 特征；

S2、根据步骤S1得到的初始特征建立初始特征矩阵，使用主成分分析法对初始特征矩阵 进行降维，进行特征提取；

S3、设有n个训练样本，输入预测数据x后，训练得到梯度提升决策树模型，通过梯度提 升决策树模型将各个决策树的决策结果线性组合起来做出预测；

S4、将步骤S2的特征数据输入步骤S3训练好的模型中进行预测，根据不同的预测目的选 择不同的预测指标。

具体的，步骤S2具体为：

S201、将初始特征矩阵记为：X＝{x₁,x₂,...,x_n}，对每一列特征进行中心化；

S202、令样本点x_i在新空间中超平面上的投影为W^Tx_i，若所有的样本点都分开，使投影后 样本点的方差最大化，确定优化目标；

S203、使用拉格朗日乘数法对优化目标件求解，对协方差矩阵XX^T进行特征分解，将求 得的特征值排序；再取前k个特征值对应的特征向量构成投影矩阵W，最终得到的特征矩阵 W^TX为一个k行n列的矩阵。

进一步的，步骤S201中，m维列向量x₁为：

其中，n是训练样本个数，m是特征维数。

进一步的，步骤S202中，优化目标为：

其中，W是转换矩阵，W^T是转换矩阵的转置矩阵，X是初始特征矩阵，X^T是初始特征矩阵的转置。

进一步的，步骤S203中，对优化目标件求解，得到

XX^TW＝λW

投影矩阵W为：

W＝(w₁,w₂,...,w_k)

其中，λ是特征值，w_i是投影矩阵的列向量，1≤i≤k，特征值的排序为：λ₁≥λ₂≥...≥λ_n。

具体的，步骤S3具体为：

S301、初始化迭代次数m＝0，设置最大迭代次数为M，初始化模型f₀(x)；

S302、每一次模型迭代在当前模型的基础上添加一颗决策树，使用残差L(y,f_m-1(x))估计 参数Θ_m；

S303、令m＝m+1，若m小于最大迭代次数，则返回步骤S302；否则，停止训练，返回以上训练的所有决策树，完成表位预测模型训练。

进一步的，步骤S301中，初始化模型f₀(x)为：

其中，N是样本数量，c是初始模型拟合的常数，L是对数似然损失函数，定义为：

其中，Y为输出变量，X为输入变量，L为损失函数，M为表位类别个数，y_ij是一个二值指标，若类别j是否是输入示例x_i的真实类别，则y_ij＝1；否则y_ij＝0，p_ij为模型预测输入实例x_i属于类别j的概率。

进一步的，步骤S302中，第m次迭代的结果为：

f_m(x)＝f_m-1(x)+β_mT(x；Θ_m)

其中，f_m-1(x)是第m-1次迭代的决策模型，使用所有的R_mi的集合i∈[1..n] 来拟合一个回归分类决策树。

更进一步的，使用残差L(y,f_m-1(x))估计参数Θ_m，决策树的参数Θ_m按照如下的优化目标 求解得到的：

损失函数在模型f_m-1上的负梯度用来近似估计残差为：

其中，i是第i个训练样本的索引。

进一步的，步骤S303中，

其中，f_M(x)为最终得到的由M颗决策树组成的集成模型，M为表位类别个数，β_m为第m颗决策树的权重，T为决策树，x为决策树T的输入，Θ_m为决策树的参数。

与现有技术相比，本发明至少具有以下有益效果：

本发明是一种基于TCR序列中联合多联体特征进行TCR表位预测的方法，逐个扫描CDR3 β序列，将多肽链解析为连续的长度为3的短肽链，并统计各三联体出现的频次。将统计结果 作为初始特征矩阵，将CDR3β序列对应的表位作为类别标签。根据生物学知识，人体内有20 种氨基酸，20种氨基酸最多能够有8000种不同的排列组合，因此初始的特征矩阵维数最高不 超过8000维，解决现有模型特征提取繁琐复杂以及模型训练耗时的问题，能够在较短时间内 完成训练，且预测性能优于现有模型。

进一步的，使用主成分分析进行特征变换，降低特征的维度。

进一步的，将特征矩阵输入到梯度提升决策树(英文名称：Gradient BoostingDecision Tree， 英文缩写：GBDT)进行训练，通过网格搜索得到模型的最佳参数，最终得到多颗决策树。

进一步的，用同样的方法对测试数据进行编码，将测试数据特征矩阵输入到模型种，取所 有决策树预测的结果之和作为最终预测结果。

综上所述，本发明仅使用三联体的统计值作为初始特征，结合梯度提升决策树模型能够在 极短的时间内完成模型的训练，且预测的准确度更高。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为对Dash et.al论文中的TCR数据进行特征选择之后得到的特征矩阵；

图2为本发明流程示意图；

图3为在Dash数据集中不同模型对比结果；

图4为在Dash数据集中做多分类ROC曲线结果。

具体实施方式

本发明提供了一种基于TCR序列中相邻氨基酸信息的TCR表位预测方法，命名为SETE (Sequence-based Ensemble learning approach for TCR Epitope bindingprediction)，训练集的数 据为CDR3β序列及对应的能够进行特异性识别的多肽链，测试集的数据为CDR3β序列。

基于以下学术界的普遍共识：

1.TCR序列中CDR3区域和MHC呈递的多肽链具有明显的相互作用，并且该区域的β链 在肽的识别中有着显著贡献；

2.目前已知的人体内组成蛋白质的氨基酸共有20种。

请参阅图2，本发明一种基于多连体特征的T细胞受体对应表位预测方法，包括以下步骤：

S1、提取初始特征

请参阅图1，由于输入数据是CDR3β链以及对应的表位，表位是模型预测的类别。由于 无法直接利用输入的氨基酸序列作为特征，需要将其解析为长度为3的碱基，并统计每种三联 体的频次作为初始特征。对得到的初始特征进行特征选择之后，可以发现不同类别表位对应的 TCR序列特征之间具有一定的相似性，对Dash et.al论文中的TCR数据进行特征选择之后得 到的特征矩阵，图中x轴表示特征，y轴表示样本；最右侧的色带表示TCR对应的表位类别； 颜色越深表示三联体数目越多。

S2、特征提取

一方面由于一共有20种氨基酸，长度为3的氨基酸短链最多会有20³种组合，因此特征 最多能够达到8000维，需要做特征筛选降低特征的维度；其次，由于相似的TCR序列之间存 在着相似性，同一类别的TCR序列的三联体之间可能存在冗余信息。因此使用主成分分析对 数据进行降维，具体包含以下步骤：

S201、将初始特征矩阵记为：X＝{x₁,x₂,...,x_n}，对每一列特征进行中心化；

m维列向量x₁为：

其中，n是训练样本个数，m是特征维数；

S202、令样本点x_i在新空间中超平面上的投影为W^Tx_i，若所有的样本点都分开，使投影后 样本点的方差最大化，确定优化目标为：

其中，W是转换矩阵，W^T是转换矩阵的转置矩阵，X是初始特征矩阵，X^T是初始特征矩阵的转置。

S203、使用拉格朗日乘数法对优化目标件求解，得到XX^TW＝λW，对协方差矩阵XX^T进行特征分解，将求得的特征值排序：λ₁≥λ₂≥...≥λ_n，再取前k个特征值对应的特征向量构成 投影矩阵W，最终得到的特征矩阵W^TX为一个k行n列的矩阵；

投影矩阵W为：

W＝(w₁,w₂,...,w_k)

其中，λ是特征值，w_i是投影矩阵的列向量，1≤i≤k。

S3、表位预测模型训练

提出一种新的基于梯度提升决策树的预测模型；假如有n个训练样本，输入预测数据x后， 梯度提升决策树模型通过将各个决策树的决策结果线性组合起来做出预测，具体包含如下步 骤：

n个训练样本为：

{(x₁,y₁),...,(x_n,y_n)}

其中，i＝1，2，...，n；

S301、模型初始化

初始化迭代次数m＝0，设置最大迭代次数为M，初始化模型f₀(x)为：

其中，N是样本数量，c是初始模型拟合的常数，L是对数似然损失函数，定义为：

其中，Y为输出变量，X为输入变量，L为损失函数，M为表位类别个数，y_ij是一个二值指标，若类别j是否是输入示例x_i的真实类别，则y_ij＝1；否则y_ij＝0，p_ij为模型预测输入实例x_i属于类别j的概率。

S302、模型迭代

模型的每一次迭代都会在当前模型的基础上再添加一颗决策树，第m次迭代的结果为：

f_m(x)＝f_m-1(x)+β_mT(x；Θ_m)

其中，f_m-1(x)是第m-1次迭代的决策模型，决策树的参数Θ_m是按照如下的优化目标求 解得到的：

由于基函数是线性相加的，目标是使用残差L(y,f_m-1(x))估计参数Θ_m。

为此，损失函数在模型f_m-1上的负梯度用来近似估计残差。

其中，i是第i个训练样本的索引。

使用所有的R_mi的集合i∈[1..n]来拟合一个回归分类决策树(英文名称： Classification and Regression Tree，英文缩写：CART)，求解其参数Θ_m。

S303、将m赋值为m+1，若m小于最大迭代次数，则返回步骤S302；否则，停止训练，返回以上训练的所有决策树；

S4、表位预测

用同样的方法提取初始特征并进行特征提取，将最终的数据输入训练好的模型中，进行预 测，根据不同的预测目的选择不同的预测指标。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图， 对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分 实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通 过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并 非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实 施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本 发明保护的范围。

本发明一种基于多连体特征的T细胞受体对应表位预测方法，解决了现有算法训练时间 长，预测结果并不理想的问题。

由于目前并没有模型直接能够进行TCR表位多分类预测问题，为了验证本发明的有效性， 首先测试其二分类的预测效果。由于现有方法TCRGP使用受试者工作特征(英文名称：Receiver Operating Characteristic，英文缩写：ROC)以及ROC曲线线下面积(英文名称：Area Under Curve， 英文缩写：AUC)作为模型的评价指标，使用AUC评价本发明的预测性能。

另外，还对比了两个模型在相同数据集上的运行时间；然后，进行了多分类预测的测试， 由于各个类别的数据样本量不均衡，而ROC受数据不均衡的影响较小，因此，仍然使用ROC 和AUC指标衡量模型的预测性能。指标名称：真阳性(TP)、假阳性(FP)、真阴性(TN)和假阴 性(FN)。

定义假阳性率FPR＝FP/(FP+TN)。

定义真阳性率TPR＝TP/(TP+FN)。

ROC曲线由不同阈值下FPR和TPR的值绘制而成，AUC是ROC曲线的线下面积。

在公共数据集VDJdb上进行了测试。经过筛选，从VDJdb数据种选择了22类表位数据。 由于目前已有的模型均只能处理二分类任务，为了和其他模型进行对比，首先进行二分类测试。 在每一个二分类任务中，使用所有正例数据，并随机从其他类别中抽样等量的TCR数据作为 负例。二分类结果如表1所示。

表1：SETE与TCRGP二分类结果比较(*指FRDYVDRFYKTLRAEQASQE)

由上表得出，在二分类任务中，本发明与现有方法TCRGP相较，预测效果相当，但耗时 显著缩短，大大的降低了训练时长，这种优势在数据量大的数据集上尤为明显。

在多分类任务中，本发明也进行了一系列实验来验证其有效性。我们使用ROC曲线作为 评价模型的指标，为了绘制出多分类的ROC曲线我们使用了OneVsRest策略，对每个类别的 数据训练一个分类器，每个分类器将一类TCR序列看做正例，将其他类TCR序列看做负例， 最终对十个分类器的输出结果进行投票得到最终的分类结果。使用五折交叉验证得到的结果如 表2所示。

表2：SETE在VDJdb数据集上的多分类预测结果

为了进一步验证本发明预测TCR对应表位的能力，在Dash et.al论文中公开的数据集中 进行测试，该数据集共收集了3类人类的表位数据和7类小鼠的表位数据。

由于模型更适合进行多分类任务，因此首先在该数据集上进行多分类测试，使用ROC曲 线和AUC结果评估模型效果，在Dash数据集中多分类结果如表3所示。

表3：SETE在Dash数据集中多分类结果

由上表可知，SETE在多分类问题中的整体表现优异，在个别表位基因上的预测结果较差， 如pp65，可能和该类型的表位基因的数据量较小有一定关系。SETE多分类和TCRGP、TCRdist 的对比结果如图3所示，x轴表示不同预测模型，y轴表示各个模型ROC曲线的线下面积；另 外，分别绘制人类和小鼠数据多分类的ROC曲线，结果如图4所示。图中x轴表示假正例率， y轴表示真正例率。

在Dash数据集中进行了二分类测试，二分类的预测结果如表4所示。

表4：SETE与TCRGP在Dash数据集中的二分类结果比较

和之前的结果一样，SETE能够在极短的时间内完成训练，且预测准确率优于TCRGP模 型。

综上所述，本发明比现有的方法TCPGP相比，能够在更短的时间内完成模型的训练，并 且在二分类任务中的表现优于现有方法。此外，本发明能够直接应用于多分类任务，且预测的 准确率较高。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明 提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之 内。