Follicular lymphoma B cells exhibit heterogeneous transcriptional states with associated somatic alterations and tumor microenvironments

Abstract

Follicular lymphoma (FL) is an indolent non-Hodgkin lymphoma of germinal center origin, which presents with significant biologic and clinical heterogeneity. Using RNA-seq on B cells sorted from 87 FL biopsies, combined with machine-learning approaches, we identify 3 transcriptional states that divide the biological ontology of FL B cells into inflamed, proliferative, and chromatin-modifying states, with relationship to prior GC B cell phenotypes. When integrated with whole-exome sequencing and immune profiling, we find that each state was associated with a combination of mutations in chromatin modifiers, copy-number alterations to TNFAIP3, and T follicular helper cells (Tfh) cell interactions, or primarily by a microenvironment rich in activated T cells. Altogether, these data define FL B cell transcriptional states across a large cohort of patients, contribute to our understanding of FL heterogeneity at the tumor cell level, and provide a foundation for guiding therapeutic intervention.

Keywords: B cell; follicular lymphoma; genomics; germinal center; transcriptome; tumor microenvironment.

Graphical abstract

Figure 1

FL B cell transcriptional groups (A) Schematic of samples used for data generation. Dotted lines indicate patient-matched data/samples. (B) Consensus clustering map result of 200 NMF runs on FL and benign B cell (n = 87) log2(TPM+1) gene expression. The heatmap indicates the sample-sample pairing frequency. Three groups were identified: group 1 (red, n = 20), group 2 (blue, n = 24), and group 3 (green, n = 43). Tissue type bar displays sample origin: FL (gray) and benign (blue). (C) Proportional association between diagnostic clinical attributes and each B cell NMF group. Fisher’s exact test ∗p < 0.05. (D) B cell group signature scores from the 23-gene high-risk signature and IR-1/IR-2 signatures on matching bulk, unsorted tumor RNA-seq samples. The y axis represents the average value of signature genes’ Z scores with each sample. Independent Wilcoxon rank-sum tests; ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001.

Figure 2

Divergent FL B cell gene expression states highlight unique biological activities (A, left) Gene contribution scores for each group state. (See STAR Methods.) Higher/more significant positive contribution (red) and lower/more significant negative contribution (blue). (Center) Relative gene expression heatmap (n = 92) of group signature-defining genes (n = 500 genes/group). Selected genes of interest are highlighted at right. (Top) Individual sample state coefficient values for each factor/state from the NMF H matrix. (B) Normalized enrichment scores (NESs) from selected significant gene sets for each group state. Gene sets are grouped vertically based on the group with which they share the most significance. (C) GSEA plots of significant gene sets for group states 1 (top row), 2 (center row), and 3 (bottom row). (D) VIPER regulator activity enrichment plot for group states. Selected top significant results (p < 0.01) are shown for positive enrichment (top proteins, red bar) and negative enrichment (bottom proteins, blue bar). Red lines depict genes positively regulated by the protein, and blue lines depict genes negatively regulated by the protein. Genes are ordered according to their group state contribution score (bottom bar plot).

Figure 3

FL B cell programs associate with independent GC phenotypes (A) Weighted voting prediction scores of GC B cell cluster gene sets for each sample. Empirically derived null distribution 95% confidence interval (CI) from imputed data, in gray. Scores range from −1 to 1, with 0 representing equal number of votes for and against a gene set in a sample. (B) Heatmap of normalized RNA-seq values of FL B cell samples (n = 92), for the top 50 genes in each gene set. Single cell clusters are listed on the left side of the map and FL B cell grouping assignments are represented on the top bar. Samples are ordered based on the highest GC B cell cluster score. (C) GC B cell cluster enrichment analysis from hypergeometric testing of B cell group assignment. Bars represent the false discovery rate-corrected p values (q-val). Positive and negative associations (residuals) are plotted directionally with –log(q) values. Gray area represents p = 0.05.

Figure 4

FL B cell states shape their local immune profile (A) Uniform manifold approximation and projection of immune metaclusters from individually profiled FL patient CyTOF samples (n = 60). Clusters identified using Phenograph and are colored by manual annotation. (B) Heatmap of mean surface protein expression in each identified cluster mentioned in (A). Cluster identities listed next to their coordinating colors from (A). (C) Boxplot of immune content as percentage of live cells in each sample by their B cell group. (D) Correlation plot between the difference in CMI vs. INFM state coefficients (positive = higher CMI, negative = higher INFM) and the non-B immune cell percentage (% of live cells). A standard linear regression line is displayed with 95% CIs, along with (r) Pearson correlation, and representative p value. (E) (Left) Heatmap of correlation matrix between cell cluster abundances in all of the samples, depicting cellular communities (black boxes). Cell values from Pearson correlation. Dendrogram from correlation distance. (Right) Heatmap of Pearson correlation between cell abundance and sample B cell group state. Pearson p values are printed where p < 0.1. (F) Heatmap of Spearman correlation between single-sample gene set scoring, using SingScore of published GC stromal gene sets, and B cell group states (n = 54). Pearson p values are printed in cells where p < 0.1. (G) CODEX multiplexed immunofluorescence analysis of CD21 (white), CD20 (blue), Ki-67 (yellow), CD4 (green), and CD8 (magenta) on 1 sample from each B cell group. (Left) Baseline cell type abundance in each sample (% of live) from CyTOF. (Right) First 3 panels, stitched images equal to 10× magnification with 500-μm scale. White box delineating the magnified view in panel 4 of follicle and interfollicular region, with 200-μm scale.

Figure 5

Landscape of somatic alterations in FL (A) GISTIC 2.0 analysis of CNVs identified from WES (n = 123). G-score plot of significant CNA peaks (q < 0.1). (B) Mutational signature trinucleotide single-base-substitution profiles derived from BayesNMF resultant signatures on all SNVs derived from WES (n = 123). Sig1 = aging, Sig2 = unknown, Sig3 = AID-like. (C) Oncoplot of individual mutation and CNV events across FL patients (n = 123). Mutation variants are filtered to >1% variant allele frequency and genes that contain significant driver scores from CHASM and/or VEST. (Left) Driver scores from CHASM and VEST gene level scores. (D) Ratios of clonal to subclonal events by gene/q-band. Analysis performed using ABSOLUTE and variant CCFs > 0.85 were considered clonal.

Figure 6

Somatic alterations associate with B cell states (A) Box and violin plots of TMB by B cell group (n = 119). Pairwise Wilcoxon test p values displayed. (B) Box and violin plots of the Sig3 mutational signature exposure by B cell group (n = 119). Pairwise Wilcoxon test p values displayed. (C) Heatmap of B cell state and mutation mutual association analysis was performed using a generalized linear model between sample B cell state values and genotype profile. Degree of association is determined from linear model coefficients and are colored from negative (purple) to positive (orange). Gray circles depict the association p value, and significant associations are marked with a thick black border (p < 0.05). Only genes with at least 1 association p < 0.1 are shown.