Identification of key modules and prognostic markers in adrenocortical carcinoma by weighted gene co-expression network analysis
YONG ZOU and LUANLIAN JING
Department of Oncology, The People’s Hospital of Hanchuan, Hanchuan, Hubei 431600, P.R. China
Received January 24, 2019; Accepted June 12, 2019
DOI: 10.3892/ol.2019.10725
Abstract. Adrenocortical carcinoma (ACC) is a rare and aggressive cancer with a high relapse rate and limited treatment options. Therefore, the identification of potential prognostic markers in patients with ACC may improve early detection, survival rates and may additionally provide novel insights into the early detection of recurrence. In the present study, clinical traits and RNA-seq data of 79 patients with ACC were obtained from The Cancer Genome Atlas (TCGA). Weighted gene co-expression network analysis was carried out and 17 distinct co-expression modules were built to examine the association between the modules and the clinical traits. Of the 17 modules, two co-expression modules, which contained 214 and 168 genes, were significantly correlated with two clinical traits, tumor stage and vital status. Functional enrichment analysis was performed on the selected modules. The results showed that one of the modules was primarily enriched in cell division and the other module was enriched in metabolic pathways, suggesting their involvement in tumor progression. Furthermore, cyclin dependent kinase 1 (CDK1) and ubiq- uitin C (UBC) were identified as hub genes in both modules. Survival analysis revealed that the high expression of the hub genes significantly correlated with the poor survival rate of patients, suggesting that CDK1 and UBC have vital roles in the progression of ACC. In the present study, a co-expression gene module of ACC was provided and the prognostic genes that may serve as new diagnostic markers in the future were defined.
Introduction
The adrenal glands reside above the kidneys and produce multiple hormones essential for development (1). Each gland consists of
an inner medulla and an outer cortex that produce catechol- amines and steroid hormones, respectively (2). Adrenocortical carcinoma (ACC) is a rare tumor of the adrenal cortex with an estimated incidence of 1 patient per million/year (3,4). There is a higher prevalence of ACC in females and an increased incidence in the first and fourth to fifth decades of life (5). There are no notable clinical phenotypic characteristics in patients with ACC during the early stage and the majority of patients are diagnosed with advanced stage ACC in the first instance (6). Patients with ACC do not respond favorably to chemotherapy and radio- therapy (7) and the patients frequently have to undergo surgical resection where the 5-year survival rate is >35%. Mitotane (o, p’-dichlorodiphe nyldichloroethane) has been used since the 60s for treating patients with ACC, despite its toxicity and narrow therapeutic index (8,9). An improved understanding of the genes associated with ACC may improve treatment options by identi- fying potential therapeutic targets. Weighted gene co-expression network analysis (WGCNA) is a frequently used method to explore the association between genes and phenotypes (10). Gene expression data are transformed into co-expression modules and provide insights into the signaling networks that may underlie certain phenotypes. WGCNA is widely used to improve understanding of various biological processes such as cancer and its progression (11,12). Yang et al (13) identified candidate biomarkers and molecular mechanisms involved in glioblastoma multiforme using WGCNA (13). WGCNA compares differentially expressed genes and identifies key interactions among different co-expression modules (12).
Next generation sequencing is used to detect genomic alterations which could be used to guide targeted therapies for treating patients with ACC and several targets have been discov- ered (14,15). However, the molecular diagnostic parameters are still not entirely known and there are only small number of studies that have cataloged relevant expression modules in patients with ACC, which has limited understanding of the disease and its mechanisms (16-18). Using WGCNA, it is possible to construct gene networks and analyze the connectivity between genes and clinical traits (19). The master regulators identified in the gene regulation network will typically exhibit important functions.
In this study, it was hypothesized that distinct ACC co-expression modules are associated with different clinical outcomes and the highly connected genes in the modules can represent the biological feature of this module and have the potential as prognosis markers. Normalized-RNA-seq
Correspondence to: Dr Luanlian Jing, Department of Oncology, The People’s Hospital of Hanchuan, 33 Huanle Street, Hanchuan, Hubei 431600, P.R. China
E-mail: luanlianjings@163.com
Key words: adrenocortical carcinoma, prognostic genes, weighted gene co-expression network analysis, co-expression modules
and clinical data was downloaded from the TCGA database of 79 patients with ACC at different stages. The candidate mRNAs related to tumor progression were identified by co-expression analysis. Furthermore, 2,472 differentially expressed genes from ACC and normal tissues were down- loaded from the Gene expression profiling interactive analysis (GEPIA) website (20), and ACC and ACC co-expression modules were constructed (20). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrich- ment analyses were performed on the selected modules, and the hub gene in each module was identified, which assisted in determining the primary function(s) of the genes in each module. These findings may play a significant role in recog- nizing the malignant potential of these genes as well as the prognosis of patients with ACC.
Materials and methods
Patients and TCGA data retrieval. The clinical and gene expression data of 79 patients with ACC were obtained from TCGA (https://tcga-data.nci.nih.gov/tcga/) (21), where the expression profile was obtained based on the IlluminaHiSeq RNA-seq platform (Illumina, Inc.). The Bioconductor (22) packages TCGAbiolinks (23) based on the R software (24) (version 3.4.0) were used to download and process the data collected from TCGA. These data included the age, sex, survival status and the cancer stage of the patients in addition to the vital status. The DEGs in normal and ACC tissues were obtained from GEPIA (http://gepia.cancer-pku.cn/) (20). The screening standards for the identification of DEGs were a log2 fold-change (log2FC)>1 and Q<0.01.
WGCNA. Co-expression networks were constructed with the identified DEGs. Pearson’s correlation coefficients were calculated for all the genes in the dataset and the correlation matrix of the entire gene dataset was obtained. The power B was used to remove weakly correlated genes, while retaining the strongly correlated ones. The process produced an adja- cent matrices weighted network that was converted into a topological overlapping matrices (TOM) network as previ- ously described (10,25). After constructing the TOM network, hierarchical clustering was used to generate a cluster dendro- gram with branches corresponding to the gene co-expression modules. The WGCNA (10) algorithm was used to identify the co-expression modules. WGCNA was implemented using R (10). The TOM representing the overlap in shared neighbors and the soft thresholding power were calculated according to a previous study (25). Colored heatmaps were used to analyze the strength of the interactions.
Construct module-trait associations of ACC. Module-trait associations were estimated using the association between the module signature and the phenotype (clinical traits) as previously described (12), thus allowing easy identification of the expression set (module) that strongly correlated with the phenotype. For each expression profile, the gene significance (GS) was calculated as the absolute value of the association between the expression profile and each clinical trait, whereas the module membership (MM) was defined as the correlation between the expression profile and each module signature.
Functional analysis of network module genes. To identify the underlying biological mechanisms responsible for the progres- sion of ACC, the DEGs derived from the brown and yellow modules were used for GO and KEGG pathway enrichment analyses, because the two modules are both correlated with the patient clinical traits. The DEGs annotated in the GO database were used to classify the GO functions in the ClusterProfiler package (26), and the DEGs for KEGG enrichment analysis were mapped to the KEGG database. After the enrichment analysis, the significant KEGG pathway and GO terms were selected according to the cut-off criterion of adjusted P<0.05.
Protein-protein-interaction (PPI) network construction. The PPI data was retrieved from a previous study, which contained protein interaction associations from 15 databases (27). The overall PPI network was based on 16,081 nodes and 231,633 interactions in these databases. The genes involved in the brown and yellow modules were mapped to the overall PPI network to obtain the module specific interaction network. Furthermore, the degree of each node was calculated using the IGRAPH package (28). The widely linked genes (hub genes) in the PPI network were more closely related to most proteins.
Determination of the receiver operating characteristic (ROC) curves and the area under the ROC curve (AUC). To validate whether the mRNA levels of the hub genes exhibit excellent diagnostic efficiency for distinguishing the tumor tissues from the normal tissues, the ROC curve analysis was performed. Specifically, the normalized the mRNA expression profile both for the normal adrenal gland tissue and the ACC tissue were downloaded from Recount2 database (29), including 159 normal adrenal gland samples and 79 ACC samples respectively. The ACC samples were labeled as ‘positive’ and the normal adrenal gland sample were labeled as ‘nega- tive’. Subsequently, the ROCR (30) package based on the R (version 3.5.2) software system was used to plot the ROC curve and calculated the AUC for the hub genes.
Survival analysis. Survival analysis for all genes in the brown and yellow modules was performed using the R survival package (https://CRAN.R-project.org/package=survival; version 2.41-3). A log-rank test was performed to determine whether the expression of a gene correlated with the overall survival. For the overall survival rate, a log-rank test was used to test for significance in the univariate analysis between each subgroup. Unless otherwise specified, P<0.05 was considered to indicate a statistically significant difference.
Results
Pre-processing of datasets and construction of ACC co-expression modules. The clinical information, including the age, sex, race, ethnicity and vital status, and the normal- ized RNA-seq data of 79 patients were obtained from TCGA. The DEGs in normal and ACC tissues were retrieved from GEPIA. A total of 2,472 DEGs were identified with a Q-value <0.01 and a log2FC>1. The expression values of the 2,472 genes in the 79 ACC samples were used to construct the co-expression module by WGCNA. The flash- Clust package (31) was used to perform the cluster analysis
(Fig. 1). A total of four abnormal samples were removed and 75 samples were used for further analysis.
The power value, which primarily affects the independence and the average connectivity degree of the co-expression module of the most critical parameters in the network, is considered an important parameter. Therefore, the power value of the data was assessed (Fig. 2A). A power value of 3, indicated the independent degree was up to 0.8 and the average connectivity degree was higher (Fig. 2B). The degree distribution of the node and the fitting relationship between the node degree and its corresponding proportion were calculated (Fig. 2C and D). The results revealed that the constructed network was a scale-free network conforming to biological characteristics. Therefore, a power value of 3 was used to construct the co-expression modules and the results generated 17 distinct gene co-expression modules in the ACC samples. These co-expression modules were constructed and depicted in different colors (Fig. 3A). They were arranged from large to small by the number of genes they included and the interac- tions of the 17 co-expression modules are shown in Fig. 3B.
Gene co-expression modules correspond to clinic traits. The data on the clinical traits were obtained from TCGA. The modules with common expression patterns that were associated with certain traits were identified based on the association between the module eigengene and the clinical traits (Fig. 3C; Tables SI and SII). The dendrogram and heat map of the eigengene were both used to identify groups of correlated eigengenes and clinical characteristics. The brown module clustered with two important clinical indexes, the tumor stage and vital status (Fig. 4). Moreover, the Pearson’s correlation coefficient (PCC) between the yellow module and the tumor stage is -0.34 (P=0.003), and the PCC between the yellow module and the vital status is -0.47 (P=2x10-05). Furthermore, the yellow module negatively correlated with these two indexes, indicating that genes in this module may be associated with the prognosis of patients (Fig. 3B). Therefore, both modules were defined as core modules for further study.
GO and KEGG enrichment analysis of genes in brown and yellow modules. The association between the genes and the clinical characteristics were calculated along with the MM and GS, where MM was the correlation of the gene expression profile with the eigengenes and GS was the absolute value of the association between the gene and the external traits. MMs and GSs with high thresholds were chosen to avoid false posi- tive prognostic genes, and the top 50% of genes in the brown and yellow modules were selected as candidate genes. GO and KEGG analyses were performed to explore the biological func- tions of the candidate genes in the brown and yellow modules. All significant terms in the annotated systems were represented as colored bars to compare the relative significance of the enriched terms, where the length and color saturation of each term were proportional to the gene count/ratio and the P-value obtained from the enrichment analyses (Figs. 5 and 6). According to the results, GO enrichment analysis indicated that the genes in the brown module were primarily involved in cell division and protein kinase activity (Fig. 6A), which was consistent with the KEGG pathway enrichment results (Fig. 5A). For the yellow module, the genes were enriched in various metabolic pathways
Sample clustering to detect outliers
1500000
Height
1000000
500000
O
A
Scale independence
B
1.0
Mean connectivity
4 5 6 7 8 910 12 14 16 18 20
350
3
1
Signed R
Mean connectivity
300
0.5
250
200
0.0
2
150
100
2
-0.5
8
3
1
0
4 5 6 7 8 910 12 14 16 18 20
5
10
15
20
5
10
15
20
Soft threshold (power)
Soft threshold (power)
C
D
Histogram of k
Check scale free topology scale R^2= 0.86 , slope= - 1.59
800
-0.5
0
Frequency
600
Log10(p(k))
-1.0
400
-1.5
200
-2.0
0
O
0
50
100
150
1.2
1.4
1.6
1.8
2.0
2.2
k
Log10(k)
in GO enrichment analysis such as valine, creatinine and isoleu- cine metabolic processes (Fig. 6B) and various tumor-associated pathways including acute myeloid leukemia and renal cell carci- noma in KEGG enrichment analysis (Fig. 5B).
PPI network-based prognostic gene identification. The brown and yellow modules contained 214 and 168 genes, respectively. To identify the prognostic genes based on the PPI network, the genes were mapped to the pre-built PPI network, and the nodes
A
Cluster Dendrogram
B
C
Network heatmap plot, all genes
Module-trait relationships
1.0
MEgreenyellow
-0.15
0.2
0.04
0.046
0.17
-0.26
(0.2)
(0.08)
(0.1)
(0.03)
MEcyan
-0.11 (0.4)
-0.14
(0.7)
0.006
(0.7)
0.15
-0.17 (0.1)
-0.46
1
(02)
-0.28
0.13
(0.4)
(0 2)
(36-05)
0.9
MEpurple
0.22 (0.06)
0.00044
-0.19
(0.02)
-0.14
(0.3)
MEred
0.18
013
9.12
0.06 (0.6)
0.066 (0.6)
0.15
(0.1) 18.1)
(0.2)
(0.1)
0.13
(0.3)
MElightcyan
0.033 (0.8)
(0.2)
(26-10)
-0.21
-0.12
0.8
(0.08)
-0.043
-0.54
(0.3)
-0.16 (0.2)
MEyellow
0.12
0.47
(0.3)
0.014 (0.9)
(0.7)
0.5
021
(0.003)
(0.3)
0.16 (0 2)
MEsalmon
-0.15 (0.2)
-0.02
(20-05)
0.027
0 18
(0.07)
0.18
0.082
Height
0.19
(0.9)
0.038
(0.8)
0.19 (0.1)
(0.1)
-0.14 (0.2)
(0.1)
0.12
(0.5)
0.7
MEpink
(0.1)
MEmagenta
0.11
(0.7)
0.026
(0.4)
0.0024
0.15 (0.2)
-0.16 (02)
(0.3)
-0.054
(0.8)
Of
0.19
-0.15
(0.6)
-0.016
10.9)
0
MEgreen
0.09
0.12
0.058
MEturquoise
(0.4)
-0.11 (0.4)
(0 4)
(0.1)
(0.2)
0.08 (0.5)
ACTA
0.0079
(02)
02
0 004
0.6
MEmidnightblue
-0.15
(0.5)
0.066 (0.6)
-0.029 (0.8)
(0.9)
0.009 (0.9)
(0.08) -0.11 (0.3)
0.094 (0.4)
(0.2)
-0.0037
MEblack
-0.043
0.061
-0.15 (0.2)
0.2
0.32 (0.005)
(1)
0.18 (0.1)
(0.00)
-0.5
(0.7)
(0.6)
0.5
MEblue
0.19
(0.1)
0.083
(0.5)
0.072
0.088
(0.5)
(0.5)
0.17 (0.1)
0.011 (0.9)
MEbrown
0.39
(50-04)
-0.11
0.62
-0.49
0.0087
0.1
(0.4)
-0.18
(40-09)
0.26
0.0081
(76-06)
0.21
(0.9)
-0.12
(0.00)
MEtan
(0.4)
0.077 (0.5)
(0.1)
-0.019 (0.9)
(0.9)
=
0.13
(0.06)
(0.3)
MEgrey
-1
Tumor stage
(03)
0.02
(09)
0.28
(0.01)
0.17
[0.1)
Module
Age at
0.15 (02)
diagnosis
Vital status
Days to
Sex
Race
colors
death
A
Eigengene dendrogram
B
Eigengene adjacency heatmap
C
Eigengene dendrogram
D
Eigengene adjacency heatmap
1.0
1
1.0
1
0.8
0.8
0.8
0.8
MEgreenyellow
MEsalmon
MEmidnightblue
MEtan
MEsalmon
MEtan
-0.6
MEgreenyellow
MEmidnightblue
0.6
0.6
MEcyan
MElightcyan
MEyellow
MEblack
MEblue
0.4
0.6
0.4
MEpurple
MEpink
MEmagenta
MEbrown
MEcyan
MElightcyan
MEyellow
MEblack
0.4
tumor_stage
MEblue
0.2
MEred
0.4
MEpink
Vital status
0.2
MEpurple
MEred
MEmagenta
0.2
0
MEbrown
vital_status
0
MEgreen
MEturquoise
Tumor stage
Tumor
stage
0.2
MEgreen
MEturquoise
Vital status
that failed to map to the PPI network were ignored. Cytoscape (V3.6.1) (32) software was used to construct the module and to calculate the intramodular connectivity. The intramodular connectivity was calculated for each gene by the connec- tion strength with other module genes and genes with high intramodular connectivity were considered intramodular hub genes. The hub genes in the two modules are represented with larger dots in Fig. 7. The brown module subnetwork contained 105 nodes and 216 edges (Fig. 7A), whereas the yellow module subnetwork contained 92 nodes and 126 edges (Fig. 7B).
To further validate the results, three genes with the highest degree, in the brown (CCNB1, CDC2 and CDK1) and the yellow (PRKCA, RAD23A and UBC) modules were selected to investigate their correlation with the overall survival of the patients (Fig. 8). The nodes with the highest degrees in the two modules were CDK1 and UBC, and their elevated expression levels were significantly associated with poor overall survival (P<0.0001; Fig. 8A and D). Other genes (CCNB1, CDC2, PRKCA and RAD23A) were also significantly associated with the overall survival and they may serve as prognostic genes in
ACC (Fig. 8B, C, E and F). The relationships between all nodes and the overall survival were calculated (data not shown). In addition, as shown in the Fig. 9, ROC curve validated that the high degree gene in the brown module (Fig. 9A) and the yellow module (Fig. 9B) exhibited good diagnostic efficiency for normal and tumor tissues.
Discussion
Compared with other methods such as differential expression analysis, WGCNA places a focus on the association between co-expression modules and clinical traits, and therefore, the results that are considered more reliable yield relevant biolog- ical significance (12). In the present study, a total of 17 distinct gene co-expression modules were identified from 79 patients with ACC by WGCNA to determine the association between the transcriptome of patients with ACC and the clinical traits. After examining the associations between the modules and the clinical traits, two modules were correlated with clinical traits. Several hub genes in the network were identified which
A
KEGG_CELL_CYCLE
KEGG_OOCYTE_MEIOSIS
KEGG_PROGESTERONE_MEDIATED_OOCYTE_MATURATION
KEGG_P53_SIGNALING_PATHWAY
GeneRatio
0.05
KEGG_PEROXISOME
0.10
KEGG_UBIQUITIN_MEDIATED_PROTEOLYSIS
0.15
KEGG_STEROID_BIOSYNTHESIS
0.20
KEGG_PYRIMIDINE_METABOLISM
KEGG_TERPENOID_BACKBONE_BIOSYNTHESIS
P-value
KEGG_CITRATE_CYCLE_TCA_CYCLE
0.04
KEGG_DNA_REPLICATION
0.08
KEGG_GLUTATHIONE_METABOLISM
0.12
KEGG_NON_SMALL_CELL_LUNG_CANCER-
KEGG_COLORECTAL_CANCER
KEGG_FOLATE_BIOSYNTHESIS
5
10
Count
B
KEGG_GAP_JUNCTION-
KEGG_CALCIUM_SIGNALING_PATHWAY
KEGG_ERBB_SIGNALING_PATHWAY
KEGG_ACUTE_MYELOID_LEUKEMIA
P-value
KEGG_CHRONIC_MYELOID_LEUKEMIA
0.025
0.050
KEGG_PHOSPHATIDYLINOSITOL_SIGNALING_SYSTEM-
0.075
KEGG_GLYCOSAMINOGLYCAN_BIOSYNTHESIS_KERATAN_SULFATE
0.100
KEGG_SPHINGOLIPID_METABOLISM
0.125
KEGG_LONG_TERM_DEPRESSION-
GeneRatio
KEGG_LONG_TERM_POTENTIATION-
0.04
KEGG_RENAL_CELL_CARCINOMA
0.06
0.08
KEGG_MELANOMA
0.10
KEGG_B_CELL_RECEPTOR_SIGNALING_PATHWAY
KEGG_VEGF_SIGNALING_PATHWAY-
KEGG_OTHER_GLYCAN_DEGRADATION
1
2
3
4
5
Count
confirmed that CDK1 and UBC serve important roles in the progression of ACC.
Early diagnosis and specific markers are important for managing and limiting tumor development. Electrochemical immunosensors have been developed to detect dehydroepian- drosterone sulfate in blood plasma samples for early diagnosis of patients with pediatric ACC (33). Mohan et al (34) demon- strated that evaluation of G0S2 hypermethylation identified a subgroup of patients with ACC with a rapidly progres- sive disease course which is feasible for clinical treatment options (34). Novel therapeutic regimens are frequently based on specific gene or protein targets present in the disease. Identification of programmed cell death protein 1 function in T cells has improved the development of cancer immuno- therapy (35,36). The identification of potential biomarkers in the present study may serve as novel targets for drugs or in diagnostic methods clinically.
The power value was the most critical parameter affecting the independence and the average connectivity degree of the co-expression modules in WGCNA (10,37). A higher average connectivity degree appeared when the power value =3. The brown module clustered with tumor stage and vital status, whereas the yellow module negatively correlated with these two indices. The Tumor-Node-Metastasis classification system was used in ACC staging and stage III ACC was character- ized by infiltration in the surrounding tissues (38,39). As discussed above tumor stage and vital status were important clinical parameters, and they may also reflect tumor prog- nosis (40-42). Functional enrichment analysis for candidate
genes in the brown module showed that these genes were primarily enriched in pathways associated with cell division. Uncontrolled self-renewal capacity and aberrant regulation of genetic material may promote tumor cell progression and recurrence (43). Tripartite motif-containing protein 3 (TRIM3) has been used as a tumor suppressor due to its ability to regu- late asymmetric cell division in glioblastoma and expression of TRIM3 additionally attenuates the stem-like quality of primary glioblastoma cultures (44). A combination of mito- chondrial division inhibitor 1 and platinum agents produced a synergistic pro-apoptotic effect in drug-resistant tumor cells (45). These results indicate that cell division is associated with tumor progression and genes in the brown module may serve an important role in the regulation of cell division in patients with ACC. Furthermore, genes in the yellow module were primarily enriched in metabolic processes and lyso- some membrane-associated pathways, which may influence the T cell-mediated immune response, tumor invasion, and malignancy (46-48).
Cytoscape was used to construct a co-expression network for the brown and yellow modules, and the genes with high intramodular connectivity were considered as hub genes. The genes with the highest degree of connectivity in the two modules were CDK1 and UBC. CDK1 is a master regulator of the cell cycle and overexpression of CDK1 increases the spheroid-forming ability of tumor cells and the tumor-initi- ating capacity, whereas the inhibition of CDK1 reduced these characteristics (49-51). CDK1 was previously identified as a biomarker in patients with ACC using PCR, western blotting
A
organelle fission
chromosome segregation
nuclear division
nuclear chromosome segregation
sister chromatid segregation
mitotic nuclear division
₹
mitotic sister chromatid segregation
microtubule cytoskeleton organization involved in mitosis
mitotic spindle organization
P-value
sister chromatid cohesion
chromosomal region
0.001
condensed chromosome
chromosome, centromeric region
0.002
condensed chromosome, centromeric region
spindle
0.003
kinetochore
8
0.004
condensed chromosome kinetochore
midbody
Count
mitotic spindle
condensed nuclear chromosome, centromeric region-
10
20
tubulin binding
microtubule binding
30
microtubule motor activity
motor activity
histone kinase activity
cyclin-dependent protein kinase activity
4
histone serine kinase activity
ubiquitin conjugating enzyme activity
ubiquitin-like protein conjugating enzyme activity
cyclin-dependent protein serine/threonine kinase activator activity
0.05
0.10
0.15
GeneRatio
B
visual learning
visual behavior
associative learning
valine metabolic process
creatinine metabolic process
allantoin metabolic process
8
isoleucine metabolic process
taurine metabolic process
Count
alkanesulfonate metabolic process
2
cellular lactam metabolic process
4
vacuolar membrane
6
lysosomal membrane
lytic vacuole membrane
8
endosome membrane
10
recycling endosome
asymmetric synapse
8
brush border
Cajal body
P-value
dystrophin-associated glycoprotein complex
0.0025
glycoprotein complex
0.0050
actin binding
0.0075
lipid transporter activity
calcium-release channel activity
0.0100
ligand-gated calcium channel activity
0.0125
intracellular ligand-gated ion channel activity
phospholipid transporter activity
류
NAADP-sensitive calcium-release channel activity
triglyceride lipase activity
estrogen receptor activity
actin-dependent ATPase activity
0.02
0.04
0.06
GeneRatio
A
SGSM3
PDCD6IP
NRCAM
CYP4F12
B
ANK2
RAD23A
FDFT1
RFX5
CDC34
TWIST1
MVK
CDK
CEP55
HSPA6
KIAA1958
TRIM54
USP13 MINK1
PHKG1
TOP2A
H2AFZ
TAFZ
FOXM1
IGSF11
DLG2
RUVBL2
PPP2R2C
EZH1
GLS
MEIS2
HDAC4
NFAT5 HIBCH
INPP4A
PNLIP POK2
CDC25C
PITPNMI
CDCA4
GARD10
PKMYT1
ORC6
TPP1
HMBOX1 ISG20
GOT1
RRM2
ENC1
TPCN1
GGA1
LACTB2
E2F1
CCNA2
NOARG GATA3
AKR7A3
ASF1B
DNALI1
DDHD2PDGFD
DIXDCTANC1
EDA
DDIT4L
TK1
DT
STMN1
NEK7ZDHHC2
TYMS
EPB41
ANAPCA
BIRC5
ÇOKN3
NPPA
LDLRAD4 ATP11A
SLC12A6
ARFGEF2
UBE2C
PTUG1
HMGCS1
BMS1P4 ||LPCAT SLFN5
CALCOCO2
GMNN
CONB
CKS2
UBE2T IDH1
NPEPPS
ARMC9
CDC20
CONB2KIF22
HSPB8
TROAP
MAGED2
UBC
FXYD6
HOOK1 FLCN
NEK2
MAD211 AURKB AURKAPRG1
KIF11
PHKA1
BIRC7
TBC1D5
KRT19 GJA1
KIF20A
TPX2
KIF4A
RBK
PRTFDC1
BBS1
ABCG1
TSPYL1
ACACA
BRE
NEU3
TMEM150ELC38A2
UBE2S
NDC8OHMMR
RACGAP1
SRCAPPARPASCPEP1
UBE2G211C12 DDR2
PRKCA
BUB
G1
NUF2 SHCBP
GGH
CKM
ITGA7
ITPR3
RFXANINNF185
SRGAP3 CALCOCO1
CBX7
FSCN1
DPP10
ATRX NAT14
TRIM68
IKBKB
ZWINT
LMNB1
PLSOR4
MCOLN3
TMOD1
CENPF
SYNPO
PMAIP1
CENPH
LRROSE
ANLN
ZNF622
ADCK3
PRSS1
STAT58
STATSA
CENPN
ANKFY1
CENPK
HNRNPH1
MELK
HP1BP3
EPHX2 SPARC
SRRM2
MYO19
AMPD3
ESR2
CENPU
KONH2
AFF1
ACTR3
ARHGEF28
PACSIN3
MYO1D
CENPM
MYH14
A
CDK1: +Low+High
B
CCNB1: +Low+High
C
CDC20: +Low+ High
1.00
1.00
1.00
Survival probability
0.75
Survival probability
0.75
Survival probability
0.75
0.50
0.50
0.50
0.25
0.25
0.25
P<0.0001
P<0.0001
P<0.0001
0.00
0.00
0.00
0
1500
3000
4500
0
1500
3000
4500
0
1500
3000
4500
Time (days)
Time (days)
Time (days)
D
UBC: +Low + High
E
PRKCA: +Low +High
F
RAD23A: +Low +High
1.00
1.00
1.00
Survival probability
0.75
Survival probability
0.75
0.75
Survival probability
0.50
0.50
0.50
0.25
0.25
0.25
P=0.046
P<0.0001
P=0.026
0.00
0.00
0.00
0
1500
3000
4500
0
1500
3000
4500
0
1500
3000
4500
Time (days)
Time (days)
Time (days)
A
1.0
B
1.0
0.8
0.8
Sensitivity
0.6
Sensitivity
0.6
0.4
0.4
0.2
CDK1 (AUC=0.88)
0.2
UBC (AUC=0.81)
CCNB1 (AUC=0.91)
PRKCA (AUC=0.78)
0.0
CDC20 (AUC=0.76)
0.0
RAD23A (AUC=0.89)
1.0
0.8
0.6
0.4
0.2
0.0
1.0
0.8
0.6
0.4
0.2
0.0
Specificity
Specificity
and immunofluorescence by Xiao et al (52). UBC is a highly conserved protein which is involved in the selective proteolysis of abnormal proteins (53). Hao et al (54) identified UBC as a differential node protein which may serve as a key regulator in the response of non-small cell lung cancer A549 cells to phycocyanin (54). Furthermore, knockdown of UBC and UBB with mixed small hairpin RNAs suppressed the growth and radio sensitivity of H1299 lung cancer cells (55). UBC was also associated with the regulation of Toxoplasma gondii rhopty protein 18, which serves a key regulatory role in cell
immunity and apoptosis (56). Survival analysis revealed that increased expression of CDK1 and UBC resulted in poorer overall survival of patients with ACC, which is consistent with previous studies (49,50).
However, there were some limitations in the present study. The finally identified hub genes in the network requires experimental verification. At present, there are no studies investigating the function of UBC in ACC, to the best of our knowledge. Therefore, experiments based on clinical samples are required to verify the results of the present study.
In conclusion, the brown and yellow modules were identi- fied as the most critical modules in the progression of ACC. The hub genes CDK1 and UBC were the most significantly expressed genes in the two modules, and they may serve as potential diagnostic and prognostic biomarkers of patients with ACC in the future. The selected candidate genes may serve as targets for the development of novel therapeutics.
Acknowledgements
Not applicable.
Funding
No funding was received.
Availability of data and materials
The datasets analyzed in the present study are available in the The Cancer Genoma Atlas repository, (https://cancergenome.nih.gov/).
Authors’ contributions
YZ designed the study. YZ and LJ developed the methods. YZ collected the sample. YZ analyzed and interpreted the data. YZ and LJ wrote, reviewed and revised the manuscript.
Ethics approval and consent to participate
Not applicable.
Patient consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
References
1. Drenthen LCA, Roerink SHPP, Mattijssen V and de Boer H: Bilaterally enlarged adrenal glands without obvious cause: Need for a multidisciplinary diagnostic work-up. Clin Case Rep 6: 729-734, 2018.
2. Mohan DR, Lerario AM and Hammer GD: Therapeutic targets for adrenocortical carcinoma in the genomics era. J Endocr Soc 2: 1259-1274, 2018.
3. Paragliola RM, Torino F, Papi G, Locantore P, Pontecorvi A and Corsello SM: Role of mitotane in adrenocortical carcinoma-review and state of the art. Eur Endocrinol 14: 62-66, 2018.
4. Fragni M, Fiorentini C, Rossini E, Fisogni S, Vezzoli S, Bonini SA, Dalmiglio C, Grisanti S, Tiberio GAM, Claps M, et al: In vitro antitumor activity of progesterone in human adrenocortical carcinoma. Endocrine 63: 592-601, 2019.
5. Breidbart E, Cameo T, Garvin JH, Hibshoosh H and Oberfield SE: Pubertal outcome in a female with virilizing adrenocortical carcinoma. J Pediatr Endocrinol Metab 29: 503-509, 2016.
6. Nikoleishvili D, Koberidze G, Kutateladze M, Zumbadze G and Mariamidze A: Bilateral adrenocortical carcinoma: Case report and review of literature. Georgian Med News 19-24, 2018.
7. Stigliano A, Cerquetti L, Lardo P, Petrangeli E and Toscano V: New insights and future perspectives in the therapeutic strategy of adrenocortical carcinoma (Review). Oncol Rep 37: 1301-1311, 2017.
8. Megerle F, Herrmann W, Schloetelburg W, Ronchi CL, Pulzer A, Quinkler M, Beuschlein F, Hahner S, Kroiss M and Fassnacht M; German ACC Study Group: Mitotane monotherapy in patients with advanced adrenocortical carcinoma. J Clin Endocrinol Metab 103: 686-1695, 2018.
9. Oddie PD, Albert BB, Hofman PL, Jefferies C, Laughton S and Carter PJ: Mitotane in the treatment of childhood adrenocortical carcinoma: A potent endocrine disruptor. Endocrinol Diabetes Metab Case Rep 2018: pii: EDM180059, 2018.
10. Langfelder P and Horvath S: WGCNA: An R package for weighted correlation network analysis. BMC Bioinformatics 9: 559, 2008.
11. Lin X, Li J, Zhao Q, Feng JR, Gao Q and Nie JY: WGCNA Reveals Key Roles of IL8 and MMP-9 in progression of involve- ment area in colon of patients with ulcerative colitis. Curr Med Sci 38: 252-258, 2018.
12. Wan Q, Tang J, Han Y and Wang D: Co-expression modules construction by WGCNA and identify potential prognostic markers of uveal melanoma. Exp Eye Res 166: 13-20, 2018.
13. Yang Q, Wang R, Wei B, Peng C, Wang L, Hu G, Kong D and Du C: Candidate biomarkers and molecular mechanism investi- gation for glioblastoma multiforme utilizing WGCNA. Biomed Res Int 2018: 4246703, 2018.
14. Ross JS, Wang K, Rand JV, Gay L, Presta MJ, Sheehan CE, Ali SM, Elvin JA, Labrecque E, Hiemstra C, et al: Next-generation sequencing of adrenocortical carcinoma reveals new routes to targeted therapies. J Clin Pathol 67: 968-973, 2014.
15. Papathomas TG, Duregon E, Korpershoek E, Restuccia DF, van Marion R, Cappellesso R, Sturm N, Rossi G, Coli A, Zucchini N, et al: Sarcomatoid adrenocortical carcinoma: A comprehensive pathological, immunohistochemical, and targeted next-generation sequencing analysis. Hum Pathol 58: 113-122, 2016.
16. Gu Y, Gu W, Dou J, Lu Z, Ba J, Li J, Wang X, Liu H, Yang G, Guo Q, et al: Diagnostic role of prostate-specific membrane antigen in adrenocortical carcinoma. Front Endocrinol (Lausanne) 10: 226, 2019.
17. Romero Arenas MA, Whitsett TG, Aronova A, Henderson SA, LoBello J, Habra MA, Grubbs EG, Lee JE, Sircar K, Zarnegar R, et al: Protein expression of PTTG1 as a diagnostic biomarker in adreno- cortical carcinoma. Ann Surg Oncol 25: 801-807, 2018.
18. Duregon E, Volante M, Giorcelli J, Terzolo M, Lalli E and Papotti M: Diagnostic and prognostic role of steroidogenic factor 1 in adrenocortical carcinoma: A validation study focusing on clinical and pathologic correlates. Hum Pathol 44: 822-828, 2013.
19. Ding M, Li F, Wang B, Chi G and Liu H: A comprehensive analysis of WGCNA and serum metabolomics manifests the lung cancer-associated disordered glucose metabolism. J Cell Biochem 120: 10855-10863, 2019.
20. Tang Z, Li C, Kang B, Gao G, Li C and Zhang Z: GEPIA: A web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res 45: W98-W102, 2017.
21. Zheng S, Cherniack AD, Dewal N, Moffitt RA, Danilova L, Murray BA, Lerario AM, Else T, Knijnenburg TA, Ciriello G, et al: Comprehensive pan-genomic characterization of adrenocortical carcinoma. Cancer cell 29: 723-736, 2016.
22. Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, et al: Bioconductor: Open software development for computational biology and bioinformatics 5: R80, 2004.
23. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I, et al: TCGAbiolinks: An R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res 44: e71, 2016.
24. Team RC, R: A language and environment for statistical computing, 2013.
25. Shi H, Zhang L, Qu Y, Hou L, Wang L and Zheng M: Prognostic genes of breast cancer revealed by gene co-expression network analysis. Oncol Lett 14: 4535-4542, 2017.
26. Yu G, Wang LG, Han Y and He QY: clusterProfiler: An R package for comparing biological themes among gene clusters. OMICS 16: 284-287, 2012.
27. Cheng F, Desai RJ, Handy DE, Wang R, Schneeweiss S, Barabási AL and Loscalzo J: Network-based approach to prediction and population-based validation of in silico drug repurposing. Nat Commun 9: 2691, 2018.
28. Csardi G and Nepusz T: The igraph software package for complex network 31. research. Inter Journal, Complex Systems 1695: 1-9, 2006.
29. Collado-Torres L, Nellore A, Kammers K, Ellis SE, Taub MA, Hansen KD, Jaffe AE, Langmead B and Leek JT: Reproducible RNA-seq analysis using recount2. Nat Biotechnol 35: 319-321, 2017.
30. Sing T, Sander O, Beerenwinkel N and Lengauer T: ROCR: Visualizing classifier performance in R. Bioinformatics 21: 3940-3941, 2005.
31. Langfelder P and Horvath S: Fast R Functions for robust correlations and hierarchical clustering. J Stat Softw 46: 1-17, 2012.
32. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B and Ideker T: Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res 13: 2498-2504, 2003.
33. Lima D, Inaba J, Clarindo Lopes L, Calaça GN, Los Weinert P, Lenzi Fogaça R, Ferreira de Moura J, Magalhães Alvarenga L, Cavalcante de Figueiredo B, Wohnrath K and Andrade Pessoa C: Label-free impedimetric immunosensor based on arginine- functionalized gold nanoparticles for detection of DHEAS, a biomarker of pediatric adrenocortical carcinoma. Biosens Bioelectron 133: 86-93, 2019.
34. Mohan DR, Lerario AM, Else T, Mukherjee B, Almeida MQ, Vinco M, Rege J, Mariani BMP, Zerbini MCN, Mendonca BB, et al: Targeted assessment of GOS2 methylation identifies a rapidly recurrent, routinely fatal molecular subtype of adrenocortical carcinoma. Clin Cancer Res 125: 3276-3288, 2019.
35. Ishida Y, Agata Y, Shibahara K and Honjo T: Induced expression of PD-1, a novel member of the immunoglobulin gene super- family, upon programmed cell death. EMBO J 11: 3887-3895, 1992.
36. Ascierto PA, Capone M, Grimaldi AM, Mallardo D, Simeone E, Madonna G, Roder H, Meyer K, Asmellash S, Oliveira C, et al: Proteomic test for anti-PD-1 checkpoint blockade treatment of metastatic melanoma with and without BRAF mutations. J Immunother Cancer 7: 91, 2019.
37. Liu X, Hu AX, Zhao JL and Chen FL: Identification of key gene modules in human osteosarcoma by co-expression analysis weighted gene co-expression network analysis (WGCNA). J Cell Biochem 118: 3953-3959, 2017.
38. Grubbs E and Lee JE: Limited prognostic value of the 2004 International Union Against Cancer staging classification for adrenocortical carcinoma: Proposal for a revised TNM classifi- cation. Cancer 115: 5848, 2009.
39. Libe R: Adrenocortical carcinoma (ACC): Diagnosis, prognosis, and treatment. Front Cell Dev Biol 3: 45, 2015.
40. Petru E, Huber C, Sampl E and Haas J: Comparison of primary tumor size in stage I and III epithelial ovarian cancer. Anticancer Res 38: 6507-6511, 2018.
41. Sylvestre E, Bouzille G, Breton M, Cuggia M and Campillo-Gimenez B: Retrieving the vital status of patients with cancer using online obituaries. Stud Health Technol Inform 247: 571-575, 2018.
42. Roseweir AK, Kong CY, Park JH, Bennett L, Powell AGMT, Quinn J, van Wyk HC, Horgan PG, McMillan DC, Edwards J and Roxburgh CS: A novel tumor-based epithelial-to-mesenchymal transition score that associates with prognosis and metastasis in patients with Stage II/III colorectal cancer. Int J Cancer 144: 150-159, 2019.
43. Palm MM, Elemans M and Beltman JB: Heritable tumor cell division rate heterogeneity induces clonal dominance. PLoS Comput Biol 14: e1005954, 2018.
44. Chen G, Kong J, Tucker-Burden C, Anand M, Rong Y, Rahman F, Moreno CS, Van Meir EG, Hadjipanayis CG and Brat DJ: Human Brat ortholog TRIM3 is a tumor suppressor that regulates asym- metric cell division in glioblastoma. Cancer Res 74: 4536-4548, 2014.
45. Qian W, Wang J, Roginskaya V, McDermott LA, Edwards RP, Stolz DB, Llambi F, Green DR and Van Houten B: Novel combination of mitochondrial division inhibitor 1 (mdivi-1) and platinum agents produces synergistic pro-apoptotic effect in drug resistant tumor cells. Oncotarget 5: 4180-4194, 2014.
46. Noël G, Langouo Fontsa M and Willard-Gallo K: The impact of tumor cell metabolism on T cell-mediated immune responses and immuno-metabolic biomarkers in cancer. Semin Cancer Biol 52: 66-74, 2018.
47. Herrero-Ruiz J, Mora-Santos M, Giráldez S, Sáez C, Japón MA, Tortolero M and Romero F: [TrCP controls the lysosome-medi- ated degradation of CDK1, whose accumulation correlates with tumor malignancy. Oncotarget 5: 7563-7574, 2014.
48. Dykes SS, Gao C, Songock WK, Bigelow RL, Woude GV, Bodily JM and Cardelli JA: Zinc finger E-box binding homeobox-1 (Zeb1) drives anterograde lysosome trafficking and tumor cell invasion via upregulation of Na+/H+ Exchanger-1 (NHE1). Mol Carcinog 56: 722-734, 2017.
49. Ravindran Menon D, Luo Y, Arcaroli JJ, Liu S, KrishnanKutty LN, Osborne DG, Li Y, Samson JM, Bagby S, Tan AC, et al: CDK1 interacts with Sox2 and promotes tumor initiation in human melanoma. Cancer Res 78: 6561-6574, 2018.
50. Warfel NA, Dolloff NG, Dicker DT, Malysz J and El-Deiry WS: CDK1 stabilizes HIF-1a via direct phosphorylation of Ser668 to promote tumor growth. Cell Cycle 12: 3689-3701, 2013.
51. Zeng Y, Stauffer S, Zhou J, Chen X, Chen Y and Dong J: Cyclin-dependent kinase 1 (CDK1)-mediated mitotic phos- phorylation of the transcriptional co-repressor Vgl14 inhibits its tumor-suppressing activity. J Biol Chem 292: 15028-15038, 2017.
52. Xiao H, Xu D, Chen P, Zeng G, Wang X and Zhang X: Identification of five genes as a potential biomarker for predicting progress and prognosis in adrenocortical carcinoma. J Cancer 9: 4484-4495, 2018.
53. Eichner R, Fernández-Sáiz V, Targosz BS and Bassermann F: Cross talk networks of mammalian target of rapamycin signaling with the ubiquitin proteasome system and their clinical implica- tions in multiple myeloma. Int Rev Cell Mol Biol 343: 219-297, 2019.
54. Hao S, Li S, Wang J, Zhao L, Yan Y, Cao Q, Wu T, Liu L and Wang C: Transcriptome analysis of phycocyanin-mediated inhibitory functions on non-small cell lung cancer A549 cell growth. Mar Drugs 16: pii: E511, 2018.
55. Tang Y, Geng Y, Luo J, Shen W, Zhu W, Meng C, Li M, Zhou X, Zhang S and Cao J: Downregulation of ubiquitin inhibits the proliferation and radioresistance of non-small cell lung cancer cells in vitro and in vivo. Sci Rep 5: 9476, 2015.
56. Xia J, Kong L, Zhou LJ, Wu SZ, Yao LJ, He C, He CY and Peng HJ: Genome-wide bimolecular fluorescence complemen- tation-based proteomic analysis of Toxoplasma gondii ROP18’s human interactome shows its key role in regulation of cell immu- nity and apoptosis. Front Immunol 9: 61, 2018.
3 BY NO ND
CC ® S
This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) License.