Immunomediated Pan-cancer Regulation Networks are Dominant Fingerprints After Treatment of Cell Lines with Demethylation.

Pan-cancer studies are particularly relevant not only for addressing the complexity of the inherently observed heterogeneity but also for identifying clinically relevant features that may be common to the cancer types. Immune system regulations usually reveal synergistic modulation with other cancer mechanisms and in combination provide insights on possible advances in cancer immunotherapies. Network inference is a powerful approach to decipher pan-cancer systems dynamics. The methodology proposed in this study elucidates the impacts of epigenetic treatment on the drivers of complex pan-cancer regulation circuits involving cell lines of five cancer types. These patterns were observed from differential gene expression measurements following demethylation with 5-azacytidine. Networks were built to establish associations of phenotypes at molecular level with cancer hallmarks through both transcriptional and post-transcriptional regulation mechanisms. The most prominent feature that emerges from our integrative network maps, linking pathway landscapes to disease and drug-target associations, refers primarily to a mosaic of immune-system crosslinked influences. Therefore, characteristics initially evidenced in single cancer maps become motifs well summarized by network cores and fingerprints.

Introduction

New patient stratifications and disease classifications are centered on profiling and re-phenotyping processes that put, at the nexus of cancer studies, the role of in-depth analysis of complex regulatory mechanisms linked to cancer hallmarks. The role of immune system is recognized as central to the identification of neoplastic aberrations and explanation of the pathogenesis of inflammatory conditions, directly or indirectly related to cancer. The effectiveness of cancer therapies also depends on the immune system response.1 In particular, RIP receptor-interacting protein (RIP) kinases play a major role in inflammatory signaling, apoptosis, and inflammation-associated cell death.2 Overall, the established links among dysregulated immune cells, autoimmune disorders, and cancer refer to the immunomodulatory compounds as the most promising treatment options.

First, emerging pan-cancer studies have emphasized first the role of mutations, for instance with 12 different cancers.3 Then, a pan-cancer regulatory network approach was presented with six different cancers.4,5 Epigenetic mechanisms influence immune-related cell functional regulation. Recently, a pan-cancer study proposed prognostic landscapes of genes and infiltrating immune cells.6 Another pan-cancer study proposed a tool to identify methylation-driven genes.7 In particular, epigenetic alterations may precede cancer transformation.8 For instance, DNA methylation can reveal regulatory functions exerting influence on both cancer onset and progression.9 Aberrant DNA methylation of the promoter region induces tumorigenesis suppression through gene inactivation, eg, silencing. Notably, the reversible nature of epigenetic alterations has motivated the development of targeted therapies.10,11 Unlike other agents targeting a single gene product, epi-drugs have chromatin as their target and act through the inhibition of histone deacetylases and DNA methyltransferases. Several epigenetic modifications occurring during cancer pathogenesis are of particular interest for our proposed pan-cancer study and have been reported for primary effusion lymphoma (PEL),12,13 mesothelioma (MES),14 breast cancer (BC),15,16 renal cell cancer (RCC),17 and melanoma (MEL).18,19

We have analyzed five cancer cell lines treated with 5-aza-2′-deoxycytidine (DAC; Dagogen), a potent demethylating agent used to correct epigenetic defects, including the reactivation of tumor suppressor genes silenced by epigenetic mechanisms in tumors.20 The pan-cancer cell lines consist of A375 (for MES, MEL, RCC, and PEL) and MCF7 for BC and A375 for RCC. After computing gene expression profiles and identifying possible markers in all cancers, we explored pan-cancer epigenetic dysregulation signatures through the analysis of pathways, regulatory circuits, and network patterns in search of commonalities and specificities robustly identified across the transcriptome landscapes. Network-driven identification of these epigenetic signatures is an inference strategy revealing synergisms of regulation and control for target genes, as seen already at a smaller scale.21 Modules and communities reflect the idea of genes acting contextually in signaling pathways and regulatory networks. Thus, the detection of significantly perturbed subnetworks can offer insight on the mechanism of action of drugs or treatment effects.22,23

Inference from regulatory networks is notoriously complex when multiple layers, from transcription factors (TFs) to microRNAs, are integrated. However, network fingerprints might improve over gene signatures in three possible ways: (i) by elucidating the relation between bioentities through connectivity, weakly (in correlative mode) or strongly (measuring dependence robustly); (ii) since cancers are highly heterogeneous diseases, single-gene or gene-signature markers can be limited in their capacity to capture mechanisms for which drugs are designed, while marker panels encapsulated in either subnetworks or modules may be relevant to some complex phenotypes; (iii) both synergistic and antagonistic dynamics can almost naturally be accounted by networks. Consequently, integrative network approaches are viable inference tools whose impacts might be clinically relevant and motivate refined translational approaches.

Differentially Expressed Genes

Figure 1 reports both cancer-specific gene sets and differentially expressed gene (DEG) intersections the pan-cancer. Overall, down- and upregulated DEG detections after DAC treatment are distributed as follows: PEL (up = 239; down = 236), RCC (up = 247; down = 215), MEL (up = 161; down = 174), MES (up = 296; down = 270), and BC (up = 265; down = 204). The innermost intersection with 44 genes is the one shared by all five cancer types. A few selected annotations appear in Figure 1. More specific annotations are shown in Table 1 (concordant and discordant DEG regulations), and complete information is reported in the Supplementary Table 1.

Figure 1

Venn diagram of DEGs in the pan-cancer. Selected examples of DEGs located at the central Intersection (44 genes). Concordant down- and upregulation when the same sign occurs in the pan-cancer, otherwise regulation is discordant.

Table 1

Concordant up (labeled) and down and also discordant (u-d) DEG regulations.

GENE_SYMBOL	LOG2 (PEL)	LOG2 (RCC)	LOG2 (MELAN)	LOG2 (MESOT)	LOG2 (BC)
CTAG1 up	5.02	2.80	1.59	3.76	7.29
ERG up	3.35	1.81	2.11	3.40	2.63
XPC up	1.99	1.80	1.57	1.65	2.20
ACTC	−3.49	−2.04	−3.31	−2.97	−3.55
ARCN1	−2.61	−6.23	−3.08	−2.68	−1.93
BAI1	−4.69	−4.42	−7.37	−7.05	−1.60
BRAF	−1.83	−1.63	−1.87	−2.40	−2.66
CASP8	−3.06	−1.56	−2.25	−3.54	−1.92
CBFA2T2	−4.02	−3.71	−3.43	−2.94	−3.35
DBCCR1	−2.06	−1.50	−1.55	−1.62	−1.80
EBI2	−3.98	−2.15	−3.83	−6.06	−3.57
FGF20	−3.72	−3.96	−4.57	−4.80	−4.56
FKSG2	−4.85	−6.08	−7.05	−6.04	−5.20
IFNA5	−1.81	−1.53	−2.08	−1.50	−1.58
IGFBP1	−1.91	−3.21	−2.16	−2.77	−2.57
MGB2	−4.73	−5.25	−5.04	−5.05	−4.56
MMP16	−2.69	−2.98	−1.74	−1.68	−4.61
MYO1A	−2.96	−2.80	−3.53	−4.52	−3.58
NDRG2	−1.78	−1.62	−2.55	−3.04	−2.97
P45	−1.58	−3.67	−3.25	−4.35	−3.53
PTPN13	−2.48	−1.92	−3.16	−3.35	−2.64
RBBP9	−1.97	−2.14	−1.80	−3.37	−1.89
RECK	−3.80	−4.28	−1.95	−2.10	−2.09
RELA	−4.05	−2.70	−1.83	−2.69	−1.97
SERPINB13	−2.37	−2.32	−2.21	−4.37	−3.58
SH3BP2	−4.07	−1.77	−3.64	−3.69	−3.03
SLC22A4	−2.70	−1.75	−2.54	−3.14	−1.91
ST5	−2.18	−2.36	−1.91	−3.14	−2.42
STAG1	−3.11	−2.63	−2.37	−3.71	−3.13
TERF1	−1.59	−1.55	−1.55	−2.31	−2.01
TL132	−2.44	−2.61	−1.95	−2.67	−1.70
TNFRSF18	−3.16	−3.49	−3.77	−2.57	−2.72
VEGF	−6.02	−3.63	−3.16	−3.62	−3.19
WRN	−2.59	−2.21	−2.58	−1.69	−1.77
XRCC3	−3.98	−2.20	−2.08	−3.00	−2.53
CYP1A2 u-d	2.26	−2.14	−1.92	−3.89	−3.97
GML u-d	−3.50	1.78	−3.15	−4.35	−2.16
HHGF u-d	1.67	−2.82	−3.26	−3.82	−2.37
IL26 u-d	2.78	3.52	1.65	−2.35	−2.38
SEL1L u-d	2.35	−1.51	1.82	2.69	2.91
DCC u-d	−1.50	1.55	2.06	3.42	2.37

Among the upregulated genes, ERG encodes a member of the erythroblast transformation-specific (ETS) TF family, known regulators of embryonic development, cell proliferation, differentiation, angiogenesis, inflammation, and apoptosis processes. Among the downregulated genes, VEGF encodes a protein that specifically acts on endothelial cells with various effects, such as mediation of increased vascular permeability, induction of angiogenesis and endothelial cell growth, promotion of cell migration, and inhibition of apoptosis. PTPN13 encodes a member of the protein tyrosine phosphatase (PTP) family, regulating various cellular processes such as cell growth, differentiation, and oncogenic transformation. Interestingly, the pseudogene FKSG2 also appears. BRAF and CASP8, well-known master regulators, appear too, the former playing a role in regulating the MAP kinase/ERK signaling pathway and affecting cell division, differentiation, and secretion, and the latter being central to the sequential activation of caspase cascades that determine the execution phase of cell apoptosis. Then, FIGF encodes a protein member of the platelet-derived growth factor/vascular endothelial growth factor (PDGF/VEGF) family, active in (lymph-) angiogenesis and endothelial cell growth. The group of discordant genes includes, in particular, GML (glycosylphosphati-dylinositol-anchored molecule like), involved in the apoptotic pathway or cell-cycle regulation induced by TP53 after DNA damage, and CYP1A2, which regulates cell growth, cell motility, and morphogenesis by activating a tyrosine kinase signaling cascade after binding to the proto-oncogenic c-Met receptor. Owing to stimulation of mitogenesis, cell motility, and matrix invasion, the latter gene plays a central role in angiogenesis, tumorigenesis, and tissue regeneration.

Of relevance for disease-gene association studies, we recall the importance of RIP kinases in inflammatory signaling, apoptosis, and inflammation-associated cell death, as we have found RIP1 upregulated in MEL and RIP3 upregulated in MES. IL12 (downregulated in MEL) is an early proinflammatory cytokine that links adaptive and innate immune system,24 inducing anticancer immunity at the tumor microenvironment level. Then, (a) STAT3, which is involved in Crohn’s disease, was found downregulated in RCC; (b) PTPRC (gene encoding surface receptor subunits) was found downregulated in BC and RCC; (c) MDM2 was found downregulated in MEL and upregulated in BC, encoding a protein (nuclear-localized E3 ubiquitin ligase) that promotes tumor formation by targeting tumor suppressors (for modification and degradation) such as with TP53, which in turn transcriptionally regulates the gene; (d) TP73, which cooperates with TP53 to induce apoptosis, was found upregulated in BC and MEL; and (e) TNFAIP3 (downregulated in RCC) encodes a cytoplasmic zinc finger protein inhibiting NFκB activity, likely player in Rheumatoid Arthritis (RA).25 The NFκB factors are known activators of inflammatory response (proinflammatory),26 and among them RELA and NFKB1 were found, respectively, downregulated and upregulated in BC and MES. Another specific activator of inflammatory genes, RUNX, was found upregulated in BC. Instead, interleukin-6 (IL-6) was found downregulated in RCC. This is known as a cytokine that regulates malignant transformation or cancer progression and can trigger cancer cell proliferation, survival, and invasion while suppressing host antitumor immunity.

Pathways

Figure 2A displays by Cytoscape the pathway landscape of the pan-cancer, with a selection operated by significance level. In association with these pathways, we also build disease, drug and druggable maps (respectively, Fig. 2B–D). First, we highlighted the most significant pan-cancer common pathways, those with a q-value ≤0.001 (Supplementary Fig. 1 shows the map of other filtered pathways, ie, the most significant pathways with q-value ≤0.01). In addition, distinctly annotated pathways appear for single cancers. In Table 2 (extracted from Supplementary Table 2), a list of pathways restricted to the top-3 terms best enriched by up- and downregulated genes from each tumor histotype. Analyses were performed by using the ORA tool in ConsensusPathDB.27 Pathways were listed by associated dysregulated gene sets, and statistical significance is represented by the hypergeometric P-value (≤0.01) with corrections for multiple testing. We detected cytokines, T-cell, and Toll-like receptors, NFκB, TGF-β, and JAK-STAT signaling, among the best enriched terms. Such intercellular regulators are crucial for cells engaged in the innate and adaptive immune and inflammatory host defense system. These in particular inhibit tumor development and progression in response to immune and inflammation related conditions and suggest a role for improving cancer immunotherapies. Among cancer-specific annotations, the main emerging annotations are: (a) RCC: upregulation of PI3k-Akt, VEGF, ErbB, Atm, Hippo, platelet, extracellular matrix (ECM) organization, DNA repair, and downregulation of MAPK signaling; (b) MEL: upregulation of Hippo and Wnt and ECM organization, and downregulation of RIG-I like receptor, Caspase 8 activation, TP53, TRIF-mediated programmed cell death; (c) MES: upregulation of NGF, RAS, MAPK, ERK signaling, and DNA repair, and downregulation of PI3k-Akt and MAPK; (d) BC: upregulation of Atm, TP53, Wnt, Hippo signaling, and apoptosis, DNA repair, senescence and autophagy, and downregulation if PI3K-Akt and TP53 signaling; (e) PEL: upregulation of ErbB, insulin signaling, and Vegf hypoxia and angiogenesis, and downregulation of DNA repair and damage response, MAPK, EGF signaling, senescence and autophagy, TSH, and B-cell receptor.

Figure 2

(A) Pathway landscape. The top enriched pathways (with the q-value ≤0.001) are reported for the pan-cancer. Node color indicates enrichment frequency, ie, shared pathways among down- or upregulated genes. Therefore, blue indicated two cancers sharing, red indicates three cancers sharing, and gray indicates more than three cancers sharing. Edge color represents pathways enriched among downregulated genes (green) and upregulated genes (red). The edge size indicates the DEG number enriching the pathway. (B) Pan-cancer KEGG-disease map. Node color indicates enrichment frequency, ie, shared diseases among down- or upregulated genes in the pan-cancer. Blue: shared between two cancers; Red: shared between three cancers; Gray: shared between more than three cancers. Edge color indicates regulation sign, up (red) and down (green). Black links indicate genes connecting directly to disease pathways. (C) Pan-cancer gene-drug map. Red or green edges stand for up- or downregulation in our pan-cancer, respectively. Black links target genes and drugs that refer to the source (http://www.mycancergenome.org/) integrated in DGIdb. Disease boxes show color according to previous annotations. (D) Pan-cancer druggable gene-disease map. Search categories available in DGIdb allow to find druggable DEGs that are associated with diseases, and without having a direct known drug (removing all genes appearing in the drug map). Solid edges stand for curated gene category interactions, and dashed edges stand for predicted interactions. Disease boxes show color according to previous annotations.

Table 2A

Top-3 pan-cancer pathway terms: BC and MEL.

PATHWAY NAME	SET SIZE	MEL DEG DOWN-REGULATED	P-VALUE	Q-VALUE	PATHWAY SOURCE
Cytokine-cytokine receptor interaction	267	CSF2RB; IL2RA; IL9R; TNFSF18; IFNA5; TNFSF15; IFNA1; IL13; IL4; IL12A; IFNB1; TNFRSF10B; BMPR1B; IFNA8; TNFRSF9; TNFRSF18; LIFR	8.31E-10	1.98E-07	KEGG
Jak-STAT signaling	158	CSF2RB; IL2RA; IL9R; IFNA8; IFNA1; IL13; IL4; IFNA5; IFNB1; AKT2; IL12A; LIFR	4.51E-08	5.37E-06	KEGG
Toll-like receptor signaling	108	IFNA8; IL12A; IFNA1; CD14; IFNA5; IFNB1; AKT2; CASP8; RELA	9.89E-07	4.71E-05	KEGG
		MEL DEG Up-regulated
DNA Repair	114	MUTYH; POLD1; XPC; TDG; FANCG; FANCF; FANCA; CDK7; DDB2; BRIP1; NTHL1	3.32E-08	9.28E-06	Reactome
Extracellular matrix organization	152	TGFB1; MMP11; TIMP1; CTSL2; COL1A2; COL1A1; COL3A1; ITGB4	2.21E-04	5.13E-03	Reactome
Hippo signaling	156	TGFB1; FZD4; TCF7L2; DVL3; LATS1; TP73; TEAD4	1.41E-03	1.55E-02	KEGG
		BC DEG Down-regulated
Cytokine-cytokine receptor interaction	267	IL2RA; TNFSF10; ACVR1B; IFNA1; IFNA5; BMPR1B; IFNA8; IL2; TNFRSF9; IL26; TNFSF8; TNFRSF18	8.98E-05	7.94E-03	KEGG
PI3K-Akt signaling	346	IL2RA; IFNA8; IFNA1; THBS2; RXRA; IFNA5; GNG7; AKT2; IL2; FGF12; RELA; TSC2; FGF20	2.69E-04	9.74E-03	KEGG
Toll-like receptor signaling	108	IFNA8; IFNA1; CD14; IFNA5; AKT2; CASP8; RELA	3.00E-04	9.74E-03	KEGG
		BC DEG Up-regulated
atm signaling	18	BRCA1; RBBP8; ATM; JUN; TP73; MDM2	3.13E-07	1.25E-004	BioCarta
DNA Repair	114	BRCA1; XPC; LIG3; ATM; TDG; ERCC6; FANCE; FANCC; RAD50	1.19E-004	5.70E-003	Reactome
p53 signaling	68	PPM1D; ATM; TP73; CDK6; MDM2; SERPINE1; APAF1	1.32E-004	5.70E-003	KEGG

Some associations are strikingly similar across cancers, for instance downregulated PI3K-Akt and Toll-like signaling and upregulated ErbB and signaling to ERK. Note that Pl3K (phosphatidylinositol 3-kinases) is a family of lipid kinases coordinating intracellular signaling in response to extracellular stimuli. In particular, they integrate signals from growth factors, cytokines etc., that control many processes, from cell proliferation to growth and survival, and others. Alterations in Pl3K are common in cancers,28 legitimating that they are key drug target for anticancer therapy. Instead, among the discordant regulations, we observed Jak-STAT and DNA repair. The Janus kinase-signal transducer and activator of transcription (Jak-STAT) pathway is also targeted by multiple drugs, but in our evidences, it appears downregulated in three cases and upregulated in one (RCC) case. Gene sets of either sign appear instead in PEL.

Regulatory Networks: Transcriptional Cores

We report the TFs that were found significantly enriched among the dysregulated genes by using the TFactS analysis tool.29,30 Then, given the target gene list identified in each tumor, we predicted their miRNAs regulators using experimentally curated interactions (complete coverage is given in other supplementary files). Figure 3A (two panels and two cancers) and 3B (two panels and three cancers) shows the “cores” of the five cancer regulatory maps, which were constructed by subnetworks involving gene sets enriching the top-3 pathway terms. Therefore, the corresponding connectivity patterns extracted from all regulatory maps give shape to such cores.

Figure 3

(A) BC and (B) Melanoma regulatory cores. (C) Mesothelioma, and (D) PEL (left panel) and RCC (right panel) regulatory cores. Symbols in use: node shape indicates epigenetically dysregulated genes (ellipse), TFs (hexagons), and miRNA (diamond); node color indicates DEG either up- (red) or down- (green) regulated. Node border color indicates the significantly enriched TFs (with P-value < 0.05), while the non-colored nodes correspond to TFs that regulate at least one target of the epigenetically modified genes. Then, the label color node reflects the presence of miRNA families.

Among the BC targets, RAF1 encodes the mitogen-activated protein kinase kinase kinase (MAP3K), which functions downstream of the Ras family, and once activated can phosphorylate to activate the dual specificity protein kinases MEK1 and MEK2. These in turn phosphorylate to activate the serine/threonine-specific protein kinases, ERK1 and ERK2. The two extracellular signal-regulated kinases are pleiotropic effectors of cell physiology and control gene expression affecting cell division cycle, apoptosis, cell differentiation, and cell migration. Another target is STAT5B, encoding a member of the STAT family of TFs, operating in response to cytokines and growth factors, ie, phosphorylating by the receptor associated kinases. It is involved in diverse biological processes, such as apoptosis, and its protein mediates the signal transduction triggered by various cell ligands, including IL2 and different growth hormones. Its target interleukin-2 receptor alpha (IL2RA) is found downregulated and jointly targeted by the upregulated JUN (also linked to the downregulated IL2), encoding a TF involved in cAMP signaling and interacting directly with specific target DNA sequences to regulate gene expression. It is mapped to a chromosomal region involved in both translocations and deletions in human malignancies. Another JUN target is the upregulated MYB, which may play a role in tumorigenesis.

Among the MEL targets, the downregulated CFLAR is a regulator of apoptosis and is structurally similar to caspase-8, but lacking caspase activity. Two downregulated TFs are involved, AR and CEBPA, the latter epigenetically involved in cancer.31 The androgen receptor (AR) was recently established as significantly mutated in a pan-cancer WGS analysis,4 and in combination with FGF family receptors, the two receptors are targets for tumorigenesis.32 Then, the downregulated TNFRSF10B, a TP53 target, encodes a protein member of the TNF-receptor superfamily, which contains an intracellular death domain and transduces an apoptosis signal. Finally, the downregulated CASP8 encodes a protein involved in the programmed cell death induced by FAS and various apoptotic stimuli. Among the MES targets, there are again the upregulated RAF1 and the upregulated MAP2K1, whose encoded protein acts as a MAP kinase. This also acts as an integration point for multiple biochemical signals, involved in many cellular processes such as proliferation, differentiation, transcription regulation, and development. Similarly, the upregulated neighbor gene DUSP6 encodes a protein member of the dual specificity protein phosphatase subfamily, which negatively regulates members of the MAP kinase superfamily (MAPK/ERK, SAPK/JNK, p38), in turn associated with cellular proliferation and differentiation and involved in cancer progression and resistance mechanisms. The master tumor suppressor gene TP53, whose diversity of oncogenic variants makes it atypical,33 appears downregulated. Also, an IL-dense hub with downregulated IL2, IL4, and IL2RA appears, which is involved in autoimmune diseases due to their causal variants.34 The upregulated IL1B is also emphasized. This is centered on the upregulated NFKB1 (known transcription regulator activated by various intra- and extracellular stimuli, such as cytokines) and can induce the expression of several growth factors and chemokines, playing a role in stromal cells. The cytokine IL2 regulates growth, proliferation, and differentiation of T-cells. Associated with IL2RA, there are Multiple Sclerosis (MS), RA, IBDs (inflammatory bowel diseases), and T1D. Note that NFκB is a generally expressed master transcriptional regulator of inflammation and that immunomediated diseases (IMD) are especially associated with REL and NFKB1 genes.35 Interestingly, the expression profiles of some of these diseases (IBD) show correlation with cancers, like colorectal cancer and the pathogenesis is regulated by the PI3K-AKT pathway through cell growth, proliferation, and cell death.36

Among the RCC targets, there is a quite dense hub centered on JUN, this time nondifferentially expressed, with IL-related genes (both downregulated, IL4, IL6 – proinflammatory cytokine activating STAT3 pathway and whose expression is regulated by NFκB- and IL18, and upregulated IL5). There are also TGF-β related genes (both upregulated TGF-β1 and downregulated TGF-βR2), regulating multiple physiological and pathological processes including control of mesenchymal cell proliferation and differentiation, ECM production, immuno suppression, and carcinogenesis. In particular, TGF-β is on one hand a tumor suppressor because of its action of contrast against epithelial cell growth, and on the other hand, it induces invasiveness and metastasis via epithelial mesenchymal transition (EMT) and is also correlated with resistance to anticancer agents. Linked to the hub, there is a tumor suppressor gene in several types of cancer, the downregulated TFPI2. Also appears downregulated PTPN12, known member of signaling molecules that regulate various cellular processes (cell growth, differentiation, mitotic cycle, and oncogenic transformation). Among the PEL targets, the connected central path includes the upregulated JUN, targeting the downregulated IFNG, encoding a protein member of the type II interferon family, with antiviral, immune-regulatory and antitumor properties. This is also targeted by the upregulated EGR1, whose encoded protein functions as a transcriptional regulator, promoting differentiation and cancer suppression. Of interest to remind that immune-suppressive cells, such as Treg cells, play a crucial role in maintaining the immune homeostasis, depending on the balance between the immune responses controlling infectious pathogens and cancer and the reciprocal immune responses preventing inflammation and autoimmune diseases. In particular, depletion of Treg cells presents a high risk of developing autoimmune diseases. CXCL12/CXCR4 (stromal-derived factor and its receptor, important markers of activated fibroblasts) acting at the interface between cancer and stroma, relevant for progression and also trafficking of cancer cells to organs such as lymph nodes (thus inducing metastasis), including Treg trafficking to the bone marrow.

Regulatory Networks: Post-Transcriptional miRNA Dynamics

The BC regulation activity of mir-7 and miR-195 involves the upregulated RAF1. By acting as an oncomiR in the epithelial cellular context, mir-7 promotes cellular transformation and tumor growth.37 Already found upregulated in RCC,38 affecting cellular migration, proliferation, and apoptosis, we do not detect it the RCC core. Then, miR-155 has important targets such as the upregulated RUNX2 involved in the vascular remodeling abnormalities in chronic kidney disease,39 and the upregulated IL8 involved in the inflammation,40 and widely associated with the miRNA in this pan-cancer (see common regulatory map of Fig. 4B). Then, miR-195 was lowly expressed in BC cells and multidrug-resistant BC tissues, in association with reduced RAF1 expression in vitro and ex vivo, and inducing miR-195 expression was shown to cause apoptosis and inhibit cell viability, but also sensitivity to treatment (Adriamycin).41 Two other interesting miRNAs are miR-21, which has different target genes, such as WNT1, and reappears in other cores and in the common map too, and miR-30a (which we see acting on Jun, the same also occurring in PEL) shown to function as a tumor suppressor gene in BC development and metastasis.42 The same mir-30a action is exerted in the MEL core by acting (jointly with TP53) on the tumor necrosis factor receptor superfamily member 10B (TNFRSF10B). The mir-7 to Raf1 activity is also visible in MES, in proximity of the upregulated target MAP2K1, target for instance of mir-128, reviewed in its associations with tumorigenesis and metastasis,43 and also indicating that DNA methylation in promoter regions may cause downregulation of the miRNA gene expression. MAP2K1 is also a target of mir-181a, part of a family of top predictors associated with T2D,44 due to the protein decarboxylase 2 also found in T1D, and to sirtuin-1 acting as positive regulator of insulin signaling,45 Then, NFKB1 is target of mir-21-5p and also of mir-155-5p, already found in the previous cancers, and predicted also in PEL with TGF-β1 as a target involved in EMT.

Figure 4

(A) Union of pan-cancer regulatory map. Panoramic view with symbols as before in other regulatory maps. (B) Common regulatory pan-cancer map. Intersection between individual regulatory networks. (C) Melanoma BeC. Red and green nodes for down- and upregulated genes, respectively, computed with Cytoscape plugin for Betweenness Centrality. JUN is the highest traffic hub, but a circled region is also relatively denser in crossings at some high-traffic nodes. (D) BC MST. Red and green nodes refer to down- and upregulated genes, respectively. MST computed from the R package igraph: http://cran.r-project.org/web/packages/igraph/. Other pan-cancer component networks are reported in SM. Two hubs are emphasized in yellow, IGF2 and CDK6.

Finally, note that the cores in MES (with upregulated targets MAPK1 and KIT, ie, the tyrosine-protein kinase regulating cell survival and proliferation, hematopoiesis, and related to the PI3k-Akt signaling cascade), and those in BC (with multiple upregulated targets, such as: IL8, PTGS2 producing prostaglandin E2 (PGE2) that modulates motility, proliferation, and resistance to apoptosis, and POMC (proopiomelanocortin) involved in a wide range of physiological functions, including pigmentation, energy homeostasis, inflammation, immunomodulation, steroidogenesis, and temperature control) share miR-335-5p, which acts as tumor suppressor blocking proliferation and invasion in renal carcinoma.46 Moreover, mir-30a plays a role in Toll-like receptors47 and acts on JUN in the PEL core, in combination with mir-155-5p, which has a protective function in human dermal fibroblasts by negatively regulating this TF.48 As in the BC core earlier, there are mir-15a and mir-16 proximal to JUN and MYB targets, while in proximity of SERPINE-like targets, there are mir-1 and mir-21 predicted as interactors in PEL and BC. In particular, as a tumor suppressor in multiple cancers, mir-1 high expression downregulates NOTCH3 and allows differentiation of myoblast cells.49 Only two miRNA appear in the RCC core, associated with the target CXCR4 (chemokine receptor 4), this last resulting highly expressed in BC too, and these are mir-146a and mir-224. The former is known as a major player in the control of hematopoiesis, immune function, and cancer,36 while the latter is involved in metastasis of human BC cells to the bone by directly suppressing the RKIP tumor suppressor (note that CXCR4 is a missing target in the BC core).

Topological Summaries

Topological network configurations offer various views at different granularity. Together with the common regulatory map previously discussed, we show for demonstrative purposes a much denser connectivity map in Figure 4A, displaying the union of regulations from all pan-cancer components. The complexity of the common map, given the reduced dimensions, is easier to visualize. It reveals regulatory hubs and patterns that are shared at pan-cancer scale, in which only some of the previously described master regulatory genes re-appear. Thus, TP53, IL4, REL A, and CASP8 are among the downregulated targets. Although PTPN12 appeared quite an isolated small hub in PEL, the common map gene PTPN13 encodes another protein of the PTP family (ie, signaling molecules that regulate cell growth, differentiation, and oncogenic transformation). Due to interactions with the FAS receptor, a possible role is in FAS-mediated programmed cell death. The TERF1 gene encodes a telomere-specific protein that is a component of the telomere nucleoprotein complex, functioning as an inhibitor of telomerase. The MMP16 product belongs to the matrix metalloproteinase (MMP) family, involved in the breakdown of ECM in metastatic disease processes, among other functions. The ECM components control cell behavior through differentiation, proliferation, cell morphology, etc. The NDRG2 contributes instead to the regulation of the Wnt signaling pathway, and since it downregulates CTNNB1-mediated transcriptional activation of target genes (say, CCND1), it may act as a tumor suppressor. Finally, the insulin growth factor (IGF) plays a role especially for stimulating neoplastic growth and in therapy (metformin).50 IGF2 is the most overexpressed gene in CRC, where it represents an actionable target for patients.51 In our study, it appears upregulated in BC. The bioactivity of IGFs is modulated by IGFBPs, stimulated by tumor suppressors (say, TP53). Both such genes are downregulated in the common map.

Moving to upregulated genes, interesting new players emerge together with the IL-family. An example is ERG, encoding a member of the ETS family of TFs, key regulators of embryonic development, cell proliferation, differentiation, angiogenesis, inflammation, and apoptosis. Interestingly, xeroderma pigmentosum C (XPC) loss has been associated with enhanced MEL photocarcinogenesis.52 The gene collagen, type 1, alpha 1 (COL1A1) has among its related pathways that are PI3K-Akt and ERK. Then, INS (downregulated in PEL) and COL1A1 (upregulated in the common map, while COL1A2 is upregulated in MEL) appear linked to T1D.53 The beta-2 microglobulin (B2M) gene has among its related pathways: the FGFR signaling, coupling to the MAPK and PI3-K/Akt intracellular signaling cascades, and crosstalking with the Notch signaling pathway. Upregulation of FGFR expression may lead to cell transformation and cancer. The FGF20 gene product is a member of the fibroblast growth factor family, with cell survival activities and involvement in various biological processes including embryonic development, cell growth, morphogenesis, tissue repair, tumor growth, and invasion. The RECK gene encodes a cysteine-rich protein (an extracellular protein with protease inhibitor-like domains) whose expression is suppressed strongly in many tumors and cells transformed by various kinds of oncogenes. In terms of miRNAs, mir-335-5p appears again to target the upregulated IL8, and mir-155-5p appears again too, but targeting the upregulated IL8 (inflammatory marker) and XPC (damage recognition). In relation with the IL8 target also mir-1 is active, while mir-21 previously associated with PEL is now directed to the downregulated RECK (thus facilitating tumor invasion and metastasis) and TP53 and to the upregulated IL1B. Notably, the target CASP8 appears exactly like seen before targeted in the MEL core. The previously unmet ERG target is surrounded by mir-133a/b and mir-145-5p. Overall, we have provided further evidence of the fact that relatively small miRNA sets regulate multiple oncogenic processes in pan-cancers, confirming what observed earlier.54

Network Fingerprints

The two selected network signatures calculated from each individual regulatory networks are BeC and MST, graphically displayed in Figure 4C and D, respectively. For reasons of space, only one representative pan-cancer component is chosen in each case for the display, leaving the rest to Supplementary files 12–15 for BeC and Supplementary files 16–19 for MST. To summarize some of the evidences from these fingerprints, the master BeC nodes are JUN in BC, RCC, and MEL (Fig. 4C), TP53 in MES, and MYC in PEL. The blue circles in Figure 4C embed a few high-traffic gene region including TFs such as AR, STAT3, CEBPA, CTNNB1, RUNX2, FOXO1, DEGs such as IL2RA and CDKN29 (both downregulated), and the other three upregulated DEGs, TGF-β1, XPC, and COL1A1. The signaling traffic between nodes is regulated primarily but not completely by such genes at network scale. As for the MST, they display maximally efficient communication between the nodes of a given network. More or less localized or centralized subnetworks are identified within each cancer map, and the reference hub genes are identified. In BC, the upregulated IGF2 and CDK6 appear as important hubs (Fig. 4D, yellow nodes). Other evidences (supplementary files) indicate that for MEL: JUN and AR appear as main hubs (respectively, up- and downregulated), followed by CEBPA and STAT3 (respectively, down- and upregulated); MES: the main hub is the downregulated TP53, while the other more peripheral hubs are almost always not detected as DEGs; PEL: the upregulated CDK6 is the main detected hub; and RCC: CDK6 is again a main hub, but together with other upregulated hubs, ITGB1, APP, and HOXA9. Recently, gene expression profiling across 12 cancers allowed to identify a supercluster significantly enriching for TP53 mutations and genomic loss of (cyclin-dependent kinase inhibitor 2A (CDKN2A).55 This provides instructions for making proteins p16 (INK4a) and p14 (ARF) function as tumor suppressors (ie, keeping cells from growing and dividing too rapidly or with no control). In particular, p16 binds to CDK4 and CDK6, regulators of the cell cycle. Thus, binding of p16 has the effect of blocking these proteins’ ability to stimulate cell cycle progression, thus controlling cell growth and division (the entire series of BeC and MST figures is reported as other supplementary files).

Morbidities

Associated with the pathway landscape map of Figure 2A, a disease map is obtained from KEGG-specific disease annotations (therefore, it is limited to such annotation domain). The most peripheral nodes in Figure 2B indicate the pan-cancer components, surrounded by the DEG chain involved in the highly significant annotations enriched by the DEG sets. Up-and downregulation directionalities are expressed as red and green links, respectively. Annotations for other different cancers appear, together with other pathological conditions directly or indirectly associated with the five cancers through the genes. T1D appears, together with other sporadic annotations for autoimmune conditions (thyroid), genetic diseases (fanconi anemia, associated with a high incidence of tumors), and chronic diseases (asthma). It is noticed the presence of viruses, and other infection terms from bacteria (tuberculosis and legionellosis) and parasites (leishmaniasis). Many infectious agents are linked with cancer through NFκB factors, as in case of human papilloma, HIV, human herpes viruses, EBV, hepatitis B and C, etc. However, also hyperglycemia is involved with NFκB activation, in association with diabetes and obesity. Notably other growth factors (EGF) and cytokines associated with tumorigenesis are linked to NFκB.56

Drug Maps

The drug map for the reference pan-cancer is illustrated in Figure 2C and also derived from the KEGG-driven diseases (Fig. 2A). The source of gene–drug interactions is MyCancer-Genome (http://www.mycancergenome.org/). A closer look at the significant annotations indicates two master genes, BRAF and PIK3CG as hubs because of incoming links (from diseases) and outcoming links (to drugs). The MAP2K1 gene (top left) is instead high in degree only for the external connectivity (drugs). Note that BRAF is targeted by colorectal, pancreatic, bladder, thyroid cancer nodes, and also by hepatitis C. The path observed from BRAF to the drug REGORAFENIB, an inhibitor of multiple tyrosine kinases responsible for turning on tumor cells, particularly relevant as treatment in metastatic cases (GIST and colorectal cancer). REGORAFENIB links also to RAF1, in turn seen upregulated in MES. Finally, again linked to the pathway landscape and its disease terms, we display in Figure 2D the druggable map. The outermost layer represents selected druggable gene categories: kinase, cell surface, clinically actionable, drug resistance, growth factor, TF binding, histone modification, druggable genome. Our gene sets have identified druggability linked to targets through the diseases. Note that many druggable proteins, ie, products of the mapped genes, function as protein-interaction network complexes whose modulation involves TFs and implies consideration as therapeutically targets.57 Of interest, the category of drug resistance, with only one predicted gene as interactor, ie, APC (downregulated in PEL), encoding a tumor suppressor protein antagonist of the Wnt signaling pathway was involved in cell migration and adhesion, transcriptional activation, and apoptosis processes. Tumor-promoting mutations seem to be involved in cell survival, cell fate, and genomic stability.58 The second type of process is influenced by mutations in APC, among other genes (NOTCH, AR, etc.), and the loss of normal functions due to mutation can contribute to chemotherapeutic resistance. Mutations in APC have been identified in early stages of cancer development making it a gatekeeper of tumor progression and therefore an ideal therapeutic target.59

Discussion

Back to the key domain of cancer microenvironment, the reduced expression of inflammatory markers (cytokines, chemokines) and growth factors could allow longer treatments and introduction of novel inhibitors targeting stromal cells against cancer growth and metastasis and toward lower risks in terms of drug resistance.60 Our better understanding of the role of tumor microenvironment (a mixture of fibroblasts, leukocytes, endothelial cells, extra cellular matrix, etc.) and the immune responses in mediating and regulating anticancer immunity has motivated the development of vaccine and immunotherapeutic approaches for treatment. Benefits of combined approaches have appeared in mice studies, targeting multiple immune system points, as in the case of gene silencing of TGF-β1 reducing the Treg cells.61 PL3K inhibitors are in clinical trials due to the ability to arrest cell proliferation and induce apoptosis and also for inhibited induction of Treg cells. This indicates that a promising direction is to enhance tumor-specific immune responses and transiently blocking immunosuppressive networks (at the local lymph node level). Among the mentioned IMD, autoimmune and inflammatory diseases are included, such as asthma, RA, IBD, and T1D. Key elements regulating the immune response are the human leukocyte antigen (HLA) genes (we found HLA-related genes downregulated in the RCC regulatory core).

Immune-system regulations are complex ones and strongly overlay other types of cancer mechanisms, despite remaining only partially known. The complexity of the synergistic modulation between cancer and the immune system justifies the fact that cancer immunotherapies are increasingly investigated and proposed.62 A requirement for inference models is to be able to deal with data generations from integrative omics approaches, some focused on pan-cancers. Our study is based on a relatively small pan-cancer, treated with demethylation experiments, and yielding as a testable hypothesis that cancer proliferation and other hallmarks of disease (cancer and noncancer) could be contrasted by targeting epigenetically driven alterations. The latter do not involve alterations of DNA sequence but induce changes in gene expression and thus represent complex regulators. Indeed, several detected DEGs provided, through their multiple annotations, evidences of cancer and immune function signatures in terms of pathways, regulatory networks, drug–target, and druggable maps. Our special focus on multilayered analysis of regulatory networks presents a methodological novelty with regards to pan-cancer inference, in particular with epigenetic effects at play. Networks’ power comes from their multiresolution, useful in biology to embrace relatively fast (protein–protein) and relatively slow (transcription) interactions, and visible through the configuration hierarchies displayed as communities, ie, by reorganizing modularly nodes and links of different connectivity densities. Various topological features of networks can be described, and this has been tremendously useful for inference purposes. Furthermore, networks are highly flexible for biological applications, crosslinking various bioentities, and specializing into ad hoc representations, some simply associative or causal, others integrating information on disease and drug–target interactions.

Notably, pan-cancers drive the attention over the current limitations of inference approaches that we can try to adapt to cancer diversity or heterogeneity. The examples of differentiated therapy effects in cancers, given the same targets, say BRAF,63 induce to reconsider the complexity of cancer systems, and explore both the commonalities and specificities of cancers emerging from comparative systems analyses rather than conducting replicated inferences on an individual basis. Asking whether cancers are more or less similar in network terms is probably not the ultimate question, because there are so many ways to look at networks that we need to make our analyses contingent on choices to some extent. But once some hypotheses are established, and for instance, they specifically refer to cancer phenotypes, networks allow to build precious fingerprints, as we have shown with regulatory cores or through some of the topological characterizations. Among them, we have used two types of molecular regulation mechanisms (TFs and miRNAs) and then two types of topological properties (MST and BeC), and we could observe particular features emerging from shared versus specific pathway landscapes, and partially concordant or discordant when looking at regulations underlying the enriched terms. With combined TF-miRNA regulatory networks, complexity assessment was dramatically reduced by building network cores, under simplified constructions. A certain differentiation in connectivity patterns and biological annotations was found in the individual maps, reflecting cancer heterogeneity that is hardly summarized by signatures or motifs. Notably, relevant gene regulations subject to direct or indirect epigenetic influences can be captured even by moderate levels of connectivity, restricted in our case to a few gene sets.64 We stress a fact that the more computationally demanding these integrative modeling approaches are expected to be, owing to the Electronic Medical Data becoming more open and interlinked, the more the personalized therapies are destined to obey to: (a) the translation of patient molecular profiles into re-phenotyped signatures with renewed functional and disease impacts; and (b) combined multitarget therapies, for instance ad hoc to adapt to cancer cell subpopulations.

Methods

Data generation and expression values

The cell lines were cultured in Dulbecco’s Modified Eagle’s medium (DMEM) supplemented with 10% fetal bovine serum, 2 mM l-glutamine, at split ratio of 1:4 twice a week. After 24 hours of split, the culture medium was changed with media containing 2,5 μM 5-aza-2-dC (DAC). The treated cells were collected after 48 hours, and the mRNA was extracted for pelleted cells. The two steps for cDNA microarray analysis were mainly: (a) total RNA samples were isolated from treated/untreated cells using TRIZOL reagent (Invitrogen); and (b) concentration of purified RNA samples were determined by A260 measurement and the quality was checked by Lab-on-a-chip analysis (total RNA nanobiosizing assay; Agilent) with the Agilent 2100 Bioanalyzer. RNAs isolated from different tumor cells, and transcribed in cDNAs, were used to carry out the analysis. The cDNAs from treated BC were labeled with cy5 red fluorescent dye and untreated BC with cy3 green fluorescent dye. Hybridization was done on a microarray chip called MWG Human Cancer Array containing 50-mer oligo probes for 1920 genes (1853 human genes associated with cancer, 27 control genes, and 40 replicated genes). Spots of fluorescence intensity were read by dual laser scanner (BioDiscovery) and analyzed with Mavi Pro-2.6.0 (MWG Biotech) to process the gene expression values. Background subtraction, normalization to a number of housekeeping genes, and comparison between treated and untreated cancer cell lines were performed. In order to select the dysregulated genes, we considered fold change calculated as a cy5/cy3 normalized ratio (NR), which was calculated for each gene and by taking the ratio of the intensity in cy5 (Ic5) and the normalized intensities in cy3 (nIc3). Then, to reduce variability, all ratio values were transformed in log2. We thus selected differentially expressed genes, both upregulated and downregulated, with their values in a range established by the following thresholds: log2 (NR) ≥1.5 and log2 (NR) ≤−1.5, respectively (Supplementary file 1).

Pathways

Pathways were sorted by enrichment scores associated with the DEG detected in each cell line. Over-representation analysis (ORA) was performed through Consensus PathDB (http://cpdb.molgen.mpg.de/). In particular, ORA allows interactive querying to compute a functional enrichment by the hypergeometric test, by assigning a P-value to each predefined gene set on the basis of the number of gene IDs present in both the predefined sets and the user-specified lists (up- or downregulated genes). The background size is the number of ConsensusPathDB entities that are annotated with a gene ID and participate in at least one pathway (depending on which of these predefined classes are considered by the user). The P-values are corrected for multiple testing using the false discovery rate, thus yielding q-values. The gene annotations are found in standard database resources, such as KEGG (http://www.genome.jp/kegg/), Biocarta (http://www.biocarta.com/genes/index.asp), Reactome (http://www.reactome.org/), and Wikipathways (http://www.wikipathways.org/).

Regulatory networks

We constructed a comprehensive set of regulatory network through Cytoscape (http://www.cytoscape.org/) for each pan-cancer component, including transcriptional and post-transcriptional interactions between all DEGs, and extracted the specific versus common signatures of the epigenetically modified genes. Transcription factor (TF) and microRNA regulators were predicted from the targets in the DEG lists, in an attempt to reveal regulation mechanisms underlying the epigenetically modified genes. The tool TFactS (http://www.tfacts.org/) contains TF-responsive genes obtained from experimental evidence reported in literature. Two data-sets are included: (i) a sign-sensitive catalog, which indicates the type (up or down) of TF-regulation exerted on its targets; and (ii) a sign-less catalog, which includes all regulatory interactions contained in sign-sensitive, and further interactions without the specific type of regulation. TFactS takes as a query the two lists of up-/downregulated genes and compares them with sign-sensitive catalog of manually curated annotated target genes, then returning the lists of activated and inhibited TFs whose annotated target genes show a significant overlap with the query genes. MiRNAs regulate protein biosynthesis at the post-transcriptional level and participate in the pathogenesis of different types of human cancers inducing or suppressing their progression. Our predicted miRNAs may affect the regulation of DEGs in different cancers and play a role in the most enriched pathways. We then merged experimentally validated miRNA-target gene databases. In particular, three major tools were used for the analyses: miRTarBase V.4.5 (http://mirtarbase.mbc.nctu.edu.tw/),65 miRecords V.3 (http://mirecords.biolead.org/),66 and miR2Disease (http://www.mir2disease.org/).67 We built a nonredundant dataset of miRNA-target genes regulatory interactions in human and used such aggregate dataset to predict the miRNA regulators among the list of significantly dysregulated genes from each cancer cell line.

The pan-cancer circuit of mixed regulations was reduced by extracting the top-3 enriched pathways and highlighting the connectivity patterns activated by the corresponding DEG sets. We called such subnetworks as “regulatory cores.” Such reduced network configurations help the interpretation of the resulting evidences, by establishing criteria to improve over single-gene signatures.68 Apart from similarity/dissimilarity of network configurations in light of the role of major regulative motifs, identifying the presence of recurrent paths in pan-cancers is an interesting goal. Common targets could be the reason, and not only miRNA biomarkers could be identified as potentially therapeutic too but also annotations could be found of complementary value and significance within pathways; further contextual information (protein interactions, wide spectrum transcriptome profiles, etc.) could help determining function diversity. Notably, cancer synergism is a crucial factor that would imply the possibility for some complex relationships to be pinpointed by targeting marker ensembles, such as complexes or modules, rather than at individual ones,69 such as single genes or proteins. This difference may be particularly true with multiple cancers, despite the hallmark of cancer heterogeneity suggests to consider the distinct characteristics of each cancer.

Integrative Network Inference

Topological maps

The topology of a network consists of the structure of connectivity patterns and the corresponding dynamic communication exchange between nodes. Topological signatures are functional to the characterization of the pan-cancer in terms of connectivity and centrality. We selected a few measures suitable to our scopes. First, the Minimum Spanning Tree (MST) offers a tightly connected view of the entire network configuration. A spanning tree of an undirected graph is a subgraph that connects all the nodes, and is a tree. One problem is how to find the minimum length spanning tree, and this can be solved in various ways according to different algorithms. A single graph can have many different spanning trees. Even if MST offers an all-connected network view, which remains redundant for some aspects, this may be a valid support to the determination of gene prioritization strategies. Then, the Betweenness Centrality (BeC), which identifies the main routing mechanisms in the network, ie, the nodes concentrating through their links a substantial traffic volume. Notably, genes that are critical to cancer mechanisms have been observed to code for high-connectivity and central proteins,70–73 thus reflecting a diffuse control of the communications between network nodes. The BC of a node v is the sum of the fraction of all-pairs shortest paths that pass through v: where V is the set of nodes, σ(s,t) is the number of shortest (s,t)-paths, and σ(s,t|v) is the number of those paths passing through some node v other than s or t.

Disease maps

Ensemble network dynamics determine complex diseases, intuitively enough. Furthermore, when associations between diseases are of interest, similarity from a pathobiological perspective correlates with interactome distance, thus leading to a better identification of disease phenotypes,74 with many still to be explained,75 and including false positives76 and also comorbidities (involving co-occurrence of multiple diseases in the same patient).77,78 The disease maps here proposed are associated with gene sets annotated from the KEGG disease db (http://www.genome.jp/kegg/disease/), in which different types of diseases were filtered according to the most enriched pathways to reduce term redundancy. Each disease entry contains a list of known genetic factors (disease genes), environmental factors, diagnostic markers, and therapeutic drugs bringing information on disease associations. The simple idea underlying network views employed in disease-omics and pharmacogenomics is that diseases can be associated with perturbed states of molecular systems, and while the known disease genes are genetically perturbing agents, drugs are therapy-induced perturbing agents. The Human Diseases category of the KEGG db is a collection of perturbed molecular networks including multifactorial diseases such as cancers, immune system, infectious, neurodegenerative, cardiovascular, and metabolic diseases.

Drug maps

Drawing network interaction maps allows linkages between drugs and bioentities such as genes (proteins) and diseases to be established, including the interrelationships with novel compounds that may explain drug mechanisms of action. Consequently, network-centered studies about drug–target interactions have significantly expanded, especially those aimed to find co-occurrences and interdependencies of drugs and their effects in patients. Another aim is to assess how multiple cancer treatments may connect, which would open for more efficient and timely targeted combined therapies, considered the preferred ones in cancer patient treatments. As a result, two critical aspects emerge: (1) targeting multiple pathways allows to investigate synergizing drugs effects, and (2) critical nodes in the cascading network paths become real identification targets, because subtly enough it is through them that one can hope to decrease the robustness of interdependent cancer subnetworks. Identification of critical genes in multiple cascades is suggesting that hierarchies could be revealed, such as targets that are key for drug response and others that enable either compensatory or antagonist behaviors, for instance contributing in the latter scenario to drug resistance. For instance, betweenness centrality is a topological measure revealing the gatekeeping checkpoints that regulate cascade signaling. In order to build gene–drug links from the pan-cancer, including the annotated diseases, we used DGIdb and two of its embedded functions designed to: (a) retrieve the interactions between gene targets and drugs, and (b) search for druggable categories. The queries were again our DEGs involved in the disease map. We built drug-target maps in association with a simple KEGG-driven disease map, ie, by annotating disease genes linked to known drug–gene interactions. Then, we used the information available on potential druggability. Among the many available tools, we exploited the drug–gene interaction database (DGIdb) (http://dgidb.genome.wustl.edu/).79 This is a mining tool retrieving information from various sources such as PharmGKB (https://www.pharmgkb.org/), Therapeutic Target (http://bidd.nus.edu.sg/group/cjttd/), DrugBank (http://www.drugbank.ca/), etc. Hypotheses can thus be generated about mutated genes, whether these are existing therapeutic targets from a library of ∼15,000 drug–gene interactions involving >2000 genes and >6000 drugs, or prioritized future drug targets, ie, >6000 druggable genes from 39 established categories.

Recently, an in silico drug prescription strategy has been reported.80 It aimed to come up with a landscape offering: (i) identification of driver events; (ii) collection of therapeutic agents targeting them; and (iii) connection of each patient to all targeted therapies. As a result, especially the identified dysregulative genes, ie, cancer drivers of cancer-related regulatory processes linked to pathways or hallmarks, can be usefully summarized in functional groups such as chromatic regulatory factors, splicing and mRNA processing, and ubiquitin-mediated proteolysis system. Interestingly, only a limited number of drivers clonally dominate the tumors. Therapeutic options can be molecules known to inhibit activated drivers (direct targeting, say Vemurafenib for V600E-activating BRAF mutation), and molecules inhibiting nonaltered proteins functionally connected to the altered drivers (indirect targeting, say temsirolimus, an MTOR inhibitor treating patients with PTEN-inactivating mutations). Gene therapies can be compensating the loss of activity of a tumor-suppressor driver. Apart from targeted therapeutic options, also drivers with a protein structure suitable to small molecule binding (drug-gable, therefore), and also others can be identified by accessing antibody and protein therapies (biopharmable). Following the above line of research, a restricted gene set bearing mutations in >5% of tumors of at least a single type could represent the best candidate set from drug development. Conversely, polypharmacology could target multiple altered drivers in a tumor, following the principle of drug combinations, thus enabling smarter anticancer therapies and strategies linked to the landscape of combination therapies. Overall, two main limitations are always applying: some genes may exhibit divergent modes of action on different cancer types and also the incompleteness of the multiple drug–target interaction DBs leads to an overall underestimation of the number of targeted drivers.

The symbols used in the regulatory networks of all types, global and cores, and later also for other signature networks, are shown in Box 1.

Supplementary Table 1. Pan-cancer DEG list.

Supplementary Table 2. Annotations for all pathway terms.

Supplementary Table 3. Gene list from intersection of pan-cancers, and expression values.

Supplementary Figure 1. Significant pathway map (q-value ≤ 0.01).

Other supplementary files Tables and files covering pan-cancer evidences and networks.

Betweenness centrality figures (75MB) - http://www.la-press.com/cr_data/large_files/31809_bec.zip.

Regulatory networks (65 MB)- http://www.la-press.com/cr_data/large_files/31809_reg_net.zip.

Table 2B

Top-3 pan-cancer pathway terms: MES, RCC, and PEL.

PATHWAY NAME	SET SIZE	MES DEG DOWN-REGULATED	P-VALUE	Q-VALUE	SOURCE
Cytokine-cytokine receptor interaction	267	CSF2RB; TGFBR2; IL9R; TNFSF15; IL2RA; IFNA1; IL13; IL4; IFNA5; TNFRSF4; BMPR1B; TNFSF8; IL2; IL26; IL11RA; TNFRSF18	7.09e-06	0.001	KEGG
Jak-STAT signaling	158	CSF2RB; IL2RA; IL9R; IFNA1; IL13; IL4; IFNA5; PIK3CG; AKT2; IL2; IL26; IL11RA	1.01e-05	0.001	KEGG
PI3K-Akt signaling	346	IL2RA; PPP2R1B; TP53; ITGB6; TSC2; IFNA1; RXRA; IFNA5; PIK3CG; ITGB3; AKT2; IL2; FGF20; IL4; FGF12; RELA; FGFR2; PPP2R5B	1.3e-05	0.0011	KEGG
		MES DEG Up-regulated
Signalling to RAS	25	MAPKAPK3; NTRK1; MAP2K1; MAPK1; RAF1; RALB	3.88E-06	4.62E-004	Reactome
DNA Repair	114	BRCA1; MRE11A; XPC; LIG3; CDK7; MUTYH; POLD1; FEN1; RAD50; DDB2	5.43E-06	4.62E-004	Reactome
Signalling to ERKs	35	MAPKAPK3; NTRK1; MAP2K1; MAPK1; RAF1; RALB	3.07E-05	1.42E-003	Reactome
		RCC DEG Down-regulated
T cell receptor signaling	108	IFNG; FYN; IL4; IL5; PTPRC; RELA; CD3D; LCK; MALT1; CD3G	7.03E-07	5.55E-05	KEGG
NF-kappa B signaling	92	ATM; TNFAIP3; CD14; RELA; LCK; MALT1	8.26E-04	7.05E-03	KEGG
PI3K-Akt signaling	346	BRCA1; IL2RA; GH1; IFNA5; INS; IL4; MCL1; FGF12; RELA; MYB; FGF20; EIF4B	9.45E-04	7.46E-03	KEGG
		RCC DEG Up-regulated
Jak-STAT signaling	158	IL13RA2; IL5; STAT3; EPO; IL21R; PIK3CG; STAT5B; IL22; IL26; CREBBP; IL9	1.38E-05	9.35E-04	KEGG
ErbB signaling	88	TGFA; SHC1; CRKL; PIK3CG; STAT5B; MAPK1; BTC; CRK	3.36E-05	1.32E-03	KEGG
Focal adhesion	206	PDGFRA; SHC1; CRKL; PPP1CA; PIK3CG; MAPK1; PDPK1; HGF; CRK; FLT4; VTN; CAV1	3.39E-05	1.32E-03	KEGG
		PEL DEG Down-regulated
Cytokine-cytokine receptor interaction	267	IL18; TGFBR2; IFNA8; TNFSF15; ACVR1B; CSF2RB; IL13; IL4; IFNA5; IL6; TNFRSF10B; BMPR1B; TNFSF8; TNFRSF9; TNFRSF18; TNFRSF11A	1.5E-06	1.70E-04	KEGG
DNA Repair	114	BRCA1; LIG4; POLD1; TDG; TCEA1; MUTYH; FANCF; PCNA; XRCC1; PRKDC	6E-06	3.85E-04	Reactome
TGF beta Signaling	131	CDKN2B; TGFBR2; ITGB4; RBL2; LIMK2; JUN; MAP2K3; SKIL; MYC; BCAR1	2.06E-05	8.37E-04	Wikipath-ways
		PEL DEG Up-regulated
Jak-STAT signaling	158	IL13RA2; IL5; STAT3; EPO; IL21R; PIK3CG; STAT5B; IL22; IL26; CREBBP; IL9	1.38E-05	9.35E-04	KEGG
ErbB signaling	88	TGFA; SHC1; CRKL; PIK3CG; STAT5B; MAPK1; BTC; CRK	3.36E-05	1.32E-03	KEGG
Signalling to ERKs	35	CRK; MAPK1; NTRK1; CRKL; SHC1	1.24E-04	2.66E-03	Reactome