Copy number and transcriptome analyses reveal conserved driver mutations in JAK-STAT pathway genes
We interrogated genomic and transcriptomic landscape of 133 JAK-STAT pathway genes in 18,484 patients across 21 cancer types (Additional file 1). To investigate the effects of JAK-STAT pathway genomic alterations on transcriptional output, we first analyzed copy number variation (CNV) of all 133 genes. CNVs were classified into four categories: low-level amplifications, deep amplifications, low-level deletions (heterozygous deletions) and deep deletions (homozygous deletions). Lung squamous cell carcinoma (LUSC) and papillary renal cell carcinoma (KIRP) had the highest and lowest fraction of samples with deleted JAK-STAT pathway genes respectively (Additional file 2). In terms of gene amplification, this was also the highest in lung squamous cell carcinoma (LUSC) and the lowest in pancreatic adenocarcinoma (PAAD) (Additional file 2). To identify pan-cancer CNV events, we focused on genes that were deleted or amplified in at least one-third of cancer types (≥ 7 cancers). We identified 71 and 49 genes that were recurrently deleted and amplified respectively in at least 20% of samples within each cancer type and at least 7 cancer types (Additional file 2). Esophageal carcinoma (ESCA) had 70 genes that were recurrently deleted while only four recurrently deleted genes were found in papillary renal cell carcinoma (KIRP). When considering recurrent gene amplifications, lung adenocarcinoma (LUAD) had the highest number of gains (48 genes), while the lowest number of gene gains was observed in glioma (GBMLGG) (3 genes) (Additional file 2). CNV events associated with transcriptional changes could represent candidate driver genes. Loss-of-function genes can be identified from genes that were recurrently deleted and downregulated at the transcript level. Similarly, genes that were concomitantly gained and upregulated could represent a gain-of-function. Differential expression profiles (tumor vs. non-tumor) were intersected with CNV profiles and we identified 40 driver genes representing potential loss- or gain-of-function (Fig. 1). In at least 7 cancer types, 18 genes were recurrently deleted and downregulated (log2 fold-change < − 0.5, P < 0.01), while a non-overlapping set of 22 genes were recurrently amplified and upregulated (log2 fold-change > 0.5, P < 0.01) (Fig. 1).
JAK-STAT driver genes predict overall survival in diverse cancer types
Focusing on the 40 driver genes identified above, we next analyzed whether expression levels of individual genes were associated with overall survival outcomes. Cox proportional hazards regression analyses demonstrated that all 40 genes harbored prognostic information in at least one cancer type. IL11, PTPN6 and CISH were significantly associated with survival outcomes in patients from 9 cancer cohorts (Fig. 2a). In contrast, IL19, CNTFR and JAK2 were some of the least prognostic genes (Fig. 2a). When deciphering the contribution of individual genes across cancer types, we observed that the glioma (GBMLGG) cohort had the highest number of prognostic genes (31/40 genes) (Fig. 2a). On the other hand, none of the 40 driver candidates were prognostic in esophageal carcinoma (ESCA) and cholangiocarcinoma (CHOL), which suggests the minimal contribution of JAK-STAT signaling in driving tumor progression in these cancer types (Fig. 2a). To identify a core set of prognostic genes denoting pan-cancer significance, we performed Spearman’s correlation analyses on hazard ratio (HR) values obtained from Cox regression and identified 28 highly-correlated genes: IL7, IFNG, MPL, IL11, IL2RA, IL21R, OSMR, IL20RA, IFNGR1, CDKN1A, CISH, SOCS1, IL10, IL10RA, STAT2, IL24, IL23A, PIAS3, IFNLR1, EPO, TSLP, BCL2, IL20RB, IL11RA, PTPN6, IL13, IL17D and IL15RA (Fig. 2b). These genes were collectively regarded as a pan-cancer JAK-STAT signature. To determine the extent of JAK-STAT pathway variation across cancer types, we calculated an activity score based on the mean expression of the 28 genes. When cancers were sorted according to their pathway activity scores, chromophobe renal cell cancer (KICH) had the lowest median score while the highest median score was observed in head and neck cancer (HNSC) (Fig. 2c). Hierarchical clustering of the 28 driver genes demonstrated that they exhibited a wide range of expression depending on the cellular context where they could serve as potential candidates for therapy. For example, IL7, IL15RA, IL21R, IL10, OSMR, IFNGR1, IL10RA, IL2RA, IFNG, IL24, SOCS1, IL20RA, IL11 and IL23A were highly expressed in gastrointestinal cancers (PAAD, STAD and STES) (Fig. 2d). When the 28-gene scores were employed for patient stratification, we observed that the JAK-STAT signature conferred prognostic information in five diverse cancer cohorts (Fig. 2e). Intriguingly, the significance of the signature in predicting overall survival was cancer type-dependent. Kaplan-Meier analyses and log-rank tests revealed that patients with high scores (4th quartile) had higher death risks in glioma (P < 0.0001), pan-kidney (consisting of chromophobe renal cell, clear cell renal cell and papillary renal cell cancers; P < 0.0001) and clear cell renal cell (P < 0.0001) cohorts (Fig. 2e). In contrast, high expression of signature genes was linked to improved survival rates in lung (P = 0.025) and endometrial (P = 0.032) cancers (Fig. 2E). These results were independently corroborated using Cox regression analyses: glioma (HR = 6.832, P < 0.0001), pan-kidney (HR = 3.335, P < 0.0001), clear cell renal cell (HR = 4.292, P < 0.0001), lung (HR = 0.624, P = 0.028) and endometrium (HR = 0.504, P = 0.027) (Additional file 3). Since high expression of signature genes was associated with poor survival outcomes in brain and renal cancers, as expected, we observed a significant increase in expression scores according to tumor stage (Fig. 2g). An opposite trend was observed in lung and endometrial cancers where more aggressive tumors had lower expression scores (Fig. 2f). Lastly, multidimensional scaling analyses of signature genes in the five cohorts revealed significant differences between tumor and non-tumor samples, implying that dysregulated JAK-STAT signaling may serve as a diagnostic marker for early detection in pre-cancerous lesions (Fig. 2F).
The JAK-STAT 28-gene signature is an independent prognostic factor
Multivariate Cox regression analyses confirmed that the signature was independent of tumor, node and metastasis (TNM) stage and age (where information is available): glioma (HR = 4.316, P < 0.0001), pan-kidney (HR = 2.379, P < 0.0001), clear cell renal cell (HR = 2.552, P = 0.00047), lung (HR = 0.636, P = 0.031) and endometrium (HR = 0.434, P = 0.033) (Additional file 3). Given that the signature was an independent predictor of overall survival, we reasoned that its predictive performance could be increased when used in conjunction with TNM staging. Employing the receiver operating characteristic (ROC) analysis, we demonstrated that a combined model uniting the signature and TNM staging could outperform the signature (higher area under the curve [AUC] values) when it was considered alone: pan-kidney (0.838 vs. 0.800), clear cell renal cell (0.836 vs. 0.789), lung (0.724 vs. 0.703) and endometrium (0.760 vs. 0.713) (Fig. 3a). Independently, Kaplan-Meier analyses and log-rank tests confirmed that the signature allowed further delineation of risk groups within similarly-staged tumors: pan-kidney (P < 0.0001), clear cell renal cell (P < 0.0001), lung (P < 0.0001) and endometrium (P < 0.0001) (Fig. 3b).
High JAK-STAT scores were associated with decreased survival rates in glioma patients. We confirmed that this was also true for histological subtypes of glioma: astrocytoma (P = 0.015) and oligoastrocytoma (P = 0.037) (Fig. 3b) but not for patients with glioblastoma multiforme (Additional file 4). Independently confirmed using Cox regression, patients within the 4th survival quartile had lower survival rates: astrocytoma (HR = 2.377, P = 0.018) and oligoastrocytoma (HR = 2.730, P = 0.038) (Additional file 3). In terms of the signature’s predictive performance, ROC analyses revealed that it performed the best in oligoastrocytoma (AUC = 0.951), followed by astrocytoma (AUC = 0.878) and glioma patients when considered as a full cohort (AUC = 0.853) (Fig. 3a).
Consequences of dysregulated JAK-STAT signaling and significant crosstalk with tumor immunity
Since dysregulated JAK-STAT signaling was associated with survival outcomes (Fig. 2 and Fig. 3), we reasoned that patients from diverse cancer types might harbor similar transcriptional defects caused by aberrant activation of JAK-STAT. Differential expression analyses performed between the 4th and 1st quartile patients revealed that a striking number of over 200 differentially expressed genes (DEGs) were shared between all five prognostic cohorts (Fig. 4a; Additional file 5). Significant overlaps were observed in DEGs where 555 genes were found in at least four cohorts, 1009 genes in at least three cohorts and 2034 genes in at least two cohorts (Fig. 4a; Additional file 5). The highest number of DEGs was observed in glioma (2847 genes), followed by pan-kidney (2810 genes), lung (1221 genes), clear cell renal cell (1161 genes) and endometrial (782 genes) cancers (Fig. 4a; Additional file 5). To determine their functional roles, the DEGs were mapped to Gene Ontology (GO) and KEGG databases. All five cohorts exhibited remarkably similar patterns of enriched biological processes (Fig. 4b). Pan-cancer enrichments of ontologies related to inflammation and immune function were observed, e.g., cytokine and chemokine signaling, T cell and B cell receptor signaling, natural killer cell-mediated processes, Toll-like receptor signaling and NOD-like receptor signaling (Fig. 4b). Additionally, genes associated with other oncogenic pathways (Wnt, MAPK, TGF-β, PPAR and VEGF) were frequently altered at both transcriptional and genomic levels (Fig. 4b and c) (Chang & Lai, 2019e; Chang & Lai, 2019f). CNV analyses performed on DEGs that were found to be in common in at least three cohorts demonstrated that transcriptional dysregulation of the aforementioned oncogenic pathways was attributed to activating or inactivating CNVs (Fig. 4c). For instance, except for THBS1 and THBS2, a vast majority of TGF-β DEGs exhibited somatic gains (Fig. 4c). To further explore which upstream transcriptional regulators were involved, we mapped the DEGs to ENCODE and ChEA databases (Fig. 4b). Interestingly, we observed significant enrichment of transcription factors (TFs) involved in modulating immune function: IRF8 (Mattei et al., 2012), RUNX1 (Ono et al., 2007), RELA (Bonizzi & Karin, 2004) and EZH2 (Vasanthakumar et al., 2017) (Fig. 4b). Taken together, our analyses revealed that pan-cancer JAK-STAT drivers underpin numerous aspects of tumor oncogenesis and immunity, which play important roles in tumor progression and ultimately patient prognosis.
IRF8 and JAK-STAT pathway synergistically influence survival outcomes in glioma and renal cancer
IRF8 was among one of the most enriched TFs implicated in the regulation of transcriptional outputs of patients with dysregulated JAK-STAT signaling (Fig. 4b). As a member of the interferon regulatory factor family, IRF8 is needed for the development of immune cells and is often regarded as a tumor suppressor gene since a loss of function is associated with increased metastatic potential (Lee et al., 2008; Yang et al., 2007). Thus, we predict that tumors with low expression of IRF8 would be more aggressive. To evaluate the combined relationship between JAK-STAT signaling and IRF8 expression, patients were categorized into four groups based on median IRF8 and JAK-STAT scores: 1) high-high, 2) low-low, 3) high IRF8 and low 28-gene score and 4) low IRF8 and high 28-gene score. The combined model encompassing JAK-STAT and IRF8 offered an additional resolution in patient stratification: glioma (full-cohort, P < 0.0001), astrocytoma (P = 0.0007), pan-kidney (P < 0.0001) and clear cell renal cell (P < 0.0001) (Fig. 5a). Indeed, patients with low IRF8 and high 28-gene scores performed the worst in cancers where hyperactive JAK-STAT signaling was linked to adverse survival outcomes: glioma (full-cohort: HR = 5.826, P < 0.0001), astrocytoma (HR = 3.424, P = 0.0032), pan-kidney (HR = 5.131, P < 0.0001) and clear cell renal cell (HR = 5.389, P < 0.0001) (Fig. 5b). Our results support a model in which IRF8 influences the behavior of tumors with aberrant JAK-STAT signaling.
Hyperactive JAK-STAT signaling attenuates regulatory T cell-mediated tumor immunity
Given the wide-ranging effects of JAK-STAT signaling on immune-related functions, we hypothesized that JAK-STAT activity would correlate with immune cell infiltration. We retrieved genes implicated in regulatory T cell (Treg) function from four studies and isolated 31 genes that were common in all four Treg signatures (De Simone et al., 2016; Plitas et al., 2016; Tirosh et al., 2016; Zheng et al., 2017). Treg scores were calculated for each patient based on the mean expression levels of the 31 genes. Remarkably, we observed strong positive correlations between Treg and JAK-STAT scores, suggesting that tumor tolerance was enhanced in patients with hyperactive JAK-STAT signaling: glioma (rho = 0.85, P < 0.0001), pan-kidney (rho = 0.88, P < 0.0001) and clear cell renal cell (rho = 0.83, P < 0.0001) (Fig. 6a). As in the previous section, patients were stratified into four categories based on median JAK-STAT and Treg scores for survival analyses. Log-rank tests and Cox regression analyses confirmed that elevated Treg activity further exacerbated disease phenotypes in patients where JAK-STAT scores were already high: glioma (full cohort, HR = 6.183, P < 0.0001), astrocytoma (HR = 3.035, P = 0.00042), pan-kidney (HR = 3.133, P < 0.0001) and clear cell renal cell (HR = 2.982, P < 0.0001) (Fig. 6b and c). Taken together, elevated JAK-STAT signaling may increase the ability of tumors to escape immunosurveillance, resulting in more aggressive tumors and increased mortality rates.