A multiple network-based bioinformatics pipeline for the study of molecular mechanisms in oncological diseases for personalized medicine.

MOTIVATION: Assessment of genetic mutations is an essential element in the modern era of personalized cancer treatment. Our strategy is focused on 'multiple network analysis' in which we try to improve cancer diagnostics by using biological networks. Genetic alterations in some important hubs or in driver genes such as BRAF and TP53 play a critical role in regulating many important molecular processes. Most of the studies are focused on the analysis of the effects of single mutations, while tumors often carry mutations of multiple driver genes. The aim of this work is to define an innovative bioinformatics pipeline focused on the design and analysis of networks (such as biomedical and molecular networks), in order to: (1) improve the disease diagnosis; (2) identify the patients that could better respond to a given drug treatment; and (3) predict what are the primary and secondary effects of gene mutations involved in human diseases.
RESULTS: By using our pipeline based on a multiple network approach, it has been possible to demonstrate and validate what are the joint effects and changes of the molecular profile that occur in patients with metastatic colorectal carcinoma (mCRC) carrying mutations in multiple genes. In this way, we can identify the most suitable drugs for the therapy for the individual patient. This information is useful to improve precision medicine in cancer patients. As an application of our pipeline, the clinically significant case studies of a cohort of mCRC patients with the BRAF V600E-TP53 I195N missense combined mutation were considered. AVAILABILITY: The procedures used in this paper are part of the Cytoscape Core, available at (www.cytoscape.org). Data used here on mCRC patients have been published in [55]. SUPPLEMENTARY INFORMATION: A supplementary file containing a more detailed discussion of this case study and other cases is available at the journal site as Supplementary Data.

Introduction

Recent studies suggest that cancer can be better understood through the study of mutated or dysregulated pathways/networks [1]. Human diseases are not caused by single molecular defects but driven by complex interactions among a variety of molecular mediators [2, 3]. The analysis of human diseases using graphs and/or biological networks plays an extremely critical role in the field of precision oncology [4, 5]. The study of these complex networks can reveal new disease-associated genes and/or pathways and identify possible targets for new drug development, as well as new uses for existing drugs [6].

Network-based approaches could be the most promising strategies for identifying the specific mediators who are responsible for altering the networks, trying to facilitate the development of new combinations of drugs based on the complex interactions involved in tumor-growth, thus improving personalized medicine [7, 8]. It is well known that not all patients respond at the same extents to the same pharmacological treatments. This is true also for targeted therapies. Although patients are selected on the basis of specific biomarkers that are often represented by genetic alterations, only a fraction of patients usually respond to treatment. To further develop precision medicine, it is essential to get an in-depth view of what are the factors that determine the resistance to a therapeutic treatment and identify which patients can better respond to specific therapies [9].

Targeted therapy is based on the identification of driver alterations that allow therapeutic intervention with specific inhibitors. For example, the presence of BRAF mutations in metastatic colorectal carcinoma (mCRC) predicts the response to combinations of BRAF inhibitors, MEK inhibitors and anti-EGFR monoclonal antibodies [10]. However, mCRC have been found to carry often multiple genomic alterations. The effects of these multiple driver alterations on the probability to respond to targeted agents are not known. In this respect, recent data suggested that personalized treatment with combination therapies would improve outcomes in patients with refractory malignancies as compared with single agents [11], thus highlighting the importance of taking into consideration the complex genomic alterations of cancer in order to improve precision medicine. In the present work, we illustrate an innovative bioinformatics pipeline based on the design and analysis of biomedical and molecular networks aimed at improving the disease diagnosis, identifying the patients that could better respond to a given drug treatment, and predicting the primary and secondary effects of gene mutations involved in human disease.

Our strategy is focused on ‘multiple network analysis’ with which we look for improving cancer diagnostics using biological networks, and on generating statistical inference predictive models to probe regulatory relationships between molecular components such as genes or proteins [12, 13]. This approach allows us to represent the human physio-pathological system and biological processes in a simple and complete way [14]. It arose great interest among scientific community [15]. Once the disease-mutated genes of the patient are known, it is possible to understand how they are involved in pathological processes and what are their combined effects on the modification of molecular mechanisms in order to evaluate the correct treatment and to help physicians for what concerns diagnosis, prognosis and therapy.

Methods

Pipeline workflow

The pipeline developed by us and carried out in this work is characterized by a multistep design, which concerns different biological-molecular networks (Fig. 1): disease–disease, gene–disease, gene-variant-disease, gene–gene, protein–protein interaction and multilayer drugs. The last three steps can be performed in both directions, i.e. from gene–gene network to multilayer drug network and vice versa. All these steps have been performed through Cytoscape v3.7 Core (with its specific plugins used for the realization and analysis phase (Table 1) [16, 17].

Figure 1

Multiple network-based workflow.

Table 1

Pipeline plugins

Step of pipeline	Plugin names	Algorithm	Description of plugin	References
Disease–disease network	autoHGPEC and NetworkAnalyzer	RWRH (random walk with restart on heterogeneous network) algorithm to compute the similarities between diseases and identify what are the main and strongest correlations within the network. Organic layout algorithm allowed us to evaluate the distance between diseases beyond the evaluation of local and global topological properties. NetworkAnalyzer algorithm performs analysis of biological networks and calculates network topology.	To identify and predict novel disease-disease associations. Using NetworkAnalyzer it is possible to compute basic properties of whole network (degree distribution, clustering coefficients, centrality, etc.)	[18–21]
Gene–disease network and Gene-variant-disease network	autoHGPEC and DisGeNet	RWRH algorithm. LHGDN machine learning algorithm to extract novel gene-disease associations	To identify and predict novel gene-disease associations. Useful to analyze the role played by hub genes and investigate human complex diseases with respect to their genetic origin by a variety of built-in functions.	[18–24]
Gene–gene network	GeneMANIA plus StringApp and iCTNet2 (Cytoscape plugin Integrated Complex Traits Networks)	It uses two algorithms: (1) a linear regression algorithm to compute the functional gene–gene association networks and (2) a Gaussian field label propagation algorithm for predicting gene functions from the composite network. It uses (1) a naive Bayesian algorithm to compute combined scores from different edge types and (2) an approach based on the closest combined scores to grow the query network. It applies two different algorithms: RWRH and PRINCE, which uses network topological characteristics in the protein interaction network to prioritize candidate genes.	They have been used to analyze and investigate different types of biomedical-molecular interactions, by crossing and verifying the results obtained with what is reported in the scientific literature in the various studies. They have been used to integrate several data sources to allow automated and systematic creation of networks with up to five layers of omics information: phenotype-SNP association, protein-protein interaction, disease-tissue, tissue-gene, and drug-gene relationships.	[25–27]
Protein–protein interaction network	GeneMANIA, FunMod and ReactomeFIPlugin (Reactome)	It uses two algorithms: (1) a linear regression algorithm to compute the functional association networks setting protein–protein interaction parameters and setting molecular mechanisms and (2) a Gaussian field label propagation algorithm for predicting gene functions from the composite network. FunMod iteratively selects all edges of the network and assigns a functional annotation to an edge when two linked nodes are annotated in the same biological group or pathway in the ConsensusPathDB (DB) database. It uses dfferent kind of algorithms: HotNet to search for network modules; MCL Clustering algorithm based on spectral partition; Algorithms for detecting significantly mutated pathways in cancers	They have been used to analyze and investigate different types of biomedical-molecular interactions. FunMod extracts all pairs of nodes annotated for the same pathway in a new sub-network. Subsequently, FunMod tests the statistical significance and calculates the topological properties of the sub-network to identify the sub-networks that are statistically enriched in biological functions and that exhibit interesting topological features. The statistical significance of the sub-network is determined by performing a hypergeometric test, a well-established method used in gene enrichment analyses It explore Reactome pathways and search for diseases related pathways and network patterns using the Reactome functional interaction network	[25;28–29]
Multilayer drug network	GeneMANIA and DrugTargetProfiler (DTP-HDR); plus Clinicaltrials.gov-Eudract trial register for drugs; CIVIC Cancer Portal and SIDER 4.1 web-server in which are implemented DGIdb3.0 and Genomics of Drug Sensitivity in Cancer (GDSC); CancerDR web server implemented with Cancer Drug network map: all cancer drugs plus ResMarkerDB	It uses two algorithms: (1) a linear regression algorithm and (2) a Gaussian field label propagation algorithm. Barnes–Hut algorithm is an approximation algorithm for performing an N-body simulation and to realize multilayer drug networks. HDR is a random walk with restart algorithm applied on a heterogeneous network of drugs and diseases, to predict novel drug–disease associations	They have been used to analyze and investigate different types of biomedical-molecular interactions and to study some important side effects of drugs selected for pharmacological therapies. Applying drug-disease centric parameters is useful to calculate what drugs can be used for a given disease and for mutated genes, highlighting which group of drugs affects specific genes. Moreover, it is possible to interpret the variant involved in disease of interest describing the therapeutic, prognostic, diagnostic and predisposing relevance of inherited and somatic variants of all types. This information is useful to redesign the drug–gene interactions, to investigate the therapeutic biomarker involved in cancer cells and to identify the most important Drug-Resistance.	[25; 30–45]

Disease–disease network

As shown in the step a of Fig. 1 the realization of this network is focused on the comprehension of the disease at molecular level. This is essential to investigate the associations/correlations between our disease and the other ones that directly or indirectly could be related to them, to try to understand the disease progression and to explain the secondary symptoms that a patient could develop over time [46]. In this network, the nodes are the pathological phenotypes and the edges represent the correlations between the pairs of pathologies. In the pipeline, the disease of interest represents the input, while other diseases related to it are the output generated by Cytoscape. The realization of this network allows us to have a useful map that will help us identify the priority genes and hubs in the next study steps.

Gene–Disease network

After the generation of the disease–disease network, it is possible to identify the priority genes associated among the disease of interest and the others, creating the gene-disease network as shown in step b of Fig. 1 [47]. The output is a network in which some nodes represent the disease of interest and some others represent the priority genes mainly involved in it. The edges represent the role played by the priority genes within the network. The aim of this process is to identify ‘hubs’ and understand how a genetic perturbation/mutation can change the molecular profile. In our pipeline, the disease represents the input, while the identified hubs are the output of our analysis.

Through this network, we can extrapolate and explore the driver genes involved in the disease to explain and justify the onset of secondary symptoms developed over time by a patient with cancer. Moreover, this network is useful to identify the hub genes important to investigate how the molecular profile will change with one or more mutations. Once this network is created, it is possible to generate the networks of the following steps and extrapolate the information useful for pharmacological purposes to identify the best therapies for individual patients.

Gene-Variant-Disease network

It is useful to identify which are the single/double/multiple mutations involved in the disease of interest, to investigate the role played by the mutations involved in the disease as well as to assess which are the priority genes linked to these mutations. Therefore, it is useful to realize the gene-variant-disease network as shown in step c of Fig. 1. In this network, there are some nodes representing the disease of interest, some nodes representing the priority genes involved in specific diseases and other nodes that are the principal genetic variants involved in our context. Furthermore, edges represent one or more associations between different nodes. This network is realized using the same pipeline of the gene-disease network, but the input is a gene, or a disease, or both. This network is useful for identifying clinically relevant and useful variants or SNPs in the molecular study in order to select the most suitable therapy.

Gene–Gene network

After the study of the priority genes and of the mutations mainly involved in our disease, it is possible to create the gene–gene network as shown in step d of Fig. 1, with the aim of identifying the gene interactions involved in the physiological and pathological field [48], evaluating their differences and analyzing how molecular behavior changes in the presence of single, double or multiple mutations [49]. The study of this network can facilitate the systematic functional annotation of genes [50] and help identify the hub genes, which may lead to potential clinical applications [51]. It is a network in which the nodes represent genes (control or mutated), while the edges represent the gene interactions (physiological or pathological). In the gene–gene network, the input is given by our gene panel to which the mutations are added if the mutated networks are to be realized. Then, control gene–gene networks (healthy subjects) and mutated gene–gene networks (subjects affected by the disease resulting from single, double or multiple mutations) are created [52].

Both networks are generated using our gene panel and adding specific single/double/multiple mutations involved in the disease in case of mutated networks. Mainly missense mutations are analyzed. The generation of the gene–gene network is very important because it allows to understand how the molecular profile of the network changes for the individual patient in the presence of one or more missense or other mutations.

In this way, we can analyze the joint effect derived from the combination of multiple mutations, creating the protein–protein interaction network using the gene–gene net as input. As a result, it is possible to identify which molecular mechanisms are changing and which could be the pharmacological therapies applicable for the individual patient.

Protein–Protein interaction network

Once the gene–gene networks are generated and analyzed, it is necessary to trace the molecular mechanisms involved both at physiological and at pathological levels, evaluating the joint effects deriving from the presence of one or more mutations. This is the result of the realization of the protein–protein networks as shown in step e of Fig. 1 [53]. These are networks in which some nodes represent proteins (involved either in control or in pathological conditions), while the others represent mechanisms (either physiological or pathological). The edges represent the molecular associations between protein–protein interactions (involved in the disease) and molecular mechanisms involved in this network. The input of the network is represented by the genes derived from the gene–gene network. The result is a network in which there are molecular mechanisms and proteins involved in the disease of interest. The applied algorithms are the same used in the case of gene–gene networks; the step related to compute the protein enrichment is added to evaluate all the intermediate proteins involved in the protein–protein interaction network, which helps us trace the molecular mechanisms involved in control and mutated networks.

Multilayer drugs network

Once all the main networks are created, it is possible to accomplish the multilayer drugs network as shown in step f of Fig. 1. It contains all the information derived from the networks built in the previous steps. It is a network in which some nodes represent genes, others represent various types of gene mutations (not only missense ones), and, finally, others represent the drugs approved by FDA or present in clinical trials (phase 1/2/3 of clinical experimentation) or in the preclinical phase. Furthermore, edges represent all the interactions between the drugs and the mutated genes and all the calculated pharmacological parameters. The input is represented by genes or mutated genes, while the output is a multilayer drugs network [54]. The user-specified network is constructed based on integrated bioactivity data, namely dose-response measurements (IC50, EC50, XC50, AC50, Kd, Ki, binding activity and potency) to extrapolate the most suitable therapeutic combinations or single therapies and to evaluate the side effects, adverse effects and contraindications of single or combined therapies for the individual patient.

Results and discussion

Case study development

As an application of this pipeline, the clinically relevant case studies belonging to a cohort of 182 metastatic colorectal cancer (mCRC) patients [55] were considered and here we show one of most interesting cases, which is related to a tumor carrying both the BRAFV600E and the TP53I195N missense variants. A more detailed discussion of this case study and of other five cases is reported in the Supplementary Data. The mCRC mutations we analyzed in this study are all pathogenic mutations [55, 80, 81]. In this case, the study of combined mutations is more useful than that performed for single mutations, as it allows us to study the ‘joint effect’ derived from the presence of two or more mutations.

In this way, it is possible to analyze other molecular mechanisms, identifying new potential pharmacological therapies, which are able to block or slow down the tumor growth [56-58]. A disease–disease network was built starting from mCRC to identify novel direct/indirect associations with other diseases and to confirm old ones. The parameters used for the realization of this network were focused on two aspects: (i) biological aspect to study the first biological interactors linked to the first and second level of bio-clinical analysis without taking into consideration the role played by the MSI genes; (ii) informatics aspect based on application of RWRH algorithm to help us to compute the similarities between diseases and to identify what are the main and strongest correlations within the network (a more detailed description of the parameters used and how they were applied has been reported in section 2 of the Supplementary Data). The results showed that there are some inflammatory diseases at intestinal level and other tumors that are directly or indirectly associated with CRC (Fig. 2).

Figure 2

Disease–disease network.

This information allowed us to identify the driver and priority genes involved in mCRC and to investigate their role played in mCRC realizing the gene-disease network (Fig. 3).

Figure 3

Gene–disease network.

Next, the study of a disease–disease network and of the gene–disease network associated with it allowed us to cross data obtained from both networks, and to accomplish the extrapolation of some important disease–disease and gene–disease associations (as shown in the Supplementary Data). This information was the basis for the realization and propagation of the gene–gene network and protein–protein interaction network.

Once the biological role played by both these hub genes is known, it is possible to investigate how their molecular profile changes due to BRAFV600E-TP53I195N combined mutation and what are the joint effects of a change at the molecular level, through the molecular analysis of the gene-gene network (reported in the Supplementary Data, Fig. 4a-b;) and the protein–protein interaction network (reported in the Supplementary Data, Fig. 5a-b).

Four gene–gene networks have been created. The gene–gene control network on healthy subjects (not affected by CRC) (Fig. 4a, shown in the Supplementary Data) was built by using the genes of the specific panel for this study. Its aim is to identify (a) the gene–gene interactions involved in a healthy subject and (b) the main sub-networks around which the entire network was organized. Then, they were used as inputs for the realization of the control protein–protein interaction network (Fig.5a, shown in the Supplementary Data), which allowed us to trace and explore the molecular mechanisms that are normally triggered in humans.

The second gene–gene network (Fig. 4b, shown in the Supplementary Data) was created taking into consideration a patient (extracted from the group of 182 patients mentioned above) with BRAFV600E-TP53I195N combined mutation, to verify how the molecular profile (Fig. 5b, shown in the Supplementary Data) of the network changed by using this approach. The other two gene–gene networks were related to single mutations BRAFV600E and TP53I195N. As in the control, the gene–gene networks were used as input for the realization of the three protein–protein interaction networks, which allowed us to trace and explore the molecular mechanisms in all the three cases.

The innovative result was that the networks obtained from the combined mutations exhibit a different behavior with respect to the two single mutation networks, not only for the presence of the first and second gene interactors and for the molecular mechanisms involved, but also for the therapeutic choices. In fact, by crossing the data obtained from the study of the mutated gene–gene network and the mutated protein–protein interaction network, it was possible to examine the joint effects derived from BRAFV600E-TP53I195N combined mutation within the network.

Therefore, it was seen that the joint effects did not derive from the sum of the effects of the individual mutations, but there was a combination of triggered events that imply strong molecular changes inside the network. This study allowed us to extrapolate for single and combined mutations the lists of the top ranked genes and molecular mechanisms involved in the disease, by calculating the EdgeBetweenness parameter (Table 2; more details are shown in Supplementary Data), to get as much information as possible to help us understand the most suitable therapy and the organization of the multilayer network.

Table 2

Top ranked priority genes with their molecular mechanisms and EdgeBetweenness parameter calculated on priority genes. Priority genes that undergo changes at the molecular profile level when we have BRAF-TP53 mutations are reported in the following list

First Mutation	Second Mutation	Priority-driver Genes: First Mutation	Edge-Betweenness First Mutation	Super-pathways First Mutation	Priority-driver Genes: Second Mutation	Edge-Betweenness Second Mutation	Super-pathways Second Mutation	Combined Mutation	Priority-driver Genes: Combined Mutation	Edge-Betweenness Combined Mutation	Super-Pathways Combined Mutation
BRAF V600E	TP53 I195N	BRAF	1.00					BRAFV600E-TP53I195N	EGFR	11.62
		NRAS	4.15				cell cycle arrest,		CHEK2	5.82
		PIK3CA	8.32	RAF/MAP kinase cascade,			apoptosis,		PARP1	2.30	Apoptosis pathways,
		MAP2K1	3.70	MAPK1/3 signaling cascade,	NEUROG1	12.53	senescence,		ERBB4	6.77	cellular differentiation signaling,
		CDKN2A	10.50	Signaling by high-kinase activity BRAF mutants,	HRAS	4.36	DNA repair and changes in metabolism,		P38	8.18	survival signaling, proliferation pathways,
		PRKACA	1.00	FLT3 Signaling,	CDKN2A	4.95	p38 signaling,		NPH1-FLT3	2.00	MAP/ERK kinase signaling pathway,
		HRAS	13.30	ERK signaling pathway,	KRAS	3.58	AMPK signaling,				G1/S cell cycle phase transition,
		TGFB1	1.00	Signaling to ERKs,	ERBB2	5.42	G1/S cell cycle phase transition,		KIT	2.00	G2/M cell cycle phase transition,
		MEK1/2	1.00	Oncogene-induced senescence,	KAT6A	4.62	G2/M cell cycle phase transition,		IGF2	10.45	RTKs signaling,
		KGFR	10.60	cell division,	PARP1	9.78	GRB2 events in ERBB2 signaling,		KAT6A	3.00	KIT signaling,
		MAPK3	7.33	differentiation and secretion,	AMPK	5.30	Regulation of TP53 Activity through Phosphorylation,		PDGFRA	9.85	p38 signaling,
		MAPK1 (KSR1)	6.60	EGFR signaling,	STK11	7.22	Activation of PPARGC1A (PGC-1alpha) by phosphorylation		CHD8	20.10	EGFR signaling,
		FGFR3	1.00	IL6 and KIT signaling,	CHEK2	16.28			SMAD2	12.01	WNT signaling,
		PDGFR	1.00	ERBB signaling,	MAPK14	9.11			WNT	30.47	GSK3 signaling,
		FAK	0.80	RTKs signaling,	FBXO11	1.10			GSK/APC	13.25	angiogenesis and invasion signaling,
		VEGFR	0.90	angiogenesis-signaling pathways	P38(ALK)	9.11			JAK	13.43
									PDGFR	9.85	ERBB2 Regulates Cell Motility
		JAK	0.80						VEGFR	6.00

By crossing the data extrapolated from the gene–gene network and from propagation phase of the protein–protein interaction network, it was possible to demonstrate that in presence of the combined mutation BRAFV600E-TP53I195N, the main molecular mechanisms altered are those that regulate the process of differentiation, survival, DNA damages and apoptosis, as shown in Table 2.

Moreover, gene–gene networks are all directed graphs and so it is possible to understand the direction followed by the priority genes within each network. In this way, the last multilayer drugs network was created and only the FDA-approved drugs as well as drugs in Phase 3 clinical trials for single and combined mutation were extrapolated. This network, generated by using the combined mutation BRAFV600E-TP53I195N, is able to identify the most suitable therapy based on some pharmacological indicators as reported in the general case in section 2.7 of the Methods. As a result of this process, we obtained a table showing the FDA-approved and Phase 3 drugs with their molecular mechanisms and drug resistance parameter calculated for single-combined mutations (Table 3; a detailed discussion is shown in Supplementary Data).

Table 3

Multilayer drug network classification.The table shows FDA-approved and Phase 3 drugs applied for BRAFV600E-TP53I195N in single and combined mutation

Cases n.	First Mutation	Second Mutation	Top of drugs classified for clinical phase first mutation	Response:Resistant or Sensitive for First Mutation	Top of Drugs classified Phase Second Mutation	Response:Resistant or Sensitive for Second Mutation	Combined Mutations	Top of Drugs classified For Clinical Phase combined Mutation	Response:Resistant or Sensitive for Combined Mutation
5	BRAF V600E	TP53 I195N	FDA-approved: Binimetinib+Encorafenib+Cetuximab,	FDA-approved: Sensitive	FDA-approved: Azacitidine	FDA-approved: Sensitive	BRAFV600E_ TP53I195N	FDA-approved: Regorafenib,	FDA-approved: sensitive
			(i.MEK+i.BRAF+i.EGFR),	Resistant				Vemurafenib,	sensitive
			Cetuximab (i.EGFR),	Resistant				Sorafenib,	sensitive
			Panitumumab	Resistant				Methotrexate,	Sensitive
			(i.EGFR),	Sensitive				Azacitidine	Sensitive
			Vemurafenib (i.BRAF),	Sensitive
			Dabrafenib (i.BRAF),	Sensitive
			Regorafenib (i.BRAF),	Sensitive
			Bevacizumab (i.BRAF),	Sensitive
			Sorafenib	Sensitive
			(i.BRAF/i.PDGFR/i.VEGFR),	Phase 3:
			Pazopanib (i.BRAF/i.VEGFR)	Sensitive
			Phase 3:	Sensitive
			Masitinib (i.PDGFR/i.FGFR3/i.FAK),	Sensitive
			Motesanib (i.PDGFR),
			Trametinib (i.MEK)

The analysis of the multilayer network highlighted some important differences between the drugs indicated for single mutations (BRAF and TP53) and the drugs suggested for the combined mutation (BRAF+TP53). For the single mutation BRAF V600E, the best treatment is the triplet drugs combination of Binimetinib (MEK) plus Encorafenib (BRAF) plus Cetuximab (EGFR) giving rise to a synergistic effect that increases overall survival in mCRC patients [59-61]. In contrast, if there is a combined mutation (i.e. BRAF V600E+TP53 I195N), it is suggested to use cytotoxic drugs as Methotrexate (MTX).

The rationale of such different therapeutic approach might reside in the complex interactions between TP53 mutations, BRAF mutations and sensitivity to MTX. MTX inhibits the synthesis of DNA, RNA, thymidylates and proteins because it is a potent inhibitor of the dihydrofolate reductase (DHFR) enzyme, which converts dihydrofolate to tetrahydrofolate (THF) [62]. TP53 mutations induce an increase of E2F1 and are associated with DHFR gene amplification to start the invasion and metastasis process [63]. These mechanisms decrease the ability of tumor cells to be responsive to MTX or other antifolate drugs involved in intracellular folate metabolism, DNA synthesis and cell growth [64, 65].

BRAF encodes a serine–threonine kinase that acts upstream of MAPKK1 and MAPKK2 in response to RAS signals, and is an essential component of the RAF/MEK/ERK/MAPK kinase cascade as ‘driver’ of tumorigenesis [66]. The RAS/BRAF/MEK/ERK cascade reduces E2F1 expression level, lowering the levels of DHFR and thus increasing the sensitivity to MTX [67]. Additional mechanisms involve the activation of WNT signaling that induces a higher production of and dependency from reactive oxygen species (ROS) as well as the interaction with pathways implicated in DNA repair. In fact, MTX can trigger ROS-associated cell apoptosis [68] The anti-inflammatory actions of MTX are critically dependent upon the production of ROS. Finally, a DNA damaging activity of MTX has been described [69, 70].

BRAFV600E-TP53I195N case study was reported only as an example, as we focused on targeted therapy that could be implemented in patients with CRC derived from the BRAF mutation. From our studies we have shown that, when the combined BRAF-TP53 mutation is present, it is possible to use MTX as an alternative therapy, while when we have the combination of BRAF mutation with other genes such as PIK3CA, we see that MTX is no longer present in the therapeutic picture. This combination suggests that MTX is suitable only in patients mutated in both BRAF and TP53 genes. This knowledge indicates that using our pipeline it is possible to make new alternative hypotheses of a therapy for patients with a complex mutational status, which must be validated in preclinical/clinical studies. A detailed explanation of biological and molecular processes has been reported in the added section 3.2 of the Supplementary Data. To identify the best therapy it has been necessary to run the whole pipeline, because the information derived from both the gene–gene network and the protein–protein interaction network is needed. Two lines of drugs are obtained: the first line of drugs (directed) acting on our mutation and the second line (indirected) of drugs, which instead acts mainly on other mutations of our genes, but which could also be used as support to the therapy if the first line does not work enough well.

Conclusions and future work

Assessment of genetic mutations is an essential element in the modern era of personalized cancer treatment. The mutations in some important hubs or priority genes as BRAF and TP53 play a critical role in regulating many important molecular processes.

In the literature, the study is focused on the analysis of single mutations, while by using our pipeline based on the multiple network approach [14], it was possible to demonstrate and validate that the joint effects and changes of the molecular profile that occur in a patient with mCRC with double mutation of the BRAF-TP53 genes are able to identify the most suitable drugs for the therapy. Through this approach, it is possible to cross-validate data and extrapolate as much information as possible, trying to address the physician to the best therapy for the individual patient. In this way, on one side, improvements could be made at the diagnostic level and, on the other hand, personalized therapy could be accomplished, reducing the research time and costs.

Furthermore, following this approach, it is possible to evaluate what are the molecular changes of a network in the presence of one or more mutations and what it implies at the pharmacological level. Copy-number variation (CNV) and focal amplifications are analyzed in clinical diagnostics, but this type of information is very limited because it depends on the request from the treating physician. In our series this information was not available. However, we could include this step as a future prospective because these data might be introduced in the pipeline by modifying the last three steps, by integrating some Cytoscape plugins (i.e. FunMod [71], ClueGO [72], CluePedia [73], Enrichment Map [74] and STRING [75]). Then, the obtained information should be integrated with the mutational data, to realize a new gene–gene network, and, consequently, the protein–protein and multilayer drugs networks. In a future work, we will add the study of gene expression profile, to further validate the results obtained and to have a more complete and detailed analysis. To include the study of the gene expression profiles, we will modify our pipeline at the level of the gene–gene network. In this way, it will be possible to work both on the mutations and on the expression profile, reanalyzing the protein–protein interaction network and performing a more complete study at the drug level, considering both the mutations and the expression profiles. The Cytoscape plugins to apply for the analysis of gene expression profile of each patient are: InsituNet [76], GeneMANIA [25], HyperModules [77], NetworkAnalyst 3.0 [78] and ExAtlas [79]. Moreover, this procedure can be applied to any type of disease and not only to oncological ones.

The tumors often carry mutations of multiple driver genes. Multiple network analysis is very useful to identify and investigate joint effects derived by combined mutations involved in diseases of interest.

The application of Multiple network-based pipeline is very important for clinical researchers to find new pharmacological therapies for sick patients in few time and low cost.

The proposed procedure is useful for improving the diagnosis identifying the patients that could better respond to a given drug treatment and predicting what are the primary and secondary effects of gene mutations involved in the human disease of interest.

Click here for additional data file.