αβT细胞受体库的综合分析揭示胸腺选择的关键特征及其在免疫治疗中的应用潜力

本研究通过系统分析T细胞受体（TCR）的理化特性与序列特征，揭示了胸腺选择过程中CDR3环长度、疏水性及残基大小等关键因素对T细胞命运的决定性作用，为自身免疫病机制研究和TCR免疫疗法优化提供了重要理论依据。

　　

胸腺选择是形成T细胞池的关键过程，这些T细胞能够通过其T细胞受体（TCR）有效区分自我与非我，从而发展适应性免疫。本研究分析了TCR的多种理化特性和序列特征如何影响成功通过选择的机会。

在全局尺度上，研究发现了基于CDR3环长度、疏水性和残基大小的选择概率差异，这些差异取决于可变基因和TCR链背景。还观察到α和β链的N-糖基化位点和其他短序列基序显著减少。在局部尺度上，使用专用统计和机器学习方法结合V(D)J重排过程的概率模型，推断出在选择过程中富集或耗竭的CDR3区域模式。虽然含多甘氨酸模式的丰度可以提高所选TCR中CDR3的灵活性，但负选择诱导的TCR库中的“空洞”可能与（N）-多样性（D）-N区域（NDN）区域中的精氨酸有关。相应的模式被存储在在线数据库中。研究证明了TCR序列组成如何影响胸腺选择过程中的谱系定向。结构建模显示，具有“平坦”和“凸起”CDR3环的TCR分别更可能将T细胞定向为CD4⁺和CD8⁺谱系。最后，强调了个体MHC单倍型对选择过程的影响，表明这些“空洞”可能是供体特异性的。结果可进一步应用于识别供体库中潜在的自身反应性TCR，并有助于免疫治疗中的TCR选择。

引言

淋巴祖细胞从骨髓迁移到胸腺，通过称为可变-（多样性）-连接基因重排的过程发育T细胞受体（TCR，α和β链的异二聚体），这些受体经过选择以成为成熟T细胞。V(D)J基因重排过程涉及两个主要步骤：1）从相应基因座中选择三个（TCR β链的V、D和J）或两个（TCR α链的V和J）基因等位基因并重组在一起；2）在基因末端随机删除碱基，并在基因之间的连接位点添加非模板（N）核苷酸以增加序列多样性。这些过程导致TCR的巨大多样性，估计约为10¹⁹，比人体T细胞总数（约10¹¹）高出几个数量级。在释放到血液之前，T细胞必须经过胸腺选择以确保抗原识别的功能（阳性选择）和对健康细胞的耐受性，避免对自身反应（阴性选择）。这是一个多阶段过程，从CD4^-CD8^-双阴性（DN）T细胞重排其TCR β链并用pre-TCRα验证其功能开始。接下来，它们重排α链，形成成熟TCR，并转变为CD4⁺CD8⁺双阳性（DP）表型，经历阳性和阴性选择。最后，DP T细胞分化为CD4⁺或CD8⁺单阳性（SP）T细胞并从胸腺迁出。T细胞受体的互补决定区3（CDR3）由V(D)J连接编码并直接与MHC呈递的肽相互作用，高度多样化，并常作为研究中TCR序列的代理。本研究专注于CDR3序列以探索TCR库中的胸腺选择过程。现代TCR库研究使用高通量测序（Rep-Seq或AIRR-seq）分析来自各种生物样本的数百万个TCR序列。鉴于已发表的人类预选择胸腺细胞库数据有限，我们和其他人利用V(D)J重排的理论模型来模拟此类库。我们还依赖由单个mRNA分子支持的“单例”T细胞克隆型主要代表初始T细胞的假设来模拟选择后的初始TCR库。本研究的一般方法是通过模型数据和常规外周血单核细胞（PBMC）AIRR-seq数据比较选择前后的TCR库，并使用分选的DP和SP胸腺细胞验证发现。通过分析可在不同水平上影响选择过程的序列特征，探索了选择前后TCR的局部和全局库结构：V/J基因使用、初级CDR3序列的氨基酸组成和理化特性、库的CDR3 k-mer谱和突出序列基序。为了完善比较分析并去除源自V(D)J重排内在随机性的噪声，我们利用基于TCRNET方法的TCR序列簇富集策略，允许我们在具有内在复杂结构的全身TCR库中检测选择基序。先前的研究已经检查了胸腺选择中的TCR特征。Lu等人报告了氨基酸使用的变化，疏水和带正电荷氨基酸以及半胱氨酸的频率减少。Stadinski等人发现特定位置具有疏水残基的TCR易于交叉反应，降低其生存机会。其他研究显示了胸腺、淋巴结和脾脏的库差异。这些特征也在Isacchini等人的研究中通过库建模重新讨论，但该研究得出结论，选择特征在局部和全局尺度上相似，声称没有禁止的TCR序列和选择基序。我们使用我们的框架回顾了上述特征并得出类似结论，但通过显示选择后局部富集和耗竭的TCR模式可以是公共和供体/HLA特异性的，扩展了先前发表的发现。谱系定向与TCR的联系被广泛研究，揭示了CD4⁺和CD8⁺、CXCR3⁺和CXCR3^-初始CD8⁺以及辅助T细胞亚群之间不同的CDR3特征。这里我们通过单细胞数据更详细地探索它，显示选择基序与某些表型相关。此外，我们显示CDR3环的结构对于CD4和CD8相关的TCR基序是不同的。我们还考虑了供体HLA单倍型在选择中的作用，并证明了等位基因特异性差异如何贡献局部库特征。这种效应先前在小鼠CD4⁺库中有报道。还显示MHC背景通常塑造T细胞库。涉及双胞胎TCR库的类似方法被Tanno等人用来揭示遗传因素，特别是MHC等位基因对TCR库的影响。

材料与方法

选择后T细胞库

使用先前发表的（PBMC）TCR库测序数据用于TCRα和TCRβ链。对于TCRβ，从Emerson等人的数据集的HIP队列中随机选择10个CMV供体库用于TCRβ分析（数据集中的样本ID从1到10）。对于TCRα，使用Heikkil?等人中所有批量TCRα PBMC数据（10个样本）。仅由单读段支持的单例用于后续分析。TCRβ样本包含1,147,250个TCR序列，TCRα样本包含1,582,774个。未控制生成偏差，因为数据集包含来自多个个体的序列，这掩盖了个体生成偏差的影响。此外，生成偏差在选择过程中被印记。

分选库

分选的CD4⁺和CD8⁺初始（选择后）T细胞测序库取自Qi等人；此数据仅适用于TCRβ链测序。TCR库DP和CD8⁺（SP）胸腺细胞（预选择）取自最近Quiniou等人的研究。该数据集中所有九名供体的初始库合并在一起，产生1,346,776个克隆用于CD8⁺初始细胞和1,599,217个克隆用于CD4⁺初始细胞。对于这些数据集，未使用读段计数信息（即所有克隆型假定为单例）以避免潜在扩增偏差并使其与其他数据集兼容。

HLA匹配和不匹配库（双胞胎研究）

包含三对同卵双胞胎PBMC TCRβ库测序的数据集取自Pogorelyy等人的“精确追踪疫苗应答T细胞克隆揭示同卵双胞胎中的收敛和个性化反应”。处理前采样的疫苗接种前和第0天库用作每名双胞胎的生物重复，仅单例包含在分析中。供体的HLA等位基因列在补充表中。包含来自三对同卵双胞胎PBMC TCRα和β库的额外数据集取自Zvyagin等人的“高通量测序揭示同卵双胞胎TCR库的独特特性”。来自两对同卵双胞胎的分选初始CD4⁺ T细胞TCRβ库也从Kasatskaya等人的“功能特化人类CD4⁺ T细胞亚群表达物理化学 distinct TCRs”获得。数据集的摘要统计报告在补充表中。

模拟预选择V(D)J重排

预选择TCRα和β库基于V(D)J重排过程的理论概率模型使用OLGA软件（v1.2.4）模拟，如先前所述。软件以默认运行时参数和模型概率执行，随机种子设置为100，为每个TCR链生成10⁷个随机重排样本。

单细胞数据分析

取自Lindeboom等人的具有指定细胞类型的单细胞数据集总计178307个细胞来自88名健康患者的PBMC样本。细胞类型注释使用CellTypist执行，如原始研究中所指定。使用Fisher精确检验测试特定细胞类型上的簇丰度。

TCR氨基酸序列特征和基序分析

使用内部脚本进行TCR序列的基本特征，如V/J基因使用、单氨基酸频率、k-mer（k=3）频率和CDR3区域的理化特性，如先前所述。k-mer集未基于其在CDR3中的相对位置过滤。使用“peptides”python包（v0.3.2）计算Kidera因子（描述多肽理化性质中最方差的關鍵氨基酸特征）、电荷和疏水性。最信息丰富的Kidera因子，Kideras 2、4、6和8，分别对应于侧链大小、疏水性、部分比容和α区域中的出现。使用T检验比较Kidera因子值；使用Cohen's d作为效应大小估计。Kidera因子未进行z分数标准化，因为与OLGA生成数据比较的数据集使用与用于训练OLGA模型的数据相同的协议获得，并且胸腺细胞数据集中的所有样本在单一批次中获得。使用VDJtools（v1.2.1）中实现的TCRNET算法执行TCR CDR3序列基序的检测，如先前所述。该方法将感兴趣的TCR序列定义为那些放置在CDR3序列相似性图中更密集区域的序列，与对照（通常产生假设V(D)J重排模型无选择压力）数据集相比：比较1-汉明距离邻居的数量与调整样本和对照大小的预期邻居数量，以产生富集分数和基于二项式近似的P值。注意，为了推断由负胸腺选择耗竭的TCR序列簇，我们简单地交换“背景”（对照）和“前景”（我们感兴趣的样本），即我们搜索在预选择数据中相比选择后富集的TCR。为了产生有代表性的TCR簇集，我们基于富集P值选择前10,000个邻居富集的CDR3序列。通过选择连接序列的图的连通组件来聚类所选序列，这些序列通过单个氨基酸替换不同。使用logomaker包（v0.8）可视化所选簇的基序。随后分析前五个最大簇并根据其簇大小排名编号。簇数量是随机选择的。额外使用SoNNia（v0.2.3）模型评估氨基酸出现概率的差异，特别是它们在选择后与预选择库中的位置。为了计算序列特征的边际概率，我们使用选择后和预选择数据集对（例如选择后TCRβ和OLGA生成的TCRβ数据）训练SoNNia模型。选择后数据集用作模型推断的数据，预选择数据集用作从生成分布采样的数据。每个数据集使用“Processing”类中的方法处理，然后传递到具有适当“pgen_model”参数的SoNNia模型。每个模型训练50个周期，批量大小等于10⁴。使用内置绘图功能评估模型，如SoNNia教程中所示。

双胞胎数据集的比较分析

为了识别双胞胎数据集中阳性和阴性选择的CDR3簇，我们将每个双胞胎样本二次采样到306,553个CDR3（最小库的大小）并将所有样本合并在一起。使用随机阈值选择显著富集（与模拟序列作为对照相比，样本池的log2倍数变化 > 2和 -log10 p > 12）和耗竭（与样本池作为对照相比，模拟序列的log2倍数变化 > 1和 -log10 p > 12）的胸腺选择后簇。使用包含超过10个序列的簇进行进一步分析。接下来，通过计算簇频率之间的Jensen-Shannon发散来估计阳性和阴性选择簇之间的相似性，定义为给定样本中存在的给定簇的克隆型数量。

结构分析

使用TCRmodel网络工具建模感兴趣TCR的CDR3环结构并使用Pymol（版本2.3.0）处理。使用内部“mir”软件包注释生成的PDB文件。由于在大多数数据集中CDR3α未知，我们使用通用CAGGSSNTGKLIF（TRAV27，TRAJ37）序列，该序列是Heikkila数据集中最常观察到的变体，作为虚拟TCRα序列。通过将PCA应用于Cα原子坐标可视化TCR CDR3骨架。从VDJdb获得TCR:pMHC复合物的实验结构。我们仅选择具有可用人类TCR PDB ID的记录。结构需要包含TCRα和TCRβ以及pMHC复合物。此外，至少一个CDR3β残基必须在pMHC复合物中与肽 within 5 ?以确保TCR β链与pMHC复合物之间的直接接触。总共分析了154个结构。使用Biopython（v. 1.85）Python包计算这些结构中的二面角。

代码可用性

本研究中使用的所有代码可在https://github.com/LuppovDaniil/Thymic_selection_notebooks获得（Python版本3.11.5.，R版本4.1.2）。

结果

比较选择前后TCR氨基酸序列

从胸腺分选的DP T细胞的AIRR-seq数据可用于探索在阳性和阴性选择之前存在的重排TCRα和β序列的初始空间，类似于Quiniou等人的工作。最近的研究证明概率模型可以准确复制V(D)J重排空间的结构，生成具有氨基酸组成和频率类似于体内产生的TCR序列。这里我们使用Quiniou等人和模型数据集作为预选择库。有几种方法可以获得尚未经历强抗原暴露并获得记忆表型的选择后T细胞的TCR库。可以分选和测序初始CD4⁺和CD8⁺ T细胞，如Qi等人所做，或使用SP胸腺细胞，如Quiniou等人所做。或者，可以从未分选的PBMC AIRR-seq数据中选择仅检测一次的单例T细胞克隆型，因为它们主要代表初始T细胞。在本研究中，我们使用Qi等人和Quiniou等人的数据集，并从Emerson等人数据集中随机选择的10个样本中选择单例作为选择后但任何抗原暴露或后续扩增之前的TCR库。这些数据集的详细描述和比较分析在方法部分给出。TCRβ CDR3氨基酸频率分析显示，与V(D)J重排模型预期的相比，特定残基在选择后显著减少。带正电荷和物理上大的氨基酸，如精氨酸、组氨酸和赖氨酸，可能由于强抗原结合或抗原识别中的空间位阻而减少，而脯氨酸和半胱氨酸可能由于它们对TCR结构的影响而被阴性选择。这些发现在胸腺细胞数据中得到普遍证实。3-mer频率分析显示，具有半胱氨酸的3-mer生存选择的机会较少，与单氨基酸水平的观察一致。与N-糖基化位点相关的NX[S,T]基序的3-mer更不可能生存选择，NPT和NPS受影响最小，与这些位点无糖基化一致。我们还分析了磺化位点基序，这些在Pospelova等人的抗体工作中描述，如DDD、DDY、YYY和EEE。仅对DDD基序观察到切实的减少效应。糖基化位点和半胱氨酸相关效应以及推定的硫酸化位点缺乏效应在胸腺细胞数据中观察到。此外，我们评估了糖基化位点在CDR3序列中的位置，并发现对于CDR3β，大多数（不同数据集中85%或更多）糖基化位点位于NDN段，突出了它们的体细胞起源。然而，对于CDR3α，仅有44-53%的这些位点位于N段。为了描述TCRβ CDR3s理化特性的变化，我们利用Kidera因子。这些因子代表通过降维获得的氨基酸的关键物理特性。我们独立比较每个VJ对以减轻由V和J基因选择引起的偏差。受选择负面影响的Kidera因子是Kidera 2，它决定侧链大小，Kidera 4，它与疏水性逆相关，和Kidera 6，它决定部分比容。这些结果表明，具有物理小疏水CDR3的TCRβs有更大机会通过选择。选择后唯一增加的因子是Kidera 8（在α-螺旋结构区域中的出现）。值得注意的是，特定Kidera因子内的选择效应在所有VJ对中共享。DP和SP胸腺细胞的相同分析反映了上述结果，然而，Kidera因子6的效应不太明显。TCRβ CDR3物理特性分析显示选择后库电荷减少，与观察到的带正电荷氨基酸频率减少一致，长度的 winsorizing（太短和太长的序列都较少机会通过选择），和疏水性增加。除长度经历缩短而非 winsorizing 外，在胸腺细胞数据中检测到相同效应。接下来，我们识别了在初始细胞中相比模型生成样本富集的功能性TCRβ CDR3s簇，反之亦然。初始和模型生成样本中最大的富集簇呈现在图中。模型富集簇可以解释为倾向于不通过选择的序列，初始样本富集簇是可能存活胸腺选择的序列。为方便起见，此处和补充材料中使用的术语（富集和耗竭预选择与选择后 versus 阳性和阴性选择）在图中解释。大多数选择后富集簇在TCRβ CDR3中间含多甘氨酸序列，这已知是多肽链中最灵活的之一。此外，k-mer分析证明GGG 3-mer以未改变频率通过选择，证实上述发现。预选择富集簇中的CDR3β经常包含精氨酸和脯氨酸（可能影响CDR3环结构的残基）并在CDR3开始处显示CASS共有序列的偏差。胸腺细胞数据的富集簇分析结果大多与生成数据中获得的结果相似。为了获得选择前后氨基酸在CDR3特定位置出现的概率，我们使用SoNNia软件。通常，我们观察到与富集簇分析几乎相同的趋势。基因使用分析揭示了选择对更频繁基因的偏好（或“富者愈富效应”）的模糊图景，无论是在生成还是实验数据上。有趣的是，虽然未观察到基因频率的“富者愈富效应”，但在DP和SP胸腺细胞的CDR3β生成概率（pgens）中检测到，通过OLGA模型计算（p < 10^-5，Mann Whitney检验，SP和DP之间中位数log2 Pgen差异为2.24）。与TCRβ相同的TCRα CDR3s分析揭示了共同的选择效应和两条链之间的显著差异。对于此分析，使用来自PBMC的CDR3单例子样本 alongside 胸腺细胞CDR3α库和模型生成数据。在观察到的差异中，以下点值得注意：（i）没有清晰的选择偏好朝向单个氨基酸的趋势；（ii）TCRα VJ对的Kidera因子变化没有一致方向；（iii）选择后库疏水性不变；（iv）选择后富集和耗竭簇的结构更复杂，DS和NY 2-mer在胸腺细胞富集簇中普遍但不在初始数据中，SS 2-mer在模型和DP胸腺细胞富集簇中更丰富。值得注意的是，α和β链共同的选择效应在TCRα case中不太明显。此外，我们将选择后富集和耗竭的最大簇的TCRs和模式存储在https://github.com/LuppovDaniil/thymic_selection_motifs_database的数据库中。

单细胞分析揭示CDR3依赖性分化

为了追踪胸腺选择过程中TCR驱动谱系定向的迹象，我们利用来自88名健康供体的单细胞测序PBMC数据，总计178307个不同谱系的T细胞。我们旨在证明富集和耗竭簇与特定T细胞表型相关。我们将这些簇中的CDR3s与单细胞TCR库相交，并测试这些簇与每个特定细胞表型的联系。通常，获得了预期结果。选择后耗竭簇中的TCRs在单细胞数据中的丰度比富集簇中的低10至20倍，支持我们识别选择不利TCRs的能力。选择后富集簇主要与初始表型相关。SP胸腺细胞簇主要由CD8⁺初始细胞代表，这是预期的，因为我们研究中的SP胸腺细胞是CD8⁺谱系。我们还发现CD8⁺细胞的TCRβ比CD4⁺细胞的TCRβ更容易邻居富集。在CD8⁺ Recent Thymic Emigrants（RTE）亚型中富集的簇与其他在CD8⁺初始细胞中富集的簇相比突出。这样的观察可能是库外周选择效应的迹象。SP胸腺细胞的簇4也特别感兴趣，因为它在CD8⁺ RTE和Tregs中都表现出丰度。此外，我们观察到来自富集SP胸腺细胞簇的TCRα在MAIT细胞亚型中的丰度，这再次预期，因为MAIT细胞以其“半不变”TCRα为特征。此外，我们评估了我们关于PBMC数据中单例主要初始起源的假设。我们使用单细胞数据推断TCRβ和TCRα数据集中单例的表型，并发现TCRβ单例确实在初始亚组中富集（p < 10^-4，Fisher精确检验）；然而，我们未能检测到TCRα的相同富集。尽管存在这种特殊的不一致，我们仍然相信我们的假设是正确的，因为单细胞数据中记忆细胞的份额比初始细胞的份额低10倍，并且大多数T细胞曾经是初始的。为了巩固我们的发现，我们使用包含单细胞基因表达和TCR测序数据的标准10X Genomics数据集分析SP胸腺细胞簇。通常，获得类似结果：（i）携带来自CD8⁺胸腺细胞富集簇的TCRβ的细胞主要是CD8⁺表型；和（ii）TCRα富集簇主要与MAIT细胞相关。这些结果证明了胸腺选择过程中TCR驱动谱系定向，无论是对于MAIT细胞的充分表征案例还是CD8⁺细胞的更复杂案例。

结构分析确认CD4⁺和CD8⁺ T细胞的CDR3特征

接下来，我们比较了取自Qi等人研究的CD4⁺和CD8⁺初始细胞TCR库。此比较最显著的结果来自CD4⁺或CD8⁺ CDR3βs中富集的簇。我们使用TCRmodel2分析每组最大簇中最普遍CDR3的结构，并发现富集CDR3β CD8⁺和CDR3β CD4⁺结构之间的凸度差异。似乎CD4⁺富集簇结构平坦，而CD8⁺簇更凸。CD8⁺ CDR3β接触区域（排除前四个和最后五个残基）的质量中心位于比CD4⁺ CDR3β质量中心离环中心更远的位置。然而，在我们18个结构的样本中，这种差异不显著。这一事实可能由MHC I类和MHC II类凹槽中肽的构象解释。MHC I类，由CD8⁺细胞识别，倾向于呈现具有中间凸起的肽，而MHC II类，由CD4⁺细胞识别，倾向于呈现平坦肽。我们使用存储在VDJdb中的数据进一步调查了CD4⁺和CD8⁺ CDR3βs之间的结构差异。我们采取具有可用TCR:pMHC复合物的人类TCRs，这些复合物包含CDR3β和肽之间的接触。总共分析了154个结构。我们通过规范Ramachandran图可视化这些结构中所有CDR3β二面角，并发现图上的一个 distinct 区域被CD8⁺ CDR3β残基 abundantly 居住。Ramachandran图上此区域的结构构象主要可用于甘氨酸。此区域具有残基的大多数结构仅由一个残基代表。具有 such 构象的CD8⁺ TCRs比CD4⁺ TCRs频繁2.83倍（p = 0.017，Fisher精确检验）。Ramachandran图上右下角的此区域在Hollingsworth和Karplus的工作中被称为P_II’。P_II’被视为P_II（更常称为Polyproline-II Helix）的镜像区域，已知最大化多肽链熵并将所有能够氢键的主链原子暴露于水。此外，此结构区域经常形成蛋白质结合基序，这与TCR识别抗原相关。因此，我们期望此CD8⁺ specific 构象相对凸起，如我们上面所述。此外，讨论区域中的大多数残基是甘氨酸，它们与肽直接接触。具有相应PDB ID的Ramachandran图数据可在补充表中获得。

HLA等位基因单倍型影响选择

接下来，我们考虑了HLA等位基因对胸腺选择的影响。HLA等位基因是否影响选择仍然是一个争论的话题。我们使用三对20至23岁双胞胎的CDR3s批量测序数据来解决这个问题。此数据包含三对双胞胎（指定为S、P和Q）在两个副本中的TCRβ库：疫苗接种日射击前（第0天）和疫苗接种日前一天（预日）。我们将他们的初始CDR3簇与模型背景库比较。一方面，我们期望他们的富集初始CDR3簇彼此密切相关，因为来自同一对的双胞胎具有相同的HLA等位基因。另一方面，我们预期来自不同对的双胞胎在其初始TCR库背景下彼此不同。我们从每个双胞胎样本中提取相同数量的CDR3s，然后将每个样本中的富集序列聚类在一起。完成此操作后，我们获得了一些由来自不同样本的序列组成的簇。据推测，每个簇包含功能上彼此相似的CDR3s（识别MHC-肽复合物中的相似肽）。因此，我们预期来自相对双胞胎对的CDR3s落入相同簇。为了验证这一假设，我们分析了特定双胞胎CDR3s在获得的簇中的比例。在这种表示中，每个样本可以被视为出现在特定簇中的离散概率