TIRTL-seq: deep, quantitative and affordable paired TCR repertoire sequencing

Abstract

The specificity of T cells is determined by T cell receptor (TCR) α and β chain sequences. While bulk TCR sequencing enables cost-effective repertoire profiling without chain pairing information, single-cell approaches provide paired data but are costly and limited in throughput. Here we present throughput-intensive rapid TCR library sequencing (TIRTL-seq), an experimental and computational methodology for paired TCR repertoire sequencing (TCR-seq). TIRTL-seq is based on the parallel generation of hundreds of TCR libraries in 384-well plates at less than US$200 per plate, allowing cohort-scale paired TCR-seq studies. We benchmarked TIRTL-seq against state-of-the-art bulk TCR-seq and 10x Genomics Chromium technologies on longitudinal samples and identified severe acute respiratory syndrome coronavirus 2- and Epstein-Barr virus-specific clonal expansions after infection with distinct dynamics. TIRTL-seq offers a universal protocol scalable from a single cell to millions of T cells per sample, simultaneously delivering both precise clonal frequency estimation and accurate TCR chain pairing, combining the strengths of bulk and single-cell TCR-seq.

Fig. 1. Development and optimization of TIRTL-seq protocol.

a, Schematic of TIRTL-seq protocol. In brief, a cell suspension is distributed into 384-well plates containing an RT–lysis master mix under a hydrophobic overlay using noncontact liquid dispensers. After the RT reaction, PCR I master mix with V-segment and C-segment primers is dispensed into the same plate. The PCR I product is then diluted and transferred to the PCR II plate for indexing PCR with well-specific unique dual indices. The PCR II products are pooled by centrifugation, purified, size-selected using magnetic beads and sequenced on an Illumina platform. Total library preparation cost is listed for one 384-well plate. b,c, The sensitivity of single sorted T cells according to CellsDirect (b) and TIRTL-seq (Maxima H-based) (c). Green, both TCRα and TCRβ identified; orange, TCRα lost; yellow, TCRβ lost; blank, no cell present. Column 24 is a negative control (no cells sorted). d, Relative fraction of cells with both TCRα and TCRβ identified (green), lost TCRα chain (orange) and lost TCRβ chain (yellow) shown for Invitrogen CellsDirect RT protocol (left) and TIRTL-seq protocol (right). e, TIRTL-seq shows robustness to an increasing number of PBMCs. Number of unique clonotypes detected from each well (y axis) plotted for different numbers of cells per well (x axis) for the TIRTL-seq protocol (n = 32 wells per group), each box representing the interquartile range (IQR), the line inside the box indicating the median, whiskers showing ±1.5× IQR and overlaying points representing individual values. Panel a created with BioRender.com.

Fig. 2. Computational TCR chain pairing with TIRTL-seq.

a, TIRTL-seq scaling with different analysis methods. The total number of unique paired αβTCRs recovered (y axis) is plotted using three distinct analytical approaches: increasing the number of individual plates analyzed separately and concatenating the results (red), increasing the number of cells per well analyzed (blue) or increasing the number of total wells analyzed jointly (green). The dotted orange line shows the number of unique αβTCRs from a single 10x Genomics scTCR-seq reaction. b, Clone-size distribution of paired and unpaired (dark gray) clones by TIRTL-seq. Paired αβTCRs from a 10x Genomics scTCR-seq experiment are ordered by the number of cells with a clonotype (y axis, log scale), while the x axis shows the clone rank. Each color indicates overlap between 10x Genomics scTCR-seq and αβTCRs called by MAD-HYPE (orange), T-SHELL (green) or both (blue). Dark gray indicates lack of pairing in TIRTL-seq compared with 10X Genomics scTCR-seq. c, T-SHELL algorithm pairs TCR chain through frequency correlation. Top: correlation of relative per-well frequencies of largest TCRα no. 1 and top 5 largest TCRβs in the repertoire; red line shows linear fit. Bottom: Manhattan plot for unadjusted P values (two-sided Pearson correlation t-test) for pairing of TCRα no. 1 to 10,000 most abundant TCRβs. Dotted line shows P value cutoff after Bonferroni multiple testing adjustment. d, T cell expansion increases pairing efficiency. T cell clone frequency (y-axis) is plotted against rank (x-axis). The dotted line shows the minimal three-well occurrence threshold for αβ chain pairing by TIRTL-seq. An increased number of clones clearing the threshold after expansion (orange curve) results in more called pairs compared with before expansion (green curve). The inset shows clonal frequency distortion after antigen-independent T cell expansion. Clonal frequency before expansion (x axis) is plotted against clonal frequency after expansion (y axis). e, The fraction of TCRβs overlapping between 10X Genomics scTCR-seq (filtered or unfiltered for clones with >1 cell) and TIRTL-seq experiments (x-axis) with matching or mismatching TCRα for MAIT and non-MAIT clones. f, The fraction of clonotypes with a given chain (α or β) paired with one or more partner chains (TIRTL-seq data from one 384-well-plate experiment).

Fig. 3. Longitudinal clonal tracking with TIRTL-seq.

a, Longitudinal sampling of a donor with SARS-CoV-2 infection. scGEX, single-cell gene expression. b, Top: TIRTL-seq identifies expansions and contractions in the CD8⁺ T cell repertoire. Bottom: CD4⁺ T cell clonal frequencies from the same time points. The dashed diagonal indicates the line of equality. c, Colored dots show clonotypes matching known TCRs specific for A*02 YLQ (cyan) and B*07 SPR SARS-CoV-2 epitopes (green) and A*02 GLC (orange) and B*07 RAK EBV epitopes (red) on pairwise time point comparisons. The dashed diagonal indicates the line of equality. d, Cumulative frequency of CD8⁺ clones specific to A*02 YLQ (cyan) and B*07 SPR SARS-CoV-2 epitopes (green) and A*02 GLC (orange) and B*07 RAK EBV epitopes (red) across time points. conv, convalescent. Panel a created with BioRender.com.

Fig. 4. Benchmarking of TIRTL-seq against state-of-the-art single-cell and bulk TCR-seq.

a, Number of paired TCR clonotypes identified in each sample by 10x Genomics on CD3⁺ cells and TIRTL-seq (split into CD8⁺ and CD4⁺ subsets). Color shows clonotypes identified only in one time point (orange), two time points (yellow) and three time points (green) by each method. b, Number of expanded and contracted clonotypes between time points identified with 10x Genomics scTCR-seq, bulk 5′RACE TCRβ-seq and TIRTL-seq between baseline/acute and acute/convalescent time point pairs. c, The number of expanded (from baseline to acute) clonotypes from an independent bulk TCRβ sequencing experiment paired using 10x Genomics scTCR-seq (blue) and TIRTL-seq (green). d, The overlap between expanded clonotypes between baseline and acute time point identified with bulk TCRβ sequencing (yellow) experiments or TIRTL-seq (green, CD4⁺ and CD8⁺ combined).

Extended Data Fig. 1. Benchmarking of TIRTL-seq with Jurkat spike-in.

a-c. Plate maps showing the number of unique TCRbeta clonotypes detected in each well. Plates a and b were processed with automated 384-well protocol, plate c was processed with manual 96-well protocol without non-contact liquid handling. Row I in the 384-well plates (a,b) and well H12 in the 96-well plate (c) were intentionally left empty (except well I13, which was loaded with Jurkat cells only). d. Distributions of adjusted p-values for clones that were independently paired using 10x Genomics scTCR-seq. Vertical dashed lines indicate the p-value obtained for the Jurkat spike-in pair in each corresponding plate. e. Comparison of the number of pairs identified in the 96-well spike-in experiments by T-SHELL, MADHYPE, or both algorithms, using two different T-SHELL adjusted p-value thresholds (10⁻³ and 10⁻¹⁰).

Extended Data Fig. 2. Functional validation of paired virus-specific TCRs identified by TIRTL-seq.

a. Gating strategy and representative flow plots for in vitro validation experiments. b. In vitro validation of predicted TCR specificity. Bars show average percentage of GFP+ cells out of TCR-transgenic Jurkat cells (mCherry + ), dots show replicates (n = 2).

Extended Data Fig. 3. Anti-EBV VCA (left) and Anti-SARS-CoV-2 Spike IgG levels across time points.

Lowest dashed line shows seronegativity cut-offs.

Extended Data Fig. 4. Clone size estimation accuracy with TIRTL-seq.

a. Correlation (Spearman’s ρ = 0.71) between clonal frequencies estimated by 10x Genomics scTCR-seq and TIRTL-seq for the most abundant clones on acute timepoint (defined as those with ≥5 cells detected in the 10x data). TIRTL-seq CD8⁺ frequencies were multiplied by 0.32 (the frequency of CD8⁺ cells among CD3⁺ cells for this donor) to estimate frequency relative to total CD3+ cells, comparable to the 10x scTCR-seq measurement. (b) Comparison of log₂ fold change (log₂FC) between acute and convalescent time points for large clones (defined as those with ≥ 5 cells detected in 10x at either of the compared time points). Log₂FC calculated from 10x scTCR-seq (x-axis) is plotted against log₂FC calculated from an independent TIRTL-seq experiment on isolated CD8⁺ cells (y-axis). The dashed diagonal indicates the line of equality. Pearson’s correlation coefficient (R) and p-value are shown. Colors indicate clone status determined by the longitudinal TIRTL-seq analysis shown in Fig. 3b: orange for expanding clones, green for contracting clones, grey for stable clones.

Extended Data Fig. 5. Longitudinal clonal tracking with 10x Genomics Chromium scTCR-seq.

10x Genomics scTCR-seq identifies expansions and contractions in the CD3 + T cell repertoire. Each dot represents TCRβ clonotype, frequency at two timepoints is plotted in log-scale. Orange and green color show significantly expanding and contracting clones, respectively.