A Systems Biology Approach Provides Deeper Insights into Differentially Expressed Genes in Taxane-Anthracycline Chemoresistant and Non-Resistant Breast Cancers

Objective: To date, numerous studies have been conducted to search for reasons for chemoresistance and differences in survival rates of patients receiving chemotherapy. We have sought to identify differentially expressed genes (DEGs) between predicted chemotherapy resistance and sensitive phenotypes by a network as well as gene enrichment approach.
Methods: Functional modules were explored with network analysis of DEGs in predicted neoadjuvant taxane-anthracycline resistance versus sensitive cases in the GSE25066 dataset, including 508 samples. A linear model was created by limma package in R to establish DEGs.
Results: A gene set related to phagocytic vesicle membrane was found to be up-regulated in chemoresistance samples. Also, we found GO_CYTOKINE_ACTIVITY and GO_GROWTH_FACTOR BINDING to be up-regulated gene sets with the chemoresistance phenotype. Growth factors and cytokines are two groups of agents that induce the immune system to recruit APCs and promote tolerogenic phagocytosis. Some hub nodes like S100A8 were found to be important in the chemoresistant tumor cell network with associated high rank genes in GSEA. Conclusions: Functional gene sets and hub nodes could be considered as potential treatment targets. Moreover, by screening and enrichment analysis of a chemoresistance network, ligands and chemical agents have been found that could modify significant gene sets like the phagocytic vesicle membrane functional gene set as a key to chemoresistance. They could also impact on down- or up-regulated hub nodes. Creative Commons Attribution License

Introduction

Neoadjuvant therapy of breast cancer is a safe and effective therapeutic approach to reduce the risk of recurrence and mortality. Although in some case, application of chemotherapy may not be effective or even lead to recurrence following treatment or treatment failures, but it is a routine option for treating breast cancer (Eatemadi et al., 2016a; Tabatabaei Mirakabad et al., 2016). Neoadjuvant therapy will apply to achieve two important goals i) reduce the size of unresectable tumor which allowing surgery to be performed and ii) in operable tumors it helps for greater conservation of the breast and decrease need for mastectomy (Thompson and Moulder-Thompson, 2012; Rami and Zarghami, 2013). Over the past few years, numbers of chemotherapy regimens have used to treatment and assessment the usefulness and mechanism of action these regimens (Ghalhar et al., 2014; Eatemadi et al., 2016b). Anthracycline-based regimens increase treatment benefits in comparison with Cyclophosphamide, Methotrexate, and Fluorouracil (CMF) combinations (Bines et al., 2014). The recent reports of incorporation the Taxane with anthracycline-based regimens and comparison of survival rates showed significant improvement in patient outcome, De et al., (2008) study on 23,000 women from 13 clinical trials showed that taxane-anthracycline chemotherapy improved DFS (Distance free survival) and OS (Overall survival) in high-risk and early-stage breast cancer patients (De Laurentiis et al., 2008). Increasing in applying anthracyclines and taxanes as a treatment option for early stage breast cancers, lead to resistance and failure in the mechanisms of action of these agents in patients or sub-populations of tumor lesions. Therefore, consideration of complications and high costs of these kinds of treatments make it necessary to select qualified case to administrate these drugs and also attention has to be paid to efficacy of applying variety combinations, optimization of dosage and different sequences of similar combinations of administrated drugs (Moreno-Aspitia and Perez, 2009; Farajzadeh et al., 2017; Maasomi et al., 2017).

Recently, gene expression profiling has become useful as a robust genomic tool that used in screening the details of cells and applied to several field of biology e.g. clarification and understanding heterogeneous feedbacks to a special drug and treatment and more recently used to drive predicting signatures that have prognostic power for survival, the effectiveness of the drugs for special type of cancer or patients (Chang et al., 2003). Hatzis et al., (2011) tried to develop a predictor to define response and survival from chemotherapy in newly diagnosed invasive breast cancer. They studied two datasets, the first one that used for developing the predictor including 310 patients and the second population used to test the predictor on an independent group including 198 patients. They developed a genomic predictor which predicted chemoresistance, chemosensitivity and predicted endocrine sensitivity identified patients with high probability of survival following taxane and anthracycline chemotherapy. In this study, we aim to explore the mechanisms underlying in chemoresistance and chemosensitivity (Ch-R vs Ch-S) via PPI network approach and enrichment analysis of DEGs to find functional gene sets related to Ch-R/Ch-S driver genes.

Materials and Methods

Data preprocessing

Gene expression profile of GSE25066 including 508 samples downloaded from InSilicoDB, a genomics data repository. The GSE25066 expression profile is based on the GPL96 (Affymetrix Human Genome U133A Array) that preprocessed by FRMA method (McCall et al., 2010).

DEG analysis of predicted Ch-R vs Ch-S samples

A linear model was created with limma package in R (Ritchie et al., 2015) to finding DEGs between samples that have been recognized as chemotherapy resistance and sensitive samples in the study.

Gene enrichment analysis

A dataset consisting more reproducible genes with adjust p-value < 0.05 from limma output selected from the original dataset, to analyze by GSEA desktop software (Subramanian et al., 2005). We select default setting of GSEA. Permutation type set to gene set to find the most affected pathways and functional gene sets in Ch-R vs Ch-S phenotypes.

Pathway analysis

Limma tap-rank genes with log fold change (lfc) > 0.5 and adjust p-value < 0.05 selected to analyze with Ipathway, a web-based service that computes pathways perturbation vs over-representation (Draghici et al., 2007).

Network analysis

Constructing networks were done by STRING and NetworkAnalyst that have different options and utility to study the networks. The network constructed by STRING from significant DEGs with adjusted p-value < 0.05 and absolute lfc > 0.5 that have PPI pairs whose protein interaction scores were > 0.4. The network constructed by NetworkAnalyst based on mentioned conditions along with lfc amounts. We used Walk Trap algorithm from NetworkAnalyst to finding subnetworks in network that created by NetworkAnalyst (Szklarczyk et al., 2014).

Results

Differential expression gene analysis

Differential expression gene analysis between predicted Ch-R and Ch-S samples resulted in 2,863 differentially expressed prob-IDs with adjust p-value less than 0.05.

Gene set Enrichment Analysis of Ch-R vs Ch-S DEGs

More reproducible prob-IDs including 2863 prob-IDs related with 2328 DEGs extracted from the original dataset. GSEA was done with selected features and hallmark gene set in addition of Gene Ontology gene sets, including Molecular Function, Cellular Component, and Biological Process gene sets from MSigDB (Liberzon et al., 2011). The most significantly enriched gene sets for the chemoresistance phenotype were the innate immune response, cytokine activity, Phagocytic vesicle membrane and allograft rejection with BP, MF, CC and Hallmark gene sets respectively (Figure 1; Table 1). Enrichment results for the chemosensitive phenotype enriched in epithelium development, growth factor binding, extracellular matrix and estrogen receptor early as the most significant gene sets from BP, MF, CC and hallmark gene sets respectively (Figure 1; Table 2).

Figure 1

GSEA Plots of Ch-R vs Ch-S DEGs. (A-D) GSEA plots of chemoresistance phenotype enriched gene set with BP, MF, CC and hallmark gene sets. (E-F) GSEA plots of chemosensitive phenotype enriched gene set with BP, MF, CC and hallmark gene sets. The score at the peak of the plot shows enrichment score.

Table 1

Three Most Significant Enriched Gene Sets for Chemoresistance Phenotype Features from BP, MF, CC and Hallmark Gene Sets

	NAME	NES	FDR q-val
BP	GO_INNATE_IMMUNE_RESPONSE	-3.9926152	0
	GO_IMMUNE_RESPONSE	-3.7706816	0
	GO_DEFENSE_RESPONSE	-3.6039486	0
MF	GO_CYTOKINE_ACTIVITY	-2.4459608	0
	GO_SIGNAL_TRANSDUCER_ACTIVITY	-2.352623	0
	GO_SIGNALING_RECEPTOR_ACTIVITY	-2.350062	0
CC	GO_PHAGOCYTIC_VESICLE_MEMBRANE	-2.4337091	0
	GO_ENDOCYTIC_VESICLE_MEMBRANE	-2.3747494	0
	GO_SIDE_OF_MEMBRANE	-2.3073456	0
Hallmarks	HALLMARK_ALLOGRAFT_REJECTION	-3.799198	0
	HALLMARK_INTERFERON_GAMMA_RESPONSE	-3.7596009	0
	HALLMARK_INTERFERON_ALPHA_RESPONSE	-3.361638	0

Table 2

Three Most Significant Enriched Gene Sets for Chemosensitive Phenotype Features from BP, MF, CC and Hallmark Gene Sets

	NAME	NES	FDR q-val
BP	GO_EPITHELIUM_DEVELOPMENT	2.5826445	0
	GO_EPITHELIAL_CELL_DIFFERENTIATION	2.465178	0
	GO_EPITHELIAL_CELL_DEVELOPMENT	2.3974626	3.35E-04
MF	GO_GROWTH_FACTOR_BINDING	2.2345507	0.002986222
	GO_CORE_PROMOTER_BINDING	2.1788	0.005907953
	GO_TRANSCRIPTIONAL_ACTIV
	ATOR_ACTIVITY_RNA_POLYMERASE_II	2.1535537	0.005258679
	_CORE_PROMOTER_PROXIMAL_REGION
	_SEQUENCE_SPECIFIC_BINDING
CC	GO_EXTRACELLULAR_MATRIX	2.163985	0.011273597
	GO_APICAL_PART_OF_CELL	2.1007788	0.014542908
	GO_RIBONUCLEOPROTEIN_COMPLEX	2.0989242	0.009695271
Hallmarks	HALLMARK_ESTROGEN_RESPONSE_EARLY	3.5579836	0
	HALLMARK_ESTROGEN_RESPONSE_LATE	3.4210422	0
	HALLMARK_ANDROGEN_RESPONSE	2.2230856	0

Pathway analysis of DEGs

The limma output was load into Ipathway to filter DEGs with adjust p-value < 0.05 and lfc > 0.5. From the filtering stage of Ipathway, 51 genes selected (Figure 2A) to analyze via Ipathway to explore for the most impacted pathways between the Ch-R/Ch-S as different phenotypes (Figure 2B and 2C). The results from pathway analysis show Chemical carcinogenesis, KEGG: 05204 (P= 0.021), Complement and coagulation cascades, KEGG: 04610 (P= 0.026) and Prolactin signaling pathway, KEGG: 04917, (P= 0.046) as the most significant impacted pathways (Figure 2).

Figure 2

(A) Volcano plot of filtered DEGs by Ipathway that represents significant DEGs with adjust p-value< 0.05 and log fold change > 0.5. Genes on the right show up-regulated genes in chemosensitive and genes on the right represent up-regulated genes in chemoresistance phenotype. (B) Summary of significant impacted pathways related with Ch-R vs Ch-S DEGs. (C) Chemical carcinogenesis (KEGG: 05204) pathway that impacted with changing in expression of genes highlighted in red, including common genes between Chemical carcinogenesis pathway and anthracyclines and taxanes induced gene.

Constructing of networks were done by mapping 54 genes with adjust p-value < 0.05 and lfc > 0.5 using STRING. Resulting network includes 54 nodes and 65 edges in comparison with 12 edges as the expected number of edges for the same number of random nodes that imply this network has significantly more connections. Average node degree is 2.41 and clustering coefficient (CC) is equal to 0.777 (Figure 3A). Screening for the functional modules in this network results in negative regulation of signal transduction (GO: Biological Process), protein binding (GO: Molecular Function) and extracellular space (GO: Molecular Function) as the most significant enriched modules from GO analysis of involved proteins (Table 3). The network from NetworkAnalyst constructed by 54 proteins as seed nodes. Obtained network from NetworkAnalyst includes 1,620 nodes and 2,111 edges (Figure 3B). Three top modules that found using Walk Trap algorithm named subnetwork1 (p= 3.52e-23, size=163), subnetwork2 (p= 6.62e-239, size= 620), and subnetwork3 (p= 2.28e-20, size= 136), (Figure 4). KEGG pathway enrichment analysis of genes that involved in each subnetwork represent in Table 4.

Figure 3

PPI Networks of DEGs in Ch-R vs Ch-S. The nodes indicate the DEGs and the edges indicate the interactions between two genes. In the network in right side of the figure red nods expressed higher in the sensitive samples and the green node represent nods expressed higher in the insensitive samples.

Table 3

GO Analysis of Differentially Expressed Genes in the Protein-Protein Interaction Network

pathway ID	pathway description	count in network	false discovery rate
Biological Process (GO)
GO:0009968	negative regulation of signal transduction	14	0.000501
GO:0030155	regulation of cell adhesion	11	0.000501
GO:0045785	positive regulation of cell adhesion	9	0.000501
GO:0065008	regulation of biological quality	23	0.000501
GO:0042127	regulation of cell proliferation	16	0.000592
GO:0042592	homeostatic process	15	0.000619
GO:0048585	negative regulation of response to stimulus	15	0.000619
GO:0048519	negative regulation of biological process	26	0.000908
GO:0048523	negative regulation of cellular process	25	0.000908
Molecular Function (GO)
GO:0005515	protein binding	27	0.00817
GO:0001948	glycoprotein binding	4	0.0336
GO:0035662	Toll-like receptor 4 binding	2	0.0336
GO:0019838	growth factor binding	4	0.0466
GO:0050544	arachidonic acid binding	2	0.0466
Cellular Component (GO)
GO:0005615	extracellular space	20	1.41E-08
GO:0044421	extracellular region part	31	1.66E-08
GO:0005576	extracellular region	32	1.68E-07
GO:0070062	extracellular exosome	24	2.78E-06
GO:0031988	membrane-bounded vesicle	23	0.000537
GO:0009986	cell surface	9	0.00892
GO:0009897	external side of plasma membrane	5	0.0362
GO:0071944	cell periphery	23	0.0362

Figure 4

Top Three Modules in the PPI Network of the DEGs in Ch-R/Ch-S.

Table 4

Three Top Pathways Obtained from KEGG Pathway Enrichment Analysis of Genes That Involved in Each Sub-Network

SubNetwork1		SubNetwork2		SubNetwork3
Pathway	p-value	Pathway	p-value	Pathway	p-value
ErbB signaling pathway	1.45E-18	Pathogenic Escherichia coli	0.0000165	Cell cycle	6.09E-16
Jak-STAT signaling pathway	3.44E-17	Regulation of actin cytoskeleton	0.0000173	Chronic myeloid leukemia	1.01E-13
Chemokine signaling pathway	1.99E-14	Glycolysis / Gluconeogenesis	0.0000251	Pathways in cancer	6.89E-11

The top ten highest degree nodes in NetworkAnalyst are ESR1, CCND1, PTN, MYB, S100A9, ERBB3, MAPT, IGF1R, KRT15, NRIP1 that represent the most important hubs in Ch-R/Ch-S network. From the above list, the S100A9 was up-regulated hub node in chemoresistance tumor cell and the others are important hub nodes that up-regulated in chemosensitive tumor cells, (Figure 3B). Since there are accumulating evidences that imply the role of the phagocytic pathway in chemoresistance, significant genes related to phagocytic vesicle membrane gene set from GSEA results, extract and mapped on STRING, PPI database (Figure 5C). The network created with significant genes in vesicle membrane gene set include 13 nodes and 8 edges in comparison with 1 edge that was expected for this network with random nodes and the same size. Average node degree and CC for this network are 1.23 and 0.872 respectively (Figure 5C). Since all genes related to this pathway in Ch-R/ Ch-S DEG test, were up-regulated in chemoresistance phenotype, enrichment analysis of these genes was done using Enrich R (Chen et al., 2013) along with chemical agents and ligands perturb down, that defined potentially agents lead to perturb down genes that were up-regulated in phagocytic vesicle membrane significant genes (Figure 5; Table 5).

Figure 5

Clustergrams of Enrichment Analysis of Significantly Enriched Phagocytic Vesicle Membrane Gene Set with Chemical Agents (A) and ligands (B) lead to perturbation down, genes that up-regulate in phagocytic vesicle membrane significant genes in chemoresistance breast cancer. (C, D) Clustergrams of enrichment analysis of hub nodes in chemoresistance network with chemical agents and ligands that lead to perturbation down, hubs that up-regulate in chemoresistance, or cells. (E) PPI network of significant genes related phagocytic vesicle membrane gene module.

Table 5

Enrichment Analysis Results of Phagocytic Vesicle Membrane Significant Gene Set and Hub Nodes in Chemoresistance Network with Chemical and Ligand Database in from Enrichr

Term	P-value	Adjusted P-value	Z-score
LINCS L1000 Chem Perturbations down phagocytic vesicle membrane significant gene set
LJP006_HME1_3H-celastrol-1.11	3.31E-07	0.000459075	-1.937454563
LJP006_SKBR3_24H-alvocidib-0.04	1.09E-07	0.000300973	-1.796769478
LJP006_SKBR3_24H-HG-6-64-01-3.33	2.12E-06	0.001957743	-2.014133824
LINCS L1000 Ligand Perturbations down phagocytic vesicle membrane significant gene set
MCSF-MCF7	0.000348385	0.003135466	-1.641055509
BNGF-BT20	0.000324264	0.003135466	-1.633430948
EGF-MDAMB231	0.004760296	0.028561776	-1.500470998
LINCS L1000 Chem Perturbations down significant hub nodes in chemoresistance network
LJP005_PC3_24H-CP466722-3.33	6.46E-07	0.000415406	-2.28816605
LJP005_MCF7_24H-linsitinib-1.11	3.18E-07	0.000415406	-2.25671746
LJP008_HA1E_24H-GSK-1059615-0.04	1.14E-06	0.000415406	-2.22757269
LINCS L1000 Ligand Perturbations down significant hub nodes in chemoresistance network
KGF-MCF7	0.000827715	0.013136852	-1.975009722
BTC-HS578T	0.000799327	0.013136852	-1.920521786
INS-SKBR3	0.000856751	0.013136852	-1.797721404

Discussion

Due to the high probability of recurrence of breast cancer after chemotherapy, understanding the mechanisms by which chemoresistance can occur is important to developing novel therapeutic approaches and finding potential therapeutic targets. Anthracyclines and Taxanes are the most widely used chemotherapeutic drugs for breast cancer. Anthracyclines acts through inhibition of DNA and RNA synthesis, Inhibition of topoisomerase II enzyme, and a variety of mechanism underlying DNA damage (Offermanns and Rosenthal, 2008). On the other hand, Taxanes disrupt microtubules formation and function and perturb the mitotic organization and cell division (Rowinsky, 1997). There are numerous evidences for the roles of immune system, especially innate immune systems in chemoresistance activity of cancer cells that is in concordance with GSEA results with GO biological process in this study, e.g. enrichment of chemoresistant driver genes significantly represent GO_INNATE_IMMUNE_RESPONSE, GO_RESPONSE_TO_INTERFERON_GAMMA, and some other gene sets related to immune systems (Schmidt et al., 2008; McDonnell et al., 2011). Also, GSEA with GO molecular function gene sets represents GO_CYTOKINE_ACTIVITY as the most significant gene set with chemoresistance DEGs. Furthermore, enrichment result of chemoresistance DEGs with hallmark gene sets represents HALLMARK_ALLOGRAFT_REJECTION that imply immune systems role in chemoresistance activity of cancer cells. On the other hand, accumulating evidence has suggested a key role for the cancer stem cells that have high self-renewal activity and form a specific population that tends to be refractory to conventional treatments (Dean et al., 2005; Abdullah and Chow, 2013). One of stem cells features that involved in chemoresistance of cancer cells is involving ATP binding cassette (ABC) transporters (Ou and Guo, 2007). Screening in DEGs results, for this family of proteins showed that ABCA4 was up-regulated in predicted chemoresistant cells. It has been demonstrated that anthracycline immunogenicity relies on their ability to induce the tumor cells phagocytosis by dendritic cells (Kopecka et al., 2011). On the other hand, Phagocytosis of tumor cells often mediates the immunosuppressive properties of tumor-associated antigen-presenting cells leading to the resistance of tumors to chemotherapy (Jinushi, 2014). Screening in our results shows that the gene set related with Phagocytic vesicle membrane is up-regulated in chemoresistance samples. Searching for chemical agents and ligands that can adjust this gene set activity enriched in potentially chemical agents and ligands that could perturb-down the up-regulated gene set related with phagocytic vesicle membrane and modify this gene set related with anthracyclines and taxanes resistance.

Moreover we found that GO_CYTOKINE_ACTIVITY and GO_GROWTH_FACTOR BINDING, are up-regulated gene sets in chemoresistance phenotype. Growth factors and cytokines are two agents that induced the immune systems for recruiting APCs and promoting tolerogenic phagocytosis (Bines et al., 2014). Another example of factors that contribute in chemoresistance is the up-regulation of TLR-2, TLR-4, or TLR-7/8 that enhances tumor cell survival through the action of immune-mediated and cell autonomous NF-κB- and Bcl2-dependent antiapoptotic programs (Cherfils-Vicini et al., 2010). Screening in our result showed the TLR2 with lfc= 0.23 in predicted chemoresistance cells. In this study screening of taxane-anthracycline chemotherapy response, show common mechanisms with Chemical carcinogenesis pathways (KEGG: 05204), (Figure 2C) which have shared up-regulated genes including NAT1 (N-Acetyltransferase 1) and GST (Glutathione S-transferases) with taxane-anthracycline chemotherapy response that have higher expression in sensitive tumor cells. Chemical carcinogenesis significantly contributes to the causation of a sizable fraction of human cancers that acts through genotoxic and non-genotoxic mechanisms that imply i) attack to biological macromolecules such as DNAs and RNAs and ii) carcinogens act through the mechanisms for example induction of inflammation, immunosuppression, etc (Fine et al., 2010). Most important hub nodes that play as driver agents in chemoresistance phenotype are including S100A9, PTPRC, S100A8, SELL, SRGN. The S100A9 involved in the regulation of numbers of cellular processes including cell cycle progression and differentiation and has important roles in regulation of inflammatory processes, immune response and regulatory effect on toll-like receptor, etc. that is in concordance with previously mentioned enriched gene sets in GSEA results (Srikrishna, 2012). In particular, numerous studies represent roles of S100 family in progression after chemotherapy and chemoresistance. Zhou et al., (2016) found that patients with lymphoma that have higher expression of S100A9 displayed unfavourable treatment outcome, this result in addition to other referenced results endorsed the role of S100 family proteins in poor outcome of cancer patients after treatment (Zhou et al., 2016). Protein tyrosine phosphatase, receptor type C (PTPRC) involve in the variety of cellular processes including cell growth, immunomodulator and oncogenic transformation, etc (Tchilian and Beverley, 2006). However less attention has been paid to roles of this protein in breast cancer chemoresistance. Potentially increasing in expression of this protein in our network analysis results, in addition of immunomodulator and oncogenic transformation activity of this protein raise the important role of this protein in chemoresistance of breast cancer. S100A8 is another member of S100 family. Although there are numerous studies that referred role of S100A8 in progression of breast cancer and gaining invasive properties of breast cancer, but contribution of this protein in breast cancer chemoresistance has not attracted much attention (Moon et al., 2008; Yin et al., 2013) Although, Yang et al., (2016) and Yang et al., (2014) suggested association of up-regulation of S100A8 protein with drug resistance by promoting autophagy and influencing the apoptosis pathway in leukemia (Yang et al., 2014; Yang et al., 2016), but relation between taxane-anthracycline Chemotherapy in breast cancer and up-regulation of S100A8 is not clear. Our results represent S100A8 as one of the most important hub nodes in chemoresistant tumor cell network and high rank genes associated with chemoresistance in GSEA feature selection result. Selectin L (SELL) is another important hub node in chemoresistance network. According GeneCards, this protein encodes a cell surface adhesion molecule that belongs to a family of adhesion/homing receptors (Safran et al., 2010). Although this protein is a prominent node in chemoresistance network but it is not known as a key agent in chemoresistant breast cancer in scientific texts. Serglycin (SRGN) encodes a protein that best known as a hematopoietic cell granule proteoglycan and expressed in some tumor cells. It promotes metastasis and protects some tumor cells from complement system attack (Korpetinou et al., 2013).

There are some functional gene sets and hub nodes suggested as a treatment target in BC in this work. Also, there are several papers (Moreira et al., 2014; Hsu et al., 2015). Benzyl butyl phthalate increases the chemoresistance to doxorubicin/cyclophosphamide by increasing breast cancer-associated dendritic cell-derived) that point to some of these driver agents as important factors in cancer chemo-resistance and also there are some routine therapies for BC that these functional gene sets and hub nodes are involved in its targets, eg from hub nodes, ERBB3 is one of the hub nodes up regulated in sensitive breast cancer that, there are some drugs like Gefitinib, a kinase inhibitor drug and Iressa, which approved in cancer treatment. ESR1 is one of the hub nodes in CHR vs CHS network that is a target for Tamoxifen and Estrone drugs that approved too. From functional gene sets, innate immune gene set is one of the most significant gene sets and includes some genes like STAT1 that encodes a protein that is a target for Fludarabine, an approved drug to inhibit cell making. EGFR is one of the genes in HALLMARK_ALLOGRAFT_REJECTION, CELL-CYCLE, and some other significant gene set in this work, and is a target for Lapatinib, Gefitinib, Erlotinib, Cetuximab, Panitumumab, that could be inhibitors of this genes in chemo-resistance breast cancer cells (Safran et al., 2010).

Since the output of the gene expression analysis showed lots of significant genes and these genes with a high possibility could be find as a target in literature and databases, it has preferred to find these genes in a significant gene sets framework. Lots of functional significant gene sets of these work are subjects of routine therapies in breast cancer eg. Cell cycle, immune system and etc.

In conclusion, the resistance of breast cancer to anthracycline and taxanes raised the need of understanding underlying mechanisms in resistance to these drugs. Events responsible to chemoresistance have not been clearly identified, but a variety of mechanisms including contribution of cancer stem cells (CSCs), immune response, microenvironment, epigenetic changes or combination of these biological phenomena can be involved in changing drug sensitive cancer to insensitive tumor cells. In this study using systems approach many of functional gene sets and hub nodes highlighted mechanisms that can be potentially involved in chemoresistance.

Conflict of Interest

The authors declare no conflict of interest.