Deconstructing Intratumoral Heterogeneity through Multiomic and Multiscale Analysis of Serial Sections

Abstract

Tumors may contain billions of cells, including distinct malignant clones and nonmalignant cell types. Clarifying the evolutionary histories, prevalence, and defining molecular features of these cells is essential for improving clinical outcomes, since intratumoral heterogeneity provides fuel for acquired resistance to targeted therapies. Here we present a statistically motivated strategy for deconstructing intratumoral heterogeneity through multiomic and multiscale analysis of serial tumor sections (MOMA). By combining deep sampling of IDH-mutant astrocytomas with integrative analysis of single-nucleotide variants, copy-number variants, and gene expression, we reconstruct and validate the phylogenies, spatial distributions, and transcriptional profiles of distinct malignant clones. By genotyping nuclei analyzed by single-nucleus RNA-seq for truncal mutations, we further show that commonly used algorithms for identifying cancer cells from single-cell transcriptomes may be inaccurate. We also demonstrate that correlating gene expression with tumor purity in bulk samples can reveal optimal markers of malignant cells and use this approach to identify a core set of genes that are consistently expressed by astrocytoma truncal clones, including AKR1C3, whose expression is associated with poor outcomes in several types of cancer. In summary, MOMA provides a robust and flexible strategy for precisely deconstructing intratumoral heterogeneity and clarifying the core molecular properties of distinct cellular populations in solid tumors.

Keywords: IDH1; clonal evolution; gene coexpression; intratumoral heterogeneity; low-grade glioma; multiomic; single-nucleus analysis; tumor microenvironment.

Figure 9

Aggregating correlations to tumor purity reveals core transcriptional features of astrocytomas. (a) Gene expression correlations (n = 15,288 genes) to malignant cell abundance in case 1 and case 2. Red and blue denote significantly correlated genes that were used for enrichment analysis (b), and the star denotes AKR1C3. (b) −Log₁₀ FDR-corrected p-values (q-values) from one-sided Fisher’s exact tests analyzing gene set enrichment in red and blue genes from (a). (c) Validated protein–protein interactions (PPI) from STRINGdb [54] for red genes from (a). The 201 proteins shown formed networks of five or more proteins, with the number of interactions equal to the number of edges. (d) −Log₁₀ FDR-corrected p-values (q-values) from one-sided Fisher’s exact tests analyzing gene set enrichment for each STRINGdb interaction cluster in (c). (e,f) AKR1C3 immunostaining in FFPE tissue adjacent to the sectioned region of case 1 (e) and non-neoplastic human brain (f). Image: 200×; scale bar: 50 μm. (g–i) Immunofluorescent co-staining of IDH1 R132H (white), AKR1C3 (green), and nuclei (blue (DAPI)) in case 1 demonstrating expression of AKR1C3 in malignant cells carrying the truncal IDH1 R132H mutation. Scale bar denotes 50 µm.

Figure 1

Overview of MOMA. (a) Schematic of a heterogeneous human brain tumor. (b) Serial sectioning introduces variation in cellular composition. (c) Section usage can be flexibly tailored for diverse multiscale and multiomic assays. (d) Correlative analysis of bulk cellular frequencies and molecular feature levels predicts optimal markers of malignant clones and nonmalignant cell types of the tumor microenvironment. (e) Predictions from bulk analysis are validated by single-nucleus analysis of interpolated sections and histology.

Figure 2

Multiomic analysis of serial tumor sections reveals the clonal composition of a primary grade 2 IDH-mutant astrocytoma (case 1). Axial T2 (a) and axial FLAIR (b) images demonstrate a round, well-defined T2 and FLAIR hyperintense intraaxial left temporoparietal mass that is non-enhancing and consistent with a low-grade glial neoplasm. (c) Image of the frozen tumor sample prior to cryosectioning and nucleic acid isolation. (d,e) Immunostaining for IDH1 R132H (d) and TP53 (e). Images: 400×. Scale bars: 50 μm. (f) Schematic of serial sectioning strategy and section usage plan. Amp-seq = deep sequencing of PCR amplicons spanning mutations identified by exome sequencing. (g) Hierarchical clustering of mutations, using 1—Pearson correlation of amp-seq variant allele frequencies (VAFs) over all tumor sections (n = 69) as a distance measure, reveals three clusters. Amp-seq was performed in two sequencing runs (denoted by bold and regular fonts). (h–k) VAF patterns comprising cluster 1 (h,i), cluster 2 (j), and cluster 3 (k). Cluster 1 was split to illustrate the effects of high (h) and low (i) coverage. (l) Clone phylogeny (with arbitrary branch lengths) derived from integrated analysis of SNVs (from amp-seq data) and CNVs (from DNA methylation data). Percentages represent the average abundance of each cellular fraction over all analyzed sections (n = 68). (m) Estimated cellular fractions for all clones and nonmalignant cells over all sections (n = 68).

Figure 3

Gene coexpression modules are highly correlated with clonal abundance (case 1). (a) Hierarchical clustering of gene coexpression modules over all tumor sections (n = 69). (b) Module eigengenes (MEs) illustrate the relative expression levels of genes in each module over all tumor sections. (c) The number of genes used to form each ME. (d–g) Top left: MEs with the strongest correlations to clonal abundance (defined cumulatively). Locally weighted smoothing (LOESS) lines are shown; correlation is based on data points. Bottom left: the 12 genes with the highest correlations to the ME (k_ME). Right: enrichment analysis of gene coexpression modules using published gene sets. FDR-corrected p-values (q-values) from one-sided Fisher’s exact tests are shown. Positive values represent enrichments of genes that were significantly positively correlated to the ME, while negative values represent enrichments of genes that were significantly negatively correlated to the ME. Gene sets representing chromosomal gains or losses include all genes within affected regions (as described in Figure 2l and Table S5). See Table S9 for descriptions and sources of featured gene sets.

Figure 4

Multiomic analysis of serial tumor sections reveals the clonal composition of a recurrent grade 2 IDH-mutant astrocytoma (case 2). Axial T2 (a) and axial FLAIR (b) images demonstrate a non-enhancing, expansile, infiltrating glioma centered in the right insula and involving the basal ganglia, inferior frontal lobe, and temporal lobe. Cystic degeneration was present in the tumor. (c) Image of the frozen tumor specimen prior to cryosectioning and nucleic acid isolation. (d) The tumor was determined to harbor the IDH1 R132H mutation based on immunostaining with an antibody specific to the mutant protein. (e) TP53 immunostaining demonstrated nuclear expression with an estimated staining index of 20%. All histological images were captured at 400×. Scale bars denote 50 μm. (f) Schematic of serial sectioning strategy and section usage plan. (g) Hierarchical clustering of mutations, using 1—Pearson correlation of amp-seq VAFs over all tumor sections (n = 85) as a distance measure, revealing five clusters. (h–l) VAF patterns comprising cluster 1 (h), cluster 2 (i), cluster 3 (j), cluster 4 (k), and cluster 5 (l). (m) Controlling for gene dosage reveals discordance of IDH1 R132H VAF with respect to truncal ATRX and TP53 mutations, which is explained by a subclonal deletion of chromosome 2q (including IDH1) that occurred after the IDH1 point mutation. (n) Heatmap of the chromosome 2q deletion event frequency (as determined by FACETS [38]), with LOESS fit line (black) and smoothed 95% confidence interval (gray envelope). (o) Clone phylogeny (with arbitrary branch lengths) derived from integrated analysis of SNVs (from amp-seq data) and CNVs (from RNA-seq data). Percentages represent the average abundance of each cellular fraction over all analyzed sections (n = 85). (p) Estimated cellular fractions for all clones and nonmalignant cells over all sections. Black vertical line denotes orthogonal sample rotation.

Figure 5

Gene coexpression modules are highly correlated with clonal abundance (case 2). (a) Hierarchical clustering of gene coexpression modules over all tumor sections (n = 90). (b) Module eigengenes (MEs) illustrate the relative expression levels of genes in each module over all tumor sections. (c) The number of genes that formed each ME. (d–g) Top left: MEs with the strongest correlations to clonal abundance (defined cumulatively). Locally weighted smoothing (LOESS) lines are shown; correlation is based on data points. Bottom left: the 12 genes with the highest correlations to the ME (k_ME). Right: enrichment analysis of gene coexpression modules using published gene sets. FDR-corrected p-values (q-values) from one-sided Fisher’s exact tests are shown. Positive values represent enrichments of genes that were significantly positively correlated to the ME, while negative values represent enrichments of genes that were significantly negatively correlated to the ME. Gene sets representing chromosomal gains or losses include all genes within affected regions (as described in Figure 4o and Table S19). See Table S9 for descriptions and sources of featured gene sets.

Figure 6

Single-nucleus RNA-seq analysis validates inferences from bulk data. (a) Heatmap of p-values (one-sided Wilcoxon rank-sum test) comparing differential expression t-values for genes comprising each bulk coexpression module (colors, x-axis) to all other genes in each SN cluster versus all other clusters. (b) UMAP plot of all nuclei (n = 809) with characterizations of clusters from (a) superimposed. (c) UMAP plot of malignant nuclei (n = 360), with results of Slingshot trajectory analysis [51] superimposed. Malignancy was determined by genotyping all nuclei via single-nucleus amplicon sequencing (snAmp-seq) of cDNA spanning mutations in the truncal clone.

Figure 7

Genotyping nuclei profiled by snRNA-seq reveals the limitations of single-cell CNV-calling algorithms. Heatmap of scaled log₂ expression vectors for the five most upregulated genes in each snRNA-seq cluster vs. all other clusters (one-sided Wilcoxon rank-sum test). Far left: malignancy vector determined by snAmp-seq of cDNA spanning mutations in the truncal clone. Left: malignancy vectors inferred from CNV analysis of snRNA-seq data using the CopyKat [49], InferCNV [25], or CaSpER [50] algorithms (blue = nonmalignant; all other colors = malignant). Right: bar plots depict the total number of unique reads (UMIs) for each nucleus and the average number of UMIs for genes comprising the Gene Ontology category ‘mitotic chromosome condensation’ (GO: 0030261). Red vertical line: max expression of mitotic genes in neurons, which presumably represents background noise.

Figure 8

Correlation to malignant cell abundance predicts single-cell differential expression analysis of malignant vs. nonmalignant cells. (a–d) Analysis schematic. An adult malignant glioma consisting of malignant cells (pink) interspersed with nonmalignant cells (a). (b) Single-cell RNA-seq (scRNA-seq) reveals a hypothetical gene (gene X) that is significantly upregulated in malignant vs. nonmalignant cells. (c) Correlating the same gene’s expression pattern with a binary vector encoding malignant cell abundance (1 = malignant, 0 = nonmalignant) produces identical results. (d) Left: scRNA-seq data from 10 adult human IDH-mutant astrocytomas [24] were randomly sampled and aggregated to create 100 pseudobulk samples. Right (top): Genome-wide differential expression (DE) was analyzed for all sampled cells. Right (bottom): Genome-wide gene coexpression was analyzed for all pseudobulk samples. Each pseudobulk module was summarized by its module eigengene (PC1), which was compared to malignant cell abundance, and the correlation between each gene and each module eigengene (module conformity, or k_ME) was calculated. (e) A pseudobulk malignant cell module featuring the top 15 genes ranked by k_ME. By correlating the module eigengene to pseudobulk tumor purity (f), we see that this module is driven by variation in malignant cell abundance among pseudobulk samples. (g) The extent of DE (t-value) identified by scRNA-seq of malignant vs. nonmalignant cells predicts the correlation between gene expression and malignant cell abundance (pseudobulk k_ME).