Identification of four hub genes associated with adrenocortical carcinoma progression by WGCNA

Data collection

We obtained gene expression transcripts per million (TPM) values (Table S1) from the UCSC Xena (http://xena.ucsc.edu/) database, which included 77 ACC samples from the cancer genome atlas (TCGA) (https://cancergenome.nih.gov/) and 128 normal samples from genotype tissue expression (GTEx) (https://www.gtexportal.org/home/). The two databases raw sequencing reads were recalculated with a unifying pipeline. Clinical data were downloaded from TCGA using the “cgdsr” package in R (v3.1.3) (R Core Team, 2015; Jacobsen, 2015).

DEG screening

Of the 60,498 genes in each sample, we removed genes with a mean TPM ≤ 2.5 (>1 is a common cutoff for determining if an isoform is expressed or not (Liu, Jing & Tu, 2016)) in the cancer and normal samples and thus retained 13,987 genes. For those genes in the samples that showed significant changes, we used analysis of variance (ANOVA) in R (R Core Team, 2013) to determine the variance in genes between the two groups. ANOVA is a collection of statistical models useful for DEG analysis (Alabi et al., 2018; Monterisi et al., 2015). We obtained 2,953 significant DEGs (Table S2) in ACC with a p < 0.001 and |log2 (fold-change)| > 1 cutoff.

Co-expression network construction by WGCNA

Weighted gene co-expression network analysis (v1.49) can be applied to identify global gene expression profiles as well as co-expressed genes. Therefore, we installed WGCNA package for co-expression analysis using Bioconductor (http://bioconductor.org/biocLite.R). We used the soft threshold method for Pearson correlation analysis of the expression profiles to determine the connection strengths between two transcripts to construct a weighted network. Average linkage hierarchical clustering was carried out to group transcripts based on topological overlap dissimilarity in network connection strengths. To obtain the correct module number and clarify gene interaction, we set the restricted minimum gene number to 30 for each module and used a threshold of 0.25 to merge the similar modules (see the detailed R script in File S1).

Identification of clinically significant modules

We used two methods to identify modules related to clinical progression traits. Module eigengenes (MEs) are the major component for principal component analysis of genes in a module with the same expression profile. Thus, we analyzed the relationship between MEs and clinical traits and identified the relevant modules. We used log10 to transform the p-value from the linear regression between gene expression and clinical stage, which was defined as gene significance. Average gene significance in a module was defined as module significance.

Functional and pathway enrichment analysis

The database for annotation visualization and integrated discovery (DAVID) (v6.8) (http://david.abcc.ncifcrf.gov/) was used for functional annotation of genes to better understand their biological functions. All genes in the blue module were uploaded for gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analyses, with cutoffs of p < 0.01 and p < 0.05 established for significant biological processes and pathways, respectively.

PPI and co-expression analysis

Genes were uploaded to the search tool for the retrieval of interacting genes (STRING) (v10.5) (https://string-db.org/) database. Confidence was set to more than 0.4 and other parameters were set to default. We visualized the gene co-expression network with Cytoscape (v2.7.0) (Shannon et al., 2003).

Gene expression correlation with stage and survival analysis

The correlation between gene expression and stage was determined using GEPIA (http://gepia.cancer-pku.cn/index.html) (Tang et al., 2017). The correlation between gene expression and overall survival (OS) was established using the Cox model. A hazard ratio p-value of <0.01 was considered significant. Each gene with higher expression in the ACC samples had corresponding lower survival expectation. The “limma” (Ritchie et al., 2015) R package was used to test for the significantly expressed gene in GSE10927.

Construction and analysis of gene co-expression network with DEGs in ACC

Genes with a mean TPM ≤ 2.5 were removed from the two groups and the remaining 13,987 genes were used for differential expression analysis with ANOVA. In total, 2,953 significant DEGs were identified with a cutoff of p < 0.001 and |log2(fold-change)| > 1 (Fig. 1A), which included 1,181 up-regulated and 1,772 down-regulated genes (Fig. 1B). The 2,953 gene expression levels in ACC and normal samples are shown in the heatmap in Fig. 1C and Table S2.

Genes with similar expression patterns may participate in similar biological processes or networks (Mao et al., 2009). To better understand the gene expression network during ACC development, the co-expression network of the 2,953 DEGs was analyzed by WGCNA. First, to determine whether all 77 ACC samples were suitable for network analysis, the sample dendrogram and corresponding clinical traits were analyzed. We found that all samples were included in the clusters and passed the cutoff thresholds (Fig. 1D). The power value is a critical parameter that can affect the independence and average connectivity degree of the co-expression modules. Thus, network topology using different soft thresholding powers was screened, with β = 6 (scale free R² = 0.928) selected for later analysis (Figs. 1E, 1F and 1G). We then constructed the gene co-expression network using WGCNA based on the hierarchical clustering of the calculated dissimilarities, and nine modules were obtained (Fig. 1H; Table S3). We used eigengenes as representative profiles and quantified module similarity by eigengene correlation (Fig. 1I).

Correlation of blue module with clinical stage and progression

We investigated whether any module was correlated with clinical stage and tested the relevance between each module and ACC clinical traits. We found that module significance of the blue module was higher than that of any other, implying it had greater correlation with ACC stage (Fig. 2A). The blue module also displayed a positive correlation with ACC clinical stage (r = 0.5, p = 6e-06) and negative correlation with OS (r = −0.56, p = 3e-07) (Fig. 2B). The eigengene dendrogram and heatmap indicated that the MEblue and MEyellow modules were highly correlated with clinical stage (Fig. 2C). Finally, gene significance and module membership were plotted for the blue module (Fig. 2D), which indicated that this module was significantly related to clinical stage.

To determine the function of the 650 genes in the blue module, GO and KEGG function and pathway enrichment analyses were performed by DAVID functional annotation (Huang, Sherman & Lempicki, 2009). For GO biological processes, genes in the module were significantly enriched in cell division (p = 1.05e-26) (Fig. 2E; Table S4), whereas for KEGG pathway analysis, the genes were mainly enriched in cell cycle (p = 2.7e-19) and DNA replication (p = 8.27e-8) (Fig. 2F; Table S5) pathways. These processes and pathways all play critical roles in cancer progression (Tachibana, Gonzalez & Coleman, 2005), implying that genes in this module may be involved in ACC progression.

PPI and co-expression networks to identify hub genes in ACC progression

To clarify high confidence hub genes, we entered the blue module genes into the STRING (Szklarczyk et al., 2015) database. The genes were ranked by the protein–protein interaction nodes and the top 5% of genes (16 genes) were chosen as candidate hub genes (Fig. 3A; Table S6). As highly connected hub genes in a module play important roles in biological processes (Liu, Jing & Tu, 2016), genes in the blue module were ranked by their degree of gene co-expression connectivity (Table S7). To identify genes that may play notable roles in ACC progression, the top 5% of genes (31 genes) (Fig. 3B) in the blue module with the highest connectivity were classified as candidate hub genes for further analysis. Finally, four common genes (i.e., TOP2A, TTK, CHEK1, and CENPA) in the two analysis were identified as hub genes in ACC (Fig. 3C). These four genes were highly expressed in ACC samples compared with normal samples (Figs. 3D–3G), indicating that they likely act as oncogenes in ACC. Further analysis of the GSE10927 dataset, which included microarray data of 10 normal samples and 33 ACC samples (Human Genome U133A 2.0 Plus; Affymetrix, Santa Clara, CA, USA) (Giordano et al., 2009), demonstrated that the four genes showed significant high expression in ACC (Figs. S1A–S1D). Furthermore, based on immunoreactivity experiments, TOP2A is reported to be highly expressed in ACC (Giordano et al., 2003).

Significant associations of hub genes with ACC stage and survival

We investigated the four hub genes to better understand their functions. We found that TOP2A, TTK, CHEK1, and CENPA play critical roles in biological processes that are highly correlated with cancer (Dominguez-Brauer et al., 2015), such as DNA topological structure, cell cycle progression, and mitosis (Hoffmann et al., 2011; Liu et al., 2000; De Resende et al., 2013; Thu et al., 2018), thereby suggesting their possible role in cancer development. Further exploration of their expression patterns during ACC clinical progression showed that the levels of these genes were significantly altered with clinical stage and markedly increased at stage III and IV (Figs. 4A–4D). This correlation between the expression levels of the four genes and ACC progression may be useful in ACC diagnosis.

Tumor prognosis is an important feature in cancer and has attracted considerable attention. To assess the utility of WGCNA at identifying hub genes indicative of ACC, we conducted survival analysis (Figs. 4E–4H). We separated the samples into two groups according to median gene expression levels and performed survival analysis using the Cox model. Survival analysis showed that the expression of all four genes was significantly correlated with OS (Figs. 4E–4H), with higher expression associated with lower patient survival time. The correlation between the hub genes and ACC prognosis suggests that these four genes likely contribute to the progression of ACC.