Large-scale microarray profiling reveals four stages of immune escape in non-Hodgkin lymphomas

ABSTRACT

Non-Hodgkin B-cell lymphoma (B-NHL) are aggressive lymphoid malignancies that develop in patients due to oncogenic activation, chemo-resistance, and immune evasion. Tumor biopsies show that B-NHL frequently uses several immune escape strategies, which has hindered the development of checkpoint blockade immunotherapies in these diseases. To gain a better understanding of B-NHL immune editing, we hypothesized that the transcriptional hallmarks of immune escape associated with these diseases could be identified from the meta-analysis of large series of microarrays from B-NHL biopsies. Thus, 1446 transcriptome microarrays from seven types of B-NHL were downloaded and assembled from 33 public Gene Expression Omnibus (GEO) datasets, and a method for scoring the transcriptional hallmarks in single samples was developed. This approach was validated by matching scores to phenotypic hallmarks of B-NHL such as proliferation, signaling, metabolic activity, and leucocyte infiltration. Through this method, we observed a significant enrichment of 33 immune escape genes in most diffuse large B-cell lymphoma (DLBCL) and follicular lymphoma (FL) samples, with fewer in mantle cell lymphoma (MCL) and marginal zone lymphoma (MZL) samples. Comparing these gene expression patterns with overall survival data evidenced four stages of cancer immune editing in B-NHL: non-immunogenic tumors (stage 1), immunogenic tumors without immune escape (stage 2), immunogenic tumors with immune escape (stage 3), and fully immuno-edited tumors (stage 4). This model complements the standard international prognostic indices for B-NHL and proposes that immune escape stages 3 and 4 (76% of the FL and DLBCL samples in this data set) identify patients relevant for checkpoint blockade immunotherapies.

KEYWORDS:

• Data mining
• immune escape
• lymphoma

Abbreviations

CLL	=	chronic lymphocytic leukemia
DLBCL	=	diffuse large B-cell lymphoma
DLBCL ABC	=	activated B-cell subtype of DLBCL
DLBCL GCB	=	germinal center B-cell subtype of DLBCL
FL	=	follicular lymphoma
KS	=	Kolmogorov–Smirnov
SES	=	sample enrichment score
ssGSEA	=	single sample gene set enrichment analysis
W	=	Wilcoxon

Introduction

Cancers must evolve strategies of immune escape to develop within immunocompetent hosts. In B-NHL, genetic and epigenetic alterations create tumor phenotypes that are protected against immune cytolysis.^1-7 Immune escape also evolves through the concerted upregulation of gene expression.⁸ For example, the B-NHL subtypes FL and DLBCL upregulate the expression of genes involved in the PD1/PDL1-2 axis, the CTLA4 ligand axis, the biosynthesis of immunosuppressive Galectin 3, and the genes IDO1, VEGFA, PGE2, IL10, and HGF, among several others. Although the immune escape strategies in B-NHL may vary between individuals, the phenotypic and transcriptomic analysis of their biopsies consistently identifies a combination of several pathways.⁸

This issue is particularly relevant when considering the ongoing trials of immune checkpoint-targeting therapeutic regimens for these aggressive lymphomas.^9-13 To maximize the effects of these innovative therapies, the determinants of response in B-NHL patients must be better understood. In particular, the various immune escape pathways might constitute a mechanism for immune resistance and beyond, to immune checkpoint therapeutic resistance. Thus, the development of tools to assess the global status of immune escape in B-NHL patients is required.

Because gene expression profiles reflect the function and activity status of a cell, transcriptome data mining on a large scale allows us to detect genes that are differentially expressed between groups of samples.^14-18 The functional significance of these differentially-expressed gene sets is classically determined by gene set enrichment analysis (GSEA).¹⁹ This method requires a priori-defined groups of samples (e.g., group A vs. group B), and therefore cannot analyze single samples. In addition, it produces excessive false discovery rates and p values since it inappropriately uses the Fischer exact test.²⁰ Single sample GSEA (ssGSEA) have also been used although carry the same statistical flaw,^21,22,23 and their long computing time is an issue for large datasets. Another method for ssGSEA is based on z scores per gene set per sample ^24,25 and avoids this pitfall but still requires cohorts of samples to compute each sample score.

Here, we report a novel and algorithmically-optimized method for scoring the enrichment of pre-defined gene sets in single samples without pitfalls such as those described above. We used this method to score the enrichment of various gene sets, including a gene set for immune escape, in ∼1500 transcriptomes from B-NHL biopsies. This was validated by comparing these results with the main phenotypic hallmarks of B-NHLs, and allowed us to detect four stages of immune escape in B-NHL.

Results

Quality control for SES

We downloaded and assembled a series of 1,446 transcriptomes from 7 distinct B-NHL and 32 samples of normal B-cell controls. Each transcriptome was rank-ordered by decreasing expression levels and scored against a range of gene sets (downloaded from the specified sources: see Methods section). The score for the enrichment of a gene set in a sample (SES) was defined as: SES = (−log₁₀ (p)), using p values from testing the alternative: “empirical distribution of the gene set in the transcriptome greater than uniform distribution”. To determine which of the Kolmogorov–Smirnov (KS), t-test (T), and Wilcoxon (W) tests yielded SES that best matched the phenotype, these tests were compared as follows.

First, the phenotypic data indicated that this data set included both highly mitotic samples (centrocytes, centroblasts, lymphomas, and lymphoid tissue-derived CLL biopsies; n = 796) and weakly mitotic samples (naive B-cells, blood-purified B-cells, and blood-derived CLL and HCL; n = 635). This dichotomy allowed us to compare the KS, T, and W tests to determine which was best for discriminating between these two classes, using SES calculated for the “mitosis”, “cell proliferation”, “cell division”, and “DNA replication” gene sets. Comparing the respective ROC curves revealed that the greatest area under the curves was with the W test for all four gene sets, indicating that its SES most accurately reflected the mitotic phenotype (Fig. 1A, Table 1).

Figure 1. Quality control of the sample scoring test. (A) Performance of scoring methods by the Kolmogorov–Smirnov (KS), T-test, and Wilcoxon tests, compared by ROC curves for the gene sets “DNA replication”, “cell division”, “cell proliferation”, and “mitosis” in the cohort of (n = 1446) B-NHL transcriptomes split according to their proliferating or non-proliferating status (determined by phenotypic analysis). (B) Dot plots of sample enrichment scores (SES) for “female” versus “male” gene sets in the (n = 553) gender-annotated samples within the dataset, using the Wilcoxon test (top), or KS-, T-test-, or Wilcoxon-based scoring (bottom). The rate of discordant cases (DC) for molecularly-inferred versus clinically-annotated gender is indicated. (C) 3D plot of the (n = 553) subtype-specified DLBCL samples using Wilcoxon-based scores for the “ABC”, “GCB”, and “PMBL” gene sets (see Fig. S2).

Table 1. AUC for ROC curves.

Gene set (Database)	AUC KS test	AUC T test	AUC W test
“mitosis” (GO)	0.98	0.97	0.97
“DNA replication” (KEGG)	0.87	0.87	0.91
“cell division” (GO)	0.91	0.96	0.96
“cell proliferation” (GO)	0.81	0.96	0.95

Second, we considered gene set-defined gender. Sexually dimorphic genes defined the “male” and “female” gene sets (Table S2) ²⁶ and their respective SES were computed for the 553 gender-annotated samples. The plots for female versus male allowed us to partition most samples correctly and showed the lowest discordance rate (<2%) with the W test, the same value as has been reported elsewhere (Fig. 1).¹⁴ Whether the discordance between the clinically-annotated and molecularly-inferred gender reflected errors of scoring or of annotation could not be determined. However, this method allowed us to infer gender for all of the other unannotated samples, yielding no discordance for the DLBCL of testicular subtype (not shown).

Third, DLBCL are comprised of distinct subtypes such as activated B-cell (ABC), germinal center B-cell (GCB), and primary mediastinal large B-cell lymphoma (PMBL), which are each defined by specific gene expression patterns.^27-29 The W test was used to score each of the 553 DLBCL samples for the three subtype-specific gene sets, and yielded plots which partitioned these three groups (Fig. 1C). On this basis, the subtype of the other DLBCL samples annotated by tissue of origin or listed as unclassified could also be unambiguously assigned (Fig. S2).

Finally, this ssGSEA method discriminated between relevant and irrelevant (simulated) sets of data and was not sensitive to different sources of RMA-normalized data (Fig. S3). Therefore, these control studies qualified the Wilcoxon test-based scoring of gene set enrichment in single transcriptomes. Although its performance cannot be assessed a priori for any new gene sets, this method scored each sample independently and yielded scores matched to phenotype.

SES detect hallmarks of B-NHLs

We then scored for enrichment of gene sets from the KEGG and GO databases that correspond to well-known phenotypic hallmarks of B-NHL. First, we assessed for proliferation: for the gene set “Mitosis”, the mean SES was 4.8 for centrocytes and centroblasts, 0.3 for naive and blood B-cells, 0.5 for blood CLL, 0.8 for HCL, 5.3 for all DLBCL samples, 7.2 for BL, 2.5 for FL, 2.3 for MCL, and 1.8 for MZL (Fig. 2A). These results matched the respective mitotic phenotypes of such samples, as B-cells from healthy donor blood are non-proliferating in contrast to the highly mitotic B-centrocytes and centroblasts from lymphoid tissues. CLL and HCL blood samples also do not usually have a proliferative phenotype,³⁰ with a higher level of mitosis occuring in nodal CLL. In contrast, lymphoma samples are usually very mitotic, with BL reported to have nearly 100% of cells proliferating.³¹ Similar results were obtained by scoring other gene sets related to cell proliferation (Fig. S4A).

Figure 2. SES of the B-NHL data set matched to sample phenotype. Sample enrichment scores (SES) were computed for each of the 1,446 samples of the B-NHL data set and for the gene sets indicated on the y-axis. Each sample is shown by a dot, red bars are the means of the specified group. CC: centroblasts, CB: centrocytes. (A) SES for gene sets according to proliferation phenotype. (B) SES for tight junctions with sample groups split according to sample preparation. (C) SES for the KEGG gene sets “BCR signaling”, “glycolysis”, and “oxydative phosphorylation”, classified according to lymphoid malignancy.

Second, we assessed the cell adhesion phenotype characteristic of lymphoma in our transcriptomic dataset. The SES of the “tight junction” gene set were higher in lymphoma biopsies than in purified B-cell samples (mean SES = 1.7 vs. 0.6, respectively, unpaired Student's p < 10⁻⁴⁰), as expected ³² (Fig. 2B).

Third, different signaling and metabolic pathways were assessed. As expected, all samples scored high for the “BCR signaling pathway” gene set, although MYC-induced BL scored lower than the others ³³ (mean SES = 9.2 for BL vs. 16.5 for all others, p = 0.003) (Fig. 2C). Glycolysis and oxidative phosphorylation are metabolic hallmarks of cancer cells, and in accordance with this the SES for “oxidative phosphorylation” were high in all samples. The SES for “glycolysis and gluconeogenesis” were low in CLL, HCL, MCL, and resting normal B-cells, whereas such scores were significantly higher in normal proliferating B-cells, BL, DLBCL, FL, and MZL, in line with results from ¹⁶FDG-PET imaging of these tumors.^34-36 Combining the SES for “Glycolysis” and “OxPhos” illustrated the distinct metabolic profiles of the different neoplasms (Fig. 2C and Fig. S4B). Scores for the additional gene sets “acetyltransferase activity”, “histone modification”, and “T-cell activation” also matched the hallmarks of B-NHL (Fig. S4C).

Fourth, leucocyte infiltration of cancer biopsies can be assessed by the gene expression-based quantification of 22 gene sets composed of differentially-expressed genes (known as TM22).^14,18,23 These 22 gene sets were scored against the data set and the results were compared for each type of B-NHL. In agreement with previous reports,^14,18 all tumors presented the signature of their B-cell-of-origin and had very low mastocyte and granulocyte scores. The presence of other leucocyte subsets differed according to the disease (Fig. 3). On average, the DLBCL samples showed similar levels of lymphoid (T and NK cells) and myeloid (monocytes, macrophages, and DC) signatures, with M1 macrophages dominating in the latter. The CLL samples comprised 4-fold lower myeloid levels than those found for DLBCL. In contrast, the FL, MCL, and MZL samples had the highest T-cell compartments, as expected.³⁷ A similar procedure comparing the various subtypes of DLBCL evidenced a significant deficit of leucocyte infiltrates in central nervous system DLBCL and showed the divergent cellular and metabolic profiles of GCB and ABC (Fig. S5).

Figure 3. Scores for B-NHL tumor-associated leucocyte subsets. Distributions of sample enrichment scores (SES) for 22 leucocyte gene subsets across the B-NHL dataset, pooled by histological type of malignancy.

Scoring immune escape in B-NHL

B-NHL express immune checkpoint PD1 and CTLA4 markers that can be targeted by new therapeutic antibodies (see refs.^9,10 for review). In B-NHL other markers of immune escape frequently co-evolve with these, as shown through immunohistochemistry analysis of tissue microarrays and the expression of these immune escape markers in DLBCL and FL biopsies.⁸ We compared the expression of these immune escape genes at the single gene level in control B-cells, FL, and DLBCL biopsies using the Oncomine resource, and found a significant upregulation of several genes in both the FL and DLBCL lymphoma samples (listed by decreasing significance relative to controls): CCL2, IL6ST, IDO1, TIMP1, LGALS3, VEGFA, HAVCR2, MRC1, TIGIT, CD163, IL10, PDCD1LG2, CTLA4, LAG3, LGALS1, CSF1, MSR1, JAK2, SOCS3, CD274, ICOS, HGF, IL23A, GDF15, FOXP3, PVR, MCL1, PDCD1, CCL22, LAIR1, CD86, IDO2, and KIR2DL1. Each of these genes was upregulated to a different extent in each sample, as illustrated by the fact that the most upregulated gene was either MRC1 or PTGS2 in two different samples of FL, and CTLA4 or LAG3 in two different samples of DLBCL (Fig. S4). Hence, despite the lack of a specific gene expression signature, these distinct patterns suggested a functional convergence toward immune escape.

We therefore defined a gene set (IEGS33, Table S3) that encompassed the above 33 most shared and upregulated genes to investigate their collective contribution to immune escape at the single sample level in B-NHL. The SES for IEGS33 was computed for the 1,446 B-NHL dataset samples. A preliminary check-up of the distribution of these scores with these actual or 500 simulated samples confirmed that IEGS33 is not a random gene set and carries a significant amount of information about these samples (Fig. 4A). In the B-NHL cohort, we found that 24.9% of samples had a significantly higher IEGS33 score, while healthy B-cells did not (mean SES = 0.3). In addition, most CLL and HCL scored nearly as low as normal B-cells for IEGS33 (group mean of SES = 0.2 for CLL and HCL), whereas few CLL blood samples were elevated, confirming the heterogeneity of immune escape in these diseases.^38,39 The IEGS33 scores of BL were slightly higher (mean SES = 0.6), suggesting that these lymphoma escape immunity through different pathways.^40,41 Low scores were also found in most but not all MCL (mean SES = 0.4 and 1 for MCL and MZL, respectively), consistent with earlier reports ^42,43 that included PD1 defects and autoimmunity in MZL.^44,45 In contrast, IEGS33 scores were significantly elevated in 73% of FL (mean SE = 1.4) and 78% of DLBCL (mean SES = 1.6) (Fig. 4B).⁸ Likewise, the DLBCL ABC, as well as PMBL samples, generally scored higher than the DLBCL GCB (mean SES = 1.6 vs. 1.3, unpaired Student p = 0.0008 for ABC vs. GCB) (Fig. 4C). In addition, the DLBCL samples that were clinically annotated as “nodal” and “extra-nodal testicular” had the highest IEGS33 scores and were of ABC subtype (Fig. S2). Although primary CNS DLBCL were also of ABC subtype, their IEGS33 scores were the lowest among the DLBCL samples (Fig. 4C). Similar results were obtained with the IEGS51 gene set depicted elsewhere.⁸

Figure 4. Scores for the immune escape gene set IEGS33 in 1,446 B-NHL samples. (A) Quality control of the IEGS33 gene set. The right skewed distribution of SES for IEGS33 in the 1446 B-NHL samples or in 500 simulated samples (shown in Fig. S2) demonstrates that this gene set carries relevant information about these samples. For B-NHL samples 24.9% of SES are above 1.3 whereas less than 5% are predicted for the simulated samples. (B–C) IEGS33 scores of the B-NHL samples classified by malignancy (B) and subtype of DLBCL (C). Red line shows group means. CC: centroblasts, CB: centrocytes.

Together, these results pinpoint qualitative and quantitative differences in the immune escape strategies of NHL, with three quarters of FL and DLBCL upregulating the collective expression of immune escape genes.

Four stages of immune escape in NHL

Since the activation of immune effectors represents the substrate of immune escape, we first analyzed the relationship between these features by plotting the SES of the IEGS33 gene set versus the “T-cell activation” gene set. We arbitrarily set gates to define the following groups of samples: group I (T-cell activation⁻ IEGS33⁻), group II (T-cell activation⁺ IEGS33⁻), group III (T-cell activation⁺ IEGS33⁺), and group IV (T-cell activation⁻ IEGS3⁺). As expected, the normal B-cell samples fell into group I. Most CLL and HCL samples were also in group I, with the others in group II. BL samples were in groups I and IV only. Of the MCL samples, 57.4% were in group I, 30% in group II, and 12% in group III, with none in group IV. In contrast, the MZL, FL, and DLBCL samples were distributed between all four groups. For MZL, 25% of samples were in group I, 16.7% in group II, 50% in group III, and 8% in group IV. For FL, these figures were 13.4%, 6.7%, 72.2%, and 0.7% in groups I–IV, respectively, and for DLBCL they were 19.6%, 2.3%, 53.3%, and 24.8%, respectively. Among the three subtypes of DLBCL, we observed 21.4% of ABC, 23.3% of GCB, and 12.5% of PMBL samples in group I, 0.7% of ABC, 5% of GCB, and 0% of PMBL samples in group II, 57.9% of ABC, 51.6% of GCB, and 66.7% of PMBL samples in group III, and 20% of ABC, 20.1% of GCB, and 20.8% of PMBL samples in group IV (Fig. 5A).

Figure 5. Stages of immune escape in B-NHL. (A) SES dot plots for IEGS33 versus T-cell activation for the 1,446 samples. CC: centroblasts, CB: centrocytes. Groups of phenotypes were arbitrarily defined according to the dotted lines: group I (IEGS33⁻ T-cell activation⁻), group II (IEGS33⁻ T-cell activation⁺), group III (IEGS33⁺ T-cell activation⁺), group IV (IEGS33⁺ T-cell activation⁻). (B) Left: four groups. Right: distribution of samples in these groups (group I in black, group II in green, group III in red, and group IV in blue) shown by malignancy. (C) Overall survival (OS) of the (n = 580) DLBCL patients according to the four immune escape groups (Log rank p = 0.04).

To further characterize these four groups of samples, we compared their respective scores for various other gene sets, restricting this analysis to the 418 DLBCL samples to avoid interpretation bias. Samples from all groups scored similarly for most metabolic and housekeeping gene sets, as well as for B-cell signature gene sets, typifying tumor cells. For mitosis gene sets, samples from groups I and IV had significantly higher scores than groups II and III (Student's p <0.05), indicating higher proliferation rates. The same pattern was found for gene sets related to genetic instability, such as “DNA repair”, “response to DNA damage stimulus”, and “DNA integrity checkpoint”, indicating higher mutation rates. The level of total leucocyte infiltration was significantly lower in group I than elsewhere (Student's p < 10⁻¹⁰), and was highest in groups II and III, with the “T-cell-specific” and “NK cell-specific” scores also highest in groups II and III. In parallel, the scores for “apoptosis” and “programmed cell death” were highest in groups II and III. However, myeloid infiltration (including the macrophage and DC subsets) scored significantly higher in groups III and IV (Student's p < 10⁻²⁷) (Table 2).

Table 2. SES in the four groups defined by T-cell activation and IEGS33.

	SES (mean of group +/−SD)
Gene set	I			II			III			IV
Kegg Glycolysis	1.9	+/−	0.9	2.2	+/−	0.7	2.1	+/−	0.9	2	+/−	1.1
Kegg OxPhos	31.6	+/−	0.9	30.8	+/−	0.7	31.2	+/−	0.9	31.8	+/−	1.1
Kegg Ribosome	51.8	+/−	0.9	51	+/−	0.7	51.8	+/−	0.9	52.2	+/−	1.1
B cells (all)	14.6	+/−	4.6	17.2	+/−	3.8	13.4	+/−	4.8	12.5	+/−	3.8
Mitosis	6.3	+/−	1.5	4.2	+/−	1.2	4.6	+/−	1.4	6	+/−	1.5
Cell Cycle	10.2	+/−	2.3	6.3	+/−	2	7.2	+/−	2.1	9.4	+/−	2.3
Kegg DNA Replication	10.4	+/−	0.9	8	+/−	0.7	8.4	+/−	0.9	10	+/−	1.1
DNA Repair	8.3	+/−	1.9	5.0	+/−	1.2	5.6	+/−	1.7	7.4	+/−	7.4
Response To DNA Damage Stimulus	8.2	+/−	1.9	5.0	+/−	1.3	5.4	+/−	1.7	7.3	+/−	7.3
DNA Integrity Checkpoint	1.3	+/−	0.3	0.6	+/−	0.2	0.9	+/−	0.3	1.2	+/−	1.2
Apoptosis	7.3	+/−	1.1	8.3	+/−	0.8	8.5	+/−	0.9	7.7	+/−	1
Programmed Cell Death	7.4	+/−	1.1	8.4	+/−	0.8	8.5	+/−	0.9	7.8	+/−	1
T cells (all)	3.8	+/−	3.3	12.1	+/−	5.2	28.2	+/−	12.3	8.0	+/−	5.4
NK cells (all)	0.4	+/−	0.5	1.6	+/−	0.9	4.2	+/−	2.7	1.1	+/−	1.2
Mφ (all)	2.2	+/−	1.3	2.8	+/−	1.1	8.5	+/−	4.3	5.7	+/−	3.0
DC (all)	0.8	+/−	0.5	0.9	+/−	0.5	3.7	+/−	2.1	2.3	+/−	1.5

Bold text indicates group values significantly higher than the other groups for the same gene set (Student's p < 0.05).

Together, these group characteristics were consistent with the following four stages of cancer immune editing:

• Stage 1: Non-immunogenic tumor. Group I scores showed that cells were proliferating and genetically unstable but with very low immune infiltrates;

• Stage 2: Immunogenic tumor without immune escape. Compared to stage 1, group II scores were increased for “T- and NK cell-specific”, “T-cell activation”, and “programmed cell death” gene sets, but decreased for “proliferation” and “DNA repair” gene sets, presumably due to tumor lysis;

• Stage 3: Immunogenic tumor with immune escape. Group III scores were as those for stage 2 but with higher scores for the “IEGS33” and “myeloid cell infiltrates (monocytes, macrophages, DC)“ gene sets;

• Stage 4: Fully immuno-edited tumor. Group IV scores were as those for stage 3, but collapsed for “T- and NK cell infiltrates”, “T-cell activation”, and “programmed cell death”, while scores for “proliferation” and “DNA repair” were increased. Stage 4 is very similar to stage 1 but with immune infiltrates.

This model of four successive stages of cancer immune editing predicts a worse clinical outcome for B-NHL patients with a tumor score characterizing either the earliest or the most advanced stages of immune escape, since both lack effective immuno-protection. Using the 418 clinically annotated DLBCL samples from the data set, we compared the clinical outcomes of DLBCL patients with tumors in the different groups described above. This clearly showed that overall survival at 5 y (OS) was longest for group II and shortest for group I. In agreement with the four stage model, OS gradually decreased from group II (OS = 0.80, n = 10) to group III patients (OS = 0.62, n = 239) and group IV patients (OS = 0.58, n = 87) (Log Rank p = 0.001), with group IV patients having the same OS rates as those in group I (OS = 0.50, n = 82). Similar results were obtained for OS at 10 y (p = 0.015) and 20 y (p = 0.04) (not shown). Together these results fully support a model with four stages of cancer immune editing in B-NHL.

Finally, we correlated our four-stage immune escape model with the standard DLBCL international prognostic index (IPI) scoring system and found that all cases classified within the favorable immune escape stage 2 had a low IPI. Reciprocally however, all of the IPI classes each had a similar distribution of immune escape stages (Fischer exact test p = 0.6) (Table S4). In addition, DLBCL of GCB subtype encompassed significantly higher rates of stage 2 cases than ABC and other subtypes (Fischer exact test p = 0.04). Hence, immune escape staging is distinct from but complements the standard prognostic scores of DLBCL.

Discussion

To conclude, we used a novel data mining method that scored the enrichment of gene sets in single microarrays, and detected four stages of cancer immune editing in humans, an issue so far deemed to be “a significant technical challenge” with predictive and prognostic significance.⁴⁶ Current ssGSEA methods rely on clinical classification or prior information on sample phenotype, and require either cohorts to z-score each sample ^24,25 or the application of statistics involving equiprobable genes as inadequate null hypotheses.^21-23,47 Our method avoids such pitfalls and its algorithmic optimization allows for the rapid analysis of large series of data, such as labtop-based scoring within seconds for 1,500 gene sets in 1,446 microarrays. Although RNASeq or proteome datasets are rich sources of material for such studies, nearly two decades of microarray-based research in cancer means that large and precious amounts of this publicly-available resource have accumulated. Collecting large series of microarray data from GEO data sets involves collating many different datasets, whose eventual inconsistencies might introduce confounding divergences.⁴⁸ Here, by associating robust statistics and selecting compatibility-based data sets, we were able to assemble 1,446 consistent transcriptomes from 33 different datasets.

Although this approach is capable of analyzing other types of large-scale data sets from NGS, the microarray-based study reported here is able to detect the hallmarks of B-NHL and their specific stages of immune editing. Indeed, immune escape by B-NHL involves not only changes in gene expression but also somatic mutation-driven pathways.⁵ Our technique cannot detect somatic mutations including those that target immune activation, the downregulation of immuno-activating genes, or the upregulation of genes not in IEGS33. Nevertheless, it can detect the first stage of immune escape in DLBCL, which is characterized by high proliferation and mutation rates but a significantly low level of immune cell infiltration. This suggests that tumors in this first stage lack immunity due to an impaired primary activation. Most HLA-mutated B-NHL tumors^1,3,5,6 presumably lie at this stage, which presents the worst clinical outcome (Fig. 5). However, this might not be the case for B-NHL tumors at stage 2, which show the best survival rates and transcriptional hallmarks of functional cytotoxic immunity: cell infiltration prominently by T- and NK cells, T-cell activation and an increase in cell death alongside decreased proliferation and genetic instability. Stage 2 DLBCL patients were also found to have low IPI scores (Table S4), supporting a progressive model. However, the low number of cases meant that this finding was not significant, hence from this data we conclude that immune escape stages and IPI scores are related but distinct prognostic criteria. The stage 2 pattern of immune activation was also found in stage 3, but was accompanied by obvious signs of immune escape, as indicated by increased scores for IEGS33 and myeloid cell infiltrates. When considering cancer immune editing in terms of Darwinian evolution, through elimination, equilibrium and escape phases,⁴⁹ stage 3 would correspond to the end equilibrium-early escape and, in line with this, DLBCL patients at this stage have shorter survival rates than those at stage 2 (Fig. 5). This progression can lead to fully immuno-edited stage 4 tumors whose immune escape collapses the host-protective immunity and allows the re-emergence of highly proliferating and mutated tumors despite the presence of a significant but non-protective immune cell infiltrate. Stage 4 clearly corresponds to the late escape phase of the cancer immuno-editing process,⁴⁶ and DLBCL patients at this stage have shorter OS than those at stage 3.

It could be argued that B-NHL tumors which are either immunologically undetected or detected by an ineffective immunity should ultimately grow in the same way, a view supported by the nearly similar OS rates of DLBCL patients at stages 1 and 4 (Fig. 5). Nevertheless, these stages represent clearly distinct contexts since in stage 4 (but not stage 1) cellular targets for the emerging immune checkpoint therapeutics are present. We report that a total of 78% of DLBCL, 73% of FL, and 58% of MZL are tumors at stages 3 and 4 of immune escape. We thus propose that the immune escape stages 3 and 4 constitute a major criterion to stratify B-NHL patients to identify those most likely to benefit from immune checkpoint blockade therapies.

Materials and methods

Datasets

We first screened the GEO ⁵⁰ repository for gene expression profiles of B-NHL obtained with the Affymetrix HG U133 plus 2.0 microarray and downloaded 47 series of RMA-normalized data sets (GSE) as txt files, together with gene sets from Gene Ontology,⁵¹ KEGG (http://www.genome.jp/kegg/pathway.html), the Molecular Signatures Data base (MSigDB 3.0) ⁵² and TM22,^14,18 and the gene sets defined in the text. Each series of transcriptomes was then collapsed to common protein-encoding genes (using the HUGO nomenclature) and the compatibility for merging of the datasets was assessed. The gene expression data of each sample (GSM) were transformed into van der Waerden's scores,^53,54 then each data set was reduced to a per gene median of samples and the pairwise correlation matrix of these 47 datasets was computed (Fig. S1).⁴⁸ The 33 data sets with correlations >0.8 were deemed compatible for merging while others less homogeneous were discarded (Table 1). Since these 33 datasets comprised malignant samples as well as various normal cell controls, compatibility screening was repeated at the sample level for the GSM groups from either all malignant samples or all non-cancer cell controls (Fig. S2). The samples deemed irrelevant for this meta-analysis (e.g., those from a microarray-based study of in vitro drug activity on cell lines) were rejected. This yielded a data set composed of 1,446 samples from 7 distinct malignant histologies. These included 489 DLBCL, 149 FL, 125 MCL, 12 MZL, 4 Burkitt's lymphoma (BL), 630 chronic lymphocytic leukemia (CLL), 5 hairy cell leukemia (HCL), and control cells. The control cells included 32 normal B-cell samples including CD20 B-cells purified from peripheral blood, naive B-cells, centrocytes, and centroblasts purified from the tonsils of non-cancer donors.

Sample enrichment scores (SES)

The RMA-normalized expression data from each collapsed sample (GSM) were rank-ordered by decreasing the expression level along the [0–1] segment (highest expression level = 0, lowest expression level = 1), to yield v, a vector of ranks. Since each series from the NHL data set had been pre-selected for a high pairwise correlation (see above), their rank-ordered gene distributions were also highly correlated, despite distinct batch RMA normalizations (Fig. S3).

We rationalized that if gene set g was a subset of v then assessing the enrichment of g would mean testing g versus v. However the human genome comprises many pleiotropic genes that are involved in various different pathways, functions and cell compartments and so are more frequently upregulated in all samples and are over-represented in gene sets and database.²⁰ Since genes are far from equiprobable in gene sets, the adequate null hypothesis for g versus v gives genes with probabilities proportional to their frequency in the gene set database ⁵⁵ and requires the KS, T-, and W-tests with frequency-corrected null hypotheses (see refs.^56-58 for review). In addition, scoring the enrichment of the ∼1500 KEGG and GO gene sets in >1,000 samples would mean computing ∼1.5 million tests with a high accuracy for the most significant results. Therefore, algorithmic optimization is necessary to i) compute scores in a reasonable time, ii) precisely compute the most significant (smallest) p values, iii) allow the gene frequency correction, and iv) allow testing of the three alternatives, namely “greater” (indicating an enriched gene set), “less” (a downregulated gene set), and “two-sided” (when the gene set is enriched or downregulated). We developed an R script for which the source codes are available at: https://sites.google.com/site/fredsoftwares/products/autocompare_ses. The current version of it computes the enrichment of ∼1500 gene sets within 0.1 sec for a single sample.

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Acknowledgments

We acknowledge all members of JJF's team for their stimulating comments on this work and Dr Kelly Thornber for English editing of the text.

Funding

This work was supported in part by institutional grants from the Institut National de la Santé et de la Recherche Médicale (INSERM), the Université Toulouse III: Paul Sabatier, the Center National de la Recherche Scientifique (CNRS), the Laboratoire d'Excellence Toulouse Cancer (TOUCAN) (contract ANR11-LABX) and the Program Hospitalo-Universitaire en Cancérologie CAPTOR (contract ANR11-PHUC0001).