Integrated Proteogenomic Characterization across Major Histological Types of Pediatric Brain Cancer

Abstract

We report a comprehensive proteogenomics analysis, including whole-genome sequencing, RNA sequencing, and proteomics and phosphoproteomics profiling, of 218 tumors across 7 histological types of childhood brain cancer: low-grade glioma (n = 93), ependymoma (32), high-grade glioma (25), medulloblastoma (22), ganglioglioma (18), craniopharyngioma (16), and atypical teratoid rhabdoid tumor (12). Proteomics data identify common biological themes that span histological boundaries, suggesting that treatments used for one histological type may be applied effectively to other tumors sharing similar proteomics features. Immune landscape characterization reveals diverse tumor microenvironments across and within diagnoses. Proteomics data further reveal functional effects of somatic mutations and copy number variations (CNVs) not evident in transcriptomics data. Kinase-substrate association and co-expression network analysis identify important biological mechanisms of tumorigenesis. This is the first large-scale proteogenomics analysis across traditional histological boundaries to uncover foundational pediatric brain tumor biology and inform rational treatment selection.

Keywords: BRAF alteration; CPTAC; CTNNB1 mutation; kinase activity score; kinase substrate regulation; pediatric brain tumor; post-translational modification; proteomic cluster; recurrent versus primary tumors; tumor microenvironment.

Figure 1.. Proteomic based clustering of pediatric brain tumors.

A. Summary of pediatric brain tumor cohort. B. Presence of omics data sets for each of the 218 tumor samples. For each sample, the clinical status at sample collection (i.e., post- mortem, post-treatment or treatment naive) is also reported. C. Kaplan Meier curves for OS of patients stratified by proteomic cluster. D. Proteomic clusters and differentially expressed proteins allocated to 14 gene clusters (top heatmap). Each row represents a proteomic cluster, while each column represents a protein. Red/blue colors denote up/down regulation patterns of different proteins in a cluster. Distributions of diagnoses, clinical outcomes, and mutation status among the 8 clusters (top left pie plots), and gene members of key pathways enriched in each gene group (bottom heatmap) are shown. For each pathway, the averaged ssGSEA score in each protein cluster based on global proteomic (Protein) and RNA-seq (RNA-seq) data are illustrated to the right. E. Heatmap of kinase activity scores for the CP tumors (n=16). Silhouette scores (top) measures the cohesiveness of tumors classified as C4 and C8 based on kinase activity score. Kinases involved in AKT1 or ERK1/2 signaling are highlighted in the heatmap. F. Diagram illustrating differences between C4 and C8 CP tumors in terms of phosphorylation abundance and kinase activity for AKT and ERK1/2 signaling members. G. MRM measurements validated different activities of proteins and phosphoproteins between C4 and C8 CPs. The numbers annotated under each pair of boxplots correspond to AUC (area under the curve) for classifying the two groups of CP using the corresponding protein/phosphosite measurement.

Figure 2.. Immune infiltration in pediatric brain tumors.

A. Heatmap illustrating cell type compositions, and activities of selected individual gene/proteins and pathways across 5 immune clusters. The heatmap in the first section illustrates the immune/stromal signatures from xCell. The heatmap in the second section illustrates signatures of microglia, neurons and oligodendrocytes derived from single cell sequencing data from Darmanis et al (2017). RNA and protein abundance of key immune-related markers, and ssGSEA scores based on global proteomic data for biological pathways upregulated in different immune groups are illustrated in the remaining sections. B. Contour plot of two-dimensional density based on Macrophage (y-axis) and Microglia scores (x-axis) for different immune clusters. For each immune cluster, key upregulated pathways significant at 10% FDR are reported based on RNAseq (R), global proteomic (P) and phospho-proteomic data (Ph) in the annotation boxes. For Cold-mixed and Cold-medullo clusters, pathways upregulated in both clusters are reported. C. Distribution of pathway scores of Signaling by WNT and Oxidative Phosphorylation based on global proteomic data and RNA stratified by immune clusters. D. Heatmap showing the comparison between immune clusters (columns) with proteomic clusters and different histologies (rows). Each row sums to one, with different entries showing the proportion of tumors allocated to different immune clusters. E. xCell immune/stromal and antigen presentation signatures in BRAF^V600E or BRAF^Fusion compared to BRAF^WT in LGG. F. Distribution of RNA levels of HLA-A, HLA-B and HLA-C in LGG tumors with different BRAF statuses. G. Distribution of macrophage and microglia polarization (M2-M1) in LGG tumors with different BRAF statuses.

Figure 3.. Impact of genomic alterations on transcriptomic, proteomic and phosphoproteomic abundances.

A. Distribution of protein abundance of BRAF, CTNNB1, and NF1 across tumor samples stratified by different mutation status and diagnoses. Symbols *, **, and *** correspond to p-values less than 0.1, 0.01 and 0.001, respectively. B. DNA copy number amplification/deletion frequencies along chromosome 1 among EP, HGG and MB samples. Genes with detected CNV-RNA/protein or CNV-RNA/protein/phospho cascade events are labelled as vertical bars in the top track. C. Distribution of DNA copy number (log ratio), RNA and protein abundance of RABGAP1L, RAB3GAP2 and FDPS stratified by their amplification statuses in EP, MB and HGG tumors. For RABGAP1L, Symbols *, ** and *** mean the same as in A. “ns” stands for “not significant” (p-value>0.1). D. Illustration of the impact of CTNNB1 mutation on RNA and protein abundance in CP samples. x-axis (y-axis) represents signed -log10 FDR for testing the association between protein abundances (RNAs) and CTNNB1 mutation. Cell-Cell Contact Zone (Coagulation) pathway is enriched in the set of proteins up (down) regulated in CTNNB1 mutant samples. A few members of the WNT Signaling pathway whose protein or phosphosites are associated with CTNNB1 mutation are highlighted in red. Phosphosites are annotated with “(P)” in their gene symbols. E. Distribution of protein and phosphosite abundances among CTNNB1 mutant and CTNNB1 wild-type CP tumors for known key members of the WNT Signaling pathway interacting with β-Catenin and transcription factors regulated by CTNNB1. Symbols *, ** and *** correspond to FDR less than 0.1, 0.01 and 0.001, respectively. “ns” stands for “not significant” (FDR >0.1). F. Illustration of the regulatory role of β-Catenin.

Figure 4.. Phosphoproteomic analysis of kinase activity

A. Heatmaps showing the global abundance (right panel) and the kinase activity score (left panel) of selected kinases across different histologies. For each kinase, the Pearson’s correlation between its global abundance and kinase activity within each histology is shown in the middle panel. B. Scatterplot showing the global abundance of a particular kinase (x-axis) versus the phospho-abundance of the targeted substrates (y-axis). First row is based on the data from the discovery cohort; while the second row displays the data based on the validation cohort. C. Heatmap showing global proteomic abundance of CDK1, CDK2 and CAMK2A as well as phosphorylation abundance of MCM2 Ser 139, GJA1 Ser 325, GJA1 Ser 314, SYN1 Ser 568 and SYN1 Ser 605 among HGG in the discovery and validation cohorts. D. Diagram showing kinase-substrate associations involved in CNS development in LGG (top-middle panel). Scatter plots showing the association between the global (or phospho) abundance of each kinase (x-axis) and the phospho-abundance of the corresponding substrate (y-axis).

Figure 5.. Insights from proteogenomic analysis of LGG

A. Heatmap illustrating ssGSEA scores of selected pathways differentially expressed between LGG tumors with different BRAF statuses based on global proteomic data. Dot-plot on the left-side summarizes ssGSEA pathway scores based on RNA data among samples with different BRAF statuses. B. Distributions of RNA, TMT protein abundance (TMT Global), and MRM protein abundance (MRM Global) of AKT1, AKT2, and AKT1S1 in samples with different BRAF alteration statuses. FDR levels of two-sample comparisons between BRAF^V600E/ BRAF^Fusion and BRAF^WT are annotated. C. The network topology representing the LGG phosphosite co-expression network module enriched in sites upregulated in BRAF^v600E compared to BRAF^WT tumors. Phosphosites mapping to genes in the HNRNP family or contained in the MYC Targets pathway are highlighted in red and blue, respectively. D. Scatterplot displaying the association between each phosphosite’s abundance with the global abundance of AKT2 (y-axis) versus the association with BRAF^V600E (x-axis). Phosphosites contained in the network module in C are highlighted in red. Boxplots illustrate the distribution of the activity scores (ssGSEA) of the network module in C based on phosphoproteomic data in samples with different BRAF status. Pie-plot shows the proportion of phosphosites contained in the network module in C whose abundances are associated at 5% FDR with the global abundance of AKT2.

Figure 6.. Insights from proteogenomic analysis of HGG

A. Scatterplot showing OS of HGG patients versus the global protein abundance of IDH1 and IDH2 in the tumors. B. Heatmap of global abundance of IDH proteins in the discovery cohort. C, D. 95% CI of hazard ratio coefficients from Cox-regression for IDH1/2 scores and other covariates based on the discovery cohort (panel C) and Data Set 2 (D). E, F. Kaplan-Meier curves of overall survival for HGG H3^Mut samples (grey), H3^WT samples with low IDH1/2 abundance (red) and H3^WT tumors with high IDH1/2 abundance (blue) for the discovery cohort (panel E) and the validation cohort (F).G. Illustration of drug target analysis result. The bottom-left heatmap illustrates the targeting genes (rows) of each detected drugs (columns). For each gene, the z-score comparing its RNA and proteomic abundances between HGG and LGG is shown in the bottom-right heatmap. Mechanism of actions are annotated on the top of the heatmap together with the resulting score from the cMAP analysis. H. Distribution of kinase activity scores of CDK1, CDK2 and MAPK1 among HGG and LGG tumors, with the latter further stratified by BRAF status.

Figure 7.. Comparison between Initial and Recurrent Tumors

A. (a) Clinical properties and genomic characterization of 18 pairs of IN vs RP tumors. The bar plot illustrates the number of non-synonymous mutations in IN and RP tumors with the number of shared mutations being represented by the shaded area. The potential driver mutation track shows the allele frequencies of somatic mutations of known oncogenes and tumor suppressor genes. Chromosome arm aberrations of each sample and the change of tumor grade from IN to RP of each patient are also shown. (b) Differences in ssGSEA score between RP and IN tumors of key molecular pathways associated with different proteomic clusters. The annotation at the bottom indicates the diagnosis and clinical event IDs of the paired samples for each patient. For example, “Epen.496.3319” refers to a pair of EP tumors with IDs: 7316–496 and 7316–3319. B. Distribution of Spearman’s correlation between the proteomic abundance of any pair of tumors within a particular histology. Correlations between the 18 paired IN-RP samples were further labeled in the violin plots. C. Distribution of kinase activity scores of MAPK1/3 among all LGG samples, LGG samples allocated to C4 and LGG samples allocated to C8. IN and RP samples of patients LGG.350.944 and LGG.173.2154 are highlighted.