ELSEVIER

International Immunopharmacology

journal homepage: www.elsevier.com/locate/intimp

E

International Immunopharmacology

Identification of immune-related biomarkers in adrenocortical carcinoma Immune-related biomarkers for ACC

Check for updates

Yun Penga,1, Yuxuan Songb,1, Jin Ding®, Nan Liª, Zheyu Zhangª, Haitao Wangd,»

ª Department of Urology, Tianjin Medical University Second Hospital, Hexi, Tianjin 300000, China

b Department of Urology, Tianjin Medical University General Hospital, Heping, Tianjin 300000, China

” Department of Metabolism and Endocrinology, The Second Xiangya Hospital, Central South University, Changsha, Hunan 410011, China d Department of Oncology, Tianjin Medical University Second Hospital, Hexi, Tianjin 300000, China

ARTICLE INFO

Keywords:

Adrenocortical carcinoma Tumor microenvironment Immune Biomarker Prognosis WGCNA

ABSTRACT

Emerging evidence has suggested that the tumor microenvironment, including immune infiltration, plays a crucially important role in tumor progression. Nevertheless, limited studies have been conducted on this topic in adrenocortical carcinoma. The present study aimed to explore the immune-related biomarkers in adrenocortical carcinoma. CIBERSORT was used to estimate the abundances of 22 kinds of immune cells, and univariable Cox analysis was performed to find survival-related immune cells with both Overall Survival (OS) and Progression- Free Interval (PFI). DESeq2 was applied to find differentially expressed genes between adrenocortical carcinoma and normal control samples; subsequently, weighted correlation network analysis and protein-protein interac- tion (PPI) network analysis were conducted to identify immune-related hub genes. xCell, TISIDB, and MsigDB were searched to validate the immune associations of hub genes. Eventually, univariable Cox and Kaplan-Meier analysis were used to assess the prognostic implications of the hub gene with the GEO database. Consequently, we identified two hub immune-related genes (ERN1, CEP55), GSEA revealed that both were mainly involved in tumor progression and immune response. ROC analysis indicated that ERN1 can accurately predict the 1-, 3-, and 5-year PFI, and CEP55 had the best performance for the prediction of both OS and PFI compared with other traits. Univariable Cox and Kaplan-Meier analysis showed that both genes have a significant effect on prognosis. Furthermore, both hub genes were validated in GEO datasets. The hub genes can provide better insights into tumor microenvironment and serve as potential biomarkers for immunotherapy in adrenocortical carcinoma.

1. Introduction

Adrenocortical carcinoma (ACC) is a complex malignant tumor with a heterogeneous clinical characteristic and poor overall survival, making it a substantial clinical challenge [1-4]. Localized ACC gen- erally has a high risk of recurrence (50%) after surgical resection [5]. Nevertheless, the most widely used tumor, lymph node, and metastasis (TNM) classification system proposed by the American Joint Committee on Cancer remains unacceptable due to its lack of crucial genomic molecular characteristics [6,7]. Thus, exposing the potential genomic properties in ACC is of vital importance for the improvement of therapies and survival prediction [1,6]. Robust and trustworthy mole- cular biomarkers, except the usual clinical characteristics, are necessary for improving the diagnostic, therapeutic, and recurrence risk predic- tion for ACC.

Emerging evidence, which suggested that the tumor microenvironment

(TME) is crucially important in tumor progression, has paved the way for unprecedented progress in the understanding of tumor growth and devel- opment [8]. Tumors exhibit diverse TME composed of various kinds of cells including immune cells, fibroblasts, and inflammatory cells, whose com- munications can promote tumor progression and metastasis [9]. Thereby, identification of special TMEs in various tumors can lead to a better un- derstanding of tumor proliferation, invasion, and metastasis. Moreover, the emerging prominence of precision medicine in oncology has increased the applications of targeted immunotherapy [10]. Nevertheless, as few studies have been conducted on ACC, the present study pointed to explore the prognostic implications of immune-related genes in ACC.

In the present study, Immune cell fractions were evaluated by CIBERSORT, and we identified immune-related genes in ACC by DESeq2 and weighted correlation network analysis (WGCNA). Subsequently, hub genes were identified by combining the coexpression network and PPI network analysis. ROC curve, univariable Cox, and

* Corresponding author at: Department of Oncology, Tianjin Medical University Second Hospital, Pingjiang Road, Hexi, Tianjin, China.

E-mail address: haitao_peterrock@outlook.com (H. Wang).

1 Yun Peng and Yuxuan Song contributed equally to this work.

Kaplan-Meier analysis were employed to explore the association with the clinical traits of the hub genes. Our results provide a better insight into the TME in ACC and present potential biomarkers for ACC im- munotherapy.

2. Materials and methods

2.1. Data collection and preprocessing

The high-throughput RNA sequencing raw count data corre- sponding to 79 ACC samples and 258 normal adrenal gland samples were retrieved from the TCGA data portal [11] and GTEx (V8 release) [12] database respectively. The FPKM value of RNA sequencing data of 79 ACC samples was obtained from the TCGA database. We also ac- quired detailed clinical information of the ACC samples, including gender, age, TNM stage, and survival information encompassing both Overall Survival (OS) and Progression-Free Interval (PFI) [13] from TCGA. After filtering samples without corresponding clinical data, a total of 77 ACC samples from TCGA and 258 normal control samples was included in subsequent analysis. We also confirmed the absence of duplicated samples in our study. GSE49278, GSE19750, and GSE76021 were derived from the GEO database to validate our findings. The re- quirement of approval from the Ethics Committee was waived as the data in the study came from the TCGA, GTEx, and GEO database Committee. The Ensembl (version 99) [14,15] database was employed to annotate genes. We used the R program [16] for most of the analysis. The subsequent process was depicted in Fig. 1.

2.2. Identification of immune cells with survival implications

The CIBERSORT algorithm can provide precise profiles of cellular heterogeneity in both microarray and RNA sequencing data obtained

from bulk tumor samples [17]. Specifically, a fraction of 22 types of TIICs (tumor-infiltrating immune cells) was evaluated using the CIBE- RSORT package in an absolute mode representing the overall abun- dance of corresponding immune cells [18] combining FPKM value of RNA sequencing from ACC. A univariable Cox analysis was applied to analyze the association between the immune cells and the OS and PFI. A threshold of P-value < 0.05 was considered significant.

2.3. Differentially expressed genes

DESeq2 was used to reveal the differentially expressed genes be- tween ACC and normal control samples with a threshold of |log2-fold change| > 1 and an adjusted P-value < 0.01 using the high- throughput RNA sequencing raw count data from TCGA and GTEx.

WGCNA [19,20] was applied with the biweight midcorrelation as- sessing coexpression similarities to explore the gene coexpression net- work in differentially expressed gene profiles. First, a similarity matrix was defined using biweight midcorrelation function S(ij) = (1 + bicor (i,j))/2 (i and j represent genes). Followed by a adjacency matrix de- termined by a(ij) = |S(ij)|ß, the eligible parameter ß was chosen to ensure the scale-free topology criterion. Next the topological overlap matrix was generated from adjacency, and finally, a dynamic tree cut method was used to identify gene modules. A threshold of 9 for soft- power, 25 for minModuleSize and 0.20 for mergeCutHeight was chosen to investigate the network. 4 types of immune cells significantly asso- ciated with both OS and PFI (T cells CD4 memory resting, NK cells activated, Macrophages M0, and Dendritic cells activated) were defined as the traits in WGCNA. Gene significance (GS), module significance (MS), and module membership (MM) were explained by biweight

Fig. 1. The workflow of the identification of immune-related hub genes.

77 ACC clinical data from TCGA

77 ACC mRNA expression (FPKM) data from TCGA

77 ACC mRNA expression raw counts data from TCGA

258 normal control mRNA expression raw counts from GTEx

immune cell abundance estimated by CIBERSORT

Defferential expression gene determined by DESeq2

survival- related immune cell

Co-expression network construction by WGCNA and identification of immune-related module

High connected gene determined by GS and MM

Construction of PPI network

Verification of Immune-correlation by xCell, TISIDB and MsigDB

Identification of important cluster by MCODE

hub gene

Biological funcitons analysed by GSEA and Exploration of the relationship with clinical characteristics

prognostic implications of hub genes

Validation in Oncomine and GEO

midcorrelation coefficients. Gene Set Enrichment Analysis (GSEA) [21] was used to explore the Kyoto Encyclopedia of Genes and Genomes (KEGG) [22] and Gene Ontology (GO) [23,24] for the most significant modules associated with the immune cells using the clusterProfiler [25] R package.

Hub genes are those highly connecting with the nodes in the co- expression network which can play crucial biological functions. We defined genes with both high gene significance (GS > 0.4) and high intramodular connectivity (MM > 0.8) as highly connected genes. We selected 2089 genes in the modules with maximal module significance corresponding to 4 survival-related immune cells in the coexpression network for the subsequent analysis. The STRING database was used to predict the interactions among proteins and to construct a PPI network with a cut-off value of 400 for combined score. The PPI network was further imported to Cytoscape. The Molecular Complex Detection (MCODE) [26] algorithm was used to investigate the important clusters in the PPI network with a default threshold of degree cutoff at 2, a node score cutoff at 0.2, a K-Core at 2, and a maximum depth from seed at 100. Genes in both the highly connected genes and important clusters were considered hub genes.

2.6. Validation of immune association with hub genes

To confirm the association between hub genes and immune cells, the xCell [27], which combines deconvolution approaches with gene set enrichment analysis and can generate a reliable score for different kinds of immune cells, was performed to verify the immune correlation of hub genes by Pearson’s Correlation analysis. The TISIDB [28] was used to further strengthen these findings. We searched for correlation of hub genes with lymphocyte, immunomodulator, and chemokine. Further- more, the GSEA method was applied to identify the most strongly as- sociated immunologic gene sets in the Molecular Signatures Database (MsigDB) [29].

2.7. Exploration of clinical associations and the prognostic implications of hub genes

For analyzing the biological function of hub genes, both GO and KEGG analysis were performed using GSEA. A heatmap depicting cor- relations of 2 hub genes with clinical traits was visualized via ComplexHeatmap [30]. The gene expression differences of hub genes in different groups according to clinical traits were also compared. The prediction precision for OS and PFI among various clinical traits and hub genes were also compared by a receiver operating characteristic curve (ROC) analysis, the area under the ROC Curve (AUC) was cal- culated. The influence of hub genes on survival for both OS and PFI was identified by univariable Cox analysis and Kaplan-Meier analysis.

2.8. Validation of the clinical implication of hub genes

The differential expression of hub genes between ACC and normal control samples was further verified using Oncomine [31]. The gene expression levels of hub genes and survival characteristics in GSE49278, GSE19750, and GSE76021 was analyzed to verify the prognostic implications of hub genes. Kaplan-Meier analysis were per- formed to determine the statistical significance between different groups defined by the median gene expression values.

3. Results

3.1. Enumeration and survival significance of immune cells

CIBERSORT algorithm was run to derive the abundance fractions of

22 types of immune cells in R with a parameter of 1000 permutation and absolute mode, followed by a univariable Cox analysis for the 77 ACC samples from TCGA, the results were shown in Supplement Table 1. Only immune cells with a significance on both OS and PFI were selected for further analysis, finding 4 types of immune cells which presented significant threshold of P < 0.05, namely Dendritic cells activated, Macrophages M0, NK cells activated and T cells CD4 memory resting.

3.2. Construction of weighted co-expression network

A total of 9496 differentially expressed genes between ACC and normal control samples including 5941 up-regulated genes and 3555 down-regulated genes, was derived by DESeq2(Fig. 2A). The differen- tially expressed genes were further applied to construct a correlation network after filtering those with a median absolute deviation of 0. Thereby, a total of 7289 genes was included in WGCNA, Fig. 2B-C de- picted the topological overlap matrix (TOM) and clustering dendrogram of genes. A biweight midcorrelation was applied to define the gene significance and module significance by incorporating immune cell abundance into the co-expression network. The relationships between the 4 survival-related immune cells and modules in coexpression net- work were explored in Fig. 2D, where T cells CD4 memory resting and NK cells activated were most associated with brown (499 genes), Macrophages M0 most correlated with turquoise (1181 genes), and Dendritic cells activated was most related to green module (409 genes). GSEA was performed to reveal the biological process in GO, which showed that the brown module was most related to cell adhesion, de- fense response and biological adhesion, the turquoise module most to reproduction, microtubule cytoskeleton organization and mitotic cell cycle, and the green module to cell cycle, cell cycle process and neuron differentiation. KEGG was also explored. The brown, turquoise, and green module were mostly involved in Cytokine-cytokine receptor in- teraction, MicroRNAs in cancer, and Cell cycle, respectively.

The correlations between MM and GS in different modules were depicted in Fig. 3A, which shown that GS in different immune cells significantly associated with corresponding MM. We first retrieved the highly connected genes in the coexpression network with a threshold of MM > 0.8 and GS > 0.4. Meanwhile, a PPI network was constructed via the STRING database with 2089 genes in the 3 significant modules (brown, turquoise, and green modules). Consequently, a PPI network containing 726 nodes and 10,592 edges was developed. We applied MCODE in Cytoscape to refine important clusters from the PPI network, resulting in 29 clusters accounting for 264 genes. The top 2 clusters of the PPI network based on cluster scores were depicted in Fig. 3B. Two genes (ERN1, and CEP55) were defined as hub genes by combining highly connected genes with cluster genes.

3.4. Validation of the relationship between hub genes and immunity

We calculated Pearson’s correlation coefficient for both hub genes with immune cell abundances estimated by CIBERSORT and xCell methods separately. The results were shown in a correlation heatmap in Supplement Fig. S1A. TISIDB was searched for the relationship between hub genes and immune cells in different tumors (Fig. 4A). The che- mokines and immunomodulators including both immunoinhibitors and immunostimulators were also explored for the associations with hub genes using TISIDB as presented in Supplement Fig. S1B and S1C. GSEA was conducted to investigate the enrichment of hub genes in MSigDB immunologic gene sets, showing that ERN1 was most enriched in M3345 (GSE14000_UNSTIM_VS_4H_LPS_DC_UP: Genes up-regulated in comparison of dendritic cells (DC) before and 4 h after LPS (TLR4 agonist) stimulation) [32] and CEP55 was most related with M3577

Fig. 2. Identification of immune-related genes. (A) the volcano of differential expression genes in ACC compared with normal controlled samples, red spots represent up-regulated genes, and blue spots represent down-regulated genes; (B) topological overlap matrix (TOM) were depicted in a heatmap, each row and column corresponds to a gene, light color denotes low topological overlap, and progressively darker red denotes higher topological overlap. Darker squares along the diagonal correspond to modules. The gene dendrogram and module assignment are shown along the left and top .; (C) the clustering dendrogram of genes and corresponding modules; (D) the relationship of genes in modules with the survival-related immune cell was investigated in a heatmap where the color representing the biweight midcorrelation coefficients. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

A

B

200

150

-log10(padj)

up-deg

100

non-deg

down-deg

50

0

-10

-5

0

5

log2FoldChange

10

E

C

D

Module-trait relationships of cibersort

Cluster Dendrogram

1.00

MElightcyan

0.202

(0.078)

0.19

(0.098)

-0.187

(0.1)

-0.257

(0.024)

MEblue

0.661

(5.9e-11)

0.31

(0.0061)

-0.0406

-0.369

(0.73)

(0.00094)

0.95

MEcyan

0.181

-0.0378

0.2

(0.082)

-0.0329

(0.11)

(0.74)

(0.78)

MEmidnightblue

-0.286

(0.012)

-0.151

0.281

(0.013)

0.0296

(0.19)

(0.8)

0.90

MEgrey60

-0.126

-0.192

0.111

(0.28)

(0.095)

(0.34)

0.181

(0.12)

MEgreen

-0.468

(1.80-05)

-0.174

(0.13)

-0.0278

(0.81)

0.592

(1.40-08)

0.85

MEsalmon

-0.436

(7.5e-05)

-0.121

(0.29)

-0.139

(0.23)

0,509

(2.3e-06)

Height

MEtan

-0.439

(6.40-05)

-0.128

(0.27)

0.0121

(0.92)

0.221

(0.053)

0.80

MEturquoise

-0,422

(0.00013)

-0.272

(0.017)

0.385

(0.00054)

0,309

(0.0062)

MEblack

-0.432 (8.90-05)

-0.245

(0.032)

0.122

(0.29)

0.432

(8.6e-05)

0.75

MEpink

-0.0544

(0.64)

-0.129 (0.27)

0.11

(0.34)

0.276

(0,015)

MEmagenta

0.047

(0.68)

0.126

(0.28)

0.0713

(0.54)

0.102

(0.38)

0.70

MEgreenyellow

-0.21

(0.067)

0.00227

0,308

(0.98)

(0.0064)

-0.0131

(0.91)

MEbrown

0.671

(2.40-11)

0.41

(0.00021)

-0.138

(0.23)

-0.318

(0.0049)

0.65

MEpurple

0.487 (7.1e-06)

0.359

(D.DD14)

:0.228

(0.046)

-0.232

(0.043)

MElightgreen

-0.147

(0.2)

0.0783

(0.5)

-0.0372

(0.75)

0.259

(0.023)

MEyellow

-0.21

(0.067)

-0.0077

(0.95)

-0.112

(0.33)

0.267

(0.019)

module color

MEred

0.135

(0.24)

0.151

-0.196

(0.087)

0.107 (0.35)

(0.19)

MEgrey

0.161

0.0182

(0.16)

0.124

-0.0375

(0.88)

(0.28)

(0.75)

T cells CD4 memory resting

NK cels activated

Macrophages MO

Dendritic cells

activated

-0.6

-0.4

0.2

0

-0.2

-0.4

(GSE15750_DAY6_VS_DAY10_EFF_CD8_TCELL_UP: Genes up-regulated in comparison of wild type CD8 effector T cells at day 6 versus those at day 10) [33]. An upset plot for the enrichment results was depicted in Fig. 4B.

3.5. Exploration of the association between hub genes and clinical characteristics

To understand the biological function of hub genes, GSEA was performed according to the Pearson correlation coefficients between all genes and hub genes, finding that ERN1 was correlated with humoral immune response mediated by circulating immunoglobulin, interferon- gamma production and regulation of acute inflammatory response and that CEP55 was enriched in cell cycle DNA replication, attachment of spindle microtubules to kinetochore and regulation of chromosome segregation in Gene Ontology (Fig. 5A). KEGG analysis showed that ERN1 was most regulated in RNA transport, Protein processing in en- doplasmic reticulum and Autophagy-animal and that CEP55 was most up-regulated with RNA transport, Ribosome, and Spliceosome. (Fig. 5B). Furthermore, we explored the correlations between clinical informations and hub genes. As shown in Fig. 6, the heatmap was split

based on K-means cluster and the higher expression values of ERN1 and CEP55 contributed to higher tumor stages and shorter OS and PFI time. The clinical characteristics of 77 ACC samples were displayed in Table 1. Moreover, we compared the expression levels of both hub genes separately in different traits (Fig. 7).

3.6. Identification of the prognostic implications of hub genes

The difference in expression abundance of hub genes between ACC and normal control samples was first compared in a boxplot (Fig. 8A). A ROC analysis was conducted to compare the precision for OS and PFI in Fig. 8B, which presented that the ERN1 (AUC for 1-, 3-, and 5-year PFI: 0.707, 0.734, and 0.721) can give a more accurate prediction for the 1-, 3-, and 5-year PFI compared with gender, age, T stage, N stage and M stage and that CEP55 (AUC for 1-, 3-, and 5-year OS: 0.839, 0.927, and 0.872; AUC for 1-, 3-, and 5-year PFI: 0.797, 0.840, and 0.826) per- formed the best in predicting both OS and PFI compared with other clinical traits. We conducted a univariable Cox analysis to identify the prognostic implications of hub genes, and both hub genes (ERN1, CEP55) showed a significant (P < 0.05) relationships with OS and PFI. Kaplan-Meier analysis giving the same results was also implemented

Fig. 3. Acquisition of immune-related hub genes. (A) Gene significance versus Module membership. The x-axis stands for the biweight midcorrelation coefficients of genes with the corresponding module, the y-axis represents the biweight midcorrelation coefficients of genes with corresponding immune cell abundance. (B) The top two clusters defined by MCODE of PPI network. The node color from yellow, orange to red, and the node size from smaller to bigger represent the degree of a node in the network from lower to higher.

bicor = 0.72, p = 4.3e-66

|bicor = 0.46, p = 1.4e-63

0.5 -bicor = 0.49, p = 7.7e31

bicor = 0.8, p = 1.3e-111

0.6

Dendritic cells activated

0.4

0.4

00

Macrophages MO

NK cells activated

T cells CD4 memory resting

0.4

0.3

0.3

0.4

0.2.

0.2

0.2

0.2

0.1

0.1

0.0

0.0

0.0

0.0

O

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.4

0.6

0.8

0.3

0.4

0.5

0.6

0.7

0.8

0.3

0.4

0.5

0.6

0.7

0.8

green Module

turquoise Module

brown Module

brown Module

(IF18/

FEN

DC25A

B

SKA

CENPM

PLK4

RIP1

WIN

PRC

MCM6

DEPD

CAP

YMS

WDHD1

COND

KIF4A

CNB

PX2

DHFR

SP CBP

DC80

DIPS

MYBL2

HRH3

ITR10

DRD4

HIS

T1H2

T1H2

BARD

CNA

CNB

DCA

CORT

HI T1H2

SSTRS

IST1H2BC

MCM5

KIF15

TK1

E2F2

HI

ST1H2

ST1H3H

CAPO

MELK

TASIR

MUR2

GTSE1

BIRCE

DC25

FAM83D

I

HIS

T1H2

ST1H3F

TROS

CEP55

SPN

GABBR

KIF2

XXIII

KI.

AO

01

CXCL3

CKS2

CASC5

A

ST1H3

ST1H2AB

AD2

SF1

SUCNR

M

LDCAB

JBE2

RXFP3

HI

MCM2

CKAP2L

ST1H3

ST1H41

URKB

OP2

BUB

NEK2

GRM8

$1PR5

KIF11

HIST1H

ST1H4D

KIFC

COCE

BINS?

POLQ

GALR

SSTR1

RBL1

H2AF

T1H2AL

CDT

DTL

CENPL

SPL

GINS1

SSTR3

GNG13

IST1HAST1H2A

ANLN

DEPDC1B

KIF1

CENPH

CDK2

Fig. 4. Estimation of immune correlation of hub genes. (A) Heatmap of correlated matrix between the immune cell and diverse cancer. (B) The upset plot for the enrichment in MSigDB immunologic gene set where the y-axis represented the correlation with the hub gene. ERN1 (M9453: Genes down-regulated in CD8 T cells during chronic infection with LCMV-Clone 13: effectors at day 15 versus exhausted at day 30. M3345: Genes up-regulated in comparison of dendritic cells (DC) before and 4 h after LPS (TLR4 agonist) stimulation. M6926: Genes down-regulated in B lymphocytes: control versus anti IgM and cell anti IgM and CpG oligodeox- ynucleotide 1826. M4817: Genes up-regulated in comparison of B cells versus monocyte macrophages. M8444: Genes down-regulated in lung innate lymphoid cells versus spleen CD4 [GeneID = 920] T cells). CEP55 (M3580: Genes up-regulated in comparison of wild type CD8 effector T cells at day 6 versus those from mice defficient for TRAF6 [GeneID = 7189] at day 10. M3577: Genes up-regulated in comparison of wild type CD8 effector T cells at day 6 versus those at day 10. M5041: Genes up-regulated in comparison of splenic primary CD8 effector T cells at day 8 post-acute infection versus splenic secondary CD8 effector T cells at day 8 post- acute infection. M2984: Genes up-regulated in B lymphocytes: control versus stimulated by anti-IgM for 12 h. M2961: Genes down-regulated in natural T reg versus T conv).

A

ERN1

B

ERN1

Act CDB

Tem CDB

1.00

Tem CD8

Act CD4

Tom CD4

Tem CD4

Tff

0.75

Tgd

Thi

Th17

Th2

Treg

0.50

Act B

Immm B

Mem B

NK

CD56bright

0.25

CD56dim

MDSC

NKT

Act DC

PDC HỌC

M9453

M3345

Macrophage

Eosinophil

M6926

Mast

Monocyte

M4817

Neutrophil

M8444

CEP55

CEP55

Act CDB

Tem CDB

1.0

Tem CDB

Act CD4

Tcm CD4

Tem CD4

Tth

0.8

Tgd

Th17

Th2

.

Treg

Act B

Imm B

0.6

Mam B

NK

CD56bright

CD56dim

MDSC

NKT

Act DC

M3580

PDC

IDC

M3577

Macrophage

Eosinophil

M5041

Mast

Monocyte

M2984

Neutrophil

M2961

Fig. 5. Exploration of hub gene function (A) Gene and interaction of the significant biological process in GO were shown. (B) GSEA plot for KEGG analysis to investigate the pathway enrichment for hub genes.

A

ERN1

CEP55

humoral immune response mediated by circulating immunoglobulin

size

nuGe &YBR ANA replicati

size

27

22

39

34

50

47

62

59

interferon-gamma production

Cor

Cor

regulation of complement activation

-0.2

0.8

regulation of acute inflammatory response

-0.3

attachment of spindle microtubules to kinetochore

-0.4

0.6

-0.5

regulation of chromosome segregation

0.4

B

Running Enrichment Score

0.5

ERN1

CEP55

Autophagy - animal

Running Enrichment Score

0.6

. Ribosome

0.4

Protein processing in endoplasmic reticulum

· RNA transport

RNA transport

0.3

0.4

Spliceosome

0.2

0.2

0.1

0.0

0.0

Ranked list metric

1.0

Ranked list metric

1.0

0.5

0.5

0.0

0.0

-0.5

-0.5

5000

10000

15000

5000

10000

15000

Rank in Ordered Dataset

Rank in Ordered Dataset

Fig. 6. The expression level of hub genes asso- ciated with clinical characteristics was also de- picted in a heatmap; OS, Overall Survival, stands for the period from the date of diagnosis until the date of death from any cause; PFI, Progression- Free Interval, stands for the period from the date of diagnosis until the date of the first occurrence of a new tumor event.

-

P

- Z

5

1

5

E

-

L

5

:

-

A

1

A

1

ERN1

CEP55

Gender

Age

T Stage

N Stage

M Stage

Tumor Stage

OS status

OS

PFI status

PFI

log2(FPKM+1)

Gender

Age

T Stage

N Stage

☐ Female

☐ Male

☐ T1

T2

☐ T3

☐ T4

☐ NO

☐ N1

012345

0 20 40 60 80

M Stage

Tumor Stage

OS status

☐ MO ☒ M1

☐ Stage I

☐ Stage II

☐ Stage

☐ Stage IV

☐ Censored ☐ Event

OS

PFI status

PFI

☐ Censored ☐ Event

0

1000

2000

3000

4000

5000

0

1000

2000

3000

4000

5000

Table 1 Clinical characteristics of 77 adrenocortical carcinoma samples from TCGA. The group of lower and higher was based on median expression levels of corresponding hub genes. Categorical variables were tested by Chi-Squared and continuous variables were tested by Analysis of Variance.
LevelOverallERN1CEP55
Lower expressionHigher expressionpLower expressionHigher expressionp
N7739383938
Gender (%)Female48 (62.3)19 (48.7)29 (76.3)0.02427 (69.2)21 (55.3)0.303
Male29 (37.7)20 (51.3)9 (23.7)12 (30.8)17 (44.7)
Age (mean (SD))46.61 (15.83)46.18 (14.83)47.05 (16.98)0.81145.77 (16.44)47.47 (15.35)0.64
T stage (%)T19 (11.7)6 (15.4)3 (7.9)0.6148 (20.5)1 (2.6)< 0.001
T242 (54.5)22 (56.4)20 (52.6)28 (71.8)14 (36.8)
T38 (10.4)3 (7.7)5 (13.2)2 (5.1)6 (15.8)
T418 (23.4)8 (20.5)10 (26.3)1 (2.6)17 (44.7)
N stage (%)N068 (88.3)35 (89.7)33 (86.8)0.96736 (92.3)32 (84.2)0.453
N19 (11.7)4 (10.3)5 (13.2)3 (7.7)6 (15.8)
M stage (%)M062 (80.5)33 (84.6)29 (76.3)0.52837 (94.9)25 (65.8)0.003
M115 (19.5)6 (15.4)9 (23.7)2 (5.1)13 (34.2)
Tumor stage (%)Stage I9 (11.7)6 (15.4)3 (7.9)0.6568 (20.5)1 (2.6)< 0.001
Stage II37 (48.1)19 (48.7)18 (47.4)25 (64.1)12 (31.6)
Stage III16 (20.8)8 (20.5)8 (21.1)4 (10.3)12 (31.6)
Stage IV15 (19.5)6 (15.4)9 (23.7)2 (5.1)13 (34.2)
OS1 status (%)Censored50 (64.9)32 (82.1)18 (47.4)0.00334 (87.2)16 (42.1)< 0.001
Event27 (35.1)7 (17.9)20 (52.6)5 (12.8)22 (57.9)
OS[days] (mean (SD))1462.16 (1052.39)1620.44 (1044.18)1299.71 (1049.59)0.1831972.03 (1057.04)938.87 (756.40)< 0.001
PFI2 status (%)Censored37 (48.1)26 (66.7)11 (28.9)0.00227 (69.2)10 (26.3)< 0.001
Event40 (51.9)13 (33.3)27 (71.1)12 (30.8)28 (73.7)
PFI[days] (mean (SD))1079.79 (1088.49)1364.62 (1013.16)787.47 (1098.26)0.0191603.51 (1179.68)542.29 (646.99)< 0.001

1: Overall Survival, the period from the date of diagnosis until the date of death from any cause.

2: Progression-Free Interval, the period from the date of diagnosis until the date of the first occurrence of a new tumor event.

based on a group defined by the median gene expression levels in Fig. 8C.

3.7. Verification of the hub genes

To validate our findings, the Oncomine database was searched for the differential expression between ACC and normal control samples, which showed that ERN1 (P < 0.0001) was down-regulated, and CEP55 (P < 0.0001) was up-regulated in ACC (Fig. 9A). Additionally, GSE49278, GSE19750, and GSE76021 were applied to prove our find- ings with a Kaplan-Meier analysis for both hub genes. GSE49278 showed the significance of OS for ERN1 (P < 0.01) and CEP55 (P < 0.001) (Fig. 9B). Nevertheless, GSE19750 (OS: ERN1 P = 0.13; CEP55 P < 0.0001) (Fig. 9C) and GSE76021 (PFS: ERN1 P = 0.33; CEP55 P < 0.001) (Fig. 9D) presented that only CEP55 was sig- nificantly associated with survival.

4. Discussion

Adrenocortical carcinoma is a scarce heterogeneous malignancy with unclear tumorigenesis mechanisms and a poor prognosis [34]. In recent years, the discovery of TME has paved the way for enormous progress in the understanding of tumor progression and im- munotherapy. Accumulating evidence has presented that im- munotherapy can serve a key role in treating diverse tumors [35]. Therefore, discernment of TME in ACC could lead to a better under- standing of ACC proliferation, relapse, and metastasis, and could pro- vide novel biomarkers for diagnosis, therapy, and recurrence risk eva- luation. In the current study, we identified two immune-related hub genes by combining coexpression network analysis and PPI network analysis using the TCGA ACC dataset and GTEx adrenal gland dataset; and further validation for the prognostic implications of hub genes was conducted with GSE49278, GSE19750, and GSE76021 data.

Previous studies have presented immunity as a crucial part of the TME [36], and immune cells in TME have been shown to have a con- siderable influence on tumor progression and clinical outcomes. In the presented study, we found 4 prognostic-related immune cells (T cells

CD4 memory resting, NK cells activated, Macrophages M0, and Den- dritic cells activated). CD4 + T cells with class II restricted and tumor specific indeed innately permeate in tumors and make up the tumor microenvironment, and CD4 + T cells can serve anticancer functions with the participation of CD8 + T cells or, in some cases, by directly secreting type 1 cytokines or killing tumor cells following the identifi- cation of processed antigens on the surface of tumors, moreover, The infiltration of memory T cells in TME is negatively related to the relapse risk [37]. NK-cell-mediated tumor killing can be prohibited by TME with the suppression of the effector NK-cell function and the selection/ editing of poorly immunogenic tumor cells [38]. M0 macrophages can be polarized to pro-inflammatory (M1) or anti-inflammatory (M2) phenotypes by disparate stimulus [39], and tumor-associated macro- phages whose phenotypes display different roles in tumorigenesis can mediate tumor progression and therapeutic resistance of multiple cancer types [8]. Dendritic cells are the basics of the immune responses against tumor cells, however, tumors can restrain normal Dendritic cells function by influencing the metabolism of Dendritic cells in these ways, Nutrient competition and hypoxia, persistent activation of un- folded protein response, and lipid uptake [40]. Therefore, these 4 im- mune cells are vital important in tumors; For example, T cells CD4 memory resting has been revealed to be associated with CMTM6, a regulator of PD-L1 protein in lung adenocarcinoma [41]. Multiple myeloma has been reported to suppress NK cells activated by promoting major histocompatibility complex class I-related chain A(MICA) shed- ding [42]; Macrophages M0 have been also identified to be associated with hepatocellular carcinoma [43], in addition, acute lymphoblastic leukaemia cells can induce dendritic cells with immunosuppressive features by up-regulating BMP5 [44]. All these reports collectively strengthen our findings.

Next, we derived 2089 immune-related genes of ACC in 3 modules by combining DESeq2 with WGCNA, which was further validated to correlate with immunity by GSEA. Furthermore, two hub genes (ERN1, CEP55) were identified by high connected genes in the co-expression network and important clusters in the PPI network. xCell (Supplement Fig. 1A), TISIDB (Fig. 4A and Supplement Fig. 1B-C) and MSigDB (Fig. 4B) database were searched to test the immune-relation of both

Fig. 7. The correlation between hub genes and various clinical traits. The categorical variable with two levels was tested by t-test, and the categorical variable with more than two levels was tested by Analysis of variance. The continuous variable was correlated with gene expression levels by Pearson correlation coefficient.

ERN1

5

P=0.033

5

P-0.333

5

P-0.494

5

P=0.147

5

P=0.344

4.

·

4

4.

4.

S

.

4.

@-3.23

@=2.95

- 3.09

3

3

=2.73

3

@-2.7

3.

@-2.65

A-3.09

3

@-2.72

@-2.71

@-2.93

@-2.4

@-2.23

@ =2.77

@7-2.23

2

2

2

2

O

2

.

1

1

1

1

1

Female (n = 48)

Male

T1

T2 (n = 42)

T3 (n = 8)

T4

NO

N1 (n = 9)

MO

M1

Stage II (n = 37)

Stage III (n = 16)

Stage IV (n = 15)

Gender

(n = 29)

(n = 9)

T Stage

(n = 18)

(n = 68)

Stage I (n =9)

N Stage

(n = 62)

M Stage

(n = 15)

Tumor Stage

5

P=0.003

5

P=5.9€-04

5

P=0.609

·

.

5

P=0,37

5

.

.

·

0.044

A

·

·

4.

4

4

4

4

.

.

mean = 46.61

mean = 1462.16

mean

1079.79

3

f-3.21

3

@-3.14

3

3

O

3

.

.

mean =2.73

mean = 2.73

mean = 2.73

.

D

2

@-2.47

2

₱=2.29

2

.

2

2

.

..

..

.

1

1

1

1

1

.

.

Censored (n = 50)

Event

Censored (n = 37)

Event (n = 40)

20

40

60

80

0

1000

2000

3000

0

2000

3000

(n = 27)

4000

1000

4000

CEP55

OS status

PFI status

Age

OS

PFI

5- P=0.797

5-

P=9.4e-06

5

P-0.390

5

P=3.8e-04

5

P=1.3e-05

4.

·

4

4

4.

4

3

3

@-2.82

3

3

3

F

=2.89

1

-2.05

=2.89

2

F=1.74

2

@=2.12

2

@-1.73

£=2.04

2

P=1.5

2

@-1.81

A-1.1

@-1.1

-1.35

1

1

1

·

1

1

·

1

1.39

0

Female (n = 48)

Male

0

T1

T2 (n = 42)

T3 (n = 8)

T4

0

(n = 18)

NO

N1 (n =9)

0

MO

M1 (n = 15)

D

(n = 29)

(n =9)

Stage I (n =9)

Stage III (n = 16)

Stage IV (n =15)

Gender

(n = 68)

Stage II (n = 37)

T Stage

N Stage

(n = 62)

M Stage

Tumor Stage

5

P=2.2c-06

5

P=4.3c-06

5

P=0.687

5.

P=4,6c-06

P=1.4c-06

4.

4

.

4

4

·

4

.

.

.

3

3

3

3

.

mean = 46.61

0

mean

1482.16

mean - 1079.79

@-2.71

.

2

2

mean = 1.77

mean = 1.77

2

2

P=2.32

2

mean

1.77

w

@=1.17

1

..

·

1.

F=1.26

1

1

·

·

0

D

0

Censored (n = 50)

Event

0

.

(n = 27)

Censored (n = 37)

Event (n = 40)

0

20

40

60

80

-1

0

1000

2000

3000

4000

0

1000

2000

3000

4000

OS status

PFI status

Age

Os

PFI

hub genes. At the same time, GSEA was applied to reveal the immune- related biological function and pathway of both hub genes (Fig. 5A-B). The relationship between the hub genes and the clinical characteristics was also explored (Fig. 6). Oncomine was applied to validate the im- portance in ACC compared with normal control samples. ROC and Kaplan-Meier analysis were used to validate the prognostics implica- tions of both hub genes in ACC with GEO database.

ERN1 encodes the transmembrane protein kinase inositol-requiring enzyme 1 with two functional catalytic domains, a serine/threonine- protein kinase domain and an endoribonuclease domain, which serves as a sensor of unfolded proteins in the endoplasmic reticulum and can activate an intracellular unfolded protein response signaling pathway [45]. Especially, ERN1 can protect KRAS mutant colon cancer from synthetic lethality by MAPK pathway [46], and could act as a ther- apeutic target for triple-negative breast cancer [47]. Moreover, ERN1 has been reported to involve in the treatment of Myc-driven cancers [48], and was demonstrated to lead to the progression of glioblastoma multiform [49]. These findings collectively verify the importance of ERN1 in tumor progression.

The high expression of CEP55 is related to genomic instability [50]; CEP55 protein, which mainly exists in tumor extracellular vesicles [51] and is a key regulator of cytokinesis, is unique and important to cancer and is a hallmark of cancer. Concretely, the high expression of CEP55

can trigger prosurvival signaling pathways, resulting in the prolifera- tion and migration of tumor cells [52]. CEP55 was found to involve in dismal prognosis and immune-related biomarker for liver cancer [53], moreover, CEP55 has been reported to related to the prognosis in ACC [54], which suggested CEP55 were able to serve a biomarker for ACC in diagnostic and prognosis.

To our knowledge, our research is the first to comprehensively analyze the prognosis of immune-related genes based on TCGA, and GTEx database, moreover, our finding that the relationship with im- munity and the prognosis implications of both hub genes was further validated by xCell, TISIDB, MSigDB, Oncomine and GEO database. Nevertheless, some limitations still existed in our study. First, larger samples are necessary to further strengthen our findings. Second, the exact mechanism of both hub genes had not been investigated in the present study. Hence, future efforts in exploring the specific molecular mechanism of ERN1 and CEP55 in ACC in vitro and in vivo based on a larger sample size are needed.

In conclusion, both immune-related hub genes were identified to play a vital role in ACC proliferation, invasion, and metastasis. ERN1 with a significant influence on PFI of ACC may be of crucial importance in ACC progress, and CEP55 which performed the best in predicting both PFI and OS compared with other clinical characteristics was more important for ACC prognosis. Both hub genes can provide a better

Fig. 8. The prognostic implications of the hub gene in ACC. (A) the difference of hub gene expression between ACC and normal control samples; (B) Time-dependent ROC curve for 1-, 3-, and 5-year for in OS and PFI among hub genes and other clinical characteristics; (C) Kaplan-Meier analysis for hub genes with median expression levels as the cut-off.

A

B

17

ERN1

OS

ROC curve for 1 year

ROC curve for 3 year

ROC curve for 5 year

1.00

1.00

1.00

15

0.75

0.75

0.75

True positive rate

13

0.50

0.50

0.50

..

ERN1(AUC=0.592)

CEP55(AUC=0.839)

ERNI(AUC=0.576)

CEP55(AUC=0.927)

ERN1(AUC=0.694)

CEP55(AUC=0.872)

0.25

gender(AUC=0.488)

0.25

gender(AUC=0.493)

0.25

gender(AUC=0.474)

11.

age(AUC=0.607)

T.stage(AUC=0.699)

age(AUC=0.549)

N.stage(AUC=0.438)

T.stage(AUC=0.857)

age(AUC=0.603)

N.stage(AUC=0.562)

T.stage(AUC=0.834)

..

M.stage(AUC=0.503)

tumor.stage(AUC=0.611)

M.stage(AUC=0.691)

N.stage(AUC=0.545)

tumor.stage(AUC=0.818)

M.stage(AUC=0.751)

0.00

0.00

tumor.stage(AUC=0.864)

NT

TP

0.00

0.25

False positive rate

0.60

0.75

0.00

1.00

0.00

0.25

False positive rate

0.50

0.75

1.00

0.DO

0.25

False positive rate ROC curve for 5 year

0.60

0.75

1.00

CEP55

PFI

ROC curve for 1 year

ROC curve for 3 year

1.00

1.00

1.00

11

0.75

0.75

0.75

True positive rate

9

0.50

0.50

0.50

ERN1(AUC=0.707)

ERN1(AUC=0.734)

ERN1(AUC=0.721)

7

CEP55(AUC=0.797)

CEP55(AUC=0.840)

CEP55(AUC=0.826)

0.25

gender(AUC=0.636)

0.25

age(AUC=0.500)

gender(AUC=0.521)

0.25

age(AUC=0.574)

gender(AUC=0.510)

T.stage(AUC=0.679)

T.stage(AUC=0.722)

age(AUC=0.546)

N.stage(AUC=0.577)

N.stage(AUC=0.568)

T.stage(AUC=0.747)

M.stage(AUC=0.618)

M.stage(AUC=0.665)

N.stage(AUC=0.585)

5

0.00

tumor.stage(AUC=0.731)

M.stage(AUC=0.648)

tumor.stage(AUC=0.804)

0.00

0.00

tumor.stage(AUC=0.752)

0.00

NT

TP

0.25

False positive rate

0.50

0.75

1.00

0.00

0.25

False positive rate

0.50

0.75

1.00

0.00

0.25

False positive rate

0.50

0.75

1.00

C

ERN1 ++ Low ++ High

ERN1 + Low ++ High

CEP55 + Low + High

CEP55 + Low + High

1.00

4

J

-

Survival probability

0.75

-

-

5

0.50

-

0.25

p = 0.0015

-

p = 0.00013

p < 0.0001

p < 0.0001

0.00

0

1000

2000

3000

0

1000

2000

3000 0

1000

2000

3000

0

1000

2000

3000

OS(days)

PFI(days)

OS(days)

PFI(days)

Fig. 9. Validation of the survival-relationships for hub genes. (A) differential expression of hub genes using the Oncomine database. Log-rank test for the relationship with survival of hub gene with median expression levels as cut-off using dataset (B) GSE49278 (C) GSE19750 (D) GSE76021.

A

B

ERN1 + Low + High

C

ERN1 + Low + High

D

ERN1 + Low + High

0.5

1.00-

1.00-

1.00-

log2 median-centered intensity

Survival probability

.

0.75

0.75

0.75

0.0

0.50

0.50

0.50

ERN1

.

Gene Rank: 979 (in top 6%)

0.25

p = 0.0026

0.25

p = 0.13

0.25

P-value: 2.45E-6

p = 0.33

t-Test: - 6.175

0.00

0.00

0.00

-0.5

0

40

80

120

160

0

5

10

15

20

0

4

8

PFS(years)

12

1

2

OS(months)

OS(years)

16

1.5

CEP55

CEP55 + Low + High

CEP55 + Low + High

CEP55 + Low + High

Gene Rank: 80 (in top 1%)

P-value: 2.87E-11

.

log2 median-centered intensity

1.0

t-Test: 9.548

1.00-

1.00-

1.00-

Survival probability

0.75

0.75

0.75

0.5

0.50

0.50

0.50

0.0

.

0.25

p = 0.00067

+

0.25

p

0.0001

0.25

p = 0.00038

.

-0.5

1

2

0.00

0.00

0.00

0

40

80

120

160

0

5

10

15

20

0

4

8

12

OS(months)

OS(years)

PFS(years)

16

insight into TME in ACC and serve as potential biomarkers for im- munotherapy in ACC.

Data sharing and data accessibility

The datasets analyzed for this study can be found in the GTEx (https://www.gtexportal.org/home/datasets), TCGA (https://portal. gdc.cancer.gov/), and GEO (https://www.ncbi.nlm.nih.gov/geo/) da- tabases.

Funding

This work was supported by the National Natural Science Foundation of China (program no. 81572543) and the Science and Technology Support Program of Tianjin, China (Grant No. 17ZXMFSY00040).

CRediT authorship contribution statement

Yun Peng: Methodology, Formal analysis, Writing - original draft. Yuxuan Song: Conceptualization, Writing - review & editing. Jin Ding: Data curation, Investigation, Writing - review & editing. Nan Li: Data curation, Visualization. Zheyu Zhang: Data curation, Software, Validation. Haitao Wang: Supervision, Project administration, Writing - review & editing.

Declaration of Competing Interest

None.

Acknowledgements

Not applicable.

Appendix A. Supplementary material

Supplementary data to this article can be found online at https:// doi.org/10.1016/j.intimp.2020.106930.

References

[1] G. Assié, et al., Integrated genomic characterization of adrenocortical carcinoma, Nat. Genet. 46 (6) (2014) 607-612, https://doi.org/10.1038/ng.2953.

[2] M.A. Habra, et al., Epidemiological risk factors for adrenocortical carcinoma: A hospital-based case-control study, Int. J. Cancer 146 (7) (2020) 1836-1840, https://doi.org/10.1002/ijc.32534.

[3] J.P. Luton, et al., Clinical features of adrenocortical carcinoma, prognostic factors, and the effect of mitotane therapy, N. Engl. J. Med. 322 (17) (1990) 1195-1201, https://doi.org/10.1056/NEJM199004263221705.

[4] T. Else, et al., Adrenocortical carcinoma, Endocr. Rev. 35 (2) (2014) 282-326, https://doi.org/10.1210/er.2013-1029.

[5] L. Ardolino, A. Hansen, S. Ackland, A. Joshua, Advanced Adrenocortical Carcinoma (ACC): a review with focus on second-line therapies, Horm. Cancer (2020), https:// doi.org/10.1007/s12672-020-00385-3.

[6] G. Assié, et al., Value of molecular classification for prognostic assessment of adrenocortical carcinoma, JAMA Oncol. 5 (10) (2019) 1440, https://doi.org/10. 1001/jamaoncol.2019.1558.

[7] M. Fassnacht, et al., Limited prognostic value of the 2004 International Union Against Cancer staging classification for adrenocortical carcinoma, Cancer 115 (2) (2009) 243-250, https://doi.org/10.1002/cncr.24030.

[8] J.A. Joyce, J.W. Pollard, Microenvironmental regulation of metastasis, Nat. Rev. Cancer 9 (4) (2009) 239-252, https://doi.org/10.1038/nrc2618.

[9] S. Brassart-Pasco, S. Brézillon, B. Brassart, L. Ramont, J .- B. Oudart, J.C. Monboisse, Tumor microenvironment: extracellular matrix alterations influence tumor pro- gression, Front. Oncol. 10 (2020), https://doi.org/10.3389/fonc.2020.00397.

[10] B. Somarouthu, S.I. Lee, T. Urban, C.A. Sadow, G.J. Harris, A. Kambadakone, Immune-related tumour response assessment criteria: a comprehensive review, Br. J. Radiol. 91 (1084) (2018), https://doi.org/10.1259/bjr.20170457.

[11] “GDC TCGA portal.” https://portal.gdc.cancer.gov/ (accessed Mar. 27, 2020).

[12] “GTEx Portal.” https://www.gtexportal.org/home/ (accessed Mar. 27, 2020).

[13] J. Liu, et al., An integrated TCGA pan-cancer clinical data resource to drive high- quality survival outcome analytics, Cell 173 (2) (2018) 400-416.e11, https://doi. org/10.1016/j.cell.2018.02.052.

[14] D.R. Zerbino, et al., Ensembl 2018, Nucleic Acids Res. 46 (D1) (2018) D754-D761, https://doi.org/10.1093/nar/gkx1098.

[15] Homo sapiens - Ensembl genome browser 99. https://asia.ensembl.org/Homo_ sapiens/Info/Index (accessed Mar. 27, 2020).

[16] R Core Team, R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing, 2020.

[17] A.M. Newman, et al., Robust enumeration of cell subsets from tissue expression profiles, Nat. Methods 12 (5) (2015), https://doi.org/10.1038/nmeth.3337.

[18] B. Chen, M.S. Khodadoust, C.L. Liu, A.M. Newman, A.A. Alizadeh, Profiling tumor infiltrating immune cells with CIBERSORT, Methods Mol. Biol. Clifton NJ 1711 (2018) 243-259, https://doi.org/10.1007/978-1-4939-7493-1_12.

[19] P. Langfelder, S. Horvath, WGCNA: an R package for weighted correlation network analysis, BMC Bioinf. 9 (1) (2008) 559, https://doi.org/10.1186/1471-2105-9-559.

[20] B. Zhang, S. Horvath, A general framework for weighted gene co-expression net- work analysis, Stat. Appl. Genet. Mol. Biol. 4 (2005), https://doi.org/10.2202/ 1544-6115.1128.

[21] A. Subramanian, et al., Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles, Proc. Natl. Acad. Sci. 102 (43) (2005) 15545-15550, https://doi.org/10.1073/pnas.0506580102.

[22] M. Kanehisa, S. Goto, KEGG: kyoto encyclopedia of genes and genomes, Nucleic Acids Res. 28 (1) (2000) 27-30, https://doi.org/10.1093/nar/28.1.27.

[23] The Gene Ontology Resource: 20 years and still GOing strong, Nucleic Acids Res. 47(D1) (2019) D330-D338, doi:10.1093/nar/gky1055.

[24] M. Ashburner, et al., Gene Ontology: tool for the unification of biology, Nat. Genet. 25 (1) (2000) 25-29, https://doi.org/10.1038/75556.

[25] G. Yu, L .- G. Wang, Y. Han, Q .- Y. He, clusterProfiler: an R package for comparing biological themes among gene clusters, OMICS J. Integr. Biol. 16 (5) (2012) 284-287, https://doi.org/10.1089/omi.2011.0118.

[26] G.D. Bader, C.W.V. Hogue, An automated method for finding molecular complexes in large protein interaction networks, BMC Bioinf. 4 (2003) 2, https://doi.org/10. 1186/1471-2105-4-2.

[27] D. Aran, xCell: digitally portraying the tissue cellular heterogeneity landscape, 2017, p. 14.

[28] TISIDB: an integrated repository portal for tumor-immune system interactions | Bioinformatics | Oxford Academic. https://academic.oup.com/bioinformatics/ article-abstract/35/20/4200/5418799?redirectedFrom =fulltext (accessed Apr. 30, 2020).

[29] A. Liberzon, C. Birger, H. Thorvaldsdóttir, M. Ghandi, J.P. Mesirov, P. Tamayo, The Molecular signatures database (MSigDB) hallmark gene set collection, Cell Syst. 1 (6) (2015) 417-425, https://doi.org/10.1016/j.cels.2015.12.004.

[30] Complex heatmaps reveal patterns and correlations in multidimensional genomic data | Bioinformatics | Oxford Academic. https://academic.oup.com/ bioinformatics/article/32/18/2847/1743594 (accessed Apr. 30, 2020).

[31] “ONCOMINE: a cancer microarray database and integrated data-mining platform. - PubMed - NCBI.” https://www.ncbi.nlm.nih.gov/pubmed/15068665 (accessed Apr. 30, 2020).

[32] M. Ceppi, et al., Ribosomal protein mRNAs are translationally-regulated during human dendritic cells activation by LPS, Immunome Res. 5 (2009) 5, https://doi. org/10.1186/1745-7580-5-5.

[33] Enhancing CD8 T-cell memory by modulating fatty acid metabolism. - PubMed - NCBI. https://www.ncbi.nlm.nih.gov/pubmed?CrntRpt =DocSum&cmd =search& term = 19494812 (accessed May 02, 2020).

[34] M. Fassnacht, R. Libé, M. Kroiss, B. Allolio, Adrenocortical carcinoma: a clinician’s update, Nat. Rev. Endocrinol. 7 (6) (2011) 323-335, https://doi.org/10.1038/ nrendo.2010.235.

[35] D.S. Chen, I. Mellman, Oncology meets immunology: the cancer-immunity cycle, Immunity 39 (1) (2013) 1-10, https://doi.org/10.1016/j.immuni.2013.07.012.

[36] Screening and Identifying Immune-Related Cells and Genes in the Tumor Microenvironment of Bladder Urothelial Carcinoma: Based on TCGA Database and … - PubMed - NCBI. https://www.ncbi.nlm.nih.gov/pubmed/32010623 (accessed May 02, 2020).

[37] S. Hadrup, M. Donia, P. thor Straten, Effector CD4 and CD8 T cells and their role in the tumor microenvironment, Cancer Microenviron 6 (2) (2012) 123-133, https:// doi.org/10.1007/s12307-012-0127-6.

[38] M. Vitale, C. Cantoni, G. Pietra, M.C. Mingari, L. Moretta, Effect of tumor cells and tumor microenvironment on NK-cell function, Eur. J. Immunol. 44 (6) (2014) 1582-1592, https://doi.org/10.1002/eji.201344272.

[39] Monocyte differentiation and macrophage polarization. https://vpjournal.net/ article/view/3016 (accessed Jul. 22, 2020).

[40] P. Giovanelli, T.A. Sandoval, J.R. Cubillos-Ruiz, Dendritic cell metabolism and function in tumors, Trends Immunol. 40 (8) (2019) 699-718, https://doi.org/10. 1016/j.it.2019.06.004.

[41] H. Wang, J. Gao, R. Zhang, M. Li, Z. Peng, H. Wang, Molecular and immune characteristics for lung adenocarcinoma patients with CMTM6 overexpression, Int. Immunopharmacol. 83 (2020) 106478, https://doi.org/10.1016/j.intimp.2020. 106478.

[42] Y. Wang, et al., BCMA-targeting bispecific antibody that simultaneously stimulates NKG2D-enhanced efficacy against multiple myeloma, J. Immunother. Hagerstown Md 1997 (2020), https://doi.org/10.1097/CJI.0000000000000320.

[43] Y. Zhang, L. Zhang, Y. Xu, X. Wu, Y. Zhou, J. Mo, Immune-related long noncoding RNA signature for predicting survival and immune checkpoint blockade in hepa- tocellular carcinoma, J. Cell. Physiol. (2020), https://doi.org/10.1002/jcp.29730.

[44] J. Valencia, et al., Acute lymphoblastic leukaemia cells impair dendritic cell and macrophage differentiation: Role of BMP4, Cells 8 (7) (2019) 14, https://doi.org/ 10.3390/cells8070722.

[45] PubChem, ERN1 - endoplasmic reticulum to nucleus signaling 1 (human). https://

pubchem.ncbi.nlm.nih.gov/gene/ERN1/human (accessed Jul. 22, 2020).

[46] T. Šuštić, et al., A role for the unfolded protein response stress sensor ERN1 in regulating the response to MEK inhibitors in KRAS mutant colon cancers, Genome Med. 10 (2018), https://doi.org/10.1186/s13073-018-0600-z.

[47] S.E. Logue, et al., Inhibition of IRE1 RNase activity modulates the tumor cell se- cretome and enhances response to chemotherapy, Nat. Commun. 9 (2018), https:// doi.org/10.1038/s41467-018-05763-8.

[48] H. Xie, et al., IRE1a RNase-dependent lipid homeostasis promotes survival in Myc- transformed cancers, J. Clin. Invest. 128 (4) (2018) 1300-1316, https://doi.org/10. 1172/JCI95864.

[49] S. Lhomond, et al., Dual IRE1 RNase functions dictate glioblastoma development, EMBO Mol. Med. 10 (3) (2018), https://doi.org/10.15252/emmm.201707929.

[50] CEP55 is a determinant of cell fate during perturbed mitosis in breast cancer_

Kalimutho et al_2018.pdf.

[51] F. Qadir, et al., Transcriptome reprogramming by cancer exosomes: identification of novel molecular targets in matrix and immune modulation, Mol. Cancer 17 (2018), https://doi.org/10.1186/s12943-018-0846-5.

[52] Y .- F. Yang, et al., SPAG5 interacts with CEP55 and exerts oncogenic activities via PI3K/AKT pathway in hepatocellular carcinoma, Mol. Cancer 17 (2018), https:// doi.org/10.1186/s12943-018-0872-3.

[53] L. Yang, Y. He, Z. Zhang, W. Wang, Upregulation of CEP55 predicts dismal prog- nosis in patients with liver cancer, BioMed Res. Int. 2020 (2020) 4139320, https:// doi.org/10.1155/2020/4139320.

[54] H. Xiao, D. Xu, P. Chen, G. Zeng, X. Wang, X. Zhang, Identification of five genes as a potential biomarker for predicting progress and prognosis in adrenocortical carci- noma, J. Cancer 9 (23) (2018) 4484-4495, https://doi.org/10.7150/jca.26698.