Landscape of B cell immunity and related immune evasion in human cancers

Abstract

Tumor-infiltrating B cells are an important component in the microenvironment but have unclear anti-tumor effects. We enhanced our previous computational algorithm TRUST to extract the B cell immunoglobulin hypervariable regions from bulk tumor RNA-sequencing data. TRUST assembled more than 30 million complementarity-determining region 3 sequences of the B cell heavy chain (IgH) from The Cancer Genome Atlas. Widespread B cell clonal expansions and immunoglobulin subclass switch events were observed in diverse human cancers. Prevalent somatic copy number alterations in the MICA and MICB genes related to antibody-dependent cell-mediated cytotoxicity were identified in tumors with elevated B cell activity. The IgG3-1 subclass switch interacts with B cell-receptor affinity maturation and defects in the antibody-dependent cell-mediated cytotoxicity pathway. Comprehensive pancancer analyses of tumor-infiltrating B cell-receptor repertoires identified novel tumor immune evasion mechanisms through genetic alterations. The IgH sequences identified here are potentially useful resources for future development of immunotherapies.

Figure 1 ∣. TRUST performance on tumor samples with matched BCR-seq data.

a, Evaluation of the TRUST reported CDR3s under different cutoffs on the minimum clonal frequency. Precision is the fraction of TRUST called CDR3s validated by BCR-seq, and sensitivity is the fraction of BCR-seq CDR3s called by TRUST. b, Evaluation of the TRUST reported immunoglobulin (Ig) isotypes with the same CDR3 in BCR-seq. Precision is the fraction of TRUST called isotypes validated by BCR-seq, and sensitivity is the fraction of BCR-seq isotypes called by TRUST. c, An example of B cell cluster where three sequences (#1, 3, 7) were identified by TRUST with the same Ig isotype as in BCR-seq and TRUST recovered a partial CDR3 for sequence #6. Bases in BCR-seq clones but missed by TRUST are highlighted in violet.

Figure 2 ∣. Immunoglobulin (Ig) isotypes and somatic hypermutations of tumor-infiltrating B cell Ig heavy chain repertoire.

a, Distribution of 9 Ig isotypes across cancer types. Fractions of each Ig isotype were calculated as described in Methods. Top 20 cancer types with most IgH CDR3 assemblies were displayed. Circle size is proportional to the total number of CDR3 assemblies in each cancer type. b, Distribution of somatic hypermutations in the CDR3 region were grouped by CDR3 length as indicated by the number of amino acids (aa). The number of mutation events is shown by the size of horizontal bars, where CDR3 with 14 aa had the largest number of mutations (~0.5 million). c, Proportion of nucleotide pairs (n = 3,808,402) before and after somatic hypermutations. The most likely mutations were within the pyrimidine (C/T) and purine (A/G) groups. d, Distributions of somatic hypermutation rates within each Ig isotypes. The somatic hypermutation rates were calculated from CDR3 pairs with the same Ig isotype in each patient sample. Statistic significances were evaluated using Wilcoxon rank-sum test (two-sided P < 1 × 10⁻⁴ labeled with **** and < 0.01 labeled with **). The boxes in violin plots show the lower, median, and upper quartile of values, and the curves show the density of values (see Supplementary Table 3 for sample sizes and P-values). e, The spectrum of base changes identified in the CDR3 region with the mutation type displayed in the title. The x-axis is the adjacent 3’ and 5’ bases of the mutation, and the y-axis shows the contribution of the mutation context. The (W)RCY motif of AID (R = A/G, Y = T/C) is highlighted as dark red.

Figure 3 ∣. Immunoglobulin (Ig) isotypes subclass switches are associated with B cell and T cell activations.

a, Correlations of Ig switch levels and highly expressed genes in immune cells and malignant cells with corrections on the tumor purity (Methods). The partial Spearman’s correlations and two-sided P-values were from 6,430 samples (see Supplementary Table 3 for sample sizes). Ig switch levels are the number of B cell clusters with any type of switches divided by the total number of unique CDR3s in a patient sample. b, Visualization of Ig isotype co-occurrence within B cell clusters by cancer types. Circle size represents the number of clusters carrying a given Ig isotype. Lines connecting two circles indicate the enrichment level of observing switches in the two corresponding Ig subclasses. The enrichment level is the ratio of observed and expected switches if assuming Ig isotypes are independently distributed among B cell clusters. Cancer types were sorted by the total number of CDR3s.

Figure 4 ∣. Features of tumor-infiltrating B cell repertoire compared between tumor and adjacent normal tissues.

Multiple aspects of infiltrating B cells are displayed on the heatmap, including B cell infiltration levels, B cell diversity, and levels of different Ig switches. B cell Ig switch levels were calculated as the number of B cell clusters with the Ig switch divided by the total number of unique CDR3s (Methods). Cancer types with over 100 samples for tumor and over 10 samples for adjacent normal tissues were presented (see Supplementary Table 3 for the numbers). Log fold change (log FC) was calculated as the log₂ ratio of the average values between tumor and adjacent normal tissues. Statistical significance was evaluated using Student’s t-test and two-sided P-values. FDR correction was performed using Benjamini-Hochberg procedure for testing on multiple cancer types.

Figure 5 ∣. Interaction between IgG3-1 switch and somatic hypermutations.

a, Violin plots showing the distribution of somatic hypermutation (SHM) rates across cancer types. The points show the median values, lines extend to 1.5 times of the interquartile range, and the curves show the density of values (see Supplementary Table 3 for samples sizes). b, Violin plot showing the SHM rates for samples with different immunoglobulin subclass switches. From left to right, the order of switches follows the germline gene positions and so BCRs with a subclass switch on the right might arise from other subclass-switch events on the left. Statistical significance was evaluated using Wilcoxon rank sum test (two-sided P < 1 × 10⁻⁴ labeled with ****, < 0.01 labeled with **). The boxes in violin plots show the lower, median, and upper quartile of values, and the curves show the density of values (see Supplementary Table 3 for sample sizes and P-values). c, Kaplan-Meier curves showing the interaction between SHM rate groups and IgG3-1 switch groups. SHM high (n = 2,338) or low (n = 2,341) were split by the median SHM rate in all patient samples. IgG3-1 switch high or low were split by the median IgG3-1 switch level in each SHM rate group. Statistical significance comparing different groups was evaluated using multivariate Cox regression corrected for tumor purity and patient age at diagnosis.

Figure 6 ∣. Potential immune evasion of ADCC through MIC shedding.

a, NK cell activity visualized by box plots of granzyme B and CD16a. Both genes were expressed at higher levels in tumors with high level of isotype switches (CSR). Statistical significance was evaluated using Wilcoxon rank sum test on 4,273 low CSR samples and 4,270 high CSR samples. Two-sided P-values are 6 × 10⁻²⁹⁵ and 5 × 10⁻¹⁵⁶ for GZMB and FCGR3A respectively labeled with ****. b, Heatmap showing the associations between the fractions of IgG isotypes and NK cell Fc receptor FCGR3A (CD16a) expression across multiple cancers. Cancer type was selected based on at least 100 patient samples with no missing values (see Supplementary Table 3 for sample sizes). Each entry in the heatmap represents the partial Spearman’s correlation corrected for tumor purity (Methods). FDR correction was performed on two-sided P-values using Benjamini-Hochberg procedure for testing on multiple cancer types. c, NK cell inactivation through MIC shedding. Interaction between IgG1/3 Fc ligand and CD16a triggers NK response. MICA amplification is observed in 20% of tumors, and metalloproteinase overexpression may lead to the production of soluble MIC in these tumors, which results in internalized NKG2D (iNKG2D) and inactivated NK cell. d, Interaction between MICA amplification and IgG1/3 levels visualized by Kaplan-Meier curves for breast invasive carcinoma (BRCA) and skin cutaneous melanoma (SKCM), respectively. Tumor samples were categorized into 4 groups based on the presence of MICA amplification and IgG1/3 levels. MICA Amp refers to tumors with MICA amplification, and MICA WT refers to the remaining. High or low IgG1/3 groups were split by the median number of IgG1 or IgG3 B cell clusters divided by the total number of unique CDR3s in each MICA groups. There were 579 MICA WT samples and 211 MICA Amp samples for BRCA, and 126 MICA WT samples and 130 MICA Amp samples for SKCM. Hazard ratios and statistical significance were evaluated using Cox proportional hazard regressions corrected for tumor purity and patient age.