Bioinformatics analysis of RNA sequencing data reveals multiple key genes in uterine corpus endometrial carcinoma.

In the present study, the RNA sequencing (RNA-seq) data of uterine corpus endometrial carcinoma (UCEC) samples were collected and analyzed using bioinformatics tools to identify potential genes associated with the development of UCEC. UCEC RNA-seq data were downloaded from The Cancer Genome Atlas database. Differential analysis was performed using edgeR software. A false discovery rate <0.01 and |log2(fold change)|>1 were set as the cut-off criteria to screen for differentially expressed genes (DEGs). Differential gene co-expression analysis was performed using R/EBcoexpress package in R. DEGs in the gene co-expression network were subjected to Gene Ontology analysis using the Database for Annotation, Visualization and Integration Discovery. Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis was also performed on the DEGs using KOBAS 2.0 software. The ConnectivityMap database was used to identify novel drug candidates. A total of 3,742 DEGs were identified among the 552 UCEC samples and 35 normal controls, and comprised 2,580 upregulated and 1,162 downregulated genes. A gene co-expression network consisting of 129 DEGs and 368 edges was constructed. Genes were associated with the cell cycle and the tumor protein p53 signaling pathway. Three modules were identified, in which genes were associated with the mitotic cell cycle, nuclear division and the M phase of the mitotic cell cycle. Multiple key hub genes were identified, including cell division cycle 20, cyclin B2, non-SMC condensin I complex subunit H, BUB1 mitotic checkpoint serine/threonine kinase, cell division cycle associated 8, maternal embryonic leucine zipper kinase, MYB proto-oncogene like 2, TPX2, microtubule nucleation factor and non-SMC condensin I complex subunit G. In addition, the small molecule drug esculetin was implicated in the suppression of UCEC progression. Overall, the present study identified multiple key genes in UCEC and clinically relevant small molecule agents, thereby improving our understanding of UCEC and expanding perspectives on targeted therapy for this type of cancer.

Introduction

Endometrial cancer is the most commonly diagnosed female genital cancer and, in 2012, was ranked the fourth most common cancer in women worldwide, after breast, lung and colorectal cancer (1). In China, an estimated 63,400 new cases and 21,800 mortalities were reported for endometrial cancer in 2015 (2). Uterine corpus endometrial carcinoma (UCEC) is a common type of endometrial cancer. The incidence of UCEC increases with age and it is most frequently diagnosed in women aged between 45 and 65 years (3). Previous studies have sought to identify tumor biomarkers and have identified multiple UCEC biomarkers, including activated leukocyte cell adhesion molecule (4), sperm-associated antigen 9 (5), L1 cell adhesion molecule (6), progestogen-associated endometrial protein (7), heat shock protein family A (Hsp70) member 5 (8) and CD151 molecule (Raph blood group) (9). To enhance our understanding of UCEC and explore more effective and targeted therapies, the present study aimed to identify gene co-expression networks, hub genes and small molecule drugs associated with the development of UCEC.

The Cancer Genome Atlas (TCGA) is a comprehensive genomic database that holds data for >20 types of cancer obtained from thousands of patients. Data in TCGA includes whole-genome measurements of multiple genomic features, including DNA copy numbers, DNA methylation, and gene and microRNA expression, thereby assisting researchers in assessing cancer mechanisms at multiple molecular and regulatory levels (10). Furthermore, TGGA data is open access and available to all researchers in individual work settings.

In the present study, UCEC RNA sequencing data were collected and analyzed using bioinformatics tools. Differentially expressed genes (DEGs) were identified and a gene co-expression network was constructed, enabling the hub genes to be identified. Relevant small molecule drugs were also assessed. The results of the present study may provide novel insights into the pathogenesis of UCEC and thereby expand perspectives on the treatment of this type of tumor.

Materials and methods

Gene expression data

The RNA expression data (level 3) of UCEC patients and normal controls without UCEC were downloaded from the TCGA data portal up until December 2016 (http://cancergenome.nih.gov/) (11). A total of 552 UCEC samples and 35 normal control samples were included in the dataset. The exclusion criteria were as follows: A histological diagnosis not of UCEC; and samples without complete data for analysis.

DEG screening

The DEGs between UCEC and normal control tissues were screened using edgeR software (v..5; http://www.bioconductor.org/packages/release/bioc/html/edgeR.html). A false discovery rate <0.01 and |log2(fold change)|>1 were set as the cut-off values to identify the DEGs.

Cluster analysis

Bidirectional hierarchical clustering according to the expression of the DEGs was performed using the pheatmap package (v.1.0.8; https://CRAN.R-project.org/package=pheatmap) in R, as previously described (12). A heat map was used to represent the data (Fig. 1).

Figure 1.

Bidirectional hierarchical clustering of the 3,742 differentially expressed genes within the 587 samples (552 uterine corpus endometrial carcinoma and 35 control) in the dataset.

Gene co-expression network construction

Correlations among the DEGs were calculated using the R/EBcoexpress package in R (version 3.5; http://www.bioconductor.org/packages/release/bioc/html/EBcoexpress.html), as previously described (13). Genes with a correlation coefficient >0.9 were included in the gene co-expression network, which was visualized using Cytoscape software (version 3.5.1; http://www.cytoscape.org/).

Functional enrichment analysis

A Gene Ontology (GO) analysis (14) was performed on the DEGs in the gene co-expression network using the Database for Annotation, Visualization and Integration Discovery (http://david.abcc.ncifcrf.gov/) (15). P<0.05 was set as the cut-off criteria. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (16) was performed using KOBAS 2.0 software (http://kobas.cbi.pku.edu.cn/) (17), with P<0.05 set as the threshold.

Module analysis

Modules were identified using the MCODE plug-in in Cytoscape with the cut-off criteria of degree ≥2 and k-core ≥3. Each module was functionally annotated using the Cytoscape plug-in BiNGO based on hypergeometric distribution (adjusted P-value <0.01).

Screening of relevant small molecule drugs

ConnectivityMap (cMap) (18) is a public database containing >7,000 expression profiles representing 1,309 compounds (www.broad.mit.edu/cmap/). Relevant small molecule drugs were predicted using the cMap tool and those with a |score|>0.6 were included.

Results

DEGs among UCEC and control samples

The differential expression analysis identified 3,742 DEGs, comprising 2,580 upregulated and 1,162 downregulated genes. The results of bidirectional hierarchical clustering of the 3,742 DEGs in the 587 samples (552 UCEC and 35 normal control samples) that were screened were provided (Fig. 1). The UCEC and control samples differed in gene expression pattern, suggesting that DEGs could distinguish between the two types of sample.

GO annotations of DEGs

The GO annotations of the DEGs were provided (Fig. 2). Cancer-associated processes, including cell adhesion, ion transport and biological adhesion, were among the significantly associated terms.

Figure 2.

Gene Ontology annotations of the differentially expressed genes.

Gene co-expression network

Correlated genes with a correlation coefficient >0.9 were included in the gene co-expression network (Fig. 3). A total of 129 DEGs (nodes) and 368 edges (lines between nodes) were included, divided into 110 upregulated and 19 downregulated genes. The sizes of the nodes were proportional to the number of genes in the gene set, and a number of hub genes (nodes with degrees ≥8) were highlighted, including cell division cycle 20 (CDC20), cyclin B2, non-SMC condensin I complex subunit H, BUB1 mitotic checkpoint serine/threonine kinase, cell division cycle-associated 8, maternal embryonic leucine zipper kinase, MYB proto-oncogene like 2, TPX2, microtubule nucleation factor and non-SMC condensin I complex subunit G.

Figure 3.

Gene co-expression network of differentially expressed genes. Red, upregulated genes; green, downregulated genes.

Co-expression network functional enrichment analysis

GO enrichment analysis identified 20 significantly overrepresented terms for gene function in the co-expression network (Table I), including M phase, nuclear division, mitosis and the M phase of the mitotic cell cycle. KEGG pathway enrichment analysis indicated that pathways such as those of the cell cycle, oocyte meiosis, progesterone-mediated oocyte maturation and tumor protein p53 signaling were significantly enriched among genes in the co-expression network (Table II).

Table I.

GO terms significantly overrepresented in the genes from the gene co-expression network.

GO category	GO term	Count	P-value	FDR
GO:0000279	M phase	45	1.27×10−46	1.93×10−43
GO:0000280	Nuclear division	40	5.98×10−46	9.07×10−43
GO:0007067	Mitosis	40	5.98×10−46	9.07×10−43
GO:0000087	M phase of the mitotic cell cycle	40	1.27×10−45	1.93×10−42
GO:0048285	Organelle fission	40	3.20×10−45	4.85×10−42
GO:0022403	Cell cycle phase	46	1.33×10−43	2.02×10−40
GO:0000278	Mitotic cell cycle	43	3.22×10−41	4.89×10−38
GO:0007049	Cell cycle	53	2.40×10−40	3.63×10−37
GO:0022402	Cell cycle process	48	3.05×10−40	4.62×10−37
GO:0051301	Cell division	35	4.17×10−33	6.33×10−30
GO:0005819	Spindle	24	8.32×10−27	9.68×10−24
GO:0015630	Microtubule cytoskeleton	35	2.86×10−26	3.33×10−23
GO:0007059	Chromosome segregation	18	3.26×10−21	4.94×10−18
GO:0044430	Cytoskeletal part	37	1.25×10−20	1.45×10−17
GO:0005856	Cytoskeleton	41	5.11×10−19	5.94×10−16
GO:0043228	Non-membrane-bounded organelle	52	5.37×10−18	6.24×10−15
GO:0043232	Intracellular non-membrane-bounded organelle	52	5.37×10−18	6.24×10−15
GO:0000775	Chromosome, centromeric region	17	2.29×10−17	2.67×10−14
GO:0000793	Condensed chromosome	17	4.39×10−17	5.10×10−14
GO:0000779	Condensed chromosome, centromeric region	14	8.86×10−17	1.33×10−13

GO, Gene Ontology; FDR, false discovery rate.

Table II.

Kyoto Encyclopedia of Genes and Genomes pathways significantly overrepresented in the genes from the gene co-expression network.

Term	Count	P-value
hsa04110: Cell cycle	13	3.65×10−16
hsa04114: Oocyte meiosis	11	3.06×10−13
hsa04914: Progesterone-mediated oocyte maturation	7	3.60×10−8
hsa04115: p53 signaling pathway	4	8.02×10−5
hsa05166: HTLV-I infection	5	1.41×10−3
hsa04270: Vascular smooth muscle contraction	3	7.17×10−3
hsa04068: FoxO signaling pathway	3	9.19×10−3

HTLV-I, human T-cell lymphotropic virus type 1; FoxO, forkhead box O; p53, tumor protein 53.

Modules and functions

Three modules (designated A-C) were identified from the gene co-expression network (Fig. 4). Module A included 11 DEGs implicated in the mitotic cell cycle and cell cycle phases. Module B contained 6 DEGs associated with nuclear division. Module C comprised 7 DEGs associated with the M phase of the mitotic cell cycle (Table III).

Figure 4.

Modules identified from the gene co-expression network.

Table III.

Functional terms of the three modules.

GO-ID	P-value	FDR	No. of genes
Module A
GO:0000279~M phase	9.60×10−7	1.12×10−3	6
GO:0000278~mitotic cell cycle	1.71×10−6	2.01×10−3	6
GO:0022403~cell cycle phase	2.98×10−6	3.50×10−3	6
GO:0007067~mitosis	8.04×10−6	9.43×10−3	5
GO:0000280~nuclear division	8.04×10−6	9.43×10−3	5
GO:0000087~M phase of mitotic cell cycle	8.64×10−6	1.01×10−2	5
GO:0048285~organelle fission	9.43×10−6	1.11×10−2	5
GO:0022402~cell cycle process	1.37×10−5	1.60×10−2	6
GO:0007059~chromosome segregation	1.69×10−5	1.98×10−2	4
GO:0000793~condensed chromosome	3.41×10−5	3.10×10−2	4
GO:0005819~spindle	5.04×10−5	4.59×10−2	4
GO:0005694~chromosome	5.32×10−5	4.83×10−2	5
Module B
GO:0000280~nuclear division	3.36×10−7	3.63×10−4	5
GO:0007067~mitosis	3.36×10−7	3.63×10−4	5
GO:0000087~M phase of mitotic cell cycle	3.61×10−7	3.90×10−4	5
GO:0048285~organelle fission	3.95×10−7	4.27×10−4	5
GO:0000279~M phase	1.69×10−6	1.82×10−3	5
GO:0000278~mitotic cell cycle	2.69×10−6	2.93×10−3	5
GO:0015630~microtubule cytoskeleton	3.37×10−6	2.94×10−3	5
GO:0022403~cell cycle phase	4.22×10−6	4.56×10−3	5
GO:0022402~cell cycle process	1.46×10−5	1.57×10−2	5
Module C
GO:0000280~nuclear division	3.36×10−7	3.63×10−4	5
GO:0007067~mitosis	3.36×10−7	3.63×10−4	5
GO:0000087~M phase of mitotic cell cycle	3.61×10−7	3.90×10−4	5
GO:0048285~organelle fission	3.95×10−7	4.27×10−4	5
GO:0000279~M phase	1.69×10−6	1.82×10−3	5
GO:0000278~mitotic cell cycle	2.69×10−6	2.91×10−3	5
GO:0015630~microtubule cytoskeleton	3.37×10−6	2.94×10−3	5
GO:0022403~cell cycle phase	4.22×10−6	4.56×10−3	5
GO:0022402~cell cycle process	1.46×10−5	1.57×10−2	5

GO-ID, Gene Ontology identification; FDR, false discovery rate.

Relevant small molecule drugs

A total of 7 small molecule drugs were identified, including esculetin, antazoline and isometheptene. Esculetin had the highest negative connectivity score compared with the 6 other small molecule drugs (−0.844; Table IV).

Table IV.

Relevant small molecule drugs.

cMap name	Correlation	P-value
Esculetin	−0.844	7.55×10−3
Antazoline	−0.729	1.07×10−2
Isometheptene	−0.722	1.20×10−2
Oxamniquine	0.717	1.31×10−2
Pyrimethamine	−0.642	1.49×10−2
Carmustine	−0.795	1.75×10−2
Cefapirin	0.636	4.28×10−2

cMap, connectivity map.

Discussion

The results of the present study indicated that multiple DEGs are associated with UCEC pathophysiology. Matrix metalloproteinases (MMPs) serve an important function in degrading the extracellular matrix and basement membrane components (19). Karahan et al identified that MMP-9 was expressed in a high percentage of primary endometrial carcinoma tissues and that MMP-9 expression was potentially associated with parameters of tumor aggressiveness (20). Furthermore, MMP-9 has been revealed to be overexpressed in UCEC and associated with UCEC progression (21). Topoisomerase 2α, a key enzyme in DNA replication, has been implicated as a predictive marker in treating endometrial cancer with taxane-containing adjuvant chemotherapy (22). Enhancer of zeste homolog 2 (EZH2), a critical component of the polycomb repressive complex 2, is associated with cell proliferation, invasion, adhesion and metastasis in several cancer types (23), and is overexpressed in numerous types of cancer, including breast cancer (24), colorectal cancer (25), ovarian cancer (26) and nasopharyngeal cancer (27). A previous study indicated that EZH2 may predict a more aggressive endometrial carcinoma and may therefore potentially represent a therapeutic target for this type of cancer (28). Baculoviral inhibitor of apoptosis (IAP) repeat containing 5 gene (BIRC5), a member of the IAP gene family, may inhibit the activation of caspase to negatively regulate apoptosis (29). In endometrial cancer, increased BIRC5 expression has been associated with poor progression-free survival and may serve as an independent prognostic factor (30).

It is well established that CDC20 is an essential developmental gene; disrupting CDC20 function in mice induces mortality in embryos (31). Previous studies have also suggested that CDC20 functions oncogenically in human tumorigenesis. CDC20 serves important functions in the cell cycle (32), cell apoptosis (33) and in targeting downstream substrates for ubiquitination and subsequent degradation (34). Furthermore, CDC20 has been revealed to be upregulated and may represent a useful prognostic biomarker in endometrial cancer (35).

In the present study, esculetin was identified as a significant small molecule drug in the development of UCEC. Esculetin (6,7-dihydroxycoumarin), a coumarin derived from natural plants, inhibits cancer cells from proliferating and induces multiple types of human cancer cells to apoptose, including breast (36), colon (37,38) and gastric cancer (39), malignant melanoma (40), hepatocellular carcinoma (41,42) and oral squamous cancer (43,44). Esculetin has also been reported to assist in inducing cancer cells to apoptose (45,46). In the present study, esculetin had a higher negative connectivity score compared with the 6 other small molecule drugs, suggesting that esculetin could inhibit UCEC from progressing. To the best of our knowledge, no study has yet reported on the association between esculetin and endometrial cancer. Verifying whether esculetin represents a novel therapeutic agent in treating UCEC requires further studies. Overall, the present study identified multiple key genes in UCEC and clinically relevant small molecule agents, thereby improving our understanding of UCEC and expanding perspectives on targeted therapy for this type of cancer.