Integrated transcriptome interactome study of oncogenes and tumor suppressor genes in breast cancer.

Breast cancer is the leading cause for mortality among women worldwide. Dysregulation of oncogenes and tumor suppressor genes is the major reason for the cause of cancer. Understanding these genes will provide clues and insights about their regulatory mechanism and their interplay in cancer. In the present study, an attempt is made to compare the functional characteristics and interactions of oncogenes and tumor suppressor genes to understand their biological role. 431 breast cancer samples from seven publicly available microarray datasets were collected and analysed using GEO2R tool. The identified 416 differentially expressed genes were classified into five gene sets as oncogenes (OG), tumor suppressor genes (TSG), druggable genes, essential genes and other genes. The gene sets were subjected to various analysis such as enrichment analysis (viz., GO, Pathways, Diseases and Drugs), network analysis, calculation of mutation frequencies and Guanine-Cytosine (GC) content. From the results, it was observed that the OG were having high GC content as well as high interactions than TSG. Moreover, the OG are found to have frequent mutations than TSG. The enrichment analysis results suggest that the oncogenes are involved in positive regulation of cellular protein metabolic process, macromolecule biosynthetic process and majorly in cell cycle and focal adhesion pathway in cancer. It was also found that these oncogenes are involved in other diseases such as skin diseases and viral infections. Collagenase, paclitaxel and docetaxel are some of the drugs found to be enriched for oncogenes.

Introduction

Breast cancer is the most common cancer among women worldwide and about two million cases and six lakh deaths were reported this year by the World Health Organization (WHO).1, 2 Nearly ninety percentage of the women diagnosed with breast cancer at the earliest stage can survive up to several years, compared to women who are diagnosed at later stages. In this scenario, there is a requirement for the awareness among women and need for the early detection and diagnosis. Understanding cancer related genes could provide better insights of the pathogenesis and lead to the identification of targets for early detection. Breast cancer progression involves various genetic events – by activating the oncogenes or disrupting the regular functions of specific tumor suppressor genes. Over expression and amplification of oncogenes often assist in the development of cancer and cooperatively regulate genetic and epigenetic changes. A number of studies have reported that the oncogenes ErbB2, PI3KCA, MYC and CCND1 are often found to be deregulated in breast cancer. Among the oncogenes, HER-2 activation is found in about 20% of all primary breast cancer cases.4, 5 On the other hand, tumor suppressor genes are the negative regulators of cellular progression and growth which regulate the invasiveness and metastatic potential. Loss of function of these genes often lead to malignancy. Till date, there are several tumor suppressor genes reported to be associated with breast cancer - BRCA1, BRCA2, P53, PTEN, ATM and CHK2 are some among those. Recent advancements in high throughput technologies like DNA microarrays, next generation sequencing had produced massive data which are analysed with modern approach to gain better insights and also to identify novel targets.7, 8 In the present study, various bioinformatics analysis has been carried out in order to understand mechanism and interaction of oncogenes and tumor suppressor genes in breast cancer. With available breast cancer microarray data, differentially expressed genes (DEG) were identified and classified into five gene sets. Further these gene sets were subjected to various in silico analysis to understand their characteristic role and mechanism based on evidences from publicly available experimental data. From the findings, it was observed that the oncogenes were having high mutation frequency rate and also enriched with guanine-cytosine content than the tumor suppressor genes. It was also found that the oncogenes are highly interconnected and the genes CDK1, FOS, CCNA2, MMP9, CDH1, CCNB1 and TOP2A were identified to be hubs. This exploration of oncogenes and tumor suppressor genes in breast cancer could aid cancer biology research for early diagnosis and treatment options.

Materials and methods

Dataset collection

Breast cancer microarray datasets were collected from Gene Expression Omnibus (GEO), a public functional genomics data repository of NCBI. The criteria for selection of dataset is that, it must have samples of both healthy and breast cancer tissue with no drug treatment or any other illness. Comparison of normal versus cancer tissues will help in identification of genes that are deregulated.

Identification of differentially expressed genes

The collected microarray datasets10, 11, 12, 13, 14, 15 were analysed individually by comparing as groups of breast cancer tissue versus healthy tissue as controls using GEO2R, a web-based tool for gene expression analysis. The tool is based on R packages GEOquery and limma for calculation of p-value, logFC, adjusted p-value, t-statistic and B-statistic. Differentially expressed genes (DEG) were filtered with the cut-off of |logFC| ≥±2 with p-value <0.05 from each dataset. The genes which are found to be differentially expressed in more than one dataset was considered as breast cancer associated genes. The identified differentially expressed breast cancer genes were validated with available genes in the International Cancer Genome Consortium (ICGC) data portal.

Enrichment analysis

The enrichment analysis of the identified breast cancer associated genes were performed using WEB-based GEne SeT AnaLysis Toolkit (WebGestalt), an online software toolkit comprising information from various public resources for biological analysis. The enrichment analysis such as Gene Ontology (GO), pathways, diseases and drugs were carried out with top 10 results as significant using hypergeometric test and Benjamini & Hochberg method.

Gene set classification & analysis

The identified differentially expressed breast cancer genes were classified into five gene sets namely, oncogenes (OG), tumor suppressor genes (TSG), druggable genes (DG), essential genes (EG) and other genes (OtG). The classification of the gene sets was performed based on the mapping of DEG with various databases and resources. For OG, the collection of all oncogenes from Bushman lab (http://www.bushmanlab.org/links/genelists) was used and TSGene (https://bioinfo.uth.edu/TSGene/) database, a web resource for tumor suppressor genes was considered for TSG classification. Cancer Gene Census list from the Catalogue Of Somatic Mutations In Cancer - COSMIC (https://cancer.sanger.ac.uk/cosmic) was also used for the classification of TSG and OG. The Drug Gene Interaction database (DGIdb), an online database (http://www.dgidb.org/) of drug–gene interactions and druggable genome data was used for the classification of druggable genes. The essential genes of humans were collected from the database of essential genes (http://www.essentialgene.org/) and mapped to identify essential genes from DEG. The genes which did not map to any of these resources were considered as other genes (OtG). The classified five gene sets were subjected to enrichment analysis (KEGG pathways, GO, diseases and drugs) using WebGestalt to understand their biological role and properties.

GC content and mutation frequency

The Guanine-Cytosine (GC) content percentage for each of the five classified gene sets were calculated using bioMart tool (https://www.ensembl.org/biomart) from Ensembl. GC content is one of the fundamental features in a genome which is widely studied for methylation profiles and binding of transcription factors in regulating gene expression. The mutation frequency of each set of genes were calculated from The Cancer Genome Atlas (TCGA) breast cancer data. The mutation frequency is the ratio of samples where the gene is mutated among the entire samples sequenced for a gene.

Protein–protein interaction and cluster analysis

The protein–protein interaction (PPI) of the DEG was constructed using STRING (https://string-db.org/) database, an online biological database for known and predicted protein–protein interactions. The network of interacting proteins was downloaded and visualized using Cytoscape v3.5.1, an open source software tool for visualizing molecular interactions. The top 10 modules of highly interacting gene clusters among the DEG were found using MCODE plugin with default parameters. For the classified five gene sets, the protein–protein interaction network was constructed and the network topological parameters such as degree, betweenness centrality, shortest path and closeness centrality were calculated using NetworkAnalyzer in Cytoscape.

Results

Identification of differentially expressed genes (DEG)

Oncogenes and Tumor suppressor genes are the major key players in the dysregulation of cancer pathways. To study their properties and biological mechanism of action in breast cancer, seven microarray datasets were analysed and the identified differentially expressed genes were subjected to various bioinformatics analysis. The methodology adopted in the present study is depicted in Fig. 1. The detailed information about the collected seven microarray datasets is given in Table 1 along with the number of identified up and down regulated genes. The gene expression analysis results of the seven microarray datasets with up and down regulated genes is provided as supplementary data.

Figure 1

Overview of methodology adopted in the present study.

Table 1

Gene expression omnibus (GEO) datasets used in the present study with number of samples, platform information along with the number of identified differentially expressed genes (DEG).

Sl. No.	Dataset accession	No. of samples	Platform	No. of DEG
				Upregulated genes	Downregulated genes
1	GSE45584	90	GPL6480	98	124
2	GSE45581	45	GPL6480	707	1279
3	GSE21422	19	GPL570	1020	1155
4	GSE6883	24	GPL96	606	495
5	GSE79058	76	GPL19956	39	35
6	GSE45827	155	GPL570	2917	1118
7	GSE1299	22	GOL96	212	346

Enrichment analysis and classification of DEG

A total of 416 DEG were identified and subjected to enrichment analysis and the results are given in Fig. 2 and Table 2. From ICGC data portal it was observed that 99% of the genes were matched to breast cancer genes. The enrichment results include enriched GO, pathways, diseases and drugs with the number of genes involved. The GO enrichment result of the DEG is summarised in Fig. 2 which shows cell cycle, cell division, mitosis, cell cycle phase are the major biological processes these genes are involved and the molecular functions includes, protein binding, extracellular matrix structural constituent, ATP binding and microtubule motor activity. The DEG cellular component were mostly found to be in nucleus, extracellular matrix, kinetochore and chromosome. The overview of GO analysis results for the DEG with adjusted p-value is given in supplementary. The KEGG pathway enrichment results indicated that the DEG are involved in pathways such as cell cycle, ECM receptor interaction, focal adhesion, p53 signalling and small cell lung cancer pathway. The enriched disease terms for the DEG were cancer, neoplasm, breast diseases, carcinoma, adenocarcinoma, breast neoplasm and collagenase, paclitaxel, heparin, urokinase, progesterone, epirubicin, doxorubicin, alteplase, podofilox, zidovudine are the top 10 drugs found in drug enrichment analysis. The DEG were classified into five gene sets as 102 oncogenes (OG), 41 tumor suppressor genes (TSG), 77 druggable genes (DG), 243 essential genes (EG) and 112 other genes (OtG) with the help of existing resources discussed in the methods section. Some of the genes were found to be overlapped in multiple gene sets, the classified gene set list is given in supplementary. Further the classified five gene sets were subjected for enrichment analysis to understand their properties and role in cancer. It was not surprising that the same pathways, diseases and drugs of DEG were enriched for the gene sets as the genes had some overlaps. The enrichment results for the gene sets are given in supplementary.

Figure 2

Comprehensive bar chart of gene ontology (GO) results for differentially expressed genes showing number of genes in various biological process, cellular components and molecular function.

Table 2

Enrichment analysis result of pathways, diseases and drugs for the identified differentially expressed genes.

KEGG Pathway Enrichment		Disease Enrichment		Drug Enrichment
Pathways	No. of Genes	Diseases	No. of Genes	Drugs	No. of Genes
Cell cycle	27	Cancer or viral infections	86	Collagenase	26
ECM-receptor interaction	20	Neoplasms	72	Paclitaxel	15
Pathways in cancer	28	Breast Diseases	50	Heparin	17
Focal adhesion	23	Breast Neoplasms	51	Urokinase	11
p53 signalling pathway	12	Neoplastic Processes	52	Progesterone	12
Small cell lung cancer	12	Carcinoma	55	Epirubicin	7
Amoebiasis	13	Neoplasm Invasiveness	40	Doxorubicin	10
Progesterone-mediated oocyte maturation	11	Neoplasm Metastasis	40	Alteplase	9
Toll-like receptor signalling pathway	11	Skin and Connective Tissue Diseases	45	Podofilox	9
Bladder cancer	8	Adenocarcinoma	40	Zidovudine	6

The GC content percentages of all the classified gene sets were imported using bioMart tool and the average was considered for this study. The list of all genes and their GC percentages are given as supplementary. From the average percentage of all gene sets, it was observed that the druggable genes (DG) were having the highest GC content of 45.51% followed by oncogenes (OG) with 44.66%. The essential genes (EG) and tumor suppressor genes (TSG) were seen to have more or less similar GC content of 44.27% and 44% respectively. On comparison, the TSG is having genes with less GC content than the OG. The probability of change of a unit length of DNA with time is referred as rate of mutation. Of the five gene sets, OG was found to have the highest mutation frequency rate of about 1.19% and DG with 0.928%. The TSG and EG sets were having 0.746% and 0.753% average mutation frequency rates. The comparison of GC content and mutation frequency percentage for the five gene sets is given in Fig. 3.

Figure 3

Comparison of guanine-cytosine (GC) content percentage (A) and mutation frequency percentage (B) of the five gene sets.

Analysis of PPI network and clusters

The protein–protein interaction (PPI) network of DEG was retrieved and analysed by considering proteins as nodes and their interactions as edges. The constructed PPI network comprises 412 nodes and 2628 edges interactions based on experiments, co-expression, text mining, neighbourhood, gene fusion and databases. To understand the interplay, network analysis was performed by calculating the network topological parameters such as degree, betweenness centrality, closeness centrality and shortest path distance. Each of these parameters determine the role and property of gene/protein in the network. The average degree of the DEG network was 12.8, showing each gene on average has about 12.8 interactions with others in the network and the average clustering coefficient is 0.605. To know the highly interconnected genes in the network, the DEG were subjected to cluster analysis and the top 10 modules of clusters were taken into consideration, which is given in Table 3. The highest interacting cluster was found to have 37 nodes with 530 edges and the density score is 29.44. DLGAP5, AURKB, ECT2, CDCA5, RAD21, CASC5, ZWINT, KIF23, RACGAP1, CCNB2, MLF1IP, AURKA, TOP2A, CENPA, CDK1, NCAPG, NDC80, NEK2, KIF11, KIF4A, BUB1, PRC1, CDCA8, CCNA2, CDC20, BIRC5, BUB1B, MAD2L1, ZWILCH, CKAP5, KIF18A, KIF2C, CENPE, CENPF, CENPK, KIF20A, CCNB1 are the genes found in the top cluster with high interaction. It was observed that genes are majorly involved pathways related to cancer, cell cycle and p53 signalling. Furthermore, the PPI network for each of five classified gene sets were constructed and the network topological measures were calculated for the comparison of their network properties. The calculated network topologies of the five gene sets is provided in Table 4.

Table 3

Top 10 clusters of highly interconnected genes among the differentially expressed genes.

Cluster	Score (Density × No. of Nodes)	No. of Nodes	No. of Edges	Node IDs
1	29.444	37	530	DLGAP5, AURKB, ECT2, CDCA5, RAD21, CASC5, ZWINT, KIF23, RACGAP1, CCNB2, MLF1IP, AURKA, TOP2A, CENPA, CDK1, NCAPG, NDC80, NEK2, KIF11, KIF4A, BUB1, PRC1, CDCA8, CCNA2, CDC20, BIRC5, BUB1B, MAD2L1, ZWILCH, CKAP5, KIF18A, KIF2C, CENPE, CENPF, CENPK, KIF20A, CCNB1
2	11	13	66	ITGA6, COL11A1, ITGB6, COL5A2, COMP, COL1A2, COL12A1, COL5A1, COL4A6, ITGB4, COL10A1, COL3A1, ITGB1
3	6.1	21	61	FOXM1, ANLN, MCM4, PLK4, MELK, CDC45, CEP55, NUF2, CDC6, PCNA, CKS2, MCM2, TACC3, TTK, CCNE2, KIFC1, UBE2C, SMC4, PBK, RRM2, NUSAP1
4	6.08	26	76	POSTN, COL1A1, FPR3, RGS1, EZH2, MMP1, MYB, MYBL1, CXCR4, PIK3CA, TBL1XR1, THBS1, LAMA3, SPP1, LAMB3, CXCL11, CCR5, SERPINE1, CCL5, LAMC2, SDC1, FN1, CXCL9, TIMP3, RGS20, CXCL10
5	5.417	25	65	VEGFA, EIF5A, HMOX1, OAS2, CUL2, SQLE, RET, CTSS, IFI30, OASL, FGFR3, SRGN, MMP9, IFI6, OAS3, HAPLN1, VCAN, MAPK13, IRF6, FOS, MMP11, PLAUR, MMP3, MAX, STAT1
6	4	4	6	FEN1, EXO1, TRIP13, RAD51
7	4	4	6	HIST1H2BD, HIST1H3H, HIST1H2BH, HIST1H2BK
8	4	4	6	KYNU, KMO, QPRT, TDO2
9	3.333	4	5	DBF4, CDC7, CCNG2, CHEK1
10	3	3	3	ISG15, RSAD2, DDX58

Table 4

Calculated average network topological parameters for the five gene sets. OG: Oncogenes; TSG: Tumor Suppressor Genes; DG: Druggable Genes; EG: Essential Genes; OtG: Other Genes.

Network topological parameters	OG	TSG	DG	EG	OtG
Degree	11.1364	4.5334	7.1045	17.5	1.7074
Average Shortest Path Length	2.3984	2.4276	2.8151	2.6036	1.3701
Betweenness Centrality	0.0163	0.051	0.028	0.0077	0.1056
Closeness Centrality	0.4291	0.4272	0.3689	0.3963	0.798
Clustering Coefficient	0.444	0.2828	0.4779	0.4716	0.2651

Discussion

The exploration of relationship between OG, TSG and DG in breast cancer can provide basic understanding of the genetic and functional association of genes, which could aid in identification of disease markers. Breast cancer is a heterogeneous disease comprising various molecular interactions among different cell types. These molecular and cellular interactions are the major drive for the expression of cancer related genes and progression of the disease. Recent explosion of biological big data with high throughput technologies has evolved the understanding of tumor progression in cancer and mere understanding the role of cancer genes is inexpedient, stressing the need for integrated analysis.26, 27 The present study was aimed to understand the integrated network of interactions among the differentially expressed cancer related genes, especially in focus to oncogenes and tumor suppressor genes. For the integrated approach, various bioinformatics analysis such as GO, pathways, disease, drugs and networks were utilised. From the gene expression profiles, the oncogenes were found to be majorly deregulated along with the tumor suppressor genes. Further, expression analysis of the identified differentially expressed OG and TSG in TCGA portal for breast cancer showed that only a certain number of genes were found to have significant variations/mutations related to breast cancer. The enrichment analysis of these genes gave an overview of possible dysregulated pathways and associated diseases according to the genes. The role of OG and TSG in cancer and cell growth is mediated by various other factors, computing the GC content percentage for these genes will extend knowledge about the physicochemical properties and expression pattern. In mammalian cells, the genes with high GC content have several fold high expression levels compared to GC poor genes. From the study it was observed that the OG were having higher GC content percentage than TSG, suggesting that OG might be evolutionarily conserved and are mediated by methylation of CpG islands in the genes. Genome wide study of GC content and methylation profiles gave evidences that expression level and methylation are correlated stating, high GC rich genes have possibilities of high methylation and suppression in gene expression. The constant mutation rate of a genome is about 10-10/bp per generation. In a cellular environment, several physiological and pathophysiological conditions can alter the rate of mutation drastically. Cancer cells are often believed to have high mutation rate, leading to resistance mechanism and progression. With the help of TCGA data, the mutation rates of the breast cancer associated genes were calculated, which clearly shows that OG have higher mutation frequency rate, compared to the other gene sets. The genes PIK3CA, GATA3 and CDHI were the top genes with high mutation frequency rate. PIK3CA was reported as one of the reasons for cervical cancer and oncogenic in nature. GATA3 is a gene which encodes transcription factor binding proteins belonging to GATA family, involved in T-cell development, endothelial cell biology and also been reported in breast cancer.

Proteins encoded by genes often interact physically/structurally to perform their functions, these interactions between proteins are referred as protein–protein interaction (PPI). By considering proteins as nodes and the interactions among them as edges, PPI networks can be studied with the help of network topological parameters to identify key proteins.33, 34 The present study was designed to understand the interplay of interactions among the gene sets and in particular to compare the network parameters of OG and TSG. Interestingly, these two gene sets had no difference in measures of their network properties except degree centrality. It must be noted that the classified number of OG and TSG from the DEG is 102 and 41 respectively. Degree centrality of a node is an important parameter which is measured by the number of interactions with other nodes in a network. In case of OG, the average degree was found to be 11.13 which is significantly higher than the degree of TSG, which is 4.53. On the other hand, the essential genes were having highest average degree of 17.5 and druggable genes with 7.10, which clearly indicates that the essential genes are major players in a network and are indispensable for an organism. The betweenness centrality is a measure to find number of shortest paths passing through a particular node, which possibly act as a intermediate for exchange of information. The TSG are found to have slightly higher betweenness centrality than OG, indicating their intermediate level in process like activation. Clustering coefficient (CC) of a node represents the interconnectivity to form triangular sub clusters in a network. The gene sets OG, EG and DG were found to have almost similar clustering coefficient of about 0.4, which shows their interconnectivity in the network is on average, yet a complete picture of network might provide more reasonable explanation. The CC of TSG and OtG were about 0.2, which is comparatively less than that of OG indicating the relatedness in a network. The path which costs least number of edges to reach from one node to other is called as shortest path in a network. On average, OG and TSG has the shortest path distance to other nodes of 2.39 and 2.42 respectively meaning they interact more closely and can be traversed to other genes in the network. Druggable genes (DG) are vital and can act as therapeutic target for a disease. The average degree of interaction of DG is 7.10 and with clustering coefficient 0.47 which clearly shows that these genes are significantly interacting to genes in the network to a certain extent (Fig. 4). From the DEG interactions, a subnetwork of OG and TSG was constructed to look at the interplay in particular, which is depicted in Fig. 5 Analysing of this network showed that the tumor suppressor genes CDH1, CDKN2A and the oncogenes CDK1, CCNB1, FOS, TOP2A and AURKA were having high degree of interaction.

Figure 4

Comparison of various network topological parameters of five sets of genes, where oncogenes are showing high centrality measures compared to tumor suppressor genes.

Figure 5

Interaction network of oncogenes and tumor suppressor genes, where high degree nodes are represented with larger node size. The nodes coloured in green and red are oncogenes and tumor suppressor genes respectively.

The TSG genes namely CDH1, CDKN2A, CUL2, E2F1, ITGB1, STAT1 and the oncogenes CDK1, E2F1, FGFR3, FN1, FOS, LEF1, MAX, MMP9, RET, TGFBR2 and VEGFA are reported to be involved in many of the major cancer pathways like WNT signalling pathway, focal adhesion pathway and p53 signalling pathway.35, 36, 37 CDH1, a calcium dependant cell adhesion protein is involved in regulation of β-catenin, loss of this protein plays a critical role in cadherin-based adhesion and it also act as a co-activator of WNT signalling pathway to initiate differentiation process. From the observations, it was noticed that the genes CDH1 and LEF1 have remote interactions suggesting that the lymphocyte enhance factor (LEF1) or the T-Cell factor with cadherin free β-catenin complex act as control to facilitate WNT mediated gene expression to promote chromatin remodelling, transcription initiation and elongation to evade apoptosis. The cell matrix proteins FN1, ITGB1 and FGFR3 are involved in focal adhesion and cytokine–cytokine interactions which were identified to be differentially expressed in breast cancer. FN1, an ECM glycoprotein is known to be a key regulator of breast cancer cell adhesion and migration by binding to interleukin in focal adhesion pathway. The integrin subunit β1 (ITGB1) and FN1 are facilitating in the focal adhesion pathway and these genes mediate the ECM interaction by dysregulating the focal adhesion pathway which is observed in transition from ductal carcinoma in situ (DCIS) to invasive breast cancer. The oncogene COL11A1, a collagen type XI alpha 1 is overexpressed in most of the samples in the study and is also known to promote tumor aggressiveness. Cyclin dependant kinase inhibitor 2A protein (CDKN2A) regulated the tumor suppression of p53 signalling pathway and loss of p16 decrease the ability to repair DNA damage. It was also reported that the thrombospondin-1 (THBS1) is a negative regulator of new blood vessel formation and metastasis. SFN and CDK1 are noticed to downstream target genes for G2 phase arrest in cell cycle. Stratifin (SFN 14-3-3 ς) binds to CDK2 and CDK4 to arrest cell cycle in eukaryotes and inactivation of SFN lead to carcinogenesis through p53 signalling.

In conclusion, the present study reported insights on oncogenes and tumor suppressor genes in breast cancer using an integrated gene expression and interactome analysis. From the observations it was noticed that the dysregulated genes in breast cancer are involved in cell cycle, neoplasm related pathways and majorly involved in protein binding, nucleotide binding and regulating cellular functions. Oncogenes were found to have high GC content, mutation frequency and interactions when compared to the tumor suppressor genes. This study explored several key genes which are majorly involved in cancer progression and development, providing complementary clues for further experimental studies to validate and identify potential biomarkers.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Compliance with ethical standards

No ethical approval was required for this study.

Consent for publication and competing interests

All of the authors declare that they have no competing interests. All authors have read and approved the final manuscript.