High-throughput sequencing of the T-cell receptor repertoire: pitfalls and opportunities

Abstract

T细胞特异性是由T细胞受体决定的，T细胞受体是一种异二聚体蛋白，由一组由不精确的体细胞基因重组产生的极为多样化的基因编码。大规模并行高通量测序允许数百万不同的T细胞受体基因从一个单一的血液或组织样本的特点。然而，免疫系统的异常异质性对随后的数据分析提出了重大挑战。我们概述了处理免疫组库（repertoire）数据的主要步骤，考虑到原始序列文件的低级处理和高级算法，寻求提取生物或病理信息。最新一代的生物信息学工具允许数以百万计的DNA序列被准确和快速地分配到各自的可变V和J基因片段，并重建一个几乎没有错误的表示非模板添加和删除发生。高级处理可以测量不同样本中全基因组的多样性，量化V和J的使用情况，并识别私有和公共t细胞受体。最后，我们讨论了将t细胞受体序列与功能，特别是与抗原识别联系起来的主要挑战。目前正在开发复杂的机器学习算法，可以将单个t细胞受体的反常退化和交叉反应性与整个t细胞免疫反应的特异性结合起来。计算分析将为解开t细胞受体的潜能提供钥匙，从而深入了解适应性免疫系统的基础生物学，并为疾病提供强有力的生物标志物。

Introduction

下颌脊椎动物的适应性免疫系统利用不精确的体细胞DNA重组，在B细胞(BCR)和t细胞(TCR)上产生丰富多样的抗原特异性受体阵列。关于可变抗原受体多样性产生的机制已经进行了非常详细的研究(如文献[1,2]，图1中的图解)。简而言之，编码每个链的位点由多个基因片段(' minigenes ')组成，其中包括V、D和J基因。在淋巴细胞发育过程中基因组DNA的重组导致每一种可用基因片段的物理连接，并切除中间的DNA。由多个V、J和(在某些链中)D基因组成的组合多样性可以“选择”，而在连接过程中，片段连接处发生的非模板核苷酸的添加和删除极大地增强了这种组合多样性。这个系统产生了大量多样的受体序列[4,5]，而每一个血液或组织的生物样本通常都有成千上万或数百万个受体。这阻止了使用常规DNA测序对B细胞或t细胞抗原受体的全部序列进行全局分析。

t细胞重组和产生多样性。个体的V和J基因是随机选择的(但不是一致的)，并在胸腺的t细胞发育过程中重新组合。在重组过程中，碱基对可以在最终结扎前被移除或添加到连接处(A)。β基因在V和J之间包含了一个额外的D区域minigene，产生了两个连接(未显示)。最后，转录alpha和beta V区域，剪接到各自的恒定区域(B)上并翻译，两种蛋白异质二聚形成一个TCR (C)。图中所示的TCR/MHC/peptide complex来自PDB结构1FYT，使用RasMol[3]显示。TCR显示为空格填充，肽/MHC复合物显示为右边的stick表示。粉色Vα;黄色Jα;蓝色Vβ;绿色Jβ;和红CDR3上。

t细胞重组和产生多样性。个体的V和J基因是随机选择的(但不是一致的)，并在胸腺的t细胞发育过程中重新组合。在重组过程中，碱基对可以在最终结扎前被移除或添加到连接处(A)。β基因在V和J之间包含了一个额外的D区域minigene，产生了两个连接(未显示)。最后，转录alpha和beta V区域，剪接到各自的恒定区域(B)上并翻译，两种蛋白异质二聚形成一个TCR (C)。图中所示的TCR/MHC/peptide complex来自PDB结构1FYT，使用RasMol[3]显示。TCR显示为空格填充，肽/MHC复合物显示为右边的stick表示。粉色Vα;黄色Jα;蓝色Vβ;绿色Jβ;和红CDR3上。

然而，在过去的几十年中，高通量DNA测序(HTS)的快速发展为日益强大和广泛的BCR和TCR全谱研究开辟了道路。免疫库的异常异质性不仅对实验室测序管道，而且对随后生成的基因组数据集的分析提出了重大挑战。在这篇综述中，我们特别关注TCR系统的计算和生物信息学分析。我们概述了处理免疫组库数据的主要步骤，考虑到原始序列文件的低级处理和试图从数据中提取生物或病理信息的高级算法。我们调查了为这类分析而开发的一些工具，并指出了保留序列分析的潜力，以便深入了解适应性免疫系统的基础生物学，并提供强有力的疾病生物标志物

TCR系统可以被看作是一个高维和密切个性化的生物标志物的最终例子，这是未来精准医疗的标志。可以预见的是，高通量TCR测序分析的第一个十年集中在基因分配和错误纠正等具有挑战性但技术性的问题上。随着这项技术变得更加成熟和强大，人们的注意力正转向开发复杂的计算工具，这些工具是理解这些新颖但往往难以理解的免疫功能指标所必需的。

The processing pipeline: an overview

研究免疫系统的主要阶段如图2所示。它们大致可分为文库准备和测序、低级处理，包括确定受体序列并将其分配给基因组V、D和J基因，以及高级处理和分析。

图2：

免疫系统研究的主要阶段。灰框(上):库准备和排序。绿框(中间):低层次的处理，包括序列装配、基因组V、D和J基因的分配、CDR3区域的提取和错误纠正。参见本综述的基因分配、序列和丰度错误纠正策略和基准及其挑战部分。蓝盒(下):高级处理和分析，包括多样性测量、克隆频率分布测定、V和J使用差异分析、TCR序列个体间共享分析(公共与私有)以及序列与抗原特异性之间的关系。参见章节《高水平全谱处理:揭示生物学和临床意义、本综述的多样性和抗原特异性的测量方法》。

虽然我们的重点是在硅分析，简要介绍目前使用的主要方法为TCR组库的准备和测序是有益的。其他地方对高通量TCR库的制备进行了详细的综述[10-13]。许多测序技术已被应用于研究TCR库，但与基因组学和转录组学一样，Illumina平台已成为事实上的标准([14]，http://www.illumina.com/)。这一地位的实现是由于每个基本测序成本的大幅降低，以及阅读长度和质量的巨大改善。然而，在测序前有几个处理生物样品的管道。所有管道的目标是尽可能完整地记录重新排列的alpha和beta TCR基因序列，并确定它们在样本中的丰度。TCR测序的主要商业服务提供商从基因组DNA或信使RNA (mRNA)扩增([5]，http://www.adaptivebiotech.com/immunoseq)。虽然这些方法的大部分细节都是专有的，但是这些方法是基于使用多重聚合酶链反应(PCR)扩增重组TCR基因，使用一组引物可以捕获V和J基因[5]的所有可能组合。由于V基因和J基因之间存在较大的内含子，非重组基因不会扩增。相比之下，一些研究小组已经开发了基于RNA/互补DNA (cDNA)的技术来进行组库分析。以RNA为起始材料，可以应用互补DNA (cDNA)端(RACE)快速扩增技术[15,16]，从而减少不同V和J基因引物效率不同可能导致的PCR偏倚量。从DNA入手，在样品采集的便捷性和储存的稳定性方面有很多优势。此外，基于dna的策略不受mRNA转录异质性或稳定性的影响。相反，基于rna的技术的一个主要好处是，它们很容易适应引入独特的分子标识符(UMIs)，这可以提供样本中TCR丰度的精确定量估计。正如下面详细讨论的，PCR固有的异质性意味着这对于实现对保留序列的稳健量化至关重要。

一旦TCR基因得到丰富、扩增和测序，原始序列数据文件(通常以FASTQ格式输出)就需要处理，以产生有意义的生物信息。与基因组或RNA-seq协议相同，处理的第一阶段是将序列匹配到已知的基因组引用，在本例中，将每个TCR分配到其种系组成基因片段(图2，低级处理)。然而,这个过程是细胞更加困难,因为识别位点都是由许多类似V, J,和(在某些情况下)D基因,必须尊敬的准确,也因为这些重组区域还包含删除和无模板添加在复合过程中引入(图1)。几个基因赋值方法讨论了下面的基因分配部分。已经投入了相当大的努力来开发算法，以纠正由此产生的核苷酸序列的PCR和测序错误，并纠正有偏差的PCR扩增。这些将在序列和丰度错误纠正策略和基准测试及其挑战部分进行讨论。

最后，一旦一组经过修正的指定序列被组装(一个汇编)，就可以探索挖掘这些数据集的方法(分类和比较不同的汇编，测量多样性，注释TCR抗原特异性，等等)，并提取生物学或病理学意义(图2，高级处理)。高水平处理TCR库的方法将在以下章节中讨论:揭示生物学和临床意义、多样性和抗原特异性的测量。

Gene assignment

在过去的5年中，已经发布了相当数量的低层处理软件程序。通常，该软件是与实验库准备管道一起开发的，主要用于处理这个特定管道产生的数据。我们编制了一份清单，列出了我们所知的所有用于底层处理TCR库序列数据的开源软件工具(表1)。分别讨论了用于分析全局转录组RNA-seq数据中的TCR序列的工具，包括单细胞RNA-seq数据。
这里是表格：

国际免疫遗传学信息系统(IMGT)于1989年首次开发，作为免疫遗传学相关数据的集中存储库。它仍然是免疫遗传学中使用最广泛的参考点，特别是用于维护IMGT/GENE-DB数据库，该数据库包含来自多个物种(包括小鼠、大鼠、兔子和人类)的基因组V、D和J的一组注释良好的基因。这组基因为大多数TCR基因分配工具提供了参考。IMGT的一个重要作用是标准化TCR基因片段的命名法，尽管应该指出，即使在今天的文献中，更古老的命名法也很常见，并且可能在免疫组库分析中造成相当大的混淆。在这里可以找到不同命名法之间的有用比较。www.imgt.org/IMGTrepertoire/LocusGenes/#J.

除了提供TCR基因的参考序列外，IMGT resource page还提供了多种序列分配工具，包括IMGT/V-QUEST[17,32]和一个高通量版本IMGT/HighV-QUEST[18,33]。这两个版本都可以通过Web门户访问。IMGT/V-QUEST允许在一个会话中提交多达15万个序列，但这通常还不足以处理典型的高温超导实验产生的数千万或数亿个读取。基因分配采用全局成对比对算法来确定最合适的V、D、J基因，然后采用Smith Waterman局部比对算法来确定V(D)J基因末端的缺失，以及它们之间的插入。这些算法相对较慢，限制了这些工具在高温超导分析中的使用。

使用与种线序列对齐相同的原则，专门为TCR高温超导系统开发了更快版本的类似全局对齐算法。IgBLAST最初是作为一种分析免疫球蛋白序列的工具开发的，但后来添加了一个用于TCR序列分析[19]的选项，并使用BLAST算法(一种局部比对方法)在IMGT和NCBI数据库中搜索针对种系序列的查询序列。类似地，IMonitor[23]也使用BLAST，另外还有第二个比对步骤，用于在非cdr3序列中找到精确匹配。recovery TCR或RTCR[27]使用Bowtie2[34]作为其默认对齐模块，该模块允许本地对齐或端到端对齐。TCR基因分配的另一层复杂性是t细胞mRNA经常包含不遵循VDJ重组经典规则的序列，但可能包括部分重组事件、保留的基因间序列或连接的J基因[35]。这些非正则序列的分析直接由TRIg[28]来处理，它测试与整个TCR轨迹的一致性，而不仅仅是外显子。

免疫遗传序列分析(IMSEQ)[36]还通过将输入对齐到生殖系来表征TCR序列，此外还使用了一个检查步骤，该步骤查找V和J区域以及CDR3区域的侧边。该算法首先在短核心序列中寻找匹配项，然后对错误得分低于某个阈值的序列扩展对齐。

另一种基因分配策略是使用部分序列或标记来识别特定的V和J区域。解码器[20]实现了一种改进的Aho-corasick[37]搜索算法来查找匹配的字符串(或单次不匹配的标签)，从而实现了更快的基因分配。一旦确定了V和J区域，与实际序列的比对就确定了基因片段区域的末端，从而确定了插入和删除。除了速度之外，标记策略对短标记本身之外的V或J区域中的错误并不敏感。因此，潜在的错误仅限于CDR3区域(见下文)。最新版本的解压缩器，连同用于人和鼠标的V序列和J序列的完整列表，以及标识符标记列表，可以在https://github.com/innate2adaptive/解压缩器中找到。Vidjil算法[21]使用了一种启发式的方法，使用种系基因中唯一的子串来分配V和J基因。它使用子字符串指定一个“窗口”，该窗口以V和J之间的序列为中心，足够大，可以同时包含这两个区域。相同的窗口(如果没有排序错误)被聚集和计数，表示克隆类型及其各自的丰度。进一步对VJ进行细化，Vidjil输出最丰富的20个克隆进行进一步分析。另一个名为LymAnalyzer[25]的分析工具实现了一种不同的策略，使用短序列标记进行基因分配。通过识别与数据库中的引用序列最密切相关的短连续序列集(或标记)，可以对引用序列进行比对(默认情况下是从IMGT数据库中获取的，尽管用户可以指定不同的引用)。该算法允许在无法在第一个标记中找到精确匹配的情况下(从参考V基因的3 '端开始分配，J基因从5 '端开始分配)移动后续标记，从而导致不匹配。然后将序列转换为CDR3序列并进行聚类。相同的序列被组合在一起，并计数来表示克隆型频率，其中组被分类为“核心”或“最小”序列。使用一种称为汉明距离的距离度量方法对“最小”组中的每个序列与“核心”序列进行检查，如果序列满足特定的阈值(默认值为2)，则将它们与最近的“核心”组合并。TCRklass[26]提供了部分序列匹配的另一种变体，它使用k -string匹配从给定的引用集中选择J或V基因，从而为查询序列提供最佳匹配。

MiTCR[38]及其继任者MiXCR[22]将标签的使用与更经典的对齐工具结合起来。首先，确定匹配V和J区域保守的开始和结束模式的子序列或种子。然后扩展子序列的对齐并给出分数。考虑到保守残基的位置，输出是得分最高的对齐。包含超过用户指定的不匹配核苷酸阈值的比对将被丢弃。MiXCR还提供了对特定D区域的赋值。然而，D区域短且彼此相似，使得精确的赋值变得困难，而大多数TCR的 repertoire 研究仅仅将D区域包含在CDR3序列中。

可靠和明确的V和J基因的分配只需要150个碱基对(bp) 5 '的顺序到恒定区域的开始。由于大多数TCR测序流程产生扩增子，这些扩增子专门针对TCR的这一区域，现在的高通量测序仪的读取长度(通常为>PE150 测序)意味着在基因分配之前很少需要序列组装。然而，也有例外，使用V和J区域baits[39]从随机破碎的DNA中富集TCR基因区域，或者从总RNA-seq数据中恢复TCR基因区域[40-42]。后者的一个特别重要的例子是单细胞RNA-seq的应用，它可以匹配α链和β链(从而可能恢复抗原特异性)，并同时表达t细胞功能状态的分析[41,43,44]。上述研究描述了从RNA-seq数据组装TCRs所需的专门生物信息学工具。

除了V (D) J比对,TCR组库分析通常需要对高变的CDR3序列进行鉴定和翻译,这一鉴定和翻译也在决定抗原特异性中起着关键作用 (图2)。从核苷酸翻译成氨基酸序列本身是非常简单而直接的过程, 在很多脚本语言中都有现成的包或者函数,例如在Python中使用Biopython[46]或R中Biostrings包 [47]。翻译后的序列可以确定CDR3区域的保守氨基酸基序。IMGT定义CDR3区域是通过从V基因3 '端保守的第二个半胱氨酸残基到J基因[48]中保守的FGXG基序中的苯丙氨酸残基。虽然第二守恒的半胱氨酸残基通常是V基因上最后一个c端半胱氨酸，但是有一些基因并不是这样的,半胱氨酸是可以在随机重组过程中产生的(尽管很少),因此需要一个更广泛的背景或参考位置来确定CDR3的真确的起始位置。Noncanonical C-terminal motifs (e.g. FXXG or XGXG) also exist （非正则c端基序(如FXXG或XGXG)也存在）。在小鼠和人类中，保守基序分别位于-11和-10位置分别在TRAJ和TRBJ的J基因末端[49]，这使得从不规则的J基因或重组过程中的部分基序缺失中识别非典型序列成为可能。一旦CDR3被确认后,生产重组可能导致一个表达细胞链的定义通常是指那些包含CDR3上序列,是在坐标系对领袖V基因序列的开始和结束的恒定区和不含任何过早停止密码子。上面讨论的所有包都提供了cdr3的翻译，以及生产性和非生产性TCRs的标识。

测序和丰度纠错方法

Sequence and abundance error-correction strategies

核苷酸序列的错误，即最终数据不能真实地反映输入的分子，是在文库准备的不同阶段产生的。这种错误的发生主要是由于酶(通常是逆转录酶和DNA聚合酶)无法结合正确的核苷酸，或者在DNA测序反应过程中调用了错误的碱基。这些错误是不可避免的:即使现代聚合酶所宣称的100万个碱基中有1个错误是可信的[50,51]，一个从1000个输入分子中扩增出短的300 bp扩增子的20个循环PCR预计也会包含超过30万个错误碱基。

这样子的错误在分析高通量测序数据时呈现出了一个严重的问题。大多数TRCs只出现一次或两次在一个样本中而且可能只在单碱基上有真正的不同。因此，错误的序列可能被误认为是真正的序列，并人为地夸大了一个库(repertoire)的多样性。例如，我们对破伤风类毒素特异性t细胞克隆进行了测序，并观察到>150 alpha基因和>200 beta基因变异序列，其中很多在实验中只出现过一次。

标准的FASTQ输出为每个核苷酸提供了一个质量(Phred)评分，该评分估计碱基被错误called的可能性。也许最直接的减少错误的策略是过滤reads或删除低质量分数的base，通常使用Q30作为阈值[53,54]，Q30是指在Phred分数上，一个base call 不正确的估计概率为千分之一[55]。这种技术只会消除在base call过程中产生的错误，因为PCR错误不太可能受到低质量的影响。此外，由于质量评估受文库准备、测序平台和DNA序列上下文的影响[56,57]，无监督质量过滤可能会消除特定TCR重排的偏见。

一个相关的方法是简单地删除所有低频序列，因为它们是最有可能由于错误而产生的。这个阈值可以预先确定，例如每个TCR[58]至少需要5次读取，或者可以根据数据估计阈值。在他们的论文测序人类TCRβ链[53],沃伦等人报道,没有任何过滤,似乎有成上千的novel Jgene在beta-chain重组时检测到;通过只保留前96%的reads (D96 cutoff)，人工序列被删除，这个数字被减少到预期的13个TRBJ基因。然而，Nguyen等人的[54]报告说，在不同的实验中，将来自单克隆TCR转基因小鼠的基因repertoires减少到单个TCRs所需的阈值有所不同。设置这种任意的cutoff，这种方法的一个主要缺点是，它毫无疑问地从数据集中删除了大量真实但罕见的序列+。

更复杂的高通量TCR序列数据错误编辑使用聚类（clustering）来识别相似的序列集，然后将每个集群中较少见的成员吸收到更常见的成员中。早期的一篇repertoire论文使用最近邻算法进行聚类，将其短readTCR数据(54 bp)分解为序列簇，序列簇之间的差异高达2 bp(即汉明距离<=2)[5]。聚类/合并步骤可以通过将较低的丰度TCRs合并到某些阈值[25]内较高的TCRs中，或者只有在V(D)J种系序列中存在不匹配时才这样做。这保留了初始TCR多样性[59]的更大比例，但仍然有可能删除真正的但罕见的TCRs，这促使一些开发人员省略所有这些步骤，以保留许多不常见的序列[26]。一些管道(例如MiTCR/MiXCR)允许用户更改阈值和集群参数，以支持删除错误或保留多样性[22]。在PCR的早期周期中，单靠频率滤波无法检测到逆转录酶或错误，因为这些错误会随着原始序列一起放大。然而，基于数据生成过程的不同统计模型的几种日益复杂的纠错算法已经发表[22,25,27,38]，这些算法声称对TCR序列几乎没有错误分析。

在纠正错误方面的一个主要进展来自于UMIs的合并(图3)。UMIs是在PCR扩增之前(通常在逆转录过程中或紧接着)合并的随机核苷酸的短片段(有时称为分子条形码)。PCR后，UMIs可以用来识别来自同一起始mRNA分子的序列，因为它们都被标记为相同的UMI，因此可以一起计数。这些信息可以用于两个不同但相关的处理目标。第一种方法是从相同的初始模板分子中通过PCR扩增得到的重复序列。这对于在混合体系中获得关于TCR丰度的准确信息至关重要。第二步是识别和纠正PCR或测序protocol本身引入的序列错误。

Figure 3. Error correction using UMIs. (A) Schematic of the error-correction process. Each TCR is associated with a UMI, which acts as a molecular barcode. TCRs are clustered based on UMI. Identical TCRs within a cluster (i.e. with the same molecular barcode) are collapsed to a count of 1. Minority variants within a cluster are similarly merged with the majority variant. The number of clusters (i.e. same TCR, different UMI) gives the corrected abundance count for that TCR. Optionally, barcodes within a specified molecular distance of each other (usually 1 or 2 Hamming units) can be clustered together. (B) The effects of error correction on sequence abundance data for a set of TCR alpha and beta sequences obtained from a sample of unfractionated peripheral blood. The number of TCRs with each abundance observed is plotted against the abundance itself (labeled TCR abundance), e.g. the leftmost point represents the number of TCRs that occur only once in the sample, the next point the number that occurs twice, etc. The figure shows the distribution obtained before (left) and after (right) error correction using UMIs.

（A）UMI的纠错原理图，
（B）