If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Adrem Data Lab, Department of Computer Science, University of Antwerp, Antwerp, BelgiumAntwerp Unit for Data Analysis and Computation in Immunology and Sequencing (AUDACIS), University of Antwerp, Antwerp, Belgium
Adrem Data Lab, Department of Computer Science, University of Antwerp, Antwerp, BelgiumAntwerp Unit for Data Analysis and Computation in Immunology and Sequencing (AUDACIS), University of Antwerp, Antwerp, BelgiumClinical Immunology Unit, Department of Clinical Sciences, Institute of Tropical Medicine, Antwerp, Belgium
Adrem Data Lab, Department of Computer Science, University of Antwerp, Antwerp, BelgiumAntwerp Unit for Data Analysis and Computation in Immunology and Sequencing (AUDACIS), University of Antwerp, Antwerp, Belgium
Antwerp Unit for Data Analysis and Computation in Immunology and Sequencing (AUDACIS), University of Antwerp, Antwerp, BelgiumCentre for Health Economics Research & Modeling Infectious Diseases (CHERMID), Vaccine & Infectious Disease Institute (VAXINFECTIO), University of Antwerp, Antwerp, BelgiumAntwerp Center for Translational Immunology and Virology (ACTIV), Vaccine and Infectious Disease Institute (VAXINFECTIO), University of Antwerp, Antwerp, BelgiumDepartment of Paediatrics, Antwerp University Hospital, Antwerp, Belgium
Adrem Data Lab, Department of Computer Science, University of Antwerp, Antwerp, BelgiumAntwerp Unit for Data Analysis and Computation in Immunology and Sequencing (AUDACIS), University of Antwerp, Antwerp, Belgium
Adrem Data Lab, Department of Computer Science, University of Antwerp, Antwerp, BelgiumAntwerp Unit for Data Analysis and Computation in Immunology and Sequencing (AUDACIS), University of Antwerp, Antwerp, Belgium
T cells exercise a multitude of functions such as cytotoxicity, secretion of immunomodulating cytokines or regulation of tolerance, collectively resulting in an effective control of immune-related disease. Through the unique mechanism of V(D)J recombination, T cells express a highly specific receptor complex known as the T-cell receptor (TCR). Single-cell sequencing technologies have paved the road for interrogating the transcriptome and the paired αβ TCR repertoire of a single T cell in tandem. In contrast, conventional bulk methods are restricted to only one layer of information. This combination of transcriptomic- and repertoire information can provide novel insight into the functional character of T cell immunity. Recently, single-cell technologies have gained in popularity due to improvements in throughput, decrease in cost and the ability for multimodal experiments that integrate different information layers. Consequently, this prompts the need for the development of novel computational tools that integrate transcriptomic profiles and corresponding features of the TCR repertoire. Here we discuss the current progress in the field of single-cell T cell sequencing, with a focus on the multimodality of new approaches that allow the paired profiling of cellular phenotype and clonotype information. In addition, this review provides detailed descriptions of recent computational developments for analyzing single-cell TCR sequencing data in an integrative manner using novel computational approaches. Finally, we present an overview of the available software tools that can be used to perform integrative analysis of gene expression and TCR profiles.
T cells are central players of the adaptive immune system, and play a crucial role in the control of immune-related diseases. In addition, T cells are indispensable mediators of vaccination and immunotherapy response. By generating a highly diverse T cell repertoire, the adaptive immune system is equipped with a powerful toolkit to protect against a broad range of pathogenic organisms and cancer. This diversity in the T cell repertoire is achieved by the generation of a multitude of different T-cell receptor (TCR) complexes through V(D)J recombination. The TCR, expressed on the cellular surface, recognizes small peptides (epitopes), derived from foreign or self-antigens, and presented by Major Histocompatibility Complex (MHC) molecules on antigen-presenting cells. Upon binding of a TCR to its cognate peptide-MHC (pMHC) complex, a T-cell-mediated immune response is triggered [
]. The TCR complex of most T cells consists of an α and β chain, and the diversity of the complex is produced by the recombination process of V and J gene segments in the α chain, and an additional D gene segment in the β chain. This is the first level of diversity, known as combinatorial diversity. During this process of recombination, non-templated nucleotides are added and deleted at the junctions of the segments, drastically increasing the potential diversity of the TCR repertoire [
]. This is known as junctional diversity. Finally, an additional level of diversity is established through the near unconstrained pairing of α and β chains [
]. The total theoretical diversity of the TCR repertoire remains largely disputed but estimations range from 1015 up to 1061, although the true diversity is limited to the total number of T cells in the human body (3•1011) and restricted by selection processes in the thymus [
TCR sequencing has become an invaluable tool for understanding the complex TCR repertoire dynamics. A clonotype is typically defined by the unique combination of consecutive V gene, CDR3 amino acid sequence and J gene. In general, the number of unique clones within a typical bulk repertoire sample may vary between 103 to 106, depending on clonality of the repertoire, sampling conditions and sequencing depth [
]. However, even though TCR sequencing has become indispensable, typical bulk approaches provide only one layer of information on the T cell, as they only capture the TCR characteristics. Considering T cells exhibit a wide array of immune phenotypes, enabling various functionalities ranging from secreting (anti-)inflammatory cytokines to releasing cytotoxic effector molecules to induce cell death, the receptor characteristics do not fully capture the functionality of T cell they originate from. Traditionally, immunologists rely on targeted techniques, such as flow cytometry, to characterize these distinct immune cell phenotypes. Flow cytometry separates cells based on the expression of certain marker proteins, using fluorescently-tagged antibodies to label these targets [
]. Although the number of measurable parameters can be increased to 50 by using altered and improved techniques, such as mass cytometry (also known as CyTOF), this still limits users to a specific number of pre-defined markers [
]. In contrast, RNA sequencing can capture the phenotype of cellular subsets in an unbiased manner as it is not restricted to a limited and targeted selection of markers. However, gene expression at the transcriptional level can be insufficient for discriminating between certain cellular subsets, whereas a protein marker might be more descriptive. For example, the expression of distinct CD45 isoforms, enabling us to discriminate between naive and memory T cells, cannot be identified on the transcriptional level [
]. Additionally, bulk RNA sequencing typically results in a composite mix of gene expression profiles derived from all cells in the sample, which does not sufficiently reflect the cellular diversity. Thus, bulk RNA sequencing often requires prior cell sorting with fluorescent-labeled antibodies to target protein markers, in order to purify cell types within a sample.
A more promising alternative to these conventional techniques is single-cell RNA sequencing, which leverages the power of combining multiple information layers such as paired sequencing of the gene expression and TCR sequences within a single-cell. This feature of multimodality is not restricted to the transcriptional level. For example, the addition of antibodies linked to specific oligonucleotide barcodes (Feature Barcoding) enables the characterization of surface proteins, similar to flow cytometry. Nonetheless, despite its promise, several aspects of single-cell sequencing remain challenging. Compared to conventional (bulk) techniques, single-cell sequencing is still costly and labor-intensive. Consequently, sample sizes are typically lower. However, as the single-cell sequencing field is rapidly growing [
], more recent developments have allowed easy sample multiplexing by including oligonucleotide-labeled antibodies that allow discrimination between samples (Cell hashing), reducing costs and enabling greater sample sizes [
]. The number of cells that can be sequenced with single-cell platforms is several orders of magnitude lower than with bulk methods as well. For example, most single-cell sequencing technologies only allow the assessment of up to 104 cells, whereas bulk approaches can typically confidently assess > 105 cells [
]. However, this number is rising for single-cell sequencing as technology advances.
There are several platforms that enable the characterization of T cells at the single-cell level, each of them differing in how the cells are prepared and how the genetic material is subsequently enriched for sequencing. These factors have substantial impact on the sequencing throughput, depth, costs, and even the ability of generating data from multiple modalities. As prior reviews have extensively compared multiple single-cell sequencing methods, we will not further discuss the respective methods [
]. Next to the previously described advantages, the sequencing of the TCR of a single cell has the benefit of conveniently pairing α and β chains. This is difficult to achieve with conventional bulk methodologies, due to the unknown origin of TCR molecules in bulk experiments. In addition, technical limitations of bulk techniques, and the larger heterogeneity introduced by the additional recombination of the D segment, have resulted in a preferential interest in the TCRβ chain. As a result, most of our understanding of TCR recognition has been based on the sequencing of the β chain alone [
Epstein-barr virus epitope–major histocompatibility complex interaction combined with convergent recombination drives selection of diverse t cell receptor α and β repertoires.
]. Moreover, the power of multimodality provided by single-cell sequencing also allows the inclusion of T cell targeting peptide-MHC dextramers, enabling the identification of antigen-specific T cells, their TCR sequences and their functional phenotypes. This information is vital not only for elucidating the immunopathology of immune-mediated disease, but can be utilized to identify potential immunotherapeutic targets or help guide immunomonitoring in clinical trials [
Collectively, the combination of single-cell RNA and TCR sequencing offers multiple benefits over conventional (bulk) techniques. In this review, we will outline several of these benefits. Furthermore, we aim to introduce experimental or computational immunologists working on bulk TCR sequencing data or bulk RNA-seq data to common high-resolution and unbiased multimodal single-cell workflows that are able to generate new biologically relevant insights. Due to emerging applications in both targeted enrichment of TCR transcripts and single-cell sequencing of TCR repertoires, a large number of different tools are being developed for the downstream analysis of this data, in an attempt to integrate gene expression with clonotype information. The final purpose of this review is to provide an overview of the current state-of-the-art methods and software tools for the integration and downstream analysis of single-cell TCR and single-cell gene expression data. Although this review will primarily focus on the analysis of T cells and their receptors, many of the discussed methods and technologies are also applicable to B cells.
To fully introduce the reader to the benefits of paired single-cell RNA and TCR sequencing, we will first outline the typical workflows for unpaired single-cell RNA sequencing and TCR repertoire analysis. Next, we discuss the advantages of integrating these two layers of information, and the available tools that enable this integration. Finally, we identify current challenges in the field of TCR repertoire analysis and provide novel perspectives on how to bridge the research gaps in this field.
2. General workflow for single-cell transcriptomics
2.1 Power calculation for single-cell transcriptomics
Power calculation is an important component of experimental design. While power calculation approaches for bulk sequencing can be applied on the single-cell level, these often fail to take single-cell specific characteristics such as data sparsity into account. Several factors determine the statistical power of a single-cell sequencing experiment, including the depth of sequencing (# of reads per cell), the number of cells per sample and the number of samples. These factors are typically influenced by budgetary restraints, technical limits of the chosen sequencing platform and sample availability. Although information on the recommended sequencing depths is typically provided by the assay manufacturer (e.g. 10x Genomics), deciding on adequate sample sizes and the required number of cells per sample are more challenging. In addition, other prior knowledge, depending on the research question, might be required to calculate the power. For example, when trying to identify a rare cell population in a sample, prior knowledge on the proportion of cells within that sample type may be required to determine how many cells need to be sequenced for adequate power [
]. As different prior knowledge is required depending on the research question, several single-cell power analysis tools have been developed with varying purposes, such as scPower [
]. While SCEED is focused on power calculation for cell type identification, scPower can be rather wheeled for the power calculation for differential gene expression testing and expression quantitative trait loci analysis. Finally, SCOPIT uses a multinomial distribution to calculate the required cell number, based on the minimum number of cells that must be sequenced per subpopulation, the desired probability of sampling that number of cells for each population, and the frequency of the rarest subpopulation. The model can be used as an intuitive web interface or as an R package.
2.2 Pre-processing of single-cell RNA sequencing data
Sequencing data generated by single-cell platforms, similar to standard (bulk) sequencing applications, requires some processing before proceeding to downstream analysis. The main pre-processing steps are summarized in Table 1. While we will briefly discuss some of these in this section, readers are referred to the supplementary text for a more detailed explanation, including popular software tools that enable pre-processing of (single-cell) RNA-seq data. Some of these steps are also further explained in an excellent review by Luecken and Theis [
]. In addition, a recent benchmark demonstrated that the choice between pre-processing tools is relatively unimportant, observing minor differences after downstream processing [
]. Nonetheless, for V(D)J profiling using the 10x Genomics platform, CellRanger is recommended as it simultaneously processes both the gene expression and paired TCR data.
Table 1Pre-processing steps for single-cell RNA sequencing data.
Pre-processing step
Goal
Method
Acquisition of initial gene expression matrix
Transformation of raw sequencing data to raw gene expression matrix.
Benjamin Kaminow, Dinar Yunusov, and Alexander Dobin. Starsolo: accurate, fast and versatile mapping/quantification of single-cell and single-nucleus rna-seq data. bioRxiv, 2021.
Dongze He, Mohsen Zakeri, Hirak Sarkar, Charlotte Soneson, Avi Srivastava, and Rob Patro. Alevin-fry unlocks rapid, accurate, and memory-frugal quantification of single-cell rna-seq data. bioRxiv, 2021.
Using combinations of 4 metrics per barcode: (1) number of UMIs (2) number of detected genes (features) (3) fraction of reads that map to mitochondrial genes, (4) fraction of reads that map to ribosomal genes
Identification and removal of multiplets
Manual inspection: Removal of multiplets when visualizing the distribution of the number of detected genes for all cells (Caution: potential loss of relevant cells) Computational doublet-detection methods (benchmarked in
Limit the amount of noise by reducing the high-dimensional gene expression matrix to a low- dimensional space, retaining the bulk of the information captured by the data.
Selecting the top n most variable genes (generally 1000 to 5000). Dimensionality reduction techniques: PCA followed by t-SNE
Leland McInnes, John Healy, and James Melville. Umap: Uniform manifold approximation and projection for dimension reduction. arXiv preprint arXiv:1802.03426, 2018.
The first step in pre-processing consists of the conversion of raw sequencing output data to a more interpretable gene expression matrix. The acquisition of the raw gene expression matrix includes read assignment to individual samples and cells based on cellular barcodes (demultiplexing), genome and/or transcriptome alignment, and read quantification. Even after assigning all reads to a cellular barcode, not every barcode will correspond to living cells of sufficient quality. Additionally, a small but non-negligible proportion of the data will consist of distinct multiple cells that have been captured as one single cell (multiplets), potentially confounding downstream analysis. Thus, quality control, involving the removal of low quality cells and multiplets, is typically performed directly after acquisition of the raw gene expression matrix. Next, the gene expression level will be normalized to capture the relative gene expression level between cells, accounting for potential technical variation, such as sequencing depth. Some unwanted variability might still exist after normalization, including both technical and biological factors. Consequently, these factors may need to be corrected for. Generally, technical variation originate from batch effects (cells sequenced in different runs or different sequencing lanes) or the integration of data from multiple experiments, among others. Batch effects can be clearly observed during downstream visualization, as similar cells will separate based on their batch or donor origin. While correcting for technical variation is often warranted, correcting for biological factors, such as cell cycle variability, is not always advisable as it can mask relevant biological information [
]. Lastly, a single-cell RNA-seq gene expression matrix is very high dimensional and prone to noise. The high dimensional matrix will be reduced to a low-dimensional space during dimensionality reduction, thereby reducing noise in the data while retaining essential information.
2.3 Common downstream analysis for single-cell transcriptomics
In this section, the most widely used analysis methods for single-cell gene expression data will be covered. However, available downstream analyses are not restricted to only these popular methods, and include, for example, the identification of gene regulatory networks or the inference of cell-cell communication as well. The techniques discussed in the upcoming section are summarized in Fig. 1.
Fig. 1Frequently used analyses for scRNA-seq data. A. Clustering of cells based on gene expression profiles and subsequent labeling of clusters with suspected cell types. B. Compositional analysis to compare the proportion of particular cell types across conditions. C. Volcano plots can be used to visualize differentially expressed genes identified through statistical testing. D. Ligand-receptor pair analysis or similar analyses based on correlative co-expression of particular genes. E. Functional enrichment analysis using publicly available databases determines up- or downregulated biological processes. F. Cell trajectory analysis to infer pseudotemporal evolution of cells in a dynamic process.
Leland McInnes, John Healy, and James Melville. Umap: Uniform manifold approximation and projection for dimension reduction. arXiv preprint arXiv:1802.03426, 2018.
] create representations where cells that are highly similar in gene expression tend to be grouped together (Fig. 1A). The lowdimensional embeddings created by these algorithms can be used to identify distinct cell clusters or communities based on similarity scores or distance metrics. To achieve this, classical machine learning clustering techniques or graph-based algorithms (community detection) can be applied on either a distance matrix (e.g. k-means clustering) or a graph-based representation, respectively. The Louvain community detection algorithm [
] is the most popular graph-based method, offering great computational performance. Several improvements of the Louvain algorithm have been suggested, contributing towards enhancements in modularity, speed and scalability. These improvements include the smart local move [
] algorithms. More recently, the Leiden algorithm was introduced as an augmentation of the Louvain community detection algorithm, by integrating these earlier improvements [
Once clusters have been established, they can be annotated with relevant biological information such as cell type identity. Generally, this is either performed manually or by using reference-based mappings. Manual cluster annotation requires a priori knowledge of marker genes, genes that are typically only (abundantly) expressed by a particular cell type. Cells are then compared for differential gene expression. Next, annotation labels are provided based on differential expression of the marker genes. In contrast, reference-based mappings make use of reference ’atlases’ to transfer annotation labels to the dataset in question. Commonly employed reference mapping tools include SingleR [
]. Aside from differences between tools in underlying methodologies for mapping, these tools use different reference atlases as well. For example, SingleR makes use of bulk RNAseq datasets of pure cell types to map the query dataset on. In sharp contrast, Azimuth uses atlases consisting of multimodal single-cell datasets that leverage the combination of transcriptome profiling with antibodies targeting known cell markers. Reference-based mapping, however, remains challenging. As these reference atlases consist of datasets that have been generated under different conditions than the query dataset of interest, they may contain an incorrect label for the cells under investigation. However, new datasets that could facilitate novel reference atlases are continuously generated, potentially providing an extensive catalog of high-quality reference atlases that better fit the dataset under investigation. In addition to this dataset generalization problem, droplet-based single-cell methods have a low recovery rate for transcription factors and cytokines [
]. As immune cells are notoriously heterogeneous, subtle differences in the expression of transcription factors or cytokines may pose additional difficulties for reference-based mapping tools to correctly assign T cell subsets [
]. It should be mentioned that this problem also affects manual annotation. However, considering these problems, it is generally recommended to combine reference-based mapping with manual annotation to confirm the validity of the annotations.
2.3.2 Differential cell type composition analysis
Differential cell type composition analysis includes the comparison of the proportions of particular cell types, relative to the total number of cells, between conditions or state (Fig. 1B). Some pathogenic organisms or diseases are known to affect the abundance of certain cell types. For example, CD4+ T-cell depletion is a general hallmark of HIV infection [
]. Cell type composition analysis can thus be a crude way to identify affected cell types in a particular disease, without a priori knowledge. However, the composition of a sample is heavily dependent on the preparation protocol, potentially confounding the compositional analysis [
]. For instance, certain cells may be more prone to stress and damage during library preparation, potentially skewing the proportions as these cells may be depleted in the sample as a result [
As a recent example, researchers applied compositional analysis to study the immune system's response against SARS-CoV-2 infections. In accordance with the observation that COVID-19 patients experience a cytokine storm induced by inflammatory monocytes and pathogenic T cells, compositional analysis was used to illustrate that proliferative T cells and CD14+ monocytes are significantly enriched in patients with severe symptoms [
2.3.3 Differential gene expression and functional enrichment
Although differential gene expression analysis and functional enrichment have long been employed for bulk gene expression profiling, the single-cell environment confers several advantages (Fig. 1C). In contrast to bulk gene expression profiling, resulting in a homogeneous average gene expression profile, single-cell data consists of gene expression profiles for each individual cell, offering greater resolution. Moreover, for each particular cell type cluster, the fraction of cells expressing a certain gene can be calculated. Within the context of differential gene expression testing in single cells, both conventional bulk methods and methods specifically developed for single-cell data are employed. Soneson & Robinson (2018) demonstrated that both techniques perform equally well [
Jordan W Squair, Matthieu Gautier, Claudi Kathe, Mark A Anderson, Nicholas D James, Thomas H Hutson, R´emi Hudelle, Taha Qaiser, Kaya JE Matson, Quentin Barraud, et al. Confronting false discoveries in single-cell differential expression. bioRxiv, 2021.
]. In a recent paper, Bassez et al. introduce the concept of applying differential expression testing on expanded and non-expanded T cells to illustrate that the expanding T cells pre-anti-PD1-treatment are tumor-reactive, showing higher expression of activation, effector and immune-checkpoint markers [
]. In another example, Zhang et al. applied differential gene expression analysis to identify differences in the transcriptional profile of different T cell types in colorectal tumor samples [
Differential gene expression testing typically outputs a large number of differentially expressed genes that require additional biological knowledge to interpret. A common aid herein are pathway analysis methods. Here, annotated genes are grouped into certain sets based on biological features, and computational algorithms test whether any sets are enriched in the differential gene list (over/under-representation analysis) or in the extremes of the ranked log-fold change list (gene set enrichment analysis) [
scRNA-seq provides a static snapshot of the cells at a particular time point. However, some of these cells will be involved in a dynamic process, such as cellular differentiation, the cell cycle or a gradual change in biological function. Thus, clusters annotated with just cell type labels do not fully capture the heterogeneity of the clusters, as they may contain a mixture of cells at different stages along a trajectory of a particular dynamic process. With trajectory analysis, cells are ordered along a path or trajectory based on transcriptional similarity (Fig. 1D). An inferred pseudotime variable represents the progression along this trajectory, starting from a particular cell type that is designated as the root cell. Trajectory analysis thus enables the interpretation of distinct dynamic processes, and the identification of gene expression profiles responsible for branching off along the trajectory. Differential gene expression along the trajectory is also possible [
]. In the context of T cells, integrating RNA velocity with cell trajectories may help to unravel the dynamics of the T cell response and reveal phenotypic transition between clonotypes. The choice of method generally depends on the dataset and trajectory topology, and interested readers are recommended to follow guidelines for method selection as proposed by Saelens et al. [
]. A very unique single-cell transcriptomic profiling study in supercentenarians used trajectory analysis to demonstrate that T cells of these supercentenarians were more terminally differentiated as compared to the T cells of healthy donors [
The TCR repertoire is the collection of clonotypes constituting an individual's T cell landscape. TCR repertoire data can be generated through targeted enrichment strategies or the computational reconstruction of RNA-seq reads. Similarly to gene expression profiling, TCR sequencing data also requires some processing before downstream analysis. In brief, raw sequencing reads are first aligned to a reference set of V, D and J gene sequences, after which identical sequences are grouped into single clonotypes. Subsequently, poor quality reads are removed, and PCR and sequencing errors are corrected, resulting in quantitative clonotype information. There is a wide collection of tools available for processing TCR sequencing reads from bulk experiments [
Imgt/highv-quest: the imgt® web portal for immunoglobulin (ig) or antibody and t cell receptor (tr) analysis from ngs high throughput and deep sequencing.
In recent years, the post-processing of TCR repertoire data to reveal biologically relevant insight has gained more attention. These analyses can be roughly categorized into three major sections: analyzing the repertoire's diversity, specificity and clonal composition. Various methodologies have been developed for analyzing each of these aspects of the TCR repertoire, which have been summarized in Table 2. Additionally, Fig. 2 provides an overview of the different techniques discussed in this section. For each of the methods listed in Table 2 an extensive description can be found in the supplementary materials to this article. Researchers have developed several software tools that cover most of the functionalities discussed in Table 2. These allow the calculation of repertoire statistics like diversity (Fig. 2A), clonal composition or gene usage (Fig. 2C). Some tools provide additional functionalities for comparing different repertoires, for example through the quantification of clonal overlap (Fig. 2C). Lastly, more specific tools exist for advanced analyses of TCR repertoire data such as network analysis (Fig. 2H), clonotype clustering, enrichment analysis or the prediction of epitope specificity (Fig. 2G). For a more comprehensive overview of the available methodologies for analyzing T-cell clonotype data we refer to two excellent reviews by Bradley & Thomas [
Table 2Analysis of TCR repertoire data. The table provides an overview of common methods used to extract knowledge from TCR sequencing profiles. The Description columns provides a brief overview of what the analysis involves. Output describes the potential insight that the analysis permits. AIRR tools gives an overview of open-source software tools that can be used to perform the analysis. The Examples column provides an overview of different studies where the analysis was used to gain novel insight.
Analysis
Description
Output
AIRR tools
Examples
Diversity analysis
Quantification of the diversity within a lymphocyte population. Diversity is typically expressed using Hill numbers. Alternatively, different estimation approaches exist to quantify true repertoire diversity (e.g. Chao2, ICE, DivE and Recon).
Combination of summary metrics that describe different aspects of the underlying frequency distribution.
Identification of a public cdr3 motif and a biased utilization of t-cell receptor v beta and j beta chains in hla-a2/melan-a-specific t-cell clonotypes of melanoma patients.
George Elias, Pieter Meysman, Esther Bartholomeus, Nicolas De Neuter, Nina Keersmaekers, Arvid Suls, Hilde Jansens, Aisha Souquette, Hans De Reu, Evelien Smits, Eva Lion, Paul G. Thomas, Geert Mortier, Pierre Van Damme, Philippe Beutels, Kris Laukens, Viggo Van Tendeloo, and Benson Ogunjimi. Preexisting memory cd4 t cells in na¨ıve individuals confer robust immunity upon hepatitis b vaccination. bioRxiv, 2021.
Bcr cdr3 length distributions differ between blood and spleen and between old and young patients, and tcr distributions can be used to detect myelodysplastic syndrome.
Stefan A. Schattgen, Kate Guion, Jeremy Chase Crawford, Aisha Souquette, Alvaro Martinez Barrio, Michael J.T. Stubbington, Paul G. Thomas, and Philip Bradley. Linking t cell receptor sequence to transcriptional profiles with clonotype neighbor graph analysis (conga). bioRxiv, 2020.
Decompose CDR3 sequences into overlapping k-mers and identify enriched subsequences that may represent important signatures contributing to the specificity of a TCR.
K-mers or sequence motifs enriched in one repertoire (or group of repertoires) versus another (group of) repertoire(s).
Single-cell approach to influenza-specific cd8+ t cell receptor repertoires across different age groups, tissues, and following influenza virus infection.
Perform statistical enrichment tests to identify clonotypes that are significantly enriched in one or a group of repertoires versus another (e.g. healthy versus disease).
A list of TCRs that are statistically enriched in a group of individuals.
Milena Pavlovi´c, Lonneke Scheffer, Keshav Motwani, Chakravarthi Kanduri, Radmila Kompova, Nikolay Vazov, Knut Waagan, Fabian L.M. Bernal, Alexandre Almeida Costa, Brian Corrie, Rahmad Akbar, Ghadi S. Al Hajj, Gabriel Balaban, Todd M. Brusko, Maria Chernigovskaya, Scott Christley, Lindsay G. Cowell, Robert Frank, Ivar Grytten, Sveinung Gundersen, Ingrid Hobæk Haff, Sepp Hochreiter, Eivind Hovig, Ping-Han Hsieh, G¨unter Klambauer, Marieke L. Kuijjer, Christin Lund-Andersen, Antonio Martini, Thomas Minotto, Johan Pensar, Knut Rand, Enrico Riccardi, Philippe A. Robert, Artur Rocha, Andrei Slabodkin, Igor Snapkov, Ludvig M. Sollid, Dmytro Titov, Cédric R. Weber, Michael Widrich, Gur Yaari, Victor Greiff, and Geir Kjetil Sandve. Immuneml: an ecosystem for machine learning analysis of adaptive immune receptor repertoires. bioRxiv, 2021.
Identify the epitope specificity of a TCR by matching against a database with known TCR-epitope interactions or predict the specificity of a TCR using machine learning models.
A list of matched or predicted TCR-epitope interactions.
Milena Pavlovi´c, Lonneke Scheffer, Keshav Motwani, Chakravarthi Kanduri, Radmila Kompova, Nikolay Vazov, Knut Waagan, Fabian L.M. Bernal, Alexandre Almeida Costa, Brian Corrie, Rahmad Akbar, Ghadi S. Al Hajj, Gabriel Balaban, Todd M. Brusko, Maria Chernigovskaya, Scott Christley, Lindsay G. Cowell, Robert Frank, Ivar Grytten, Sveinung Gundersen, Ingrid Hobæk Haff, Sepp Hochreiter, Eivind Hovig, Ping-Han Hsieh, G¨unter Klambauer, Marieke L. Kuijjer, Christin Lund-Andersen, Antonio Martini, Thomas Minotto, Johan Pensar, Knut Rand, Enrico Riccardi, Philippe A. Robert, Artur Rocha, Andrei Slabodkin, Igor Snapkov, Ludvig M. Sollid, Dmytro Titov, Cédric R. Weber, Michael Widrich, Gur Yaari, Victor Greiff, and Geir Kjetil Sandve. Immuneml: an ecosystem for machine learning analysis of adaptive immune receptor repertoires. bioRxiv, 2021.
Anna Weber, Jannis Born, and María Rodríguez Martínez. Titan: T cell receptor specificity prediction with bimodal attention networks. arXiv preprint arXiv:2105.03323, 2021.
Represent the TCR repertoire as a network where nodes represent TCRs and the edges between them indicate similarity (typically hamming distance = 1). The repertoire architecture can be analyzed using various graph theoretic metrics.
Visualization of the repertoire architecture. Quantitative metrics of the repertoire architecture.
Fig. 2Frequently used analyses for TCR sequencing data. A. Repertoire diversity compared across different groups or time points. B. Repertoire overlap describes the pairwise similarity of a range of TCR repertoires. A heatmap is a common type of visualization for repertoire overlap. C. Using chord diagrams, the co-occurrence of V and J genes can be visualized. D. Clonotype tracking is often used to evaluate the abundance of a limited set of TCRs across time points. E. Sequence logos can be constructed from a set of TCR sequences, revealing shared motifs that may contribute to the recognition of common antigens. F. CDR3 spectratyping is used to identify biases in the distribution of CDR3 lengths, which may disclose expansions of oligoclonal populations with similar lengths of the CDR3 region. G. Based on data of known TCR-pMHC interactions, epitope-specific machine learning models can be built for predicting whether an unknown TCR will recognize a specific epitope or not. H. Repertoire architecture can be presented by a network where the nodes represent unique TCR sequences, and the edges connecting them describe the similarity between two TCRs (typically defined as HD = 1).
There is a plethora of software tools that can be used to for exploratory analysis of TCR repertoires. The Immcantation Portal1 hosts a range of different Python and R packages [
], leveraging an ecosystem for end-to-end analysis of TCR-seq data, from mapping raw sequencing reads up to advanced post-analysis (e.g. clustering of clonotypes). Moreover, the Immcantation framework is certified as being compliant with the adaptive immune receptor repertoire (AIRR) Standards2 guidelines for software tools [
], provides an extensive suite of tools for the analysis of TCR data, which includes quantitation of clonotype abundance, repertoire diversity (Fig. 2A), repertoire overlap (Fig. 2B), gene usage estimation (Fig. 2C), clonotype tracking (Fig. 2D), CDR3 spectratyping (Fig. 2F), calculation of k-mer distributions and annotation of clonotypes with database information from VDJdb [
]. Finally, another popular package is VDJtools. This command line tool provides similar functionality to Immcantation and immunarch. VDJtools integrates a TCR neighborhood enrichment test (TCRNET) [
], which can be used to identify enriched clonotypes within a single repertoire as compared to a background distribution.
3.2 Generation probability
One of the most crucial advances in the field of immunoinformatics has been the development of probabilistic models for the V(D)J recombination process [
]. This is known to be a stochastic process that favours the generation of certain TCR sequence conformations over others. These models provide the opportunity to assign a probability of generation (Pgen) to any specific TCR sequence [
]. This Pgen is calculated by explicitly modeling the probabilities of selecting a V, J or D gene (in case of TRB), and the potential nucleotide insertions and deletions at the junctions of these gene segments. This value provides an indication of whether a specific TCR sequence is rare or common. For example, longer TCR sequences tend to have a lower Pgen (i.e. they are rarer) due larger numbers of insertions. In addition, probabilistic models of V(D)J rearrangement permit the generation of large synthetic repertoires that mimic the TCR repertoire of healthy individuals. Based on this concept, Pogorelyy et al. developed ALICE, an approach similar to TCRNET that can be used to identify enriched clones from single repertoire snapshots, using a synthetic repertoire as background distribution [
Understanding which TCRs target which epitopes is arguably the most important challenge in repertoire analysis. This allows the identification of T cells responsible for the neutralization of pathogens. Consequently, this knowledge helps us understand why certain individuals may be susceptible to infection or cancer, while others are able to mount an effective immune response. In the context of autoimmune diseases, identifying TCRs that target self-antigens may leverage potential therapeutic targets. As indicated, immunarch and VDJtools provide functionalities for annotating clonotypes with epitope-specificity using databases of experimentally verified TCR-epitope interactions or associations, such as VDJdb [
] predict the specificity of any TCR towards a finite number of epitopes, based on epitope-specific machine learning models. For this application, TCR sequences are typically transformed into a numerical encoding. Popular types of encoding include the use of physicochemical properties or one-hot-encoding. The recent DeepTCR provides a deep learning framework for generating numerical representations of TCR sequences, which can be used for downstream machine learning applications such as the prediction of TCR-epitope specificity [
Milena Pavlovi´c, Lonneke Scheffer, Keshav Motwani, Chakravarthi Kanduri, Radmila Kompova, Nikolay Vazov, Knut Waagan, Fabian L.M. Bernal, Alexandre Almeida Costa, Brian Corrie, Rahmad Akbar, Ghadi S. Al Hajj, Gabriel Balaban, Todd M. Brusko, Maria Chernigovskaya, Scott Christley, Lindsay G. Cowell, Robert Frank, Ivar Grytten, Sveinung Gundersen, Ingrid Hobæk Haff, Sepp Hochreiter, Eivind Hovig, Ping-Han Hsieh, G¨unter Klambauer, Marieke L. Kuijjer, Christin Lund-Andersen, Antonio Martini, Thomas Minotto, Johan Pensar, Knut Rand, Enrico Riccardi, Philippe A. Robert, Artur Rocha, Andrei Slabodkin, Igor Snapkov, Ludvig M. Sollid, Dmytro Titov, Cédric R. Weber, Michael Widrich, Gur Yaari, Victor Greiff, and Geir Kjetil Sandve. Immuneml: an ecosystem for machine learning analysis of adaptive immune receptor repertoires. bioRxiv, 2021.
] also provides functionalities to train and assess receptor-level machine learning classifiers using various encodings. ImmuneML offers models like K-Nearest Neighbours (KNN), logistic regression, random forests, TCRDist classifier, etc. Additional methods for predicting the epitope specificity of a TCR are presented in Table 2.
4. Generating matched single-cell gene expression and TCR data
4.1 Targeted enrichment of the V(D)J locus
Combined single-cell transcriptomics and adaptive immune profiling data is typically acquired through a targeted enrichment of the V(D)J region in conjunction with gene expression profiling. Amplifying the TCR locus can be performed using three main strategies. The first involves a multiplex PCR amplification using a set of primers that target all V and J gene segments. Alternatively, V(D)J sequences can be purified by targeting them with tagged TCR-specific oligos. These baits will anneal to the target regions, and can therefore easily be captured once the sample has been fragmented. Lastly, the most popular method for V(D)J amplification of cDNA samples is the 5′ RACE strategy. To effectively pair enriched V(D)J sequences and the rest of the gene expression profile, two main methods can be distinguished. Droplet-based methods using microfluidic devices are among the most popular strategies. Commercial examples of droplet-based approaches for isolating and barcoding individual cells are the Chromium device offered by 10x Genomics, ddSEQ by Bio-Rad, Nadia by Dolomite Bio, and inDrop by Illumina. There are also approaches that apply flow cytometric-cell sorting in 96 or 384 well plates to isolate individual cells. However, this approach limits the analysis to only one cell per well, per run. A commercial example of this approach is the C1 Single-Cell Auto Prep system by Fluidigm. A comprehensive overview of library preparation methods and sequencing strategies for paired sequencing are discussed in [
]. Generally, however, these sequencing protocols differ only in the method of amplification.
4.2 Computational reconstruction of TCR sequences
Aside from targeted enrichment, it is also possible to reconstruct TCRs from scRNA-seq data using computational methods. In contrast to targeted approaches, computational reconstruction methods provide lower coverage of the TCR repertoire, but allow re-analysis of existing scRNA-seq datasets, potentially providing additional insights. Moreover, conventional immunoprofiling kits typically contain only α/β amplification primers, resulting in minimal recovery of γδ TCRs. However, the reconstruction of γδ TCRs from the gene expression profile is possible, provided that the data was amplified from the 5′ end. There is a broad range of tools designed for recovering TCR sequences from scRNA-seq data, summarized in Table 3. For a more detailed description of each of the tools listed in Table 3, we refer to the supplementary materials of this article. TCR reconstruction tools generally use a combination of reference-based and de novo assembly, enabling the reconstruction of a considerable subset of V(D)J sequences from transcriptomic data. Although not competitive to targeted amplification methods, recent developments of TCR reconstruction tools have shown significant recovery of TCR sequences from scRNA-seq profiles. For example, the TRUST4 software was able to recover about 70% of all V(D)J sequences from scRNA-seq data [
]. The authors of MiXCR illustrated the recovery of around 3000 TRBs from lymph node metastasis samples, 1700–3000 TRBs from CD4 T cells isolated from spleen, and around 400–1000 TRBs from central nervous system tissue [
]. However, the effectiveness of TCR recovery from scRNA-seq data is highly dependent on the sequencing depth and expression level of the TCR locus, which may vary considerably between cells [
]. In conclusion, reconstruction of TCRs from scRNA-seq samples may be desirable if the aim of the experiment is the identification of expanded or dominant clones in the repertoire Table 3.
Table 3Overview of software tools for reconstructing TCR and BCR clonotypes from single-cell RNA sequencing data. A more detailed description of each method can be found in the supplementary information.
5. When to opt for single-cell over bulk approaches: features of single-cell T cell profiling
Whether the TCR data accompanying gene expression profiles is generated through dedicated enrichment of the V(D)J region or reconstructed from scRNA-seq data, having both layers of information can provide multiple advantages over conventional bulk technologies. Table 4 provides a brief comparison of the main features of bulk and single-cell sequencing approaches.
Table 4Advantages and disadvantages of bulk and single-cell approaches to TCR and gene expression profiling. 1: Repertoire coverage here refers to the total number of unique TCR sequences that can be identified. Depending on the scale of the experiment, single-cell approaches can reach similar repertoire coverage as bulk methods, but this will drastically increase the cost of the experiment. 2: It is possible to study various modalities (e.g. TCR profile, gene expression profile, antigen-specificity, etc.) using bulk approaches, but they cannot be integrated. 3: Generally, bulk approaches are more appropriate for large sample sizes, mainly due to the lower cost, efficiency and duration of the protocol.
5.1 Single-cell sequencing allows integration of functional and immune receptor characteristics
While bulk sequencing of TCRs provides a clear representation of the breadth of the antigenic response, it does not offer information about the functional characteristics of the T cells it originates from. Such information, provided by scRNA-seq, may help elucidate the mechanism of action for T cell subsets associated with pathology. For example, straight-forward analyses include overlaying clonal properties on a UMAP or t-SNE generated from the gene expression matrix of the cells (Fig. 3A). Such analyses are especially interesting when performed in parallel to the cell type annotations described in the subsection on Clustering and cluster annotation. This may reveal certain biases in specific groups of cells, such as hyperexpansion of distinct phenotypic subsets [
]. However, these analyses are not restricted to just overlaying clonality, and may also include the application of previously described TCR-specific analyses (Fig. 2) on distinct T cell subsets. Such analyses may include, but are not limited to, calculating subset-specific TCR diversity (Fig. 3B), annotation with (imputed) antigen-specificity or analyzing bias in the use of germline genes in specific cell subpopulations. For example, Bilate et al. identified a clonally restricted CD4+CD8αα+ population whose differentiation was dependent on local antigen challenge [
]. As another example, Lee and colleagues identified reduced T cell clonality in patients with hereditary chronic pancreatitis (CP), as a result of CD4+ Th cells replacing tissue-resident CD8+ T cells [
Bomi Lee, Hong Namkoong, Yan Yang, Huang Huang, David Heller, Greg L Szot, Mark M Davis, Stephen J Pandol, Melena D Bellin, and Aida Habtezion. Single-cell sequencing unveils distinct immune microenvironment with ccr6-ccl20 crosstalk in human chronic pancreatitis. bioRxiv, 2021.
]. In addition they studied clonotype sharing between different phenotypic subsets, and discovered interactions between CCR6+ Th and Th1 populations, in combination with significant upregulation of CCR6 ligand. These findings may indicate a role of a CCR6-CCL20 signaling pathway in hereditary CP. Tools such as immunarch, Immcantation or VDJtools can be used to evaluate the clonal distribution and diversity of T cells within a specific subset or to compare the overlap between multiple cell subsets. Conversely, the combination of gene expression and TCR repertoire profiles enables direct interrogation of the functional response of specific T cells using conventional transcriptomic analyses. For example, differential gene expression analysis can also be applied on distinct T cell subsets, on expanding versus non-expanding T cells, or on specific T cell clones or clusters (Fig. 3C). This functionality is offered by various software packages, including Scirpy [
]. The practical execution of this analysis is specific for each individual package. We therefore instruct the reader to carefully study the documentation provided with each of the described packages.
Fig. 3Integrative methods for studying TCR and gene expression profiles in tandem. A. Clonotype information, such as clonal expansion, can be mapped onto a gene expression-based UMAP. B. TCR-specific metrics, such as diversity, can be evaluated on the level of different celltypes. C. Evaluating the gene expression profile of clonotypes within a TCR cluster. D. Celltype information can be projected onto TCR ype similarity networks in order to identify clonotype clusters with convergent or divergent celltypes.
Vice versa, information obtained from the gene expression profile may be mapped onto a TCR similarity network (Fig. 3D), which is something that has been explored to a far lesser extent with existing tools. This type of analysis may reveal clusters of highly similar clones (therefore likely targeting identical epitopes) that belong to the same or related cell subsets, revealing expansions of T cell sets both on the phenotypic as well as the clonotype level.
5.2 Power of multimodality: antigen-specificity profiling
Novel modalities have been developed for single-cell sequencing that allow researchers to unambiguously determine the antigen specificity of T cells. In these approaches, scTCR-seq and scRNA-seq are combined with epitopeloaded MHC multimers that epitope-specific T cells will interact with [
]. This introduces a third, and very important layer of information that enables the complete characterization of T cell functionality, providing information about its cellular phenotype, receptor sequence and the peptide-MHC complex that it can recognize. In cancer research, for example, tumor-specific T cells can be identified, and subsequently used for adoptive T cell therapies, by capturing them using barcoded MHC-multimers loaded with the tumor epitope of interest [
]. Moreover, single-cell methodologies allow for the pairing of α and β chain. This provides additional resolution, by including information about both TRA and TRB. In contrast, bulk approaches typically only offer single chain information.
6. Software packages for analyzing T cells at the single-cell level
With the emerging availability of analytical techniques for studying T cells at the single cell level, there has been a need for the development of tools to analyze the growing amount of data that accompanies this technological revolution. There is a plethora of tools available for analyzing either transcriptomics [
Milena Pavlovi´c, Lonneke Scheffer, Keshav Motwani, Chakravarthi Kanduri, Radmila Kompova, Nikolay Vazov, Knut Waagan, Fabian L.M. Bernal, Alexandre Almeida Costa, Brian Corrie, Rahmad Akbar, Ghadi S. Al Hajj, Gabriel Balaban, Todd M. Brusko, Maria Chernigovskaya, Scott Christley, Lindsay G. Cowell, Robert Frank, Ivar Grytten, Sveinung Gundersen, Ingrid Hobæk Haff, Sepp Hochreiter, Eivind Hovig, Ping-Han Hsieh, G¨unter Klambauer, Marieke L. Kuijjer, Christin Lund-Andersen, Antonio Martini, Thomas Minotto, Johan Pensar, Knut Rand, Enrico Riccardi, Philippe A. Robert, Artur Rocha, Andrei Slabodkin, Igor Snapkov, Ludvig M. Sollid, Dmytro Titov, Cédric R. Weber, Michael Widrich, Gur Yaari, Victor Greiff, and Geir Kjetil Sandve. Immuneml: an ecosystem for machine learning analysis of adaptive immune receptor repertoires. bioRxiv, 2021.
] data individually, but less attention has been spent on the combination of both layers of information. Recently, researchers have shown increased interest in the development of such tools aiming at the integrative analysis of TCR and gene expression profiles. In the upcoming section, we will discuss the current landscape of computational tools specific for the analysis of scTCR-seq data. We included all tools, to the best of our knowledge, that had an associated peer-reviewed publication or pre-print article before the 1st of October, 2021. These tools build on the foundations of the rapidly progressing field of TCR repertoire analysis and provide a giant leap towards the system-based analysis of T cell immunity, thereby providing a deeper mechanistic understanding of T cell biology. Table 5 provides an overview of the different functionalities provided by the packages discussed in this review.
Table 5Tools for analyzing single-cell TCR profiles. Asterisks indicate the availability of multiple metrics. A single asterisk (*) corresponds to a single metric (e.g. only Shannon index for measuring diversity), while double asterisks (**) reflect the availability of multiple diversity or clonality metrics. Advanced visualizations may include graph representations, UMAP projections, circos plots... The Clustering column only accounts for receptor-based clustering. Clustering of samples is covered by the Repertoire overlap column. The Integration with GE column additionally indicates the single-cell RNA-seq analysis environment each tool interacts with. GE, gene expression; AIRR, adaptive immune receptor repertoire; B, BCR; T, TCR; Se, Seurat; Sc, Scanpy; N, native.
Stefan A. Schattgen, Kate Guion, Jeremy Chase Crawford, Aisha Souquette, Alvaro Martinez Barrio, Michael J.T. Stubbington, Paul G. Thomas, and Philip Bradley. Linking t cell receptor sequence to transcriptional profiles with clonotype neighbor graph analysis (conga). bioRxiv, 2020.
Integration of gene expression and TCR data typically involves mapping TCR sequence properties to cell subsets, the latter being defined by these cells’ gene expression profile. This approach impedes the identification of new cell subsets by defining them up front. Schattgen and colleagues developed a graph-theoretic approach, Clonotype Neighbor Graph Analysis (CoNGA)3 [
Stefan A. Schattgen, Kate Guion, Jeremy Chase Crawford, Aisha Souquette, Alvaro Martinez Barrio, Michael J.T. Stubbington, Paul G. Thomas, and Philip Bradley. Linking t cell receptor sequence to transcriptional profiles with clonotype neighbor graph analysis (conga). bioRxiv, 2020.
], which aims to identify correlations between gene expression and TCR profiles in an unbiased way. CoNGA builds a similarity graph based on TCR sequence similarity (as defined by the TCRdist measure) and one based on the gene expression data. As a Python package, CoNGA is built on top of the scanpy package [
] and therefore makes use of the AnnData object to store integrated gene expression and TCR sequence data. Additionally, the package integrates an implementation of TCRdist for distance calculations between TCRs. CoNGA provides a graph-vs-graph and a graph-vs-feature analysis. The first analysis involves correlating the gene expression with the TCR sequence similarity graph by identifying clonotypes whose neighbours significantly overlap in both graphs. For each clonotype, CoNGA evaluates all components that are directly connected to this clonotype (graph neighbours) in both the TCR and gene expression graphs. A score is assigned to each clonotype, reflecting the probability that observing this degree of overlap between both graphs is greater than or equal to the overlap that would be expected by chance. To limit the number of false positives, this score is multiplied by the total number of clonotypes. In the second type of analysis, graph-vs-feature, numerical features from either property are mapped onto the similarity graph of the complementary property, thereby aiming to identify graph neighborhoods with a bias in the score distribution.
By applying CoNGA to a collection of publicly available T cell datasets, the authors identified a population of HOBIT+ expressing T cells with long CDR3s enriched for hydrophobic residues. Moreover, they observed strong correlation between the usage of the TRBV30 gene segment and the expression of the conserved EPHB6 gene.
6.2 mvTCR
Supplementing gene expression profiles with functional information from the TCR provides a more detailed understanding of the behavior of different T cell subsets. Typically, this data is processed and analyzed parallel from one another, thereby impeding the identification of novel T cell phenotypes. An et al [
Yang An, Felix Drost, Fabian Theis, Benjamin Schubert, and Mohammad Lotfollahi. Jointly learning t-cell receptor and transcriptomic information to decipher the immune response. bioRxiv, 2021.
] developed a multiview variational autoencoder, termed mvTCR, that creates a joint embedding of gene expression and TCR sequence data at the level of individual cells. By integrating both modalities, it is possible to capture groups of T cells that are correlated on both a phenotypic and functional level. mvTCR applies two types of mixture models to integrate the transcriptomic and TCR embedding into a joint latent distribution. The authors show that this joint embedding improves separation of epitopespecific clusters in the UMAP, as compared to either gene expression or TCR embeddings individually. Consequently, mvTCR-generated embeddings for multimodal single-cell data may be used to improve existing models for predicting TCR-epitope specificity, by integrating an additional layer of phenotypic information. Alternatively, sub-clustering of epitope-specific clusters may reveal epitope-specific expansions in certain T cell subsets.
] is an R-based software package dedicated to the analysis of single-cell immune repertoires. Platypus is optimized for data generated through the 10x Genomics Chromium platform, but it is also compatible with other barcode-based scRNA-seq approaches like RAGE-seq [
]. By default, scaling and normalization of the gene expression data is performed using the default Seurat parameters, although the software also supports alternative normalization methods like SCTransform [
]. Platypus provides a method for extracting V(D)J sequences from Cell Ranger output and it contains a range of functions for pre-processing and calculating basic repertoire statistics. The latter include calculating the number of isotypes per clones (BCRs), CDR3 length distributions and constructing sequence logos. An interesting feature of the Platypus package is its ability to automate Seurat workflows. The results of this gene expression analysis can be subsequently integrated with clonotype information using a custom function. This allows users to project clonotype information onto the UMAP plots generated by clustering the gene expression profiles. For example, the visualize_clones_gene expression can be used to highlight expanded clones within the gene expression clusters. Finally, Platypus provides a feature for evaluating repertoire topology through the construction of sequence similarity networks.
] is a Python library built on top of the Scanpy toolkit for analyzing scRNA-seq data in Python. Data can be directly imported from various sources, including Cell Ranger, TraCeR and standardized AIRR formats. Similar to Scanpy and CoNGA, Scirpy leverages the AnnData format [
], which stores a matrix of values along with annotations of observations and variables. This data structure also keeps track of additional unstructured annotations. Moreover, Scirpy follows the practices of the intuitive API of Scanpy. To integrate V(D)J and gene expression profiles, Scirpy provides functions for merging AIRR and gene expression data into a single AnnData object. Scirpy offers tools for pre-processing and analyzing TCR repertoire and gene expression data in tandem. The pre-processing procedure allows up to two α and β chains per T cell, flagging any cell containing more than two of either chains as potential doublets and discarding them during the process. Analysis tools include calculation of clonotype abundance within a specific group of samples, clonal expansion, diversity, imbalance, as well as repertoire overlap. However, the only available diversity metric is Shannon entropy. Next to this, the package provides graph visualizations of clonotype clusters with high sequence similarity, using igraph [
] is an R package designed for the post-analysis of filtered contigs generated from the Cell Ranger pipeline. The package interacts with Seurat and SingleCellExperiment (SCE), allowing for integration of gene expression data. Various functions are provided for the visualization of T cell contigs, which includes abundance, length, gene usage and clonotype sharing plots. scRepertoire also offers more advanced types of analysis, such as clonal homeostasis (visualization of different levels of expansion) or clonal proportions (the proportions of clone sizes). Other analysis include the calculation of repertoire overlap, sample diversity and clustering of clonotypes based on the amino acid edit-distance (the number of mismatched amino acids between two sequences).
Previously described features can also be calculated for gene expression clusters. Integration with Seurat also allows projection of clonotype information on the UMAP plots. Additional advanced visualizations include alluvial plots displaying clonotypes shared across different categories. Lastly, shared clonotype gene usage patterns across cell type clusters can be analyzed using chord diagrams.
6.6 Tessa
Tessa is a tool that generates numerical embeddings for TCR sequence and integrates it with the gene expression profiles of T cells [
] of the amino acids in the CDR3β region. Tessa uses a stacked auto-encoder to reduce the size of the numeric vector, while maintaining its intrinsic structural features. For the gene expression matrix, only the top 10% genes with the highest variation in expression are kept. Tessa then uses a parametric Bayesian model to identify the influence of the TCR on the gene expression profile of matched clones. In addition, tessa uses weighted TCR embeddings to cluster clones into groups that represent their antigen specificity. The algorithm alternates between both processes (correlating TCR and gene expression matrix, and antigen-specific grouping), updating the weights of the embeddings, until the model reaches convergence. Using tessa, Zhang and colleagues firstly showed that clonotypes sharing similar TCRs are more likely to also share similar gene expression profiles, as determined from the correlation between embeddings from the TCR and transcriptomic profiles. Moreover, the correlation was stronger in PBMCs from healthy donors than tumor samples of different cancer types. This may indicate a proportionally smaller influence of the TCR on the gene expression profile in tumor samples, which may be a consequence of high cyto- and chemokine secretion in the tumor microenvironment, influencing T cells transcriptionally [
] integrates various R packages for analyzing scRNA (Scater, Seurat, SC3, Monocle & MAST) and V(D)J sequencing data (immunarch) into an easy-to-use R Shiny web application. As input, the software allows 3′- as well as 5′-generated scRNA-seq data (both 10x and SmartSeq2). Moreover, TCR sequences can be directly reconstructed from the input scRNA-seq data, using the VDJPuzzle software [
]. The tool offers various features for analyzing clonotype abundance, CDR3 length distributions, V(D)J gene usage and clonotype sharing. For the analysis of gene expression levels, the tool includes common dimensionality reduction techniques such as PCA, t-SNE and UMAP. Additionally, cell clustering (both supervised and unsupervised) is provided based on gene expression values. Finally, the software offers pseudo-time analysis for determining single-cell state trajectories based on the Monocle package [
]. In addition, the cell population may also influence the potential number of identified clonotypes, as some cell types may be rarer than others. Therefore, researchers must carefully evaluate the choice of single-cell approaches versus bulk, depending on the research question to be answered. Deep sampling approaches that allow the capture of large numbers of cells (e.g. leukapheresis) have revealed up to 2•107 unique clonotypes within a single sample [
]. From a practical perspective, the analysis of this amount of cells is only possible using bulk sequencing approaches. For single-cell experiments, the number of uniquely identified clonotypes is typically lower. Consequently, bulk sequencing approaches may be more appropriate when the goal of the study is to characterize the full repertoire from whole blood samples. However, when interested in the functional characteristics and the phenotypes of specific (sub)populations, one may opt for single-cell technologies. These may include situations where the number of profiled clonotypes is less relevant. For example, when studying certain epitope-specific T cells and the immune responses they elicit.
8. Perspective
Single-cell technologies have opened up novel opportunities for identifying specific αβTCR pairs along with the functional profile of the cells they originate from. The information obtained through the use of these technologies simultaneously delivers gene expression profiles, TCR sequence information and optionally other modalities (such as peptide specificity, epigenetic modifications, chromatin accessibility, etc.). While a plethora of techniques have been established to individually analyze these information layers, the use of single-cell technologies has offered a new way to integrate this information at the level of individual cells. This imposes a tremendous challenge from the perspective of data analysis. In this review we have discussed several excellent software modules that offer tools aimed at the integrative analysis of paired single-cell gene expression and TCR repertoire data. Although these packages provide a comprehensive toolkit for exploring and analyzing gene expression and TCR profiles, several issues remain. scTCRseq allows chain pairing, granting information about both α and β chain. While this is considered a major advantage, there remains an unresolved problem with the pairing of α and β chains, even in the case of single-cell sequencing. Occasionally, a single cell may express multiple productive α and/or β chains [
]. In this case it is not possible to know which of the αβ pairs is functional. It has long been known that post-translational silencing mechanisms exist that result in allelic exclusion [
]. Nonetheless, one question AIRR researchers should address is: what determines functional chain pairing in TCRs? Additionally, although rare, there is always the possibility of sequencing only one TRA and one TRB, while the cell may in fact express multiple TRAs and/or TRBs. Moreover, the identified TRA and TRB may not even match as they might only pair the other, unidentified chains. Consequently, this raises the question: whether identified αβ pairs in single-cell experiments are truly functional rearrangements? Another consideration is the fact that gene expression and TCR data obtained during single-cell experiments are typically processed and analyzed individually using established approaches from the fields of scRNA-seq and TCR-seq analysis. Integration is often limited to the projection of clonotype features on gene expression-based UMAP. Therefore, we advocate for the development of novel methodologies that integrate information from both sources into a reciprocal metric. Several approaches already adopt this idea, including CoNGA [
Stefan A. Schattgen, Kate Guion, Jeremy Chase Crawford, Aisha Souquette, Alvaro Martinez Barrio, Michael J.T. Stubbington, Paul G. Thomas, and Philip Bradley. Linking t cell receptor sequence to transcriptional profiles with clonotype neighbor graph analysis (conga). bioRxiv, 2020.
Yang An, Felix Drost, Fabian Theis, Benjamin Schubert, and Mohammad Lotfollahi. Jointly learning t-cell receptor and transcriptomic information to decipher the immune response. bioRxiv, 2021.
]. Integrative approaches like this may reveal distinct subpopulations of T cells that display similar gene expression and TCR sequence characteristics. Such observations may be explained by the expansion of certain T cell subpopulations, triggered by an immunogenic peptide. Similarly, UMAP is commonly applied to gene expression matrices to project distinct cell subpopulations based on a set of highly variable genes. Little attention has been focused on applying UMAP to a combination of gene expression and TCR features. Such an approach may reveal distinct clusters of epitope-specific cells that could not be identified from gene expression or TCR profile features individually. This idea has been outlined by An et al. [
Yang An, Felix Drost, Fabian Theis, Benjamin Schubert, and Mohammad Lotfollahi. Jointly learning t-cell receptor and transcriptomic information to decipher the immune response. bioRxiv, 2021.
] who developed a variational autoencoder, mvTCR, to generate a joint embedding for gene expression and TCR sequence information, thereby improving separation of epitope-specific clusters in the UMAP.
Another major challenge for scTCR-seq is the development of improved visualizations. Currently, the most common approach for visualizing scRNAseq data is a UMAP. UMAPs can be annotated with additional layers of information such as clonal expansion amongst others. For TCR sequences, the similarity network representation is one of the most common visualizations. Although this representation provides a general overview of the repertoire architecture and highlights clonal expansions [
]. It is imperative to extract relevant subsets of the TCR network (e.g. expanded clonotype clusters with low generation probability), thus enabling visualization. In addition, features from gene expression space can be mapped onto clonotype similarity networks. Such network representations enable the identification of clonotype clusters with similar expression profiles, potentially indicating common origin or preferential differentiation of a cellular subtype. Conversely, the observation of clonotype clusters with distinct transcriptomic profiles may indicate phenotypic plasticity between celltypes. These visualization strategies map features from one modality (TCR or gene expression) onto the other, but do not truly integrate both layers. Hence, there is a need for improved visualization techniques that captures both gene expression and TCR features by integrating them.
With the advent of novel experimental and computational approaches for determining the specificity of T cells, scTCR-seq combined with scRNA-seq will be an indispensable tool for fully characterizing the complete molecular profile of a T cell. Several methods exist that accurately predict the binding of any TCR with a known epitope, based on an epitope-specific model [
]. These are often referred to as seen epitopes. A major downside to these models is that they require sufficient data for a single epitope in order to make accurate prediction about which TCRs bind to it. In addition, these models are typically trained using only β chain information, thereby neglecting potential contribution of the less diverse α chain. Predicting the binding of a TCR with an unseen epitope is a problem that is considerably more difficult. Nonetheless, multiple studies have illustrated the possibility of solving this problem using deep neural networks [
Anna Weber, Jannis Born, and María Rodríguez Martínez. Titan: T cell receptor specificity prediction with bimodal attention networks. arXiv preprint arXiv:2105.03323, 2021.
]. A common conclusion is that predictions for epitopes similar to known epitopes are superior to vastly different ones. One of the current limitations is the low amount of known high-quality TCR-epitope pairs. However, due to the introduction of high-throughput methods for TCR-antigen screening, more data will become available which will allow the construction of more accurate models for predicting the specificity of any TCR sequence. Finally, we encourage the use of standardized pipelines for processing and analyzing scTCR-seq data, which will result in more transparency and improved comparability between scTCR-seq studies.
BO, KL and PM hold shares in ImmuneWatch BV, an immunoinformatics company.
Acknowledgments
This work was supported by the Research Foundation Flanders (FWO: 1S40321N to SV 1S71721N to NDV and 1S48819N to SG), the iBOF Modulating Immunity and the Microbiome for Effective CRC Immunotherapy (MIMICRY) Project, and the Flemish Government under the ”Onderzoeksprogramma Artifici¨ele Intelligentie (AI) Vlaanderen” Program.
Epstein-barr virus epitope–major histocompatibility complex interaction combined with convergent recombination drives selection of diverse t cell receptor α and β repertoires.
Leland McInnes, John Healy, and James Melville. Umap: Uniform manifold approximation and projection for dimension reduction. arXiv preprint arXiv:1802.03426, 2018.
Jordan W Squair, Matthieu Gautier, Claudi Kathe, Mark A Anderson, Nicholas D James, Thomas H Hutson, R´emi Hudelle, Taha Qaiser, Kaya JE Matson, Quentin Barraud, et al. Confronting false discoveries in single-cell differential expression. bioRxiv, 2021.
Imgt/highv-quest: the imgt® web portal for immunoglobulin (ig) or antibody and t cell receptor (tr) analysis from ngs high throughput and deep sequencing.
Milena Pavlovi´c, Lonneke Scheffer, Keshav Motwani, Chakravarthi Kanduri, Radmila Kompova, Nikolay Vazov, Knut Waagan, Fabian L.M. Bernal, Alexandre Almeida Costa, Brian Corrie, Rahmad Akbar, Ghadi S. Al Hajj, Gabriel Balaban, Todd M. Brusko, Maria Chernigovskaya, Scott Christley, Lindsay G. Cowell, Robert Frank, Ivar Grytten, Sveinung Gundersen, Ingrid Hobæk Haff, Sepp Hochreiter, Eivind Hovig, Ping-Han Hsieh, G¨unter Klambauer, Marieke L. Kuijjer, Christin Lund-Andersen, Antonio Martini, Thomas Minotto, Johan Pensar, Knut Rand, Enrico Riccardi, Philippe A. Robert, Artur Rocha, Andrei Slabodkin, Igor Snapkov, Ludvig M. Sollid, Dmytro Titov, Cédric R. Weber, Michael Widrich, Gur Yaari, Victor Greiff, and Geir Kjetil Sandve. Immuneml: an ecosystem for machine learning analysis of adaptive immune receptor repertoires. bioRxiv, 2021.
00001-5&id=ga1.jpg){kind=link}
00001-5&id=gr1.jpg){kind=link}
00001-5&id=gr2.jpg){kind=link}
00001-5&id=gr3.jpg){kind=link}