SYSMut: decoding the functional significance of rare somatic mutations in cancer

Abstract

Current tailored-therapy efforts in cancer are largely focused on a small number of highly recurrently mutated driver genes but therapeutic targeting of these oncogenes remains challenging. However, the vast number of genes mutated infrequently across cancers has received less attention, in part, due to a lack of understanding of their biological significance. We present SYSMut, an extendable systems biology platform that can robustly infer the biologic consequences of somatic mutations by integrating routine multiomics profiles in primary tumors. We establish SYSMut's improved performance vis-à-vis state-of-the-art driver gene identification methodologies by recapitulating the functional impact of known driver genes, while additionally identifying novel functionally impactful mutated genes across 29 cancers. Subsequent application of SYSMut on low-frequency gene mutations in head and neck squamous cell (HNSC) cancers, followed by molecular and pharmacogenetic validation, revealed the lipidogenic network as a novel therapeutic vulnerability in aggressive HNSC cancers. SYSMut is thus a robust scalable framework that enables the discovery of new targetable avenues in cancer.

Keywords: drug sensitivity; functional genomics; glucocorticoid receptor; lipid metabolism; multiomics integrative analysis.

Figure 1

Fraction of known drivers among recurrently Mut genes, and overview of SYSMut algorithmic framework. (A) Somatic mutations across 30 types of solid tumors were obtained from the Pan-Cancer Atlas dataset, comprising ~10 K tumor samples. Stacked bar plot showing the fraction of the previously known driver genes to the unknown Mut genes, where the Y-axis indicating the cancer types, and the X-axis denotes the number of Mut genes in each cancer. Cancer types are sorted from top to bottom based on total number of Mut genes. Recurrent Mut genes (mutation rates > %0.5 on the right) depicted a higher fraction of known driver genes, while only a very small fraction of rarely Mut genes (mutation rates <0.5% on the left) are known yet. (B) Module 1: For each Mut Gene (G₀), this module identifies its transcriptional TG by traversing a predefined transcriptional influence network. Additionally, this module prepares multiomics profiles of tumor samples for downstream analysis, including mutation status of the GI, gene expression levels of all genes in the network, along with DNAMeth and sCNA of all target genes. (C) Module 2: This module employs a hierarchical Bayesian regression framework that integrates multiomics tumor profiles to estimate the relative influence of each GI on the expression levels of a specific transcriptional TG in both Mut and non-Mut tumor samples, while also accounting for cis-regulatory influences of other confounder genes (PG), and trans-regulatory influences on the TG expression. (D) Module 3: This module uses the output of the Bayesian regression framework to determine the statistical significance of differences in the regulatory influence of the GI on the expression levels of a specific TG in the Mut and non-Mut contexts, thus enabling quantification of the transcriptional impact of mutations in GI on the expression levels of TG. (E) Module 4: Predicts probabilities of the mutational impact of each GI by summarizing the impact of GI across all downstream TGs.

Figure 2

Computational benchmarking of SYSMut using simulation and PanCancer datasets. (A–C) Simulation-based benchmarking of the performance of SYSMut in detecting low-, medium- and high-mutational impacts of the GI across conditions, including (a) Number of samples available for training the model; (b) mutation rates and (c) number of trans-regulatory confounders. Note the significant performance degradation in mutational impact detection when the trans-regulatory confounders are not considered in the model (dashed line). (D) Example of SYSMut-based quantification of the functional impact of Mut genes in HNSC. The scatter plot denotes the fraction of impacted target genes associated with each indicated GI (Y-axis), plotted against the GI’s mutation rate (X-axis). The horizontal dashed line denotes the threshold of statistical significance (P < 0.05). Shown in purple are Mut genes with significant functional impact identified by SYSMut, with red nodes denoting previously identified driver genes in HNSC. (E) The spider plot indicates the percentage of previously identified driver genes within each of the 29 cancer types whose impact was recapitulated by SYSMut, Hit’nDrive, Xseq and Moonlight. Cancers where SYSMut detected impact in <75% of the previously identified driver genes were due to either low sample size () or very low number of driver genes (). (F) Bar graph indicates the mean ± SD of the fraction of PanSoftware driver genes detected by each of the algorithms (SYSMut, Hit’nDrive, Xseq and Moonlight) across cancer types. (G) Bar graph indicates the mean ± SD of the number of genes identified as significant by each of the algorithms across cancer types. (H) Bar graphs (left) denote the percentage of Mut genes with significant transcriptional impact identified by SYSMut across cancers. Also shown in right, are the distributions of mutation rates of genes whose mutations were identified by SYSMut as exhibiting significant downstream impact.

Figure 3

SYSMut-based discovery of gene subnetworks harboring mutations with significant functional impact across cancers. (A) Violin plots denote the joint-mutation rate (Y-axis) of subnetworks of genes identified by SYSMut as harboring mutations with significant downstream transcriptional impact across cancer types. Bar graphs in dark red (top) denote the number of cancer types where SYSMut identified mutations in the subnetwork as exhibiting significant impact on downstream transcriptional programs. Note that SYSMut identified the lipid metabolism subnetwork as exhibiting significant mutational impact across the most number of cancer types. (B) Lipid metabolism gene subnetwork detailing the convergent transcriptional impact of mutations in GIs (green hexagons) across cancer types in the PanCancer dataset. The size of the TG nodes (circles) corresponds to the number of cancer types where the TG was deemed as being impacted by mutation in the GIs. Red TG nodes (circles) are transcriptional targets uniquely impacted by GIs within this subnetwork, whereas pink TG nodes are also impacted by other previously identified cancer driver genes. Directionality of the connection is indicated by arrows and the thickness of the arrows indicates the number of cancers where SYSMut identified a mutational impact of GI on the TG. The color of the arrow corresponds to the number of cancer types where SYSMut detected an impact of mutations in GI on the corresponding TG (black, corresponds to ≥4 cancer types, whereas light gray corresponds to <4 cancer types). (C) Lolli plot showing somatic mutations in the PanCancer dataset across the length of the amino-acid sequences for each of the GIs in the lipid metabolism network. Note the broad distribution of mutations throughout the length of the amino-acid sequences for each of the GIs in this network, a pattern suggestive of loss of function as opposed to activating mutations. (D) Horizontal bar plots indicating the impact of mutations on a transcriptional readout of lipid metabolism in the Cancer Cell Line Encyclopedia. Gray bars (left) indicate average lipid metabolism index across all cell lines in the cancer type or lineage. Bars on the right detail increased (blue) or decreased (red) lipid metabolism in cell lines harboring mutations in the lipid metabolism network as compared to non-Mut cell lines within each lineage. Stars indicate lineages where SYSMut detected significant transcriptional impact of mutations in the lipid metabolism network in cancers exhibiting at least 5% subnetwork mutation rate within the PanCancer dataset. ESCC, esophageal squamous cell carcinoma; LUAD Ade, NSCLC adenocarcinoma; LUAD Seq, NSCLC squamous; CSCC, cervical squamous cell carcinoma; BLCA, urothelial bladder carcinoma.

Figure 4

Impact of mutations in the lipid metabolism network in HNSC. (A and B) DSS in HNSC patients (A) or HPV-negative HNSC patients (B). Shown are Kaplan–Meier estimates of DSS for patients with tumors where one or more GIs in the lipid metabolism subnetwork are Mut (orange), as compared to their non-Mut (blue) counterparts. The statistical significance of differences in survival rates between Mut and non-Mut categories was determined using the LogRank test (LogRank P); also shown is the statistical significance estimated using the Wald test in a multivariable Cox regression model incorporating cancer stage and age (Adj. Cox P). (C) Lipid metabolism network detailing the convergent transcriptional impact of Mut GIs (green hexagon) in HNSC as identified by SYSMut. Red TG nodes are transcriptional targets uniquely impacted by GIs (green hexagons) in this subnetwork, whereas pink TG nodes are also impacted by other previously identified driver genes in HNSC. (D–G) SYSMut-based evaluation of the transcriptional impact of siRNA-based knockdown of MED1 and NCOA6 on their respective downstream target genes in distinct HNSC cell lines (BICR22 and CAL26). Shown are the distributions (magenta) of Bhattacharyya distances representing SYSMut’s assessments of the differential regulation of the downstream transcriptional targets of MED1/NCOA6 in HNSC cell lines treated with either siRNA targeting MED1 (siMED1) or NCOA6 (siNCOA6) as compared to nontargeting control (siSCRAM). Also shown, as an internal control, are the distributions of Bhattacharyya distances (light blue) representing SYSMut’s assessment of impact of random untargeted GIs on their respective TGs. (H and I) Depict the impact of MED1 and NCOA6 knockdown on target genes previously identified by SYSMut as being impacted by MED1/NCOA6 mutations in HNSC primary tumors (as shown in panel C). Individual stems denote the Bhattacharyya distance (Y-axis) corresponding to SYSMut’s assessment of the differential regulation of individual target genes (X-axis) upon MED1/NCOA6 knockdown in BICR22 (dark cyan) and CAL27 (yellow) HNSC cell lines. The box plots depict the distribution of Bhattacharya distances of randomly selected GIs on their respective TGs, thus serving as an internal control within each cell line model. TGs that are not expressed in a given cell line are denoted as NE.

Figure 5

Differential drug sensitivity of HNSC cell lines harboring mutations in the lipid metabolism network as compared to their non-Mut counterparts. (A) Box plots detailing lipid metabolic activity in HNSC cell lines harboring mutations in the lipid metabolism network (Mut, N = 9) as compared to their non-Mut, (N = 36) counterparts. Statistical significance of differences in lipid metabolic activity was estimated using a Wilcoxon signed rank test. (B) Volcano plot depicting significant differential sensitivity between Mut and non-Mut HNSC cell lines across drugs in the PRISM primary screen. X-axis denotes the difference (log2fold-change) in average cell viability between Mut and non-Mut HNSC cell lines upon drug treatment, whereas the Y-axis denotes the statistical significance (−log10FDR) of the difference assessed using a Wilcoxon test. Points above the horizontal red-dashed line denote drugs inducing a significant shift (log10FDR < 0.01) in average cell viability between Mut and non-Mut HNSC cell lines. Drugs were additionally denoted as exhibiting marked (triangles) shift in sensitivity if the Mut cells exhibited, on average, lower cell viability on drug treatment as compared to DMSO control (indicated by negative log2 cell viability in the PRISM primary screen) while the non-Mut cells were nonresponsive (indicated by positive log2 cell viability values in the PRISM primary screen). GRA are highlighted with dark green stars. (C) Average clobetasol-proprionate dose–response curves for Mut and non-Mut (blue) HNSC cell lines derived from the PRISM secondary screen dataset. Y-axis denotes the difference in cell viability (Log2fold change) upon treatment with different doses (X-axis) of clobetasol-proprionate as compared to DMSO control. (D) Box plots detailing the comparison of area under the curve (AUC) estimates for clobetasol−propionate dose–response curves in Mut versus non-Mut HNSC cell lines obtained from the PRISM secondary screen dataset. Statistical significance was assessed using a Wilcoxon test. (E) IncuCyte-based cell growth assessments in Mut (SCC9, BICR18) and non-Mut (BICR22, CAL27) HNSC cells treated with clobetasol−propionate (2.5 μM) as compared to DMSO control treatment. Y-axis represents the average fold change in cell confluence values at the final time point (120 h for SCC9; 96 h for BICR22, CAL27 and BICR18) normalized to DMSO control, plotted as mean ± SD, obtained from at least three replicate experiments. The statistical significance of differences in normalized cell confluency between the respective test versus control groups were estimated using a Student t-test assuming unequal variances.