Integrative bioinformatics approaches to map key biological markers and therapeutic drugs in Extramammary Paget's disease of the scrotum.

Extramammary Paget's disease (EMPD) is an intra-epidermal adenocarcinoma. Till now, the mechanisms underlying the pathogenesis of scrotal EMPD is poorly known. This present study aims to explore the knowledge of molecular mechanism of scrotal EMPD by identifying the hub genes and candidate drugs using integrated bioinformatics approaches. Firstly, the microarray datasets (GSE117285) were downloaded from the GEO database and then analyzed using GEO2R in order to obtain differentially expressed genes (DEGs). Moreover, hub genes were identified on the basis of their degree of connectivity using Cytohubba plugin of cytoscape tool. Finally, GEPIA and DGIdb were used for the survival analysis and selection of therapeutic candidates, respectively. A total of 786 DEGs were identified, of which 10 genes were considered as hub genes on the basis of the highest degree of connectivity. After the survival analysis of ten hub genes, a total of 5 genes were found to be altered in EMPD patients. Furthermore, 14 drugs of CHEK1, CCNA2, and CDK1 were found to have therapeutic potential against EMPD. This study updates the information and yields a new perspective in the context of understanding the pathogenesis of EMPD. In future, hub genes and candidate drugs might be capable of improving the personalized detection and therapies for EMPD.

1. Introduction

Extramammary Paget disease (EMPD) arises as rare cutaneous adenocarcinoma, accounting small proportion (about 6.5%) of all Paget’s disease cases. The very first time, EMPD was uncovered by Radcliffe Crocker about half a century ago [1-4]. The frequency of occurrence is higher at various anatomical locations such as the scrotum, vulva, and perianal region [5, 6]. In women, the vulvar EMPD accounts for 81.3% of all EMPD cases while scrotal EMPD accounts for 43.2% of all EMPD cases in men [7]. The risk of getting EMPD is increasing 3.2% per year. However, it is unknown why EMPD cases of the scrotum and other anatomic sites are increasing [8, 9]. The overall frequency of hospitalization is common under the age of 65 years [10]. It has a great diversity of clinical symptoms that vary considerably in the context of disease status and prognosis [10, 11]. However, owing to insufficient diagnostic methods, EMPD patients are diagnosed, on average, at the middle or late disease stage, which consequently, leads to a poor prognosis [12, 13]. Hence, the understanding of molecular mechanisms that contributes to the pathogenesis and prognosis of EMPD has become increasingly important for the development of multiple therapeutics and diagnostics approaches.

Until now, EMPD is incurable and dramatically increases the morbidity as well as the mortality rate in developed countries. As the information pileup, the same scenario is predicted in the developing countries [14]. Multiple remedies such as imiquimod, 5% cream, or cytotoxic agent with 1% 5-fluorouracil cream [15, 16] are available but all these revolve around the management of disease, however, the actual treatment is yet to be discovered.

The discovery of potential biomarkers that can halt the pathophysiology of the disease and can act as a virtual shortcut, will be considered as the miracle of the current era. The mind-boggling potential benefits of molecular biomarkers offer multiple innovative perspectives to improve diagnostic as well as treatment options. Several studies were conducted to identify the potential biomarkers for monitoring and diagnosis of EMPD, but all these studies might conflict due to the complicated molecular mechanism of EMPD [17, 18].

Nowadays, the use of bioinformatics is getting popular across all facets of life sciences. Recently, it has been seen as an outbreak of emerging sequencing technologies that enable researchers to make ground-breaking discoveries in the domain of computational biology. In recent decades, bioinformatics along with microarray technologies has paved the way for researchers to identify disease-related genes involved in the pathogenesis of EMPD. Various bioinformatics-related researches on plenty of human diseases had proven reliable and persuasive, so it implies that integrated bioinformatics analysis can contribute to evaluating the underlying complex molecular mechanism of EMPD of the scrotum.

In spite of the numerous studies on EMPD, no sufficient evidence is present yet to prove the existence of disease-related genes, and their involvement in the pathogenesis of EMPD. To tackle this issue, we used integrated bioinformatics approaches to figure out the disease-related functional genes as problem-solving negotiators to switch off the progression of EMPD. Moreover, investigation of drug-genes interactions in the present work can contribute to the discovery of therapeutic candidates for drug repurposing.

2. Materials and methods

2.1. Retrieval of microarray

Gene Expression Omnibus (GEO) database in National Center for Biotechnology Information (NCBI) is a freely available public database, enclosing the gene profiles [19]. Expression profile of GSE117285 based on GPL17077 platform (Agilent-039494 SurePrint G3 Human GE v2 8x60K array) was retrieved from NCBI-GEO database [20]. Dataset include 4 controls and 4 EMPD cases.

2.2. Identification of DEGs

GEO2R is considered as a useful interactive web tool for comparison and analysis of sample of multiple groups [21]. In present study, GEO2R tool was employed to explore the Differential Expressed Genes (DEGs) among controls and EMPD patients of GSE117285 dataset. Genes that satisfy the criteria of |log2 (fold change) | > 1.0 and adjusted P-value < 0.05 were distinguished as DEGs. Volcano plot was constructed using ggplot2 package of R, to visualize the significant and non-significant DEGs.

2.3. Analysis of DEGs at functional level

At the functional level, a database for annotation, visualization, and integrated discovery (DAVID) was used to perform GO enrichment analysis and KEGG pathway analysis [22]. The DEGs were subjected to DAVID for the prediction of DEGs function at three levels: Biological process (BP), Molecular function (MF), and Cellular component (CC). The top 10 significant items of GO and KEGG pathways were demonstrated in form of bubble maps which were plotted using ggplot2 package of R. In this regard, P<0.05 was believed to be statistically-significant.

2.4. Construction of protein-protein interaction network

Search Tool for the Retrieval of Interacting Genes/Proteins(STRING) was used for the functional interactions among DEGs at the combined score of more than 0.5 [23]. Molecular Complex Detection (MCODE) plugin from Cytoscape was utilized for distinguishing the module that best represents the clusters of DEGs [24]. The modules which met the cutoff criteria of node > 10 and the score >10 were selected. Further, the resulted modules were subjected to DAVID for pathway analysis.

2.5. Selection of hub genes

CytoHubba plugin was used for distinguishing hub genes among up and down expressed DEGs. A total of twelve topological analysis methods are available in the cytoHubba. Among these 12 methods, the degree method was chosen for the identification of hub genes.

2.6. Survival analysis of hub genes

EMPD is a malignant skin cancer characterized by the appearance of Paget cells. GEPIA server was employed to evaluate the expression level of hub genes among Skin Cutaneous Melanoma (SKCM) and normal tissues [25]. In this regard, P-value< 0.05 was believed to be statistically significant.

2.7. Drug-gene interaction analysis

Drug-Gene Interaction Database was used to search hub genes against known and potentially small drug-like molecules [26]. The final list was included only FDA-approved drugs. Finally, the STITCH tool was used for the visualization of interaction network among gene and their associated drug candidates [27].

3. Results

3.1. Identification of DEGs

In the present study, microarray dataset GSE117285 was obtained from Agilent-039494 SurePrint G3 Human GE v2 8x60K Microarray, which consists of 8 samples (4 healthy controls and 4 EMPD patients). A total of 786 DEGs was identified, of which 619 were upregulated genes and 167 were downregulated genes (. Moreover, gene differential expressions were presented in form of a volcano map ().

Representation of differential expressed genes in form of Volcano plot with red dots showing upregulated genes while blue dots showing downregulated genes.

Black dots represent non-significant genes.

3.2. Functional enrichment analysis of DEGs

GO enrichment and KEGG pathways analysis of DEGs were performed to analyze the gene function in terms of biological processes, cellular components, and molecular function as well as their associated pathways. GO enrichment analysis of the top 10 significantly enriched terms showed that in the BP category, the genes involved are concerned with cell division, mitotic nuclear division, cell proliferation, sister chromatid cohesion, organ regeneration, and inflammatory responses (). In terms of CC, the genes were enriched nucleus, cytosol, nucleoplasm, extracellular space, and centrosome (). For MF, a category the genes were mainly concentrated in the protein binding, ATP binding, protein kinase activity, protein serine/threonine kinase activity, and chromatin binding (). KEGG enrichment pathway analysis revealed that genes were significantly enriched in the cell cycle, viral carcinogenesis, Hepatitis B, TNF signaling pathway, and oocyte meiosis ().

Representation of GO enrichment and KEGG pathway analysis.

(A) Gene ontology in terms of BP (B) Gene ontology in terms of CC (C) Gene ontology in terms of MF(D) KEGG pathway analysis.

3.3. Construction of PPI network and the analysis of DEGs

PPI network of DEGs obtained from STRING were subjected to the MCODE plugin of Cytoscape which provided significant 22 modules. From these modules, the top two functional clusters of modules were selected based on the cutoff criteria of node > 10 and the score is ≥10 (). KEGG pathway analysis of the selected modules revealed that the genes belong to these clusters are enriched in the TNF signaling pathway, toll-like receptor signaling pathway, chemokine signaling pathway, cytokine and cytokine receptor interaction, and complement and coagulation cascades ().

Two modules were selected having cutoff criteria node > 10 and the score is ≥.

10 (A) The first module constructed from MCDOE comprised of 55 genes; (B) Pathways associated with the first module; (C) The second module constructed from MCDOE comprised of 14 genes; (D) Pathways associated with the second modules.

3.4. Selection of hub genes

Using 12 methods available in the cytoHubba, the topmost ten genes were selected and ranked by degree method. These ten genes named CDK1, MAD2L1, BUB1, CCNA2, BUB1B, NCAPG, NDC80, KIF2C, CENPF, and CHEK1 were considered as the hub genes (). Moreover, the interaction network of hub genes to their related neighboring genes is shown in

(A) Construction of PPI network of 6 hub genes (B) Observed coexpression of 10 hub genes in Homo sapiens.

3.5. Survival analysis of hub genes

GEPIA server was used for the survival analysis of 10 hub genes. The survival analysis of 10 hub genes revealed 5 significant genes (). These 5 genes named CHEK1, KIF2C, CENPF, CCNA2, and CDK1 were demonstrated as significant because of the p-value <0.05. Moreover, the expression level of 5 genes varied significantly among SKCM and normal tissues.

3.6. Drug-gene interaction

A total of 14 drugs were explored using DGIdb that might have the potential to treat PD patients. CHEK1, CCNA2, and CDK1 were chosen as possible targets of 14 drugs based on the DrugBank source (). Furthermore, using the STITCH tool, downstream interaction networks of CHEK1, CCNA2, and CDK1 were generated ().

(A) Targetable CHEK1 subnetwork. (B) Targetable CCNA2 subnetwork. (C) Targetable CDK1 subnetwork.

4. Discussion

EMPD is considered a rare malignant dermatosis with poorly defined outcomes [28, 29]. Clinically, the more frequent indication of EMPD is eczematous [30]. There has been substantial heterogeneity regarding various EMPD cases that lead to the earlier detection, a thought-provoking question [31-33]. Recently, a comprehensive study on patients suffering from EMPD has made it clear that progress in the overall pathogenesis of EMPD varies from patient to patient which turned the diagnosis as well as treatment of EMPD into a challenging task [34]. The current work planned to identify the disease-related functional genes involved in the progression of EMPD in the scrotum. This whole research revolves around the analysis of gene ontology, gene enrichment pathways, PPI, and hub genes. The gene expression microarray profiles of the GSE117285 dataset were analyzed using integrated bioinformatics analysis. Through KEGG pathway analysis, the DEGs were significantly found to be enriched in the cell cycle, viral carcinogenesis, Hepatitis B, TNF signaling pathway, oocyte meiosis, RIG I like receptor signaling pathway, Legionellosis, and NOD-like receptor signaling pathway. Our functional annotation of target genes might be helpful in understanding this targeted slicing on the pathogenesis of EMPD. In the current work, after the survival analysis of ten hub genes, a total of 5 genes were found to be altered in EMPD patients involving CHEK1, KIF2C, CENPF, CCNA2, and CDK1. Hence, it represents that these genes plays important role in the pathogenesis of EMPD. Furthermore, 14 drugs of CHEK1, CCNA2, and CDK1 were found to have therapeutic potential against EMPD.

CHK1 kinase is considered a key element of DNA damage response. There is no study available regarding the role of CHK1 in EMPD. Apoptosis is caused when DNA damage is extremist. Recently a study has shown that Chk1 inhibition encourages cancerous cells to undergo apoptosis and inappropriately arrest damaged DNA [35]. CHK1 plays important role in tumor suppression hence might prove fruitful to combat the disease condition by preventing EMPD from becoming malignant. By targeting CHK1, the pathogenic mechanisms of EMPD can be controlled, hence might serve as a molecular biomarker for the diagnosis and treatment of EMPD.

Cyclin-dependent kinase 1(CDK1) is a key member of the CDKs family, crucial for the cell cycle progression [36, 37]. During the last ten years, a slew of studies has made it clear that CDK1 encourages tumor growth as well as enables cancerous cells to proliferate spontaneously [38, 39]. According to previous reports, CDK1 has been implicated in the development of various cancers, such as liver and colorectal cancer [40-42]. EMPD is considered a rare cancer-related disease and found to have an association with skin cancer [43]. Hence all this evidence provides a precious clue that upregulation of CDK1 might control the intricate molecular mechanism behind the progression of EMPD.

CCNA2 is considered a therapeutic target in skin cancer [44]. Multiple studies on CCNA2 has made it clear that upregulation of CCNA2 interfere the normal process and leads to proliferation of cancerous cell spontaneously [45]. During the process of the cell cycle, CCNA2 regulates G1/S as well as G2/M transitions. It was discovered to be a prognostic biomarker for melanoma survival [44, 46]. Hence by considering CCNA2 as a molecular biomarker, the diagnosis and treatment of EMPD might become an easy task. Considering our analysis, we propose that the upregulation of CCNA2 might induce EMPD.

Identification of aberrant pathways in EMPD patients might reveal the molecular mechanism underlying and uncover more enthralling and promising molecular candidates with effective diagnostic and prognostic value. KEGG pathway analysis of CHEK1, CCNA2, and CDK1 revealed that this gene disturbs the cellular pathways through disruption in cell cycle and p53 signaling pathway which unfortunately leads to the disease condition. It is noteworthy that CHEK1, CCNA2, and CDK1 are involved in cell cycle regulation. Cancer is an unregulated cell division caused by the failure in cell cycle regulation. EMPD is malignant skin cancer, hence these findings conclude that EMPD interferes with the regulatory mechanism of the cell cycle and contributes to the development of skin cancer. These findings shed light on the pathogenesis of EMPD and facilitate the development of personalized treatment for EMPD. The disturbed pathways identified using integrated bioinformatics analysis may have an important role to play in the pathogenesis of EMPD. Additional studies are required to investigate the molecular mechanisms behind these aberrant pathways and EMPD development.

5. Conclusion

In the present work, a mechanism was proposed to explain the progression in the pathogenesis of EMPD. It might due to the genes that disturb the cellular pathways which ultimately leads to the disease condition. This research discerned hub genes as key biological markers and their associated pathways involved in the development of EMPD. Based on our knowledge, CHEK1, CCNA2, and CDK1 have not been previously reported to be related to EMPD. It was found that these genes might act as potential biomarkers for the diagnosis of EMPD at an early stage. Our findings reveal that EMPD causes disruption in cellular pathways which unfortunately makes the disease condition much worse. In the near future, further study and clinical trials are required for the identification of genes and small drug-like molecules having effective diagnostic and prognostic value, respectively. Our research will serve as a significant pioneer for the researchers who want to identify the associated pathways involved in the pathogenesis of EMPD. This research relies on various freely available databases to shed light on EMPD pathogenesis and treatment. In vivo and in vitro investigations of genes and pathways interaction are essential to delineate the specific roles of the identified genes, which could confirm gene functions and reveal the mechanisms underlying EMPD. In the future, additional experimental research on these hub genes could lead to an increase in our knowledge to fight against EMPD by means of novel therapeutic approaches. Based on the hub genes, experimental models may be designed in terms of the detection of pathogenesis, evaluation of risk, and determining the targeted therapies of EMPD.

List of upregulated and downregulated DEGs.

(DOCX)

Click here for additional data file.

Table 1

Top 2 modules were selected having cutoff criteria node > 10 and the score is ≥ 10.

Clusters	Score	Nodes	Edges	Nodes IDs
1	49.296	55	1331	MELK, CCNA2, PLK4, SMC4, KIF23, NCAPG, DLGAP5, MND1, FBXO5, TYMS, RAD51AP1, AURKA, TRIP13, ARHGAP11A, ATAD2, CENPK, CENPN, CENPA, HJURP, ASPM, CEP55, CENPF, PBK, FANCI, KIF2C, CDK1, KIF4A, ECT2, EXO1, KNTC1, HMMR, ERCC6L, HELLS, CDC45, SKA3, FAM83D, BUB1, CHEK1, BIRC5, MAD2L1, CDKN3, BUB1B, STIL, CDCA2, SKA1, SPC25, SPAG5, CENPU, CKS2, CDCA5, CENPE, KIF18A, NDC80, KIAA0101, RAD54L
2	10.615	14	69	MRPL13, RPS4X, MAGOHB, RPL22L1, RPS2, MRPL1, RPS11, PSMA1, RPL7A, EIF2S2, EIF4E, RPL13A, RPS6, RSL24D1

Table 2

List of FDA approved drugs.

Search_ term	Match_ term	Match_ type	Gene	Drug	Interaction_types	Sources	Pmids
CHEK1	CHEK1	Definite	CHEK1	PALBOCICLIB		DTC
CHEK1	CHEK1	Definite	CHEK1	CISPLATIN		CIViC	28490518
CHEK1	CHEK1	Definite	CHEK1	OLAPARIB		CIViC	28490518
CHEK1	CHEK1	Definite	CHEK1	GEMCITABINE		NCI	17245119
CHEK1	CHEK1	Definite	CHEK1	FOSTAMATINIB	inhibitor	DrugBank	26516587
CHEK1	CHEK1	Definite	CHEK1	ETOPOSIDE		DTC	22364746
CCNA2	CCNA2	Definite	CCNA2	TAMOXIFEN		PharmGKB	24622579
CDK1	CDK1	Definite	CDK1	SERTRALINE		DTC	16680159
CDK1	CDK1	Definite	CDK1	RUCAPARIB		DTC
CDK1	CDK1	Definite	CDK1	FENOFIBRATE		DTC	16680159
CDK1	CDK1	Definite	CDK1	CINNARIZINE		DTC	16680159
CDK1	CDK1	Definite	CDK1	CLOFIBRATE		DTC	16680159
CDK1	CDK1	Definite	CDK1	FOSTAMATINIB	inhibitor	DrugBank	26516587
CDK1	CDK1	Definite	CDK1	CLOTRIMAZOLE		DTC	16680159