A transcription factor signature predicts the survival of patients with adrenocortical carcinoma

Data collection and preprocessing

A total of 1,639 human TFs were identified from a previously published study (Lambert et al., 2018). TCGA ACC level 3 RNAseqv2 data (RSEM_genes_normalized file) and corresponding clinical information were downloaded from the TCGA database (https://tcga-data.nci.nih.gov/). A total of 79 patients with ACC from TCGA were included after excluding those lacking complete clinical and survival data, which served as training sets. The gene expression profiling of 22 patients with ACC was selected via the accession number GSE19776 from GEO (http://www.ncbi.nlm.nih.gov/geo/) using the GPL570 platform [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, which served as validation set. The clinical data included age, gender, and tumor stage. Firstly, we obtained 1,554 overlapped TFs between TFs in literature and TCGA dataset, and 1,508 overlapped TFs between TFs in literature and GEO dataset. A total of 118 TFs were excluded from the TCGA dataset by removing lower expression genes. Next, we acquired the common TFs between the TCGA and GEO datasets.

Survival analysis

After obtained common TFs between TCGA and GEO datasets, univariate Cox survival analysis was performed using Cox proportional hazards regression model. Only TFs with Cox P < 0.001 and Log Rank P < 0.001 on univariate analysis were incorporated into the Lasso Cox regression analysis. Kaplan-Meier method was used to analyze the correlation between overall survival (OS) and TF expression, and the OS of different patient groups were compared using the Log-Rank test. The “survival” package (Therneau & Lumley, 2015) in R software was used for survival analysis and the time-dependent ROC (Receiver Operating Characteristic) curve analysis was performed using the “survival ROC” package (Robin et al., 2011).

Construction of ACC prognostic signature

LASSO Cox regression model was widely used for high-dimensional predictor identification. In the present study, OS-associated TFs were used to select the significant TFs associated with the OS of ACC patients according to the coefficient value. These factors were incorporated in the multivariate Cox regression model to construct the ACC prognostic signature. The risk score for each TF-encoding gene was calculated as follows: $Riskscore={\sum }_{i=1}^{n}ex{p}_i\ast {\beta }_i$, where n is the number of selected genes, exp_i is the expression level of gene i, and β_i is the coefficient of gene i.

mRNA expression of 13 TFs in ACC cell line

Gene expression data of all 13 TFs in the prognostic signature from the SW13 cell line were downloaded directly from the Cancer cell line encyclopedia (CCLE) database (https://portals.broadinstitute.org/ccle/home), which is an open-access database covering large-scale deep sequencing raw data of more than 900 human cancer cell lines.

Gene set enrichment analysis (GSEA)

GSEA was performed to analyze the significance of the 13 TFs constituting the prognostic signature using GSEA v2.0.12 (http://www.broadinstitute.org/gsea) by computing the enrichment score for each gene set (Subramanian et al., 2005). The distributions of these TFs against the rank-ordered gene ontology (GO) hallmarks were characterized using GSEA with the default settings. Positive and negative normalized scores indicated enrichment in the high-risk and low-risk groups respectively.

Construction of the TF prognostic signature

The scheme for developing the TF signature is outlined in Fig. 1. After initially identifying 1,639 TFs by literature search, the down-regulated genes were removed and 1,304 common TFs were screened from TCGA and GEO datasets. Using the TCGA training set (n = 79), univariate regression analysis showed that 23 TFs were correlated with OS (Cox-P < 0.001 and Log-Rank P < 0.001), of which 13 were identified by the Lasso regression analysis, and the risk score was calculated by multivariate cox analysis.

The λ value was selected in the LASSO Cox regression analysis when the median of the sum of squared residuals was the smallest (Fig. 2A). The following survival-related TFs with non-zero coefficients were then screened: CREB3L3 (cAMP-responsive element-binding protein 3-like 3), NR0B1 (nuclear receptor subfamily 0, group B, member 1), CENPA (centromere protein-A), FOXM1 (Forkhead Box M1), E2F2 (E2F transcription factor 2), MYBL2 (v-myb myeloblastosis viral oncogene homolog (avian)-like 2), HOXC11 (homeobox C11), ZIC2 (Zic family member 2), ZNF282 (zinc finger protein 282), DNMT1 (DNA methyltransferase 1), TCF3 (transcription factor 3), ELK4 (ETS transcription factor ELK4) and KLF6 (Krüppel-like factor 6) (Fig. 2B). Only CREB3L3 and NR0B1 were negatively correlated with the remaining TFs (Fig. 2C), while CENPA, FOXM1, and E2F2 displayed a strong correlation (Table 1). As shown in Fig. 3A, a total of 79 patients from the TCGA database were divided into a high-expression group (n = 40) and a low-expression group (n = 39). Significantly, high expression of CREB3L3 (HR = 0.663, Cox P = 5e−05, Log-Rank P = 1.97e−07) and NR0B1 (HR = 0.799, Cox P = 6.93e−05, Log-Rank P = 3.18e−06) were associated with good prognosis, and that of other TFs with poor prognosis (HR > 1, P < 0.001).

TF signature risk score can predict prognosis of ACC patients

The clinical relevance of these TFs was further assessed by multivariate Cox regression analysis based on the TCGA training set, and the risk score based on their expression levels and coefficients was calculated. The 13-TF risk score classified the patients from TCGA training set into the high- (n = 39) and low-risk (n = 40) groups (Fig. 4A). Except for CREB3L3 and NR0B1, all TFs were overexpressed in the high-risk group (Fig. 4B). Furthermore, ACC patients in the high-risk group had significantly shorter survival compared to the low-risk patients (HR = 16.95 (5.02–57.2); Cox P = 5.11e−06; Log Rank P = 2.09e−09) (Fig. 4C). The sensitivity and specificity of the 13-TF signature were determined using time-dependent receiver operating characteristic (ROC) analysis, and the area under curves (AUCs) at all follow-up time points were greater than 0.9 (Fig. 4D). The predictive model was then validated in a GEO dataset, and as shown in Fig. 5A, the high-risk group (n = 11) had worse survival compared to the low-risk group (n = 11). In addition, the AUC values of the signature were greater than 0.75 (Fig. 5B). Next, we investigated the mRNA levels of all 13 TFs in the ACC cell line through the CCLE database (Fig. 5C). The results showed that the mRNA levels of CREB3L3, NR0B1, and HOXC11 were almost undetectable in SW13 cells. In contrast, the other 10 TFs were highly expressed in ACC cells. Taken together, the 13-TF signature can predict the survival of ACC patients with high sensitivity and specificity.

The risk score is an independent prognostic factor of ACC

The prognostic relevance of the 13-TF signature was further validated by multivariate Cox regression analysis after normalizing for variables including age, gender, and pathological stage. In both the training and validation ACC cohorts, the 13-TF risk score was an independent prognostic factor (Table 2). However, no significant correlation was seen between OS and age, gender, or pathological stage. The high- and low-risk groups of both training and validation datasets were further divided into subgroups based on age (≤50 vs >50 years), gender (male vs female), and pathological stage (I–II vs III–IV). As shown in the Figs. 6A and 6B, patients in the high-risk group had poor prognosis and significantly shorter survival compared to those in the low-risk group regardless of other variables. Thus, the 13-TF signature is an independent prognostic predictor of ACC.

The 13-TF signature is associated with cancer-related functions

GSEA results showed that four hallmarks including G2M_CHECKPOINT (P = 0.021), E2F_TARGETS (P = 0.023), SPERMATOGENESIS (P = 0.046), and MITOTIC_SPINDLE (P = 0.048) were significantly enriched in high-risk patients, suggesting a mechanistic basis of the prognostic role of the 13-TF signature in ACC (Fig. 7).