Comprehensive Characterizations of Immune Receptor Repertoire in Tumors and Cancer Immunotherapy Studies

T cells and B cells play central roles in adaptive immunity by recognizing viral, bacterial, or parasitic pathogens and tumor antigens through their diverse receptors, T-cell receptors (TCR) and B-cell receptors (BCR). Upon antigen recognition, B cells further undergo somatic hypermutation (SHM) and class switch recombination (1) to evolve BCRs or antibodies with better affinities to the antigen. TCR and BCR analyses have been widely adopted in studies of infectious disease (2), allergy (3), autoimmune disorders (4), tumor immunity (5), and immuno-oncology (6). Though TCR sequencing (TCR-seq)/BCR sequencing (BCR-seq) allows detailed investigation of the immune repertoire in tumors (7), these technologies are expensive and sometimes infeasible with limited tissue biopsies. In comparison, TCR and BCR sequences are contained in RNA sequencing (RNA-seq) data as expressed transcripts, especially the abundant and clonally expanded receptors, but they are usually ignored in transcriptome studies. Traditional analysis methods could not process the most diverse and critical region of antigen recognition on TCR/BCRs, complementary-determining region 3 (CDR3), which is inherently different from the reference genome sequence. We and others have developed computational methods to de novo assemble immune receptor CDR3 sequences from RNA-seq data, such as TRUST (8–11), V'DJer (12), MiXCR (13), CATT (14), and ImRep (15). Our previous work systematically characterized the TCR (8) and BCR (10) repertoires using earlier versions of the TRUST algorithm and thousands of tumor RNA-seq samples from The Cancer Genome Atlas (TCGA). These studies reveal a correlation between the immune repertoire and other tumor molecular features, such as immunogenetic somatic mutations and antibody-dependent cellular cytotoxicity (ADCC).

We previously published TRUST4 (11) with significantly improved performance over previous versions of TRUST in terms of computational efficiency and the number of complete CDR3s recovered. This not only led to more accurate estimations of receptor diversities, but also allows a more detailed characterization of SHM in the CDR3s. TRUST4 assembles longer receptor sequences, which improves antibody isotype identification and association with downstream immune signals. Together, these improved features of TRUST4 motivated us to reexamine the immune receptor repertoires from TCGA tumors.

Immune checkpoint blockade (ICB) therapy has revolutionized tumor treatment by reactivating T cells and the immune system to eliminate cancer cells (16). Studies also showed that B cells could promote immunotherapy response by forming tertiary lymphoid structures in tumors (17–19). Therefore, we also extracted the immune receptor repertoires in 915 RNA-seq samples from seven ICB studies across four cancer types and investigated the repertoire difference between immunotherapy responders and nonresponders.

To facilitate the characterization of immune repertoires under various oncology and immuno-oncology settings, we also developed an open-access web server VisualizIRR (Visualized Immune Repertoire Report). Existing immune repertoire databases, such as TCRdb (20) and VDJdb (21), focus on the receptor sequences without the ability to associate them with clinical information. In VisualizIRR, we processed immune repertoire data from various resources including TCGA and immunoACCESS, integrating metadata from immunotherapy trials. We curated deidentified clinical information for each sample from their original manuscripts, which allows the users to interactively associate immune repertoires with different clinical metadata.

The 10,970 RNA-seq samples from TCGA were downloaded through Genomic Data Commons, along with the gene expression file, of which 10,096 were from tumor tissue, 724 were from adjacent tumor sites, and the 150 acute myeloid leukemia (LAML) samples were from blood (more details in Supplementary Table S1). We collected seven immunotherapy studies from Riaz (22), Gide (23), McDermott (24), Hugo (25), Mariathasan (26), Van Allen (27), and Kim (28). In the Riaz study, the authors split the patients into naïve (treatment naïve) and progressive (progressed after ipilimumab treatment) groups, so we indicated them as Riaz_ipiNaive and Riaz_ipiProg in this work. The 915 RNA-seq samples from these seven immunotherapy studies were all from tumor tissues (Supplementary Table S1). For the immune repertoire reconstruction, TRUST4 was applied on the alignment BAM files, which were available in TCGA and were generated by STAR (29) on the human hg38 reference genome for the seven immunotherapy cohorts.

We applied TRUST4 v1.0.0 to obtain TCRs and BCRs, along with their abundances, for each sample. We explored the diversity of TCR and BCR in TCGA and immunotherapy cohorts, where diversity was defined as the number of clonotypes per thousand CDR3 reads (CPK) as in our previous study (8). The advantage of CPK is that it normalizes the effect of sequencing depths and immune cell infiltration in different samples. Lower CPK means fewer clonotypes given the same number of total CDR3 abundance and may suggest clonal expansion. In addition to CPK, we also measured Shannon entropy ( | ${\smash e = - \smash ∑_{i = 1}^N \ p_i \ log \ p_i}$ |⁠) and clonality ( | ${\smash 1 - \smash\frac {e} {logN}}$ |⁠) for each TCGA sample, where N is the number of distinct clonotypes. When computing the correlation of CPK and other genes’ expression, partial Spearman rank correlations were performed by controlling for the tumor purity of each sample. P values were adjusted by the Benjamini–Hochberg method. Gene Ontology analysis of genes highly correlated with CPK was performed by the clusterProfiler package v3.10.1 (30), with significance thresholds set at P < 0.01 and q-value < 0.05. TCGA RNA-seq sample immune cell deconvolution results based on MCPcounter and CIBERSORT-ABS were downloaded from TIMER2 (31).

We clustered highly similar IGH CDR3 clusters within each sample and investigated the SHM patterns based on these clusters. To establish the threshold for the cluster similarity, we computed the similarities for two CDR3s from different samples of the same cancer type as the null distribution. Because of the large number of IGH CDR3s found in each sample, we randomly selected at most 1 million pairs from each cancer type for the similarity analysis. Because small IGH CDR3 clusters contained too few sequence variations, we only considered the 1.3 million clusters with more than 30 CDR3s each to compute the substitute matrix. For each CDR3, we further exclude the first three residues, as they are highly conserved. Then, we followed the 5.0 version of the BLOSUM Matrices to produce the amino acid substitution matrix. To study the coevolution pattern in SHMs, we aligned TRUST4’s full-length assemblies to their corresponding IGHV genes and utilized AbRSA (32) to index the variations. We applied CCMpred v0.3.2, a well-established pseudolikelihood maximization (PLM) tool, to discover the direct couplings between pairs of columns in the alignment results. The direct couplings indicated the coevolution patterns in IGHs. We considered the 30 most confident coevolution pairs from each IGHV gene analysis and computed the fraction for a pair across 67 IGHV genes with enough TRUST4 full-length assemblies as the consensus coevolution pair confidence. To illustrate the residues in coevolution, a three-dimensional (3D) model of a heavy chain was modeled by Modeller (version 9.0; ref. 33) using the consensus sequence of the IGHV genes, and colored with the variation ratio at each position. The figure was plotted with VMD (version 1.9; ref. 34).

To explore the association between immune features and immunotherapy outcomes, the fold change of immune features (including diversity, richness of IGH and T-cell receptor β (TRB), expression of CTLA4, PD-1, PD-L1, etc.) were calculated between responders and nonresponders for each cohort. P value was given by the two-sided Wilcoxson signed-rank test. The P values are combined with Fisher method or adjusted by the Benjamini–Hochberg method depending on the application. To increase the statistical power from all the samples, we used mixed-effect logistic regression to test the effect of each immune feature on immunotherapy response and regarded the cohort as random effects.

To curate a comprehensive collection of immune repertoire information for VisualizIRR, we collected data from multiple assays, including RNA-seq, TCR-seq, and BCR-seq. The RNA-seq immune repertoire cohorts were sourced from published immunotherapy studies and TCGA, and the repertoires were reconstructed with TRUST4. We included the results for the receptor chains of the human, being IGH, IGK, IGL, T-cell receptor α (TRA), TRB, T-cell receptor delta, and T-cell receptor gamma. The TCR-seq and BCR-seq samples were obtained from the cancer immunotherapy section of the immuneACCESS database. We included both human and mouse data. Human data include the IGH and TRB loci and mouse data include the TRB locus. Meta information was obtained by manually cleaning the sample tags available in immuneACCESS projects associated with the cohorts. If the cohort had no or incomplete clinical information from immuneACCESS, we extracted meta information from associated studies. The meta sheets we utilized when constructing the VisualizIRR cohort collection were manually constructed amalgamations of these two sources.

We then developed an R script that combs through each cohort and performs statistical analysis on the locus level. The R script largely automated the analysis of each cohort via the inclusion of a cohort-specific configuration file. For RNA-seq–derived data, partial CDR3s were not included in analysis. This script produces sample-level and cohort-level analysis of CDR3 nucleotide sequence length distribution. Out-of-frame nucleotide sequences were also removed and compressed on the CDR3 amino acid sequence level to produce amino acid sequence length distribution analysis. Top amino acid sequence clonotypes were collected along with their associated frequencies within the repertoire. Clonotypes with out-of-frame nucleotide sequences were filtered out from the remainder of the analysis. Top V-gene and J-gene frequencies were collected, as well as combinational top V-gene and J-gene frequencies. Top C-gene frequencies were collected for the IGH locus and top D-gene frequencies were collected for the IGH and TRB loci.

To associate immune repertoire diversity and clinical conditions, we calculated commonly used diversity-associated statistics on the sample locus level, with the majority of calculations utilizing the TCR and BCR repertoire analysis R package immunarch v0.6.5 (35). Calculated statistics include Shannon entropy measure, effective number of types, Gini coefficient, Gini-Simpson index, inverse Simpson index, Chao1 index, clonal proportion, cumulative proportion clonotypes, CPK, and clonality. We also included average CDR3 length and unique CDR3 count. For IGH analysis which includes C-gene information, compartmentalization by isotype is also available. Pooled IGK and IGL analysis is also performed. The analysis for each cohort is exported to its own directory with a data structure that allows VisualizIRR to parse it appropriately.

We utilized existing TCR annotation databases VDJdb and McPAS-TCR to perform CDR3 annotation analysis. We merged the entries of the two databases by CDR3 amino acid sequence and associated species or pathology, with the names of these categories being unified between the databases. We then collected CDR3 overlap between samples and the merged database, with each overlapping CDR3 being weighted by its frequency in the sample repertoire. These overlap values were summed by associated species or pathology in a table. A special column was generated for the sum of all cancer-associated CDR3s. Top noncancer and top cancer columns were kept intact, while the rest were summed in another special miscellaneous column.

Additional processing occurred for TCGA analysis. Cancer type and subtype-level CDR3 distribution and segment usage was calculated, as well as survival analysis based on samples containing above- and below-median statistical values. Overview analysis, containing tumor versus normal sample statistical values from different cancer types and subtypes, is precalculated.

We developed an HTML web portal that utilizes widely used web development libraries, including jQuery, Plotly.js, and dataTables. The portal contains repertoire analysis modules which dynamically display the currently user-selected cohort and clinical conditions. Figures are displayed using Plotly.js so that they can be dynamically adjusted and exported. VisualizIRR is freely available at http://cistrome.org/visualizirr without the need for registration or login. It has been tested in Firefox, Google Chrome, and Apple Safari browsers as well as MacOS, Windows, and Linux operating systems.

TCGA data were downloaded through Genomic Data Commons. The deidentified immune repertoire from TCGA is available at Zenodo: https://doi.org/10.5281/zenodo.6326136. All immunotherapy data are available at http://cistrome.org/visualizirr.

The TRUST4 code is available at https://github.com/liulab-dfci/TRUST4. The code for the analysis is available at https://github.com/liulab-dfci/TRUST4_analysis. VisualizIRR web server is at http://cistrome.org/visualizirr and the code for VisualizIRR is available at https://github.com/liulab-dfci/visualizirr.

We previously published the redesigned TRUST algorithm (version 4) with substantial improvements in efficiency and sensitivity (11). When applied to TCGA, TRUST4 recovered over 35 million complete CDR3 sequences across all TCR and BCR chains from the 10,970 TCGA samples (Supplementary Table S2). Over 90% of high-quality IGH CDR3s found by TRUST3 were also recovered by TRUST4. This ratio was much higher than the 20% overlap ratio between the IGH CDR3s from the paired tumor and adjacent normal samples. In addition, TRUST4 consistently produced about seven times more complete CDR3s than previously published TRUST versions across TCGA samples (Fig. 1A), which could lead to more accurate downstream analysis.

We analyzed the V-gene and J-gene pairing pattern on different chains and noticed that V-J pairing in TRA showed similar distances to the center (Fig. 1B). In other words, the V-J recombination begins on the closest pair of V-J genes, and, if earlier attempts are unproductive, synchronously progresses to V-J genes further away from the center until successful recombination. This joining pattern was observed in mice (36) and used in a previous study to model α/β chain generation (37). Because of the limited number of J genes on other chains (6, 5, 7, 14 for IGH, IGK, IGL, and TRB, respectively, as compared with 60 for TRA) and the involvement of D genes in TRB and IGH, we could not observe the joining pattern on other chains (Supplementary Fig. S1).

B cells improve BCR binding affinity against antigens by evolving their receptor sequences through SHM. TRUST4 reported full-length assemblies, whose sequences contain all the CDRs (CDR1/2/3) and framework regions. For BCR, the alignment similarity between the assembled V gene and the reference V gene could reflect the accumulated SHMs. We observed a higher accumulated SHM rate in the downstream isotypes, in agreement with the chronological splicing of constant genes in isotype switch (Fig. 1C; ref. 38). SHM rates across different isotypes seem to be higher in stomach adenocarcinoma, head and neck squamous cell carcinoma, and colon adenocarcinoma (COAD; Fig. 1C). This may arise from the exposure to diverse food antigens in the gastrointestinal system or the gut microbiome, both of which could stimulate more activated and diverse B cells.

We explored the correlation between patient age and immune repertoire diversity [represented as clonotypes per thousand T-cell or B-cell CDR3 reads or CPK as we proposed previously (8)] in tumors. Immunosenescence, a process of immune dysfunction which occurs with age, was reported in the lymphoid tissues (39) and in peripheral blood mononuclear cells (PBMC; 40, 41). We observed a negative Spearman correlation between BCR IGH diversity and age in 26 of the 33 cancer types in TCGA, of which eight had statistical significance after multiple hypothesis correction (Fig. 1D). There was a similar negative association between TRB diversity and age in TCGA (Supplementary Fig. S2), although less statistically significant due to the lower number of TCRs extracted from the tumors. Together, our observations suggest that immunosenescence also occurs in nonlymphoid tissues.

To identify genes associated with TCR and BCR diversity, we computed the correlation between immune repertoire (TCR β chain TRB and BCR heavy chain IGH) diversity and the expression of each protein-coding gene (Fig. 2A and B). We found that most genes significantly associated with diversity are negatively correlated with diversity. Clonal expansion can decrease the immune repertoire diversity measure CPK. Therefore, we hypothesized that some of the genes negatively correlated with CPKs might be involved in or associated with T-cell and B-cell clonal expansion or activation. Indeed, pathway analysis for the top 50 correlated genes revealed these genes to be significantly enriched in T-cell and B-cell activation pathways (Fig. 2C and D), respectively. In addition, we calculated the correlation between CPK and clonality (the latter defined as the normalized Shannon entropy), and found they were significantly negatively correlated (TRB, IGH spearman correlations ρ = −0.82, −0.54, respectively, P < 1e-10 both; Supplementary Fig. S3a). The correlation between CPK and Shannon entropy is weaker (TRB, IGH Spearman correlation = −0.05, −0.29 respectively; Supplementary Fig. S3b), suggesting lower CPK corresponds to clonal expansion in tissue samples.

Detailed examination of the top 20 genes correlated to diversity revealed their potential function in lymphocyte activation. CCL5 (median Spearman correlation ρ = −0.565 across cancer types) was the top negatively correlated gene with TRB diversity. CCL5 encodes the chemokine guiding T cells to the tumor site and was reported to play an important role in orchestrating tumor-infiltrating lymphocytes (42). The next two most negatively correlated genes GZMH and GZMA (median Spearman correlations ρ = −0.556, −0.548, respectively) are granzymes expressed in cytotoxic effector cells, including CD8⁺ T cells. Moreover, CD8A and CD8B, the signature genes for the coreceptor of CD8⁺ T cells, were ranked 5 and 19 in the list (median ρ = −0.536, −0.457, respectively). On the other hand, the CD4 gene's highest rank was 154 in BRCA samples, and its median rank was 568 across all cancer types (median ρ = −0.274). This suggests that CD8⁺ T cells might have a higher clonal expansion magnitude in the tumor environment than CD4⁺ T cells.

The genes negatively correlated with TCR diversity were consistent across different cancers with some notable exceptions (Supplementary Fig. S4a). First, thymoma (THYM) showed weak or even positive correlations between these genes and diversity, which might be explained by the role of the thymus as the T-cell maturation organ, so TRB diversity is not driven by clonal expansion. In addition, many of the top 20 diversity-associated genes also had weak to no correlations with TRB diversity in pancreatic adenocarcinoma (PAAD) tumors. This might be due to tumor immune evasion through MHC-I complex protein degradation in pancreatic cancer (43), thus making antigen presentation insufficient to stimulate TCR signaling for clonal expansion. EPHB6 and HMGA2 were the two most negatively correlated genes with TRB CPK in pancreatic cancer (PAAD ρ = −0.537, −0.450, respectively; median ρ = −0.022, −0.026, respectively; Supplementary Table S3). EPHB6 and HMGA2 have been shown to be overexpressed in T-cell acute lymphoblastic leukemia (44, 45) where T cells are malignantly expanded. This may suggest abnormally high T-cell clonal expansion and dysfunction in pancreatic cancer.

A similar analysis with BCR diversity identified the genes most correlated to diversity to be related to plasma cells upon clonal expansion of memory or naïve B cells. For example, MZB1 and TNFRSF17 (B-cell maturation antigen; median ρ = −0.451, −0.433, respectively) are known plasma cell marker genes (46). Another top-ranked gene, FCRL5 (median ρ = −0.427), is expressed on plasma cells (47). It is the target of an antibody-based treatment for multiple myeloma where B cells show abnormal clonal expansion (48). In addition to plasma cell marker genes, the top-ranked CD79A is a key gene in the MHC-II signaling pathway, which occurs during the T cell–dependent activation of B cells (49). Although the correlation was observed across many cancer types, it is weak in diffuse large B-cell lymphoma where B cells exhibit malignant expansion (Supplementary Fig. S4b).

Highly similar IGH CDR3s at the nucleotide level from the same sample reflect SHMs in clonally expanded B cells of the same lineage. To find an appropriate threshold to cluster the highly similar IGH CDR3s, we compared the pairwise CDR3 similarities within the same patient tumor (intrasample) and between different patient tumors (intersample). Similarly to the observations of our previous analysis on COAD (11), the intrasample CDR3 similarities follow a bimodal distribution (Fig. 3A). The left peak follows the intersample distribution which provides the null distribution of similarity among unrelated CDR3. The right peak is likely to arise from CDR3s of the same lineage, hence we used a conservative similarity threshold of 0.8 for CDR3 clustering.

Next, we sought to investigate amino acid residue composition in naïve and expanded IGH CDR3s in the tumors. We found that residue usage in the CDR3 to be different from that of proteins in the UniProt database (Fig. 3B; Supplementary Table S4; ref. 50), potentially supporting the unique function of IGH CDR3s in antigen recognition. From each CDR3 cluster, the relatively naïve (because we might miss the original naïve B-cell clone) and clonally expanded CDR3s can be readily distinguished. For example, in each cluster, the most abundant CDR3 can be regarded as the expanded CDR3, and the CDR3 with the earliest isotype and the fewest mutations from the germline can be annotated as the naïve CDR3. We observed that Tyrosine (Tyr, Y) not only was among the most frequently used residues in CDR3s, but also was the favored residue (6% more frequent) in expanded CDR3s compared with naïve CDR3s (Fig. 3B; Supplementary Table S4). Our observation is consistent with a theoretical proposal of Tyr as a favorable residue for antibody-antigen binding (51). This is because Tyr is a hydrogen bond donor and can contribute a comparatively large hydrophobic surface to protein–protein interaction.

The CDR3 clustering also allowed us to evaluate residue substitution distributions in the CDR3s. We compared the BLOSUM matrix (52) derived from our CDR3 clusters with the public BLOSUM100 matrix which contained the general probabilities of residue substitution for evolutionarily divergent proteins. We found that CDR3 BLOSUM and BLOSUM100 to be highly consistent with each other (Fig. 3C), where residues with similar physicochemical properties tended to mutate to each other. For example, a hydrophobic residue has a substantially higher probability to mutate into another hydrophobic residue than to hydrophilic ones. The pattern implied that although SHMs occur randomly, there is a strong negative selection pressure in the tumors to eliminate the CDR3s deviating from the baseline physicochemical properties of the naïve antibodies with initial binding.

The above analysis revealed the patterns of single-residue mutations in CDR3s, and we further investigated SHMs constrained by protein structure. Whereas CDR3 is a single loop on the IGH chain, the V gene constitutes the major part of the IGH protein structure by providing two additional binding loops and seven beta sheets, thus it was the focus of our structural analysis. To this end, we employed the PLM method CCMPred (53) to compute the coevolution of residue substitution in SHM based on the BCR sequences with shared V genes. The most significant coevolutions between the residues in the framework regions are in proximity according to the 3D protein structure (Fig. 3D; Supplementary Table S5). For example, positions 82 and 67 of an IGHV gene are 24 residues apart but the estimated spatial distance is only 6.93 Å, comparable with the average amino acid contour length of 4 Å (54). This observation suggests a structural restraint of SHMs of the global BCR sequences to maintain the original folding and interactions. The coevolution information in BCR sequences can be particularly helpful to advance antibody design and engineering.

Because the antibody constant region is important for downstream immune signaling, we next investigated the association between different Ig isotypes and tumor immune signals in TCGA data. Our previous work (10) showed that IgG1 and IgG3 (denoted as IgG1+IgG3) antibody fractions are positively correlated with CD16a expression level, a surrogate of natural killer (NK)-cell activation. This observation suggested that IgG1+IgG3 antibodies could promote ADCC to control tumor growth. Besides ADCC, complement-dependent cytotoxicity (CDC) is another essential immune mechanism to eliminate malignant cells (55). The classical CDC pathway is initiated by the binding of the C1Q molecule to IgM or IgG antibodies which are bound to surface tumor-associated antigens on the target cells. A cascade of events follows, which results in the formation of a membrane attack complex and target cell lysis (Fig. 4A). We observed a positive correlation between IgM, IgG3, and IgG1 fractions and the expression of C1QA, C1QB, and C1QC across TCGA cancer types (Fig. 4B; Supplementary Fig. S5a), consistent with these three isotypes being the most effective in activating the CDC pathway (55). The correlations between C2–C9’ expression and isotype fractions were weaker and sometimes opposite from that between C1Q and isotypes (Supplementary Fig. S5b). This may be due to the fact that soluble C2–C9 proteins are mostly expressed in the liver (56).

Neonatal Fc receptor (FcRn) binds to the IgG antibodies and prolongs their half-life (Fig. 4C; ref. 57). We hypothesized that higher FcRn abundance would maintain more IgG antibodies in the tumor microenvironment and result in stronger ADCC and CDC signals. Indeed, tumor samples with higher expression of the FcRn gene (FCGRT) had higher expression of CD16a and C1Q regardless of IgG fractions (Fig. 4D). Furthermore, we noticed that FcRn expression had little or even negative correlation with IgG fractions across cancer types (Supplementary Fig. S6), suggesting that FcRn and IgG1+IgG3 antibodies could promote ADCC and CDC independently. In addition to individual gene expression, we applied MCPcounter (58) and CIBERSORT-ABS (46) to computationally deconvolve the tumor bulk RNA-seq samples. The inferred abundances of NK cells and macrophages can be regarded as the surrogates for ADCC and antibody-dependent cellular phagocytosis (ADCP). We observed the same pattern that higher IgG1+IgG3 fraction and FcRn expression related to stronger ADCC and ADCP (Fig. 4D; Supplementary Fig. S7). Together, our findings support the antitumor role of IgG1+IgG3 antibodies and FcRn and their potential association with CDC, ADCC, and ADCP.

Previous studies using TCR-seq analyses have shown that ICB treatment can reshape T-cell receptor repertoires (6, 22). In this study, we explored the tumor-infiltrating immune repertoires extracted from RNA-seq of tumors undergoing anti–PD-1, anti–PD-L1, or anti-CTLA4 ICB treatments from seven clinical studies (22–28). We first examined the two melanoma studies, Riaz (22) and Gide (23), both with pretreatment and on-treatment/posttreatment samples. We further split the Riaz study into ipiNaïve (treatment-naïve) and ipiProg (progressed after ipilimumab treatment) cohorts, as was done in the original study. We have established that the diversity measure CPK is significantly negatively correlated to clonality in TCGA analysis. By comparing the paired samples from the two timepoints, we observed trends of TRB clonal expansion and higher clonality upon ICB treatment regardless of therapy outcome in the three cohorts (Fisher method combined P = 0.095 for nonresponders, 0.029 for responders; Fig. 5A). For BCR, the responders’ diversity tends to decrease during the treatment in the patients from Gide and Riaz_naïve (Fisher method combined P = 0.238 for nonresponders, and 0.045 for responders), suggesting B-cell clonal expansion could be beneficial in melanoma immunotherapy. The BCR diversity remained similar in Riaz_ipiProg, which suggested preexisting BCR clonal expansion upon the initial ipilimumab treatment.

We extended the analysis to compare the immune repertoires between responders and nonresponders across all the cohorts. We compared the fold change of various markers between responders and nonresponders, before and after treatment, respectively (Fig. 5B). Higher expression of immune checkpoint genes, PDCD1 (PD-1), CD274 (PD-L1), and CTLA4, were associated with better immunotherapy outcomes in Gide (23) and Kim (28). Though responders had higher T-cell infiltration based on MCPcounter in Riaz_ipiProg, Gide, and Kim cohorts, we could not associate the TRUST4-derived TCR abundance and diversity with therapy outcome. This might be due to the limited number of samples and assembled TCRs. In BCR analysis, we observed a significant difference between responders and nonresponders in B-cell infiltration, with higher infiltration in responders (Fig. 5B). In addition, TRUST4-estimated IGH abundance, CDR3 richness and diversities were higher in responders in the Riaz_ipiProg, Gide, McDermott (24), and Mariathasan (26) cohorts. Responders had higher B-cell signals in these cohorts, suggesting that the B-cell immune repertoire measures could be useful for ICB treatment response prediction, as in previous studies (17–19). Furthermore, a higher IgG1+IgG3 fraction was associated with a better ICB outcome (Fig. 5B) in line with our earlier observations that IgG1+IgG3 antibodies could promote antitumor immunity. We repeated the analysis for the on-treatment samples, and found that T-cell infiltration, B-cell infiltration, and IGH richness were also higher in responders. To combine the statistical power across different cohorts, we applied mixed-effect logistic regression using the cohort as the random effect. The analysis confirmed that the responders had significantly higher B-cell abundance, IGH CDR3 richness, clonal expansions, and IgG1+IgG3 fractions (FDR < 0.1; Fig. 5B right, coefficient values and statistical test results are in Supplementary Table S6).

To help researchers analyze and explore the data from this study, we organized the TRUST4 results on TCGA and immunotherapy RNA-seq samples into a website VisualizIRR (http://www.cistrome.org/visualizirr). VisualizIRR provides a user-friendly interface for users to visualize immune repertoire features under different clinical conditions on either a computer or a smartphone (Supplementary Fig. S8). As an example, we explored TRB diversity in TCGA and confirmed our previous observations (8) that low-grade glioma had low TRB diversity, whereas colorectal cancers (rectum adenocarcinoma, colon adenocarcinoma), non–small cell lung carcinomas (lung squamous cell carcinoma, lung adenocarcinoma), mesothelioma, and melanoma (skin cutaneous melanoma) had relatively higher diversity (Fig. 6A). With TRUST4’s highly sensitive assembly results, VisualizIRR identified THYM samples as having the highest TRB diversity, which was consistent with the role of the thymus in T-cell maturation. Furthermore, VisualizIRR shows diversity comparisons between tumor samples and adjacent normal samples for each cancer type. We found that TRBs had significantly lower diversity in COAD and LUSC tumor sites than the adjacent normal sites, reflecting the lower T-cell infiltration in these tumor samples based on MCPcounter (Supplementary Fig. S9). In IGH isotype analysis, VisualizIRR showed that the IgG1 fractions are positively correlated with longer patient survival in SKCM, ovarian serous cystadenocarcinoma, and LAML (Fig. 6B). This is consistent with the beneficial role of IgG1 antibodies in tumor control in these cancer types (41, 59).

In addition to RNA-seq immune repertoires, we collected TCR-seq and BCR-seq data from an additional 28 immunotherapy cohorts and included them in the VisualizIRR website. The website allows users to examine the immune repertoire diversity under various clinical conditions. For example, TCR-seq of PBMC from healthy donors and patients with newly diagnosed classical Hodgkin lymphoma (cHL) had greater TCR diversity than from relapsed/refractory patients (6)(Supplementary Fig. S10a). VisualizIRR also has a time series feature to reproduce the result that cHL patients with complete response have increased CD4⁺ T-cell diversity upon PD-1 blockade treatment (Supplementary Fig. S10b). In contrast, CD8⁺ T cells were clonally expanded in cHL patients with progressive disease (Fig. 6C). Because cHL tumors are known to have MHC-I deficiency, this result suggests that CD4⁺ instead of CD8⁺ T cells might be the effector cells in anti–PD-1 response in cHL.

To identify pathogens stimulating T cells in patients, we further annotated the TCRs in VisualizIRR with the antigen information from the VDJdb and McPAS-TCR (60) databases. In each sample, VisualizIRR computes the fraction for different types of pathogens, such as cytomegalovirus, or tumor antigens, such as lung cancer antigen. The TCR antigen annotation provides novel insights to explore the relationship between immune repertoire and other molecular or clinical information. For example, in a PBMC TCR-seq study of patients with metastatic prostate cancer undergoing anti-CTLA4 treatment (61), VisualizIRR found that responders and nonrapid progressors had significantly more TCRs targeting tumor antigens than rapid progressors (Fig. 6D). This result suggests that the ability of ICB in expanding tumor antigen–specific TCRs might be the cause of the better patient response.

In this work, we applied TRUST4, an improved de novo TCR and BCR assembler, on TCGA and immunotherapy cohort RNA-seq data to study the immune receptor repertoire in tumors. The exploration of gene expression and immune repertoire information from RNA-seq data helped us identify the relationships of T cells and B cells with other immune pathways or processes. Our study further supports the positive correlation of IgG1+IgG3 antibodies with antitumor signals and immunotherapy outcomes and could inspire future research on establishing B cells as a biomarker for cancer patient survival or treatment predictions.

From TCGA data, we identified genes correlated with T-cell and B-cell clonal expansions such as CCL5, MZB1, and FCRL5 through correlation analysis between repertoire diversity and gene expression. Previously, there have been conflicting reports regarding the function of FCRL5 in promoting or inhibiting B-cell activation and clonal expansion (47, 62). Our analysis suggests that FCRL5 may play a more activating role in the tumor environment. We found that CD8⁺ T cells had more clonal expansion signals based on gene expression patterns, and this could be validated through CD4⁺ T-cell and CD8⁺ T-cell repertoire diversity analysis. CD4⁺ and CD8⁺ TCR information can be obtained after cell sorting, which is beyond the scope of the typical RNA-seq data in this study. We note that clonal expansion was measured by the CPK of TRB and IGH CDR3s, which might be affected by immune cell infiltration. For example, the lower CPK in elder patients could be due to clonal expansion from past infections or immune cell apoptosis in earlier life. Systematic evaluation of diversity measures or development of new metrics that are robust to the impact of T-cell or B-cell infiltrations are still needed.

In the SHM analysis, we found that Tyr was favored in BCR CDR3 SHMs, and the retained SHMs maintained the initial CDR3’s physicochemical properties. These discoveries could help researchers understand the binding affinity improvement during SHM and improve tumor-targeting antibody design (63). BCR is formed by heavy and light chains, where the light chain also undergoes SHM. We observed that SHMs were constrained by the protein structure within a chain but could not analyze the impacts from the paired chain. With the development of single-cell sequencing technology, future studies can include profiling of coevolution between heavy and light chains.

The function of B cells in cancer immunotherapies has received increasing attention in clinical studies (17–19). The observed elevated clonal expansion of B cells and higher IgG1+IgG3 fraction in immunotherapy responders could provide additional evidence that B cells may participate in antitumor immune responses to ICB treatment. Furthermore, it is likely that the B cells in the tumor are activated by cancer-associated antigens, as evidenced by isotype class switch. Characterization of these antigens may lead to novel therapeutic targets, and thus is of high clinical interest. However, due to the limited epitope annotations of BCRs, we could not computationally infer the antigens based on sequences of the clonally expanded BCRs. With the growing size of annotation databases, such as IEDB (64), or the advances of computational approaches, such as DeepCAT (65) for TCRs, we may identify the dynamics of tumor-targeted BCRs during ICB treatment in the future.

We developed VisualizIRR to consolidate various immune repertoire analyses into an interactive open-access website. We have demonstrated how an end user could reproduce previously published immune repertoire analyses and make new observations. The strength of VisualizIRR lies in the curated profiling of numerous immunotherapy and TCGA samples, combined with a variety of commonly useful repertoire analysis techniques. Going forward, we aim to not only build upon the functionality of VisualizIRR but also expand the number of collected datasets. Some potential functionality enhancements include single-cell specific modules and cross-cohort analysis. We envision continued maintenance and development of VisualizIRR to benefit the cancer immunology and immunotherapy research communities.

L. Song reports grants from NCI during the conduct of the study. D. Cohen reports grants from NCI during the conduct of the study; personal fees from ProFound Therapeutics outside the submitted work. X. Hu reports other support from GV20 Therapeutics outside the submitted work. K.J. Livak reports grants from NIH/NCI during the conduct of the study; other support from Fluidigm Corporation outside the submitted work. X.S. Liu reports grants from NIH during the conduct of the study; and X.S. Liu conducted the work while being a faculty at DFCI, and is currently a board member and shareholder of GV20 Oncotherapy, and the full-time CEO of its subsidiary GV20 Therapeutics. X. Hu conducted the work while being a postdoctoral fellow at DFCI, and is currently a full-time employee of GV20 Therapeutics. No disclosures were reported by the other authors.

L. Song: Conceptualization, resources, data curation, software, formal analysis, supervision, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. Z. Ouyang: Resources, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. D. Cohen: Resources, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. Y. Cao: Data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. J. Altreuter: Writing–original draft, writing–review and editing. G. Bai: Writing–original draft, writing–review and editing. X. Hu: Conceptualization, supervision, writing–original draft. K.J. Livak: Conceptualization, supervision, funding acquisition, writing–original draft. H. Li: Funding acquisition, writing–original draft. M. Tang: Writing–original draft, writing–review and editing. B. Li: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. X.S. Liu: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.

We acknowledge the following funding sources for supporting this work: NCI grants U01CA226196 (H. Li), U24CA224316 (CIDC at DFCI), 1R01CA245318 (B. Li), 1R01CA258524-01A1 (B. Li), and China Scholarship Council (Z. Ouyang and Y. Cao).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.