Differential expression analysis in ovarian cancer: A functional genomics and systems biology approach.

BACKGROUND: Ovarian cancer is one of the rarest lethal oncologic diseases that have hardly any specific biomarkers. The availability of high-throughput genomic data and advancement in bioinformatics tools allow us to predict gene biomarkers and apply systems biology approaches to get better diagnosis, and prognosis of the disease with a tentative drug that may be repurposed.
OBJECTIVE: To perform genome-wide association studies using microarray gene expression of ovarian cancer and identify gene biomarkers, construction and analyze networks, perform survival analysis, and drug interaction studies for better diagnosis, prognosis, and treatment of ovarian cancer.
METHOD: The gene expression profiles of both healthy and serous ovarian cancer epithelial samples were considered. We applied a series of bioinformatics methods and tools, including fold-change statistics for differential expression analysis, DisGeNET and NCBI-Gene databases for gene-disease association mapping, DAVID 6.8 for GO enrichment analysis, GeneMANIA for network construction, Cytoscape 3.8 with its plugins for network visualization, analysis, and module detection, the UALCAN for patient survival analysis, and PubChem, DrugBank and DGIdb for gene-drug interaction.
RESULTS: We identified 8 seed genes that were subjected for drug-gene interaction studies. Because of over-expression in all the four stages of ovarian cancer, we discern that genes HMGA1 and PSAT1 are potential therapeutic biomarkers for its diagnosis at an early stage (stage I). Our analysis suggests that there are 11 drugs common in the seed genes. However, hypermethylated seed genes HMGA1 and PSAT1 showcased a good interaction affinity with drugs cisplatin, cyclosporin, bisphenol A, progesterone, and sunitinib, and are crucial in the proliferation of ovarian cancer.
CONCLUSION: Our study reveals that HMGA1 and PSAT1 can be deployed for initial screening of ovarian cancer and drugs cisplatin, bisphenol A, cyclosporin, progesterone, and sunitinib are effective in curbing the epigenetic alteration.

Introduction

Ovarian cancer is one of the lethal gynaecological diseases endured by most women around the globe. It is discerned that “Ovarian cancer is more than a woman's disease”. A 2018 GLOBOCAN prediction report reveals that by the year 2040, the incidence and mortality of ovarian cancer may go up to 150,000 cases making the disease a big reason to worry about (https://gco.iarc.fr/tomorrow/). The symptomatic signs are vague and can be misinterpreted for other diseases. According to the National Ovarian Cancer Coalition (NOCC) and Cancer Treatment Centres of America, there are more than 30 different types of ovarian cancer, usually classified by originating cell-type (http://ovarian.org/). It commonly originates from three common cell-types, namely surface epithelium cells, germ cells, and stromal cells. However, epithelial ovarian cancers, developed from the cells of the outer surface of the ovary, are more prevalent and account for 85–90% of all ovarian cancers (Romero and Bast, 2012). Unfortunately, ~70% of the patients having epithelial ovarian cancers remain undiagnosed until it reaches an advanced stage. The major subtypes of epithelial ovarian cancer are namely - a) High-grade serous ovarian cancer (HGSOC), b) low-grade serous ovarian cancer (LGSOC), c) clear cell tumors, d) mucinous, and e) endometrioid (Qazi, 2018, Qazi et al., 2021).

The median age of women getting ovarian cancer diagnosed is 63, while the calculated median age of death is 70 (Jemal and Bray, 2011, Howlader and Noore, 2017). Ovarian cancer is disease-laden with paradigms, and it is a serious health problem. To know about its source of origin, tracking the natural history of the disease becomes mandatory, which is very difficult and challenging for oncologists. Ovarian cancer tends to metastasize in a very unique manner, wherein cells are exfoliated from the primary tumor as a single cell or many cell aggregates circulating via the peritoneal fluid. Since the past decade, a drastic paradigm shift in ovarian cancer was observed from being a singular disease to a nest of diseases. Another biggest risk factor apart from age is the “family history”. Even though germline mutations that indoctrinate ovarian cancer are janitors for 10–20% of cases. Genes responsible for inheriting the chances of enduring ovarian cancer in families are mainly – BRCA1 and BRACA2 (Assis and Pereira, 2018). Therefore, the major factors leading to the progression of ovarian cancer are mainly- a) age, b) family history, c) infertility issues, d) obesity, e) sedentary lifestyle, f) hormonal replacement therapy (HRT), g) in vitro fertilization (IVF), and h) endometriosis (https://www.mskcc.org/). These factors are mainly categorized as, a) genetic factors, and b) environmental factors. A 2009 study performed by Bowen and his colleagues suggest that epithelial ovarian cells in humans are dominant and may serve as the leading factor for the proliferation of ovarian adenocarcinoma (Bowen et al., 2009).

As per the National Cancer Institute (NCI) Surveillance, Epidemiology and End Results Program (SEER) reports, as of 2020 there are 21,750 (1.2%) incidence cases registered with ovarian cancer while the estimated number of mortalities is 13,940 (2.3%). The five-year survival rate (2010–2016) is 48.6% (https://seer.cancer.gov/statfacts/html/ovary.html). Fig. 1 represents ovarian cancer estimated statistics.

Fig 1

Estimated ovarian cancer cases reported in 2020 with five year survival percentage.

There are mainly four stages of cancer that describe to what extent the cancer is present in the body, thus, determines the specific treatment strategy and survival duration. Treatment of ovarian cancer in its initial stages is much easier and the patient survival is much longer when compared to advanced stages. Stage 1 or simply known as “localized” is when the cancer is found in a restricted region of the body. Stage 2 is referred to as “regional” as cancer gets spread to surrounding regions from its source, while advanced stages (stages 3 & 4) are distant and unknown and are a condition where the cancer cells get metastasized in the entire body. As far as ovarian cancer is concerned, 15.7% are diagnosed at the local stage. The 5-year relative survival for localized ovarian cancer is 92.6% (https://seer.cancer.gov/statfacts/html/ovary.html). Fig. 2 represents the stage-wise survival percentages in the case of ovarian cancer.

Fig 2

Relative survival by stages at diagnosis in ovarian cancer.

The common diagnostic practices for the screening of ovarian cancer are mainly categorized into two approaches, a) Blood-based approaches encapsulating carbohydrate antigen-125 (CA-125) and human epididymis-4 (HE-4), and b) Imaging investigations which comprise of transvaginal sonography (TVS), Doppler ultrasonography, computerized tomography (CT) scans and ultrasounds (Kurjak et al., 1991). These approaches calculate the risk of malignancy index (RMI) in ovarian cancer patients which is composed of a score given to the transvaginal ultrasound outcomes, menopausal condition, and CA-125 level. The RMI value higher than 200 indicates a greater risk of malignancy (Jayson et al., 2014). HE-4 is utilized for ovarian cancer diagnosis as its expression is observed in many organs but not in the ovary. Higher levels of HE-4 are seen in serous and endometrioid subtypes of ovarian cancer which makes it very crucial and sensitive in its diagnosis. HE-4 is not increased in benign forms, as it is in the case of CA-125, making it a specific prognostic indicator for lethal ovarian cancer. With this, we know that CA-125 has a higher sensitivity while HE-4 has a high specificity, and these two studies are amalgamated and implied as a mathematical formula named the Risk of Malignancy Algorithm (ROMA). ROMA has been successful in achieving better screening tests for ovarian cancer. The FDA in 2009 ascertained the usage of the serum OVA-1 for analyzing these prognostic indicators: CA-125, II-microglobulin, apolipoprotein A1, prealbumin, and transferrin. ROMA takes in the biomarker estimate levels to give output between 0 and 10 for menopausal patients. Patients with greater scores refer to the complexity of the malignancy and thus, are assessed by the gynaecologist or oncologist. Currently, there are not many comparative studies of ROMA and OVA-1 (Novak and Lukasz, 2015). Imaging investigations for ovarian cancer are mainly Trans Vaginal Sonography (TVS), Doppler Ultrasonography, CT Scan, and Biopsy (Kobayashi, et al., 2012, McFarlane et al., 1956, Oglat et al., 2018). Initially, trans-abdominal ultrasonography was majorly used for screening of the disease, which was later, replaced by TVS that gives more accurate, and lucid images of the ovary (Bourne et al., 1993).

Previous gene expression studies on ovarian cancer have identified various genes that are overexpressed such as - CLDN3, WFDC2, FOLR1, COL18A1, CCND1, FLJ12988, group 3 POTEs especially - C, E, and F, CLASP1, IGFBP2, and TRAIL (Peters and Kudla, 2005, Sharma et al., 2019, Barger et al., 2018, Lisowska et al., 2014, Lancaster et al., 2004, Qazi and Raza, 2021). However, the majority of these studies were executed using microarrays. Since ovarian cancer expression datasets are humongous in size and the network biology of the disease is also complex, expression analysis must be performed using high-throughput techniques.

In this proposed study, we aim to identify potential gene biomarkers using microarray gene expression profiles, apply systems biology approaches to construct regulatory networks, and perform its topological and other interaction analysis, detect modules in the network, and perform survival analysis of identified seed genes. Further, we also perform a drug-gene interaction analysis of the identified seed genes. Our in silico analysis results will help to get better disease insights, diagnosis, and prognosis of the disease with a tentative drug that can be repurposed, thus paving a way for person-centric healthcare management for ovarian cancer.

Materials and method

Datasets

We considered the gene expression profile of healthy epithelial samples (control) and samples of ovarian cancer patients (case). A dataset of 24 samples was considered, including 12 healthy ovarian surface epithelial and 12 laser capture micro-dissected serous ovarian cancer epithelial samples available at GEO with Accession No. GSE14407 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE14407).

Methods & tools

The complete methodological pipeline used in this study is shown in Fig. 3, and explained as follows:

Fig. 3

Methodological pipeline.

Gene biomarkers identification

Gene expression profiling was carried out on the GSE14407 dataset to identify the crucial genes indulged in epithelial ovarian cancer. The obtained dataset was pre-processed and log normalized using the GEO2R web-based tool of NCBI-GEO (https://www.ncbi.nlm.nih.gov/geo/geo2r/). To identify key regulators having differential expression, fold change (FC) statistics and p-value measures were considered. The FC, one of the widely used methods used for differential expression analysis, is a statistical measure that describes how much the level of expression of a gene changes over the two different conditions such as epithelial ovarian cancer and control samples (i.e., case-control analysis) (Raza and Hasan, 2015, Raza, 2014). The FC is calculated as a ratio of averages from control and disease samples (Raza, 2016), and quantified as a log of FC (logFC). Usually, logFC ≤ −1.0 is considered as down-regulated, while logFC ≥ 1.0 is treated as up-regulated (Raza, 2014). The list of short-listed genes as DEGs is forwarded for further down-stream analysis.

Gene Ontology enrichment analysis

Gene Ontology (GO) enrichment analysis is a method to interpret a set of genes and assign them to a set of predefined classes based on their functional characteristics. Given a set of genes that are either up-regulated or down-regulated under some conditions, GO enrichment analysis will answer which GO terms are overrepresented or underrepresented based on pre-defined and pre-stored annotations (Raza, 2016). Some of the GO enrichment analysis tools are DAVID (Huang et al., 2009), PANTHER (Mi et al., 2019); BiNGO (Maere et al., 2005); GeneWeaver (Bubier et al., 2015), and so on. In this study, we used the recent version of DAVID 6.8.

Reconstruction and analysis of network

Reconstruction of gene regulatory networks (GRNs) is important to further perform topological analysis, network module identification, and highly influential gene identification (seed genes) that can be further utilized for drug target identification. For the reconstruction of GRNs, we utilized GeneMANIA (Warde-Farley et al., 2010, Franz et al., 2018) database. GeneMANIA is a regularly updated, publicly available database that stores genetic interaction data such as protein-protein, protein-DNA, pathways, etc. It stores interaction data which includes physical interactions, genetic interactions, co-expression, shared protein domains, co-localization, pathway, and computationally predicted. To retrieve the interactions, we uploaded the list of identified gene biomarkers (DEGs) to the GeneMANIA and considered all the default parameter settings.

Once the interaction network is retrieved from the GeneMANIA, we used the Cytoscape 3.8 software tool (Shannon, 2003, Smoot et al., 2011) for network visualization and topological analysis. We used the NetworkAnalyzer (Assenov et al., 2008) plugin of Cytoscape to compute various topological properties of the constructed network.

Network module identification and survival analysis

A key problem in a complex and large biological network is the identification of functional units, called network modules or biological pathways. In other words, biological networks are generally very complex and large to be examined as a whole, therefore, module identification is an important step towards understanding the biological insight. Module identification is also known as community detection or graph clustering. A plethora of software tools and packages are available for module identification in a given network. We utilize a widely used and highly cited network clustering algorithm, called Molecular Complex Detection (MCODE) (Bader and Hogue, 2003) which is available as a plugin in Cytoscape software. We identified highly interconnected regions in the network (clusters or modules), which are often protein complexes and parts of pathways. To identify modules in a given network, MCODE uses a vertex-weighting approach based on the clustering coefficient C = 2n/k(k − 1), where k is the vertex size of the neighborhood of vertex i and n represent the number of edges in the neighborhood.

To analyze the roles of identified ovarian cancer gene biomarkers on patient survival, we used the UALCAN (Chandrashekar et al., 2017) tool to generate Kaplan-Meier (KM) survival plots. The UALCAN tool is based on genomics data from The Cancer Genome Atlas (TCGA). The TCGA (Tomczak et al., 2015) project makes available a large number of samples of various cancer types, allowing us to perform comprehensive molecular characterization and answer the questions related to tumor heterogeneity. The UALCAN utilized TCGA level-3 RNA-seq and clinical data from more than 30 different cancer types along with samples of various cancer stages, tumor grades, and other clinico-pathological features. It is publicly available at http://ualcan.path.uab.edu.

Gene-Drug interaction analysis

To understand the mechanism of how oncologic medications function when treated on cancers, it is important to check for the interaction affinity of the drugs being used to treat the disease and the target genes of ovarian cancer network biology. In our case, once the seed genes were identified using the previously mentioned methodologies, major synthetic drugs were identified using various Bioinformatics repositories and databases namely - PubChem (https://pubchem.ncbi.nlm.nih.gov/), DrugBank (Wishart et al., 2008) and DGIdb (Cotto et al., 2018) that either increase or decrease the genetic expression of these seed genes in the ovarian cancer network dynamics. The preliminary criteria were to search for drugs that a) have a strong interaction with the target seed genes, b) the interaction score should be >=5.0, and c) these drugs must be reported as a treatment for ovarian cancer in the literature. This was executed in a keyword-based strategy modus operandi, wherein each of the identified seed genes was separately searched in PubChem, DrugBank, and DBIdb databases.

Results

In this study, gene expression profiles of 12 healthy ovarian surface epithelial samples and 12 laser capture micro-dissected serous ovarian cancer epithelial samples are considered, whose boxplot is depicted in Fig. 4. It is observed from the boxplot (Fig. 4) that the median of Ovarian adenocarcinoma samples is greater than those of healthy ovarian surface epithelial samples. Although, quartile Q1 are almost similar in both the sample types, but their quartile Q3 and maximum range greatly vary. Further, healthy ovarian surface epithelial samples data are more skewed than that of ovarian adenocarcinoma samples. Hence, these expression profile statistics indicate that there are a significant number of genes that must be differentially expressed over the two sample types.

Fig. 4

Boxplot of gene expression samples (normal and ovarian adenocarcinoma samples).

Gene biomarkers

To identify gene biomarkers (i.e., differentially expressed genes (DEGs)), we executed the NCBI-GEO2R tool and computed fold-change (FC) statistics between the two groups of samples at genome-scale, which is usually converted to the log of FC (logFC). To control the false discovery rate (FDR) in the tests, we applied Benjamini & Hochberg procedure and FDR is set to 0.05. The scatter-plot of logFC of all the genes is depicted in Fig. 5. It is observed that more than a thousand genes are differentially expressed over the two sample types having score −2.0 ≥ logFC ≥ 2.0 with p-value ≤ 0.05, most of them are up-regulated while few as down-regulated genes.

Fig. 5

Scatter-plot of logFC of all the genes.

To narrow-down our analysis to potential gene biomarkers, we considered the top 500 genes having the best logFC scores and lesser p-value, which is mapped to gene-disease association databases DisGeNET and NCBI-Gene (with query string=’(“ovarian cancer”) AND biomarker)’. For further down-stream analysis, we considered only those genes which have an association with ovarian cancer, as reported in DisGeNET and obtained in the NCBI-Gene database after executing the said query. In this way, we obtained a list of 149 significant genes that may have a potential role in ovarian cancer. The list of identified 149 significant genes with their statistical scores such as adjusted p-value, p-value, moderated t-statistics, B-statistics, logFC, etc. is available in the supplementary Table S1.

GO enrichment analysis

We performed a GO enrichment analysis of the identified potential gene biomarkers of ovarian cancer for its enrichment with the disease term (GAD-Disease), pathways (KEGG-pathway), and tissue enrichment. It is observed that the set of supplied genes are significantly enriched in various cancer diseases (including ovarian cancer), significantly enriched in KEGG pathways, and various tissues, a snapshot of these results with higher gene-count and lesser FDR are presented in Table 1, Table 2 and Table 3, respectively.

Table 1

Gene-disease association (GAD) enrichment results.

Term	Gene-count	%	P-Value	Fold Enrichment	Bonferroni	Benjamini	FDR
Breast cancer	25	17.0	0.0000	4.2364	0.0000	0.0000	0.0000
Ovarian cancer	24	16.32	0.0000	5.4776	0.0000	0.0000	0.0000
Breast cancer	20	13.60	0.0000	3.9125	0.0008	0.0002	0.0011
Colorectal cancer	19	12.92	0.0000	4.5196	0.0002	0.0001	0.0003
Lung cancer	19	12.92	0.0000	3.7313	0.0032	0.0005	0.0044
Bladder cancer	17	11.56	0.0001	3.2565	0.0673	0.0077	0.0935
Prostate cancer	16	10.88	0.0001	3.1791	0.1515	0.0163	0.2204

Table 2

KEGG pathways enrichment analysis results.

Term	Gene count	%	P-value	Fold Enrichment	Bonferroni	Benjamini	FDR
hsa05200:Pathways in cancer	21	14.09	0.0000	4.8365	0.0000	0.0000	0.0000
hsa04151:PI3K-Akt signaling pathway	11	7.38	0.0040	2.8859	0.4505	0.0950	4.6553
hsa04014:Ras signaling pathway	9	6.04	0.0030	3.6045	0.3664	0.0872	3.5673
hsa04110:Cell cycle	8	5.37	0.0004	5.8395	0.0547	0.0277	0.4470
hsa04310:Wnt signaling pathway	8	5.37	0.0007	5.2471	0.1017	0.0351	0.8506
hsa04015:Rap1 signaling pathway	8	5.37	0.0077	3.4481	0.6840	0.1517	8.7631
hsa05166:HTLV-I infection	8	5.37	0.0202	2.8508	0.9531	0.3178	21.6204
hsa05206:MicroRNAs in cancer	7	4.70	0.0908	2.2153	1.0000	0.6140	67.9138
hsa05218:Melanoma	6	4.03	0.0010	7.6489	0.1395	0.0369	1.1893
hsa05205:Proteoglycans in cancer	6	4.03	0.0665	2.7153	1.0000	0.5479	56.0270
hsa04114:Oocyte meiosis	5	3.36	0.0325	4.0771	0.9930	0.3911	32.6281
hsa04390:Hippo signaling pathway	5	3.36	0.0821	2.9971	1.0000	0.6006	64.0428
hsa05217:Basal cell carcinoma	4	2.68	0.0207	6.7046	0.9569	0.2949	22.1504
hsa04115:p53 signaling pathway	4	2.68	0.0363	5.4037	0.9961	0.3957	35.6682

Table 3

Tissue enrichment analysis results.

Term	Gene count	%	P-value	Fold Enrichment	Bonferroni	Benjamini	FDR
Epithelium	34	22.8187	0.0048	1.6080	0.4639	0.1876	5.4309
Colon	18	12.0805	0.0151	1.8656	0.8621	0.3271	16.2602
Plasma	10	6.7114	0.0004	4.4843	0.0500	0.0253	0.4584
Mammary carcinoma	5	3.3557	0.0132	5.4583	0.8225	0.3509	14.3465
Pancreatic carcinoma	4	2.6845	0.0003	31.3372	0.0327	0.0327	0.2972
Aorta endothelial cell	4	2.6845	0.0242	6.4185	0.9585	0.4115	24.7947
Lymphocyte	4	2.6845	0.0324	5.7283	0.9861	0.4572	31.8300
Umbilical vein endothelial cell	3	2.0134	0.0461	8.6859	0.9978	0.5359	42.3169
Fetal lung	3	2.0134	0.0915	5.8757	1.0000	0.7499	67.2844

Table 1 represents the gene-disease association (GAD) results for different cancers namely, breast cancer, ovarian cancer, colorectal cancer, lung cancer, prostate cancer, and bladder cancer. The gene-disease association analysis yielded many significant signals with a p-value of 0.00 for disease breast cancer that had a gene-count of 25 and 20 significant genes, ovarian cancer had 24 prominent genes while colorectal cancer and lung cancer with both had 19 significant genes. After applying FDR control on the p-values from the gene-disease association analysis, only ovarian cancer achieved a significant level with 5.4776-fold enrichment and 0.00 Bonferroni and Benjamini FDR control.

Reconstruction and analysis of cancer network

For the reconstruction of the regulatory network, we utilized the GeneMania interactome database (Warde-Farley et al., 2010, Franz et al., 2018). All the interactions of the identified 149 gene biomarkers were retrieved from the GeneMania database by supplying their gene IDs and saved to a CSV file for further visualization and network analysis. We found around 8,987 interactions among 149 genes without any self-loop. Further, we used the Cytoscape software tool (Shannon, 2003) for network visualization and its topological analysis. Various topological properties of the reconstructed network, as computed by the NetworkAnalyzer (Assenov et al., 2008) plugin of Cytoscape, are depicted in Table 4.

Table 4

Various topological properties of the reconstruction network.

Topological parameters	Value
Cluster coefficient	0.386
Network diameter	3
Network centralization	0.223
Characteristics path length	1.825
Average number of neighbors	34.873
Network density	0.213
Network heterogeneity	0.427

The other topological properties of the reconstructed network such as node degree distribution, shared neighbor distribution, average clustering coefficient distribution, topological coefficient, and betweenness and closeness centralities are depicted in Fig. 6. The node degree distribution provides us the structure of the network and hubs of the network. The degree of most of the genes lies between 10 and 100, while few genes are having a degree between 300 and 500 (see Fig. 6(A)). The shared neighbor distribution provides the number of node partners shared between the nodes which may help us to detect network motifs in the network if any (see Fig. 6(B)). In Fig. 6(C), as the number of neighbors is increased, its average clustering coefficients are also increased, stating that the network has a modular structure. The topological coefficient, a relative measure to assess the extent to which a gene shares its neighbors with other genes, shows an increasing trend as neighbors are increased (see Fig. 6(D)). The betweenness centrality, a measure of relevance of a gene as functionally capable of holding together interacting genes, shows an increasing trend indicating that genes are functionally capable of holding together interacting genes (see Fig. 6(E)). Fig. 6(F) depicts the closeness centrality which measures how fast information spreads from a given gene to other reachable genes within the network. It shows an increasing trend in closeness centrality, varying from 0.425 to 0.64.

Fig. 6

Topological properties of the reconstructed network. The red line fits the power law.

Module identification and survival analysis

We executed the MCODE (Bader and Hogue, 2003) network clustering algorithm on the reconstructed ovarian cancer network with its default parameters (degree cutoff = 2, node score cutoff = 0.2, haircut = true, fluff = false, k-core = 2, max. depth from seed = 100), and found a total of eight clusters with one seed gene per cluster, as shown in Fig. 7. The seed genes identified in these eight clusters are NDC80 (MCODE Score = 38.85), PDGFD (MCODE Score = 3.66), SNCA (MCODE Score = 2.91), AGR2 (MCODE Score = 6.15), PSAT1 (MCODE Score = 4.61), EFEMP1 (MCODE Score = 5.16), HOXD8 (MCODE Score = 2.32), and HMGA1 (MCODE Score = 30.94), respectively shown as yellow color nodes in Fig. 7. The seed genes are the highest scoring node and most influential node in the cluster from which the cluster was derived.

Fig. 7

Identified modules in the ovarian cancer network (yellow color node represents the seed genes).

The list of seed genes identified by module detection algorithms is then subjected for expression pattern analysis (Fig. 8), cancer stage expression analysis (Fig. 9), and survival analysis (Fig. 10) using the UALCAN (Chandrashekar et al., 2017) resource and analysis tool. It can be discerned from Fig. 8 that the expression pattern of the seed genes, including HMGA1, PSAT1, HOXD8, AGR2, and EFEMP1, are higher in Ovarian serous cystadenocarcinoma. Further, Fig. 9 states that these seed genes are either over-expressed or under-expressed in various stages of ovarian cancer.

Fig. 8

Heatmap diagram showing expression pattern of seed genes in ovarian serous cystadenocarcinoma. The expression pattern is represented in the log2 transform of total transcript per million (TPM) [log2(TPM + 1)] in the TCGA samples.

Fig. 9

Box-whisker plots showing the expression of seed genes in different stages of ovarian serous cystadenocarcinoma in TCGA samples.

Fig. 10

Kaplan–Meier plots showing the association of identified seed genes expression in the survival of ovarian serous cystadenocarcinoma.

The list of seed genes is considered to further perform survival analysis. Survival analysis results for each of the seed genes are depicted in Fig. 10 with the survival probability of patients over a period of 15 years against the abnormal expression of the gene of interest. In Fig. 10, the survival probability of the patients gradually drops over time, and over a period of 10–15 years, the survival probability drops to less than 0.10.

Gene-Drug interaction analysis

We have identified 8 seed genes that are highly influential in the reconstructed network of ovarian cancer that includes NDC80, PDGFD, SNCA, AGR2, PSAT1, EFEMP1, HOXD8, and HMGA1. These high influential gene biomarkers of ovarian cancer were taken for further gene-drug interaction analysis for the drug repurposing.

When each of the genes was searched in DGIdb, no significant or common interactions were identified. It retrieved interactions for only three genes namely - a) PSAT1 with two drugs pyridoxal phosphate, which acts as a co-factor, and L-glutamate whose interaction type was not available; b) PDGFD with sunitinib, which inhibits the expression of the gene and tandutinib interaction was not mentioned and c) SNCA with cinpanemab, wherein it acts as an inhibitor for SNCA gene (Fig. 11). There were no drug interactors available in DGIdb for the remaining five genes - HMGA1, NDC80, AGR2, EFEMP1, and HOXD8. However, PubChem and DrugBank on the other hand provided drug interactors for all the genes. There were more than 30 drug interactions provided for the 8 seed genes, but common drug interactions were only 11. These common drugs are: a) Sunitinib, b) Cisplatin, c) Cyclosporin, d), Copper Sulfate, e) Bisphenol A, f) Estradiol, g) Valproic Acid, h) Progesterone, i) Dexamethasone, and j) Benzo(a)pyrene (Table 5). The pharmacokinetics (PK), pharmacodynamics (PD) and physicochemical properties of these drugs are available in Supplementary Table S2.

Fig. 11

Drug interactions with hypermethylated and hypomethylated seed genes.

Table 5

Highest scoring and prominent gene-drug associations.

Seed gene	Drug	Interaction	PMID(s)/ChEMBL Id
HMGA1	Cyclosporin	Decreased expression of TSPAN14 mRNA	25562108|27989131
	Sunitinib	Increased expression of ABCC5 mRNA and FZD9 mRNA	31533062|31533062
	Cisplatin	Increased expression of MIR188 mRNA	24,880,025
	Progesterone	Increased expression of OSR2 mRNA	20,106,945
	Valproic Acid	Increased expression of KLKB1 mRNA	23179753|24383497|24935251|26272509|28001369
	Estradiol	Decreased expression of GLYAT mRNA	20,106,945
	Bisphenol A	Increased expression of JADE3 mRNA	29,275,510
AGR2	Cyclosporin	Increased expression of TBC1D15 mRNA	20106945|25562108|27989131
	Doxorubicin	Decreased expression of phosphorylation of BAD protein	20,884,855
	Estradiol	Increased expression of EGFR mRNA	18,692,832
	Dexamethasone	Increased expression of NFIKBIB mRNA	28,628,672
	Benzo(a)pyrene	Decreased expression of FUOM, FURIN and FUS mRNA	20106945|22316170|20106945|21632981
EFEMP1	Cyclosporin	Increased expression of TNC mRNA	20,106,945
	Sunitinib	Increased expression of AHCYL2 mRNA	31,533,062
	Cisplatin	Co-treated with Jinfukang results in increased expression of TMCC1-AS1 mRNA	27,392,435
	Copper Sulfate	Decreased and increased expression of FOXRED1 mRNA	19,549,813
	Doxorubicin	Decreased expression of MAP4 mRNA	29,803,840
	Bisphenol A	Increased expression of CSRNP3 mRNA	29,275,510
	Estradiol	Decreased expression of AKT1 mRNA; Increased expression of TLL2 mRNA	23373633|23019147
	Valproic Acid	Increased and decreased expression of MXRA5 mRNA	23179753|24383497|26272509
	Dexamethasone	Decreased expression of PRPS1L1 mRNA	28,628,672
	Progesterone	Increased expression of MAPK3 protein	15358673|16175315
HOXD8	Valproic Acid	Increased methylation of HOXD8 gene	29,154,799
	Copper Sulfate	Decreased expression of HOXD8 mRNA	19,549,813
	Bisphenol A	Increased expression of HOXD8 mRNA	30,951,980
	Benzo(a)pyrene	Decreased expression of HOXD8 mRNA	22,228,805
NDC80	Cyclosporin	Decreased expression of NDC80 mRNA	20106945|21632981|25562108
	Sunitinib	Decreased expression of NDC80 mRNA	31,533,062
	Cisplatin	Increased expression of NDC80 mRNA	27,594,783
	Copper Sulfate	Decreased expression of NDC80 mRNA	19,549,813
	Doxorubicin	Decreased expression of NDC80 mRNA	30,031,762
	Bisphenol A	Increased and decreased expression of NDC80 mRNA	16474171|15223131|27685785|29275510
	Estradiol	Increased expression of NDC80 mRNA	16474171|18310284|19167446|20106945
	Dexamethasone	Decreased expression of NDC80 mRNA	28,628,672
	Valproic Acid	Decreased expression of NDC80 mRNA	19101580|23179753|27188386|28001369
	Progesterone	Increased expression of NDC80 mRNA	18,070,364
	Benzo(a)pyrene	Decreased expression of NDC80 mRNA	20064835|22316170|22579512
PSAT1	Cyclosporin	Increased and decreased expression of PSAT1 mRNA	20106945|21632981|22147139|25596134|27989131|25562108
	Benzo(a)pyrene	Increased expression of PSAT1 mRNA	20,106,945
	Progesterone	Increased expression of PSAT1 mRNA	18,037,150
	Bisphenol A	Increased and decreased expression of PSAT1 mRNA	29275510|20678512
	Copper Sulfate	Increased expression of PSAT1 mRNA	19,549,813
	Doxorubicin	Decreased expression of PSAT1 mRNA	29,803,840
	Cisplatin	Increased expression of PSAT1 mRNA	25596134|27392435|27594783
	Sunitinib	Increased expression of PSAT1 mRNA	31,533,062
SNCA	Cinpanemab	Decreased expression of SNCA mRNA	CHEMBL3833330
PDGFD	Sunitinib	Decreased expression of PDGFD mRNA	CHEMBL535

Discussion

Our in silico analysis identified a list of 149 significant genes that may have a potential role in ovarian cancer. The KEGG pathways and tissue enrichment suggests that there are 110 DEGs that are involved in the oncologic pathways and are mainly found in the epithelium cells with a gene presence percentage of 22.8%, followed by a colon with 18 genes (12.0%), and 10 genes in the plasma (6.71%). A scanty number of genes with low gene presence percentage (>5.0%) were also recorded to be present in mammary and pancreatic carcinoma, endothelial cells, fetal lung cells, and lymphocytes (refer Table 3). The bifurcation of the DEGs as observed in GO analysis is as follows (refer to Table 2). A total of 21 genes (hsa05200) play a pivotal role in cancer pathways which is evident as it holds a 4.8365-fold enrichment with 0.00p-value. The other 89 genes are crucial in cancer proliferation and progression dynamics. Out of these genes, there were 8 genes (hsa04110) involved in the cell cycle that had a prominent fold change of 5.8395, 8 genes (hsa04310) were engaged in the Wnt signaling pathway with a fold change of 5.2471. Six genes were found to be imperative for melanoma (hsa05218) with a fold change of 7.6489. Four genes are significant for basal cell carcinoma formation (hsa05217) with a fold change of 6.7046, while four genes with a fold change of 5.4037 are crucial in the p53 signaling pathway (hsa04115).

The module detection algorithm, expression pattern analysis, cancer stage expression analysis, and survival analysis identifies DNA methylation as a pivotal player in the formation and progression of ovarian cancer and in several studies, many tumor suppressor genes have been shown to be hypermethylated (Hentze et al., 2019). It is evident that genes HMGA1 and PSAT1 are majorly hypermethylated referring to the fact that there is an increase in the epigenetic methylation of cytosine and adenosine residues, indicating their crucial role in the proliferation of ovarian cancer. Genes EFEMP1, AGR2, HOXD8, and NDC80 are partially hypermethylated. This suggests that hypermethylated genes are ubiquitous in all subtypes of ovarian cancers irrespective of the type and stage of cancer. Hypomethylation, on the other hand, is not prominently observed in ovarian cancer. PDGFD and SNCA genes show clear signs of hypomethylation, referring to the detachment of methyl groups in the 5-methylcytosine nucleotide. Hypomethylation of genes leads to abnormal functioning of the DNA that further turns into incremented alterations and cancer development. This is mainly because of decreased methylation in regulatory regions of oncogenes that directly hampers the transcription leading to its greater processing. Out of 8 seeds genes, six genes (HMGA1, PSAT1, EFEMP1, AGR2, HOXD8, and NDC80) showcased hypermethylation and two (PDGFD and SNCA) displayed hypomethylation were further bifurcated to cross-validate their expression according to the four stages of ovarian cancer (OC). The HMGA1 gene is over-expressed and further increases its expression as the cancer shifts from the local region to other regions of the body (metastasis). In stage I, HMGA1 expresses itself at 500 transcripts per million (TPM) while, at stage II, it extends itself and expresses at 750 TPM. In stages III and IV (advanced stages), it can reach up to 1200 TPM. This validates that HMGA1 expression is higher in ovarian cancer than other significant genes identified. Besides, gene PSAT1 that also showcased an extreme overexpression in ovarian cancer expresses over the threshold value of 150 TPM in stage I. In stage II, it jumps over to 300 TPM and ranges more or less the same in advanced stages (stage III and IV). Apart from these, AGR2 significantly expressed at a range of 85–90 TPM in stage II, however, declined to 40 TPM in stage III and 60 TPM in stage IV, respectively. The remaining partially hypermethylated genes namely - HOXD8, expressed itself at a satisfactory level of 45–80 TPM in stages I to IV; EFEMP1 expressed itself at 10–60 TPM while NDC80 expressed itself in between 25 and 45 TPM in all the four stages. As far as hypermethylated genes are concerned, only two genes HMGA1 and PSAT1 overpower their expression in ovarian cancer (OC). Hypomethylated genes PDGFD displayed a comparatively low expression of 20–25 TPM in advanced stages III and IV while 45 TPM in stage II. It can be discerned that PDGFD is up-regulated in the localized stage (stage II) however, it becomes down-regulated in metastasis stages of III-IV. SNCA gene, on the other hand, displays a very low expression in all the four stages of ovarian cancer (I-IV). These findings suggest that HMGA1 and PSAT1 are up-regulated (overexpressed) in all the four stages of cancer and that they definitely can be deployed as specific biomarkers and as a potential therapeutic for ovarian cancer as they express themselves at an outrageously higher value in stages I itself.

The identified 8 seed genes which are highly influential in the reconstructed network of ovarian cancer that includes NDC80, PDGFD, SNCA, AGR2, PSAT1, EFEMP1, HOXD8, and HMGA1. These high influential gene biomarkers of ovarian cancer were taken for further gene-drug interaction analysis for the drug repurposing. The hypermethylated seed genes namely - HMGA1 and PSAT1 showcased a good interaction affinity with drugs cisplatin, cyclosporin, bisphenol A, progesterone, and sunitinib where their mRNA expression levels are increased when they are treated with these five drugs respectively, while hypomethylated seed genes - SNCA and PDGFD showcased an inhibited expression with drugs cinpanemab and sunitinib. Our results are further backed up by other essential findings wherein it has been discerned that HMGA1 and PSAT1 can be used as therapeutic biomarkers for ovarian cancer. A study by Kim et al. (Kim et al., 2016) showcased that HMGA1 is an important regulator for monitoring CSC-like features in ovarian cancer, thus can be deployed as a novel therapeutic marker for highly metastatic and drug resistant ovarian cancer (Kim et al., 2016). Also, in a recent study, researchers claim that their data shows that measuring urine HMGA1 can be employed as diagnostic tool for assessment of ovarian cancer patients (Zhou et al., 2015). Another study (Dai et al., 2019) suggests that PSAT1 is expressed in greater amount in ovarian cancer tissues than in adjacent normal tissues, thus hinting at it being used as a biomarker for the same. Similarly, another study indicates PSAT1 inhibitors as therapeutic biomarker option for patients with epithelial ovarian cancer (Zhang and Li, 2020). The partially hypermethylated drugs showcased an ambivert interaction with these 11 common drugs. This lucidly indicates that cisplatin; cyclosporin, bisphenol A, progesterone, and sunitinib are effective in preventing the hypermethylation of these seed genes in ovarian cancer.

Conclusion

Because of the over-expression in all the four stages of ovarian cancer, therefore, we can discern that genes HMGA1 and PSAT1 are potential therapeutic biomarkers for diagnosing ovarian cancer at an early stage (stage I). Hypermethylated seed genes, HMGA1, and PSAT1 are crucial in the proliferation of ovarian cancer, and when treated with drugs cisplatin, cyclosporin, bisphenol A, progesterone, and sunitinib, the mRNA expression levels of these genes increments which in turn helps in stabilizing the epigenetic methylation profiles in ovarian cancer subjects. Hypomethylated genes such as SNCA and PDGFD are not so common in the ovarian cancer network, however, when drugs cinpanemab and sunitinib are deployed, their expression levels are drastically decreased resulting in the alteration of de-methylation status in ovarian cancer network biology. Thus, our study reveals that HMGA1 and PSAT1 can be deployed for initial screening of ovarian cancer and drugs cisplatin, bisphenol A, cyclosporin, progesterone, and sunitinib are effective in curbing the epigenetic alteration of these genes in case of ovarian cancer.

Availability of data and materials

The microarray gene expression datasets used in this research are available on NCBI-GEO public repository.

Ethics approval and consent to participate

No participation of human is involved in this study. Gene expression data of healthy and epithelial ovarian cancer samples were used from the NCBI-GEO repository.

Patient consent for publication

No participation of patients directly used in this study. Hence, patient consent is not applicable.

Funding

No funding was received to carry out this research.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.