摘要
单细胞转录组测序（scRNA-seq）与T细胞受体免疫组库测序（scTCR-seq）是解析免疫细胞特性的两大关键技术，分别从基因表达和抗原识别维度揭示免疫系统的复杂性。然而，传统分析方法多局限于单一模态，难以有效整合两个组学所提供的互补信息。为突破这一局限，实现跨组学数据的高效融合，提出一种新型的数据整合架构scRTIA（single cell RNA and TCR integrative analysis）。该模型基于深度学习理论，以多模态变分自编码器与Transformer为核心架构，将scRNA-seq基因表达矩阵与scTCR-seq的TCR序列特征协同嵌入统一的低维潜在空间，从而构建出能同时保留转录组特征与免疫组库信息的融合细胞表征。在真实数据集上的实验验证表明，scRTIA所构建的细胞表征在细胞亚群识别方面表现出显著优越的分辨能力，能够发现传统方法难以识别的、具有特定功能状态的稀有T细胞群体。本研究通过有效深度融合转录组与免疫组库信息，突破了单模态分析的瓶颈，实现了对T细胞身份和功能的多维度刻画，在免疫相关疾病研究和精准医疗领域具有应用价值。
关键词
Abstract
Single-cell RNA sequencing (scRNA-seq) and single-cell T-cell receptor sequencing (scTCR-seq) are pivotal for deciphering immune cell characteristics, offering complementary insights into the immune system′s complexity through gene expression and antigen recognition, respectively. However, conventional analytical methods are often confined to a single modality, hindering the effective integration of complementary information from the two omics. To overcome this limitation and achieve efficient integration of cross-omics data, this study proposes a novel data integration framework named scRTIA (single-cell RNA and TCR integrative analysis). Based on deep learning theory, the model employs a multimodal variational autoencoder and a Transformer as its core architecture, jointly embedding the scRNA-seq gene expression matrix and scTCR-seq TCR sequence features into a unified low-dimensional latent space, thereby constructing a fused cell representation that simultaneously preserves transcriptomic features and immune repertoire information. Experimental validation on real datasets demonstrates that the cell representations generated by scRTIA exhibit significantly superior resolution in identifying cell subpopulations, enabling the discovery of rare T-cell populations with specific functional states that are difficult to detect using traditional methods. By effectively integrating transcriptomic and immunome data, our work transcends the limitations of single-modal analysis and enables a multidimensional characterization of T-cell identity and function, offering valuable insights for immune-related diseases research and precision medicine.
单细胞组学兴起于2009年提出的单细胞转录组测序技术（single-cell RNA sequencing，scRNA-seq），改变了之前批量转录组测序技术只能对整个批次细胞进行低分辨率转录组测序的状况，实现了对单个细胞的全转录组测序。但是，早期的单细胞转录组技术通量低，每次只能测序几十到几百个细胞。经过近几年的发展，在2014年逐渐出现了高通量单细胞转录组测序技术，实现了单细胞测序技术的降本增效。时至今日，单细胞转录组测序的通量已经可以达到数万甚至数十万个细胞[1]。
随着单细胞转录组的出现，越来越多的其他单细胞组学技术也开始蓬勃发展，如单细胞染色质可及性，关注单细胞蛋白质表达情况的数据等，均为研究人员从不同视角研究单细胞组学提供了技术支持[2]。在这些多组学技术中，T细胞免疫组库测序技术（single-cell T cell receptor sequencing，scTCR-seq）是一项重要的技术。这项技术主要关注生物体中一种重要的免疫细胞——T细胞上的免疫受体（TCR）信息，因为这部分信息直接体现了T细胞的特异性[3]。相比scRNA-seq与染色质可及性数据，蛋白组数据的整合研究，scRNA-seq与scTCR-seq的整合相关研究较少，目前为止，大致可以分为基于统计学的方法与基于深度学习的方法，这些方法各有特点且均有一些不足之处。
在基于统计学的方法中，CoNGA模型[4]使用克隆型邻域图分析，识别转录组文件与TCR之间的相关性，并且通过定义“CoNGA评分”进行量化；Tessa模型[5]使用参数贝叶斯层次模型，整合TCR数据和T细胞的基因表达数据，评估TCR序列对T细胞表型的影响。这两篇早期文章奠定了scTCR-seq与scRNA-seq整合研究的基础，阐述了该研究方向的重要性与可行性，但是，这两篇文章均假设具有相似基因表达情况的T细胞具有相似的TCR序列信息，而T细胞分化特异性复杂，此生物学假设不完全准确。自此，通过某种方法实现弱生物学假设下的双模态整合成为了研究重点。随着深度学习技术的发展，逐渐有研究转向使用神经网络对这两个模态数据进行整合。scNAT模型[6]基于联合卷积神经网络的变分自编码器架构对双模态数据进行整合，将双模态数据映射入一个联合潜在空间；mvTCR模型[7]使用Transformer模型对字符型免疫组库数据进行编码，之后利用混合变分自编码器进行整合。这两种方法需要较少的生物学假设，但是，前者只考虑了双模态整合之后的潜在表示，对T细胞的免疫组数据分布考虑不足，而后者更关心TCR的CDR3序列信息，对编码v，d，j基因的重排情况考虑较少。近日提出的MIST方法[8]充分考虑了多模态信息，但是，其免疫组数据只精确到了单细胞程度，没有对TCR双链进行分析，并且，没有对T细胞免疫组库数据分布高度离散的特性进行处理。
因此，本研究提出了scRTIA（single cell RNA and TCR integrative analysis）模型，利用Transformer[9]，基于高斯混合模型的变分自编码器（GMM-VAE）等神经网络架构，结合概率生成模型进行整合。不同于之前方法将免疫组看作一个整体模态的思路，本方法更多地从生物学角度出发，利用TCR由双链构成这个生物学事实，实现从免疫双链开始的对scRNA-seq和scTCR-seq数据的整合分析，从而使得最终的细胞潜在表示可以对T细胞双链的差异进行区分，也可以对T细胞群体进行更精细的刻画，使得方法具有发现新生物学现象的潜力。并且，本研究使用T细胞转录组和免疫组的真实配对数据集，验证方法的可行性。
1 单细胞免疫组库与转录组
1.1 T细胞免疫组库
T细胞是一类重要的特异性免疫细胞，主要执行免疫系统的适应性免疫反应，可以识别攻击被感染的细胞或者癌细胞，并且调控其他免疫细胞参与免疫活动。TCR是镶嵌在T细胞表面的一个蛋白质复合物，可以特异性识别抗原，从而使得T细胞识别目标，被激活以及发挥其免疫作用。
TCR一般由双链组成，大部分T细胞的双链为α链和β链，这类T细胞被称为αβ T细胞，少部分T细胞是γ链和δ链组成的，这类T细胞被称为γδ T细胞。对于这些链来说，都是多基因编码的，并且均分为V区（可变区）和C区（恒定区），3类不同的基因片段。一般地，v，d，j基因，编码β链和δ链的V区，而v和j编码α链和γ链的V区。T细胞成熟过程中，在胸腺中进行V区基因的重排，从而导致了T细胞的高度多样性[10]。
在TCR上，CDR3区（互补决定区3）是最核心和多变的部分，直接决定了免疫细胞可以识别哪一种特定的抗原。CDR3区位于V区，由v，d，j基因的重排进行编码，并且伴随核苷酸的随机插入或缺失，这就导致不同T细胞CDR3区的氨基酸序列几乎都是不同的。由此，T细胞免疫组库的主要内容就有每一个细胞的TCR双链，具体信息包含两部分，一部分是每一条链上的CDR3区氨基酸序列，另一部分是编码这些氨基酸序列对应使用的v，d，j基因的类别。其存储使用通用的国际AIRR社区指定的标准化格式进行[11]，如10x等商业化平台均采用这种存储格式进行数据整理。T细胞免疫组库均为字符型序列以及类别型基因信息，而非传统的欧几里得型数据，这也给其处理带来困难。如何对这种数据进行统计建模，或者使用神经网络进行处理，是此类数据研究的重点内容。
1.2 单细胞转录组
由中心法则，遗传信息由DNA转录到RNA，再由RNA翻译成蛋白质去执行细胞的生物学功能。细胞内部的mRNA既代表了细胞内的基因表达信息，又可以用于推测细胞的蛋白质表达，因此，单细胞转录组处于单细胞组学的核心地位。目前，标准的单细胞转录组数据格式为基因表达矩阵，即行代表细胞，列代表基因，而矩阵中的每一项代表某个细胞某个基因的表达量。
2 方法构建
本方法从TCR的结构出发，按照α链和β链（γ链和δ链）组成TCR免疫组完整信息，免疫组数据与转录组数据整合得到最终潜在表示的流程，模块化地进行整合，最终可以得到α单链（γ单链）的分布、β单链（δ单链）的分布、免疫组信息分布以及免疫组与转录组的整合数据分布。整个方法流程见图1，以下按照模块介绍方法各部分的细节。
图1方法流程
Fig.1Methodology flowchart
2.1 基于Transformer的单链嵌入模块
Transformer是处理序列型数据的重要神经网络架构，不同于循环神经网络（RNN），其对前后序列的关联性主要基于自注意力机制。本研究使用Transformer架构搭建编码器，对字符型氨基酸序列进行连续数值嵌入。具体地，首先，建立氨基酸字母与整数的词汇表映射，使用0将序列扩充到最大序列长度。之后，通过可学习嵌入层，获取氨基酸序列的对应数值语义向量。为了将氨基酸序列的顺序信息注入语义向量，使用正弦余弦编码生成位置向量，即
(1)
(2)
式中：PE表示位置编码; pos代表元素在序列中的位置；i为维度索引，满足0≤i≤dmodel/2，dmodel为嵌入的维数；10 000为位置编码常用参数，来自Transformer论文中的默认参数。
对于处理过后的向量，将其通过若干个Transformer编码层。在每一层Transformer编码层中，向量首先利用自注意力机制
(3)
得到氨基酸数值向量，其中，Q为查询矩阵，K为键矩阵，V为值矩阵，dK为每个键向量的维度，用于缩放注意力，使得注意力计算更加稳定。之后通过前馈神经网络，并且将输出与原始向量通过残差连接得到此层的输出，在通过若干层之后，得到编码的潜在向量。
为了充分利用免疫组库信息，得到的潜在向量分别通过两个解码部分，一部分是两个类别解码器，用于预测潜在向量对应细胞使用的重排v，j基因类型，确保潜在表示蕴含基因重排信息；另一部分是Transformer的解码器，通过该解码器重构细胞的原始氨基酸序列，确保潜在表示蕴含氨基酸序列信息。其中，类别解码器均使用单层全连接层，输出为预测的v，j基因各种类别的嵌入表示概率，对应的损失函数为
(4)
式中：vi为真实的v基因类别信息，为预测的基因类别概率，M1为v基因类别总数。类似地，
(5)
式中：ji为真实j基因的类别信息，为预测的基因类别概率，M2为j基因类别总数。对于氨基酸序列重构部分，使用交叉熵损失
(6)
式中：N为批次大小，T为序列长度。因此，总损失函数为Lchain=α1Lrecon+β1Lv+γ1Lj，其中，α1，β1，γ1为可调的权重参数，本文初始设置为α1=β1=γ1=1。
2.2 基于高斯混合潜在分布的免疫组双链整合变分自编码器
对于基本的变分自编码器架构（variational autoencoder，VAE），一般假定其潜在空间为标准高斯分布，然而，对于一些复杂的数据分布，一个普通的高斯分布难以拟合复杂的数据分布情况。对于免疫组数据而言，由于免疫克隆型的高度特异性，其类别众多，分布离散。因此，采用高斯混合模型（Gaussian Mixture Model，GMM）修正的变分自编码器，对2.1节得到的双链表示进行整合，进而得到免疫组库的潜在表示。
以αβ T细胞的处理为例，γδ T细胞的处理是完全类似的。设在2.1节中得到的双链数值嵌入分别为z1，z2，将二者拼接起来得zconcat=（z1，z2），目的是将其通过GMM变分自编码器得到双链的融合潜在表示。由于高斯混合分布具有混合系数、均值向量和协方差矩阵3个参数，其的组合决定了数据分布。虽然在网络学习中，这3个参数可以通过数据学习得到更新，但是，一个合理的初始化可以帮助参数更好地逼近分布。因此，本研究采用预训练与高斯混合训练通过EM算法结合的方法进行网络训练。
具体地，将拼接向量zconcat首先通过一个以标准高斯分布为潜在空间的变分自编码器，其损失函数为常见的向量重构损失和标准高斯KL散度，经过预训练之后，可以得到一个初始的数据潜在空间分布。在得到潜在表示后，给定聚类中心数目，通过EM算法计算每个点属于各类中心的后验概率，之后更新聚类中心，逐步迭代下去得到潜在空间作为高斯混合分布的初始参数μk，Σk，πk。
得到预训练的初始参数之后，利用其作为GMM-VAE的训练起点，进行双链整合训练。具体地，将zconcat=（z1，z2）作为编码器原始输入，映射到高斯混合模型的潜在空间，之后再利用解码器重构输入向量，得到重构向量，，对于重构项，采用均方误差作为损失函数，即
(7)
(8)
其中N为批次大小。
对于变分自编码器的KL散度项，使用编码分布与高斯混合先验的KL散度进行训练。具体地，
其中而，由其显示表达有
(9)
而另一项没有闭解，故采用蒙特卡洛估计近似逼近这个积分式，即在大量采样后，用求和式逼近这个积分。其中，d为空间维数，总损失函数为L=α2LKLGMM+β2L1+γ2L2， α2，β2，γ2为可调的参数。一般地，由于大部分T细胞的特异性体现在β链上，初始设置为α2=0.1，β2=0.2，γ2=0.7。
2.3 基于概率生成模型的双模态整合
在本模块中，采用变分自编码器的思路，通过编码器将数据映射到潜在空间，之后通过重构进行训练。这里设qφ（z|x）是近似后验分布，pθ（x|z）是解码器重构，为了考虑双模态，设基因表达信息为Xgex，免疫组信息为Ximm，潜在表示为Z，则全概率模型为
(10)
此模块最终的优化目标函数为
(11)
其中x=（xgex，ximm）。
对于基因表达信息，为了拟合基因表达的分布以及考虑单细胞转录组数据具有技术零值的特点，使用零膨胀负二项分布（ZINB）进行建模生成，ZINB分布的概率密度函数为
(12)
式中：π为零膨胀概率，NB为负二项分布，μ为负二项分布均值，φ为离散度参数，这些参数通过解码器得到，则LZINB=-log PZINB（x|μ，φ）。
对于免疫组数据，采用高斯分布的连续解码，对应的损失函数为
(13)
KL散度项为。最终的此模块损失函数为
(14)
其中ω=（ωgex，ωimm，ωKL）是可调的权重参数，本研究着重体现免疫组库数据对转录组数据分布的扰动效果，故初始设置为ω=（1，0.2，0.01）。
3 实验及结果分析
3.1 数据说明
配对的人类T细胞样本可以通过网站直接获取（https: //www.10xgenomics.com/），其中的免疫组数据集均按照AIRR标准格式进行存储，并且可以通过开源python模块scirpy直接进行读取。不同的生物条件下小鼠的单细胞转录组测序与免疫组库测序技术来源于NCBI数据库GSE275982号的公开测序数据集。本研究使用方法代码公开于https: //github.com/Jinsl-lab/scRTIA。
由于测序技术限制，得到的原始基因表达数据与免疫组库数据均受到无效数据与噪声等因素的干扰。在使用数据进行算法验证之前，需要首先对原始数据进行预处理。针对免疫组数据而言，其测序数据为基因类别与CDR3区序列，其中的每个信息均在本研究的整合方法中起到重要作用，因此，本研究只保留那些具有完整CDR3区序列、α链v，j基因类别、β链v，d，j基因类别的单链，将其他有缺失信息的单链删除。除此之外，只保留全长（fulllength=TRUE），可生成单链（productive=TRUE）的链。在将TCR拼成一个细胞时，有时一个细胞只能找到一个单链，而有的细胞可能会在某些链上出现冗余，这两种情况对训练均有影响，并且这些细胞在数据集中均有一定比例的分布，本研究删除掉这些链型，只保留可以正好组成一个单细胞的匹配双链，及只保留占数据集中多数的单链配对T细胞，之后通过神经网络进行方法训练。
针对单细胞转录组数据，其格式可以是转录组测序的.h5文件或者.mtx格式的基因表达矩阵。在预处理部分，首先，通过控制单细胞中的线粒体比例，删去线粒体含量过高、极有可能衰老凋亡的细胞，之后删除在小于3个细胞中表达的基因，并且去除掉表达基因数小于200的细胞，因为其可能是测序得到的空液滴或者死亡细胞。
3.2 实验结果
3.2.1 scRTIA可以充分提取免疫组库信息
由于本方法首先对免疫组库信息进行了提取与整合，对于免疫组库信息的数值嵌入应该整体包含大部分免疫组库信息，否则在方法的下游处理时，准确率会被大幅影响。
使用来自4个不同捐献者的CD8+T细胞的10x数据集进行测试，如图2（a）~（c）所示，随着训练轮次的增加直到神经网络触发早停，方法在免疫组数据的CDR3序列重构，v基因、j基因类别预测的准确率方面，均随训练轮次的增加而增加，说明通过Transformer架构得到的潜在表示在一定程度上保留了主要的免疫组库信息。以β链为例，由图2（d）、（e）可知，含量较多的v，j基因在其潜在空间上的分布。在潜在空间上，由于v基因的高变性，其分布除了某些基因型较为聚集之外，其他的基因型相对混合分布，而j基因相对保守。因此，在潜在空间上，按照不同基因型相对更聚集。
图2scRTIA可以提取免疫组库信息
Fig.2scRTIA extraction of immune-repertoire features
此外，本研究发现，通过添加v，j基因的类别判别损失函数，在训练结果大体相似的情况下，可以显著地减少神经网络的训练时间。如表1所示，在不添加v，j基因的判别损失时，训练时间均在11~13 min，而在添加了v，j基因的判别损失之后，训练时间均显著缩短至2~6 min。对此，在生物学上可以进行一些解释，即CDR3区的氨基酸由v，d，j基因编码产生，CDR3区序列由密码子表翻译得到，所不同的只是随机插入和缺失，由此，v，j基因的信息在一定程度上均体现了CDR3区的序列信息，故加速了神经网络的训练速度。
表1是否添加判别项的训练时间比较
Tab.1Comparison of training time with and without discriminative item
3.2.2 scRTIA可以在潜在空间中保持克隆型的类别聚集
在T细胞中，克隆型指具有相同或者相似CDR3氨基酸序列的T细胞，或者具有相同v，d，j重排方法的T细胞群体。在生物学中，分析T细胞的克隆型对分析免疫应答[12]、针对不同疾病的免疫特点[13]、设计肿瘤免疫治疗[14]、分析免疫变化[15]等重要的生物学问题起到关键作用。因此，在整合后的潜在表示中，保持T细胞的克隆型分布，对生物学分析十分重要。
事实上，由于T细胞的免疫组库数据，与其对应的转录组数据并不是强相关的组学数据，即通过其中某一个组学数据去较为准确地推断另一个组学数据的分布情况是基本不可能的，所以，如图3（a）所示，考虑了含量最多的前十位克隆型的分布情况，在转录组数据的潜在分布上，不同的克隆型虽然大体上体现出某种聚集效果，但是，总的来看还是难以体现免疫组的特异性。
在将数据集中T细胞的免疫组库信息进行提取整合，得到数据集在免疫组库下的潜在表示空间之后，由图3（b）可以看出，相同克隆型的细胞在此潜在空间上出现了部分聚集效应，由于T细胞的分布高度离散，其分布大多聚成小但密集的类别。
再按照方法流程，将T细胞的免疫组库信息与转录组信息结合，在之后的整合潜在空间中，由图3（c）可以看出，相同的T细胞克隆型基本被分配到一个聚集的类别中，而不是像免疫组库数据分布那样离散。事实上，T细胞的分布受到克隆型种类的影响，在潜在空间中比较离散，而转录组数据分布通常比较连续。因此，二者整合后的数据分布会体现出两个模态的共同特点，即有的部分呈现连续的分布状态，而有的离散小类别在潜在空间各自独立分布。
3.2.3 scRTIA可以部分利用免疫特性划分新细胞类别
对于单细胞转录组，其细胞聚类过程完全依靠于每一个细胞的基因表达信息，这种分类虽然可以体现细胞的表达特性，但对于其他模态的信息利用不足。通过整合免疫组信息，可以在体现转录组特性的同时，依据其免疫组特性进行分类，进而综合研究细胞的特性。
由图4（a）可以看出，在潜在空间聚类之后聚出了不同类别，这些类别在原始转录组空间中虽有某些对应的分类趋势，但是无法观察到明显的分类。其中，不同类别细胞的类型混杂现象比较严重，这是因为转录组数据几乎不能体现过多的免疫组信息，或者说这两个模态是弱关联模态。
为了探究不同类别的免疫特性，选择了CDR3 区氨基酸链长度这个特性。对于TCR，CDR3区氨基酸链的长度是一个重要的指标，而通过山脊图对潜在空间划分的不同类别细胞对应的免疫组库CDR3链的长度进行分布可视化，如图4（b）所示，不同的类别中，其CDR3序列长度分布差别明显。
图3scRTIA可以在潜在空间中保持克隆型的类别聚集
Fig.3scRTIA preserves clonal-type clustering in the latent space
图4scRTIA可以部分利用免疫特性划分新细胞类别
Fig.4scRTIA enables partial segregation of novel cell types using immune-repertoire features
3.2.4 scRTIA可以检测异质性数据集中的特异性细胞群体
通过前两节发现，scRTIA可以对单一T细胞群体中的免疫信息进行充分提取和整合。对于不同实验，不同条件下得到的生物样本数据集，识别其中具有某些免疫和转录组特性的细胞亚群也具有重要的生物学意义。本研究使用年老脊髓未损伤小鼠、年老脊髓损伤小鼠和年轻脊髓损伤28 d后的小鼠3种不同生物学样本构成的数据集进行分析，其中，每组样本进行了4次生物学重复，总共是12个不同的数据集。本研究将其组合到一起进行综合分析。
由图5（a）可以看出，只利用原始转录组数据，数据按照其样本分组以及多次生物学重复分成了较分散并且无联系的孤立聚集分布。在进行整合后，如图5（b）所示，不同条件样本细胞分布更加聚集，并且同一条件下的细胞没有严格按照其生物学重复进行聚集，而是融合了免疫组库信息，其独立出来的分组带有某些免疫特征。
由图5（c）~（e）中圈出的细胞亚群显示，这些细胞在T细胞的v，，j基因重排中出现了不同基因型的重排富集，对于出现较多的TRBV19、TRAV3N-3、TRBV10等基因在某离散类别中提取出来。
图5scRTIA可以检测异质性数据集中的特异性细胞群体
Fig.5scRTIA detection of specific cell populations in heterogeneous datasets
4 结论
1）采用Transformer架构对免疫组库数据进行数值嵌入，在保留关键生物学特征的基础上，使用对TCR双链的表示学习。通过基于高斯混合分布的变分自编码器将两条单链信息融合成免疫特征，最终使用概率模型进行双模态整合。
2）在方法构建上，对之前相关研究中将T细胞整体看作研究起点的思路进行细化，将分析和整合的起点放在每一个TCR单链上，通过对配对的TCR单链进行整合得到免疫组表示，可以从更精细的角度对不同的T细胞亚群进行分类与刻画。
3）在实验验证上，分别使用同质性和异质性T细胞数据集进行方法验证，均得到相关的生物学结论，实现有意义的双模态整合，得到某些既保持转录组特性，又融合免疫特性的T细胞亚群。
4）本方法也有一定的局限性，使用时需要注意数据输入的格式。免疫组格式应该具有完整的配对双链信息，而转录组数据应具有.h5格式或者.mtx格式，其他格式数据应进行格式转换之后进行输入。本方法在数据预处理方面具有一定的局限性，选择删除了全部的非配对细胞，对于数据的充分使用上具有一定可提升空间，此外，本方法未针对性地处理多数据集之间的批次效应，对于部分多样本数据集的处理有一定的局限性。
