Identification of prognostic genes in adrenocortical carcinoma microenvironment based on bioinformatic methods
Xiao Li1 D Yang Gao2 Zicheng Xu1 1 Zheng Zhang3 3 Yuxiao Zheng1
Feng Qi1,4 İD
1Department of Urology, Jiangsu Cancer Hospital & Jiangsu Institute of Cancer Research & Affiliated Cancer Hospital of Nanjing Medical University, Nanjing, China
2Department of Radiology, Jiangsu Cancer Hospital & Jiangsu Institute of Cancer Research & Affiliated Cancer Hospital of Nanjing Medical University, Nanjing, China
3Hepato-Pancreato-Biliary Center, Zhongda Hospital, School of Medicine, Southeast University, Nanjing, China
4Department of Urology, The First Affiliated Hospital of Nanjing Medical University, Nanjing, China
Correspondence
Yuxiao Zheng, Department of Urology, Jiangsu Cancer Hospital & Jiangsu Institute of Cancer Research & Affiliated Cancer Hospital of Nanjing Medical University, Nanjing, 210009, China. Email: zheng_yuxiao@163.com
Feng Qi, Department of Urology, Jiangsu Cancer Hospital & Jiangsu Institute of Cancer Research & Affiliated Cancer Hospital of Nanjing Medical University, Nanjing, 210009, China; Department of Urology, the First Affiliated Hospital of Nanjing Medical University, Nanjing, 210029, China.
Email: qf199408@163.com
Funding information
Six Talents Peaks Program of Jiangsu Province, Grant/Award Number: 2016- WSW-021; National Natural Science Foundation of China, Grant/Award Number: 81702520; Medical Research Project of Jiangsu Provincial Health and Family Planning Commission, Grant/ Award Number: H2018052; Research Project of Jiangsu Cancer Hospital, Grant/ Award Number: ZM201801 and ZN201602; Young Talents Program of Jiangsu Cancer Hospital, Grant/Award Number: 2017YQL and -04
Abstract
Background: To identify prognostic genes which were associated with adrenocorti- cal carcinoma (ACC) tumor microenvironment (TME).
Methods and materials: Transcriptome profiles and clinical data of ACC samples were collected from The Cancer Genome Atlas (TCGA) database. We use ESTIMATE (estimation of stromal and Immune cells in malignant tumor tissues using expres- sion data) algorithm to calculate immune scores, stromal scores and estimate scores. Heatmap and volcano plots were applied for differential analysis. Venn plots were used for intersect genes selection. We used protein-protein interaction (PPI) networks and functional analysis to explore underlying pathways. After performing stepwise regression method and multivariate Cox analysis, we finally screened hub genes as- sociated with ACC TME. We calculated risk scores (RS) for ACC cases based on multivariate Cox results and evaluated the prognostic value of RS shown by receiver operating characteristic curve (ROC). We investigated the association between hub genes with immune infiltrates supported by algorithm from online TIMER database. Results: Gene expression profiles and clinical data were downloaded from TCGA. Lower immune scores were observed in disease with distant metastasis (DM) and locoregional recurrence (LR) than other cases (P = . 0204). Kaplan-Meier analysis revealed that lower immune scores were significantly associated with poor overall survival (OS) (P = . 0495). We screened 1649 differentially expressed genes (DEGs) and 1521 DEGs based on immune scores and stromal scores, respectively. Venn plots helped us find 1122 intersect genes. After analysing by cytoHubba from Cytoscape software, 18 hub genes were found. We calculated RS and ROC showed significantly predictive accuracy (area under curve (AUC) = 0.887). ACC patients with higher
Xiao Li, Yang Gao and Zicheng Xu contributed equally to this article.
RS had worse survival outcomes (P < . 0001). Results from TIMER (tumor immune estimation resource) database revealed that HLA-DOA was significantly related with immune cells infiltration.
Conclusion: We screened a list of TME-related genes which predict poor survival outcomes in ACC patients from TCGA database.
1 INTRODUCTION |
Adrenal tumors are very common, approximately 3%-10% of the human population suffers from the disease, and most of them are nonfunctional adrenocortical adenomas (benign lesions).1 However, adrenocortical carcinoma (ACC) is ex- tremely rare with an incidence of 0.7-2 cases per million peo- ple annually.2-4 Surgery is the first-line therapy for patients with localized ACC (stage I, II diseases and most of stage III diseases), while different therapy methods (molecularly targeted therapy, immunotherapy, adjuvant therapy, and so on) or multidisciplinary approaches are applied in advanced patients in order to achieve improvement in the prognosis.” The prognosis of patients with ACC is poor because most of them are present with metastatic or advanced diseases, in which surgical resection is not amenable. Additionally, about two-thirds of localized patients suffer recurrence and sys- temic therapy is always required.6-8 Hence, tumor burden is one of the most important prognostic factors in ACC patients. As reported, 5-year survival rate in patients with stage I dis- ease was about 80% while only 13% in patients with stage IV diseases.9
Tumor microenvironment (TME) includes not only the tumor cells themselves, but also the surrounding fibroblasts, immune and inflammatory cells, glial cells and other cells, as well as the intercellular substance, microvessels, and bio- molecules infiltrated in the nearby area.10-12 Nowadays, nu- merous studies focused on the relationship between TME and survival outcomes or recurrence. Ren et al reported that interaction between stromal cells and epithelial/cancer cells affected the cancer progress in pancreatic cancer.13 Jian et al14 discovered that higher TH1 was an important predictor of poor prognosis in patients with nonsmall cell lung cancer (NSCLC). However, studies on the role of TME in prognosis in adrenal tumors were few and most of them focused on ad- renocortical adenomas. Kitawaki et al15 firstly histologically investigated the TME of autonomous hormone-secreting ad- renocortical adenomas. He demonstrated that immune cell infiltration uniquely existed in cortisol-producing adrenocor- tical adenomas and immune cells enter the tumor site through the CXCL12-CXCR4 signaling and eliminate the aging tumor cells.
With the advent of the era of big data, transcriptome profiles and clinical data of cancer patients can be easily obtained from some public database such as The Cancer
Genome Atlas (TCGA) database, which includes abun- dant cancer-causing genomic alterations data of various malignancies,16 and many algorithms were created to ex- plore the tumor purity based on the TCGA database. 17,18 Additionally, Estimation of Stromal and Immune cells in Malignant Tumours using Expression data (ESTIMATE), a novel algorithm proposed by Yoshihara et al17 was uti- lized to calculate the fraction of immune and stromal cells in malignancy. ESTIMATE algorithm generates stromal score (that captures the presence of stroma in tumor tis- sue), immune score (that represents the infiltration of immune cells in tumor tissue), and estimate score (that infers tumor purity and equal in number to stromal score and immune score). Previous studies had reported that scores calculated by ESTIMATE algorithm could be uti- lized as a predictor of nontumor cell infiltration in TMEs of colon cancer, prostate cancer, and breast cancer. 19-21 However, no research discussed the TME of ACC using ESTIMATE algorithm. Hence, we conducted this study to identify prognostic genes which were associated with ACC TME.
2 METHODS AND MATERIALS |
2.1 Clinical characteristics and ESTIMATE data collection |
Level 3 gene transcriptome profiles were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc. cancer.gov/) and 79 ACC patients was enrolled in our study. RNA expression for ACC Multiforme was collected using IlluminaHiSeq (version: 2017-10-13). The following cor- responding clinical characteristics were extracted: sex, age, cancer type, status, topography, lymph node and metasta- sis (TNM) based on American Joint Committee on Cancer (AJCC), pathological stage and new event of the disease. Then, we download the ACC survival outcomes through the Genomic Data Commons (GDC) tool from TCGA portal. Normalization process was performed by utilizing package limma.22 Immune scores, stromal scores, and estimate scores were calculated by applying ESTIMATE algorithm which R script was downloaded from the website (https://sourceforge. net/projects/estimateproject/). The statistical results were shown in Table S1.
2.2 Survival analysis and correlation analysis
Unpaired t test was used to analysis the relationship between im- mune scores, stromal scores, estimate scores and AJCC-N and AJCC-M. Ordinary one-way Analysis of Variance (ANOVA) was used to analysis the relationship between immune scores, stromal scores, estimate scores and AJCC-T, pathological stage, new event of the disease. P < . 05 was considered as statistically significant. Kaplan-Meier (K-M) analysis which based on the immune scores, stromal scores, and estimate scores of ACC patients was performed through survival package.23,24 In K-M analysis we used dichotomous variable. Log-rank test was per- formed and P < . 05 represented a statistically significant result.
| |
2.3 Heatmaps and clustering analysis |
Enrolled samples were divided as high- and low-level groups according to the median of three scores from the ESTIMATE method. In immune score groups, the tran- scriptome data were analyzed by limma package under the standard of | log(FC) | > 1 and False Discovery Rate (FDR) < 0.05. After that, we used pheatmap package and clustering analysis to explore the differential high and low expression genes among the different immune score lev- els. We used volcano plot performed by ggplot2 package to exhibit significantly differentially expressed genes. The cut |log2FC| was 1 and cut P value was .05. Besides, we performed the same analysis procedure in two stromal score
| Male (%) | Female (%) | Total (%) | |
|---|---|---|---|
| Sex | 31(39.2%) | 48(60.8%) | 79(100%) |
| Age (years) | 48.74 ± 15.67 | 45.38 ± 15.69 | 46.70 ± 15.67 |
| Cancer type | |||
| Usual type | 29(36.7%) | 46(58.2%) | 75(94.9%) |
| Oncocytic type | 2(2.5%) | 1(1.3%) | 3(3.8%) |
| Myxoid type | 0(0.0%) | 1(1.3%) | 1(1.3%) |
| Status | |||
| Dead | 11(13.9%) | 17(21.5%) | 28(35.4%) |
| Alive | 20(25.3%) | 31(39.2%) | 51(64.6%) |
| T* | |||
| 1 | 3(3.8%) | 6(7.6%) | 9(11.4%) |
| 2 | 15(19.0%) | 27(34.2) | 42(53.2%) |
| 3 | 4(5.0%) | 4(5.0%) | 8(10.1%) |
| 4 | 7(8.9%) | 11(13.9%) | 18(22.8%) |
| N* | |||
| 0 | 27(34.2%) | 41(51.9%) | 68(86.1%) |
| 1 | 2(2.5%) | 7(8.9%) | 9(11.4%) |
| M* | |||
| 0 | 24(30.4%) | 38(48.1%) | 62(78.5%) |
| 1 | 5(6.3%) | 10(12.7%) | 15(19.0%) |
| Stage* | |||
| I | 3(3.8%) | 6(7.6%) | 9(11.4) |
| II | 15(19.0%) | 22(27.8%) | 37(46.8%) |
| III | 6(7.6%) | 10(12.7%) | 16(20.3%) |
| IV | 5(6.3%) | 10(12.7%) | 15(19.0%) |
| New event | |||
| Distant metastasis | 6(7.6%) | 20(25.3%) | 26(32.9%) |
| Locoregional recurrence | 5(6.3%) | 5(6.3%) | 10(12.7%) |
| New primary tumor | 1(1.3%) | 0(0.0%) | 1(1.3%) |
| No new event | 19(24.1%) | 23(29.1%) | 42(53.2%) |
*There existed a partial deletion of this clinical variable when it was collected.
|
A
2000
B
2000
F
Immue score
100
Immune score
1000
Immune score
1000
Percent survival
0
0
50
-1000
-1000
L
P =. 1573
P =. 0495
-2000
-2000
P =. 0787
H
0
T1
T2
T3
T4
NO
N1
0
1000
2000
3000
4000
5000
OS (d)
C
D
2000
G
2000
Stromal score
100
Immune score
1000
Immune score
1000
Percent survival
0
0
50
-1000
-1000
L
P =. 3064
-2000
P= . 8842
-2000
P= . 0664
H
0
MO
M1
S1
S2
S3
S4
0
1000
2000
3000
4000
5000
OS (d)
E
2000
H
100
ESTIMATE score
Immune score
1000
Percent survival
0
50
-1000
P =. 0204
L
-2000
P =. 2210
No LR or DM
LR
DM
H
0
0
1000
2000
3000
4000
5000
OS (d)
groups. In this way, we obtained intersect genes among im- mune scores and stromal scores. The result was exhibited by VennDiagram package.25
2.4 Construction of PPI network, GO analysis, KEGG, and GSEA |
The protein-protein interaction (PPI) network was con- structed from (search tool for recurring instances of neigh- bouring genes) STRING database26 and the result was visualized with Cytoscape software (version 3.7.1).27 The inside-software plugin cytoHubba was helped to identified
hub genes among interested genes with 12 different topologi- cal analysis methods.28 We calculated the network nodes and excluded the networks less than 10 nodes. After that, online tool DAVID (https://david.ncifcrf.gov/) was exploited to construct the gene ontology (GO) analysis.29 What’s more, Kyoto Encyclopedia of Genes and Genomes (KEGG) analy- sis was performed for pathway analysis using org.Hs.eg.db package, clusterProfiler, org.Hs.eg.db, enrichplot, and gg- plot2 packages. Bar charts were used to reveal the pathway annotation and enrichment. Top 20 pathway enrichment was analyzed and shown by bubble chart. The q-value < 0.05 was considered as significant. The results of Cellular Component (CC), Molecular Function (MF), and Biological Process (BP)
|
A
B
Volcano Plot
TOGA-PK-ASH8-01
TOGA-08-ASK1-01
TOGA-AB-ASLG-01
TOGA-08-ASLA-01
TOGA-AB-ASL9-01
TOGA-08-ASLA-01
TOGA-OB-ASL5-01
TOGA-OB=ASLP-01
10-LISY-80-YOUI
TOGA-AB-ASTY-01
TOGA-DUASJP-01
TOGA-P6-ASOF-01
TOGA-OB-A5 12-01
TOGA-18-ASJ5-01
TOGA-08-1510-01 TOGA-08-ASJS-01
TOGA-08-4513-01
TOGA-OB-ASL3-01
TOGA-08-ASJF-01
TOGA-08-AS 19-01
16.94
13.55-
Samplellame
Ch33
-log10(P-value)
GPE132
30. 16-
VCF1
THEP1
SPB65
IL12BB1
LILRE2
6.78-
1CF18
KLERI
GLASS
POM282
3.39-
RIM2
BTK
1HLEP3
ETE
ASGB2
0
AMICAL
-2.89
-1.39
0.1
1.6
3.1
4.6
PGS18
-4.39
TUFATESLO
log2(FoldChange)
PIPPC
Ch3E
CLECIÓA
CYLIP
MARSE
ISTI
C
MACKAPIL
ADAH
FRI3
0082
im-up
CD69
ILEZ
TRAERIP3
PORY13
002
PLEK
SASH3
TECADIOC
RASGERA
ARHGAP9
CLECZ4
I PAPS
THCS
302
907
273
CSEZPE
IFPAUT
100606524
CSE2BA
GauK
LIR
13/13
FL J30679
SPMA3
PISA
DAS234E
IOSES
SEÇİMİL
MCIH3 MC28
stro-up
Clearfil
CCACIS
ARNA
POPA
IMKIR
D
KISSC
MATSL
PUTZ
SLC30A10
im-down
PGSORP
PLUS
173-2
GPBEI
Clar 150
MC1019
ECCHC12
GREBIL
PHILIPEP2
18304
PLIVI
EPMAA
SLC26A10
DOPPIO
225
215
126
PERM12
ACCI3
M13
AATK
STIL2
CITED1
212536
KAR
SEX
MIP
MIPATCA
TRAC36
17km
CHEMLA
POXDA
stro-down
10XC10
DEM
were all collected and shown as barplot. We set “c2.cp.kegg. v6.2.symbols.gmt gene sets” as gene set database, “Illumina Human.chip” as chip platform when setting parameters, FDR < 0.25, |enriched score| > 0.35, and gene size ≥ 35 in Gene Set Enrichment Analysis (GESA) software.30
2.5 Overall survival curve and risk score |
K-M plots were generated to select prognostic hub genes. In K-M analysis we used a dichotomous variable. We used P < . 05 of log-rank test to evaluate the relationship. Besides, survival package and “for cycle” R script was performed to select prognostic hub genes. Then, Multivariate Cox regres- sion analysis was used to evaluate the risk score in the formula of risk score = E (B; x Expi) (i = the number of prognostic
hub genes). After calculating risk scores of ACC patients, we separated all patients in high- and low-risk groups according to the median score. The survival ROC package was used to conduct the receiver operating characteristic curve (ROC).31 In addition, we used K-M plot to visually present the dif- ferences in survival outcome between high-risk and low-risk groups. In K-M analysis we used dichotomous variable.
2.6 Timer database analysis |
To explore the specific association of hub genes with immune cells, we used the deconvolution algorithm provided by the TIMER (tumor immune estimation resource) database (https ://cistrome.shinyapps.io/timer/) for analysis. The immune cells involved in the analysis include: B cell, CD4+ T cell, CD8+ T
Open Access
A
B
CD33
C
D
CYBA
SLORA5
VAMP8
TRPM2
CD33
TMC6
HLA-DQB1
LA-DPA1
CD33
HLA-DRA
PTPRB
HVONT
ITGB2
CD3D
CD4
CD93
C3AR1
TMC6
CD53
VAMP8
CD4
HLA-DRA
HLA-DRB5
CD53
GI
32
RAB37
TMC6
LCK
ATP8B4
VAMP8
CYBA
LCK
CD3D
CCL4
ITGB2
C3XR1
GNG2
FICK
FOR
MCC
DMNC
MNC
Radiality
E
F
G
H
I
RAC2
ITG82
HLA-DRA
TMO6
ITGB2
CYBA
LAIR1
GNĘ2
CYBA
VAMP8
CD53
C3
RT
R1
GIÁP2
CYBA
CD3D
LCK
VAMP8
ITGB2
VAMP8
VIAMARE
C
RT
004
LOX
GN
32
CCL4
CD4
CD3D
GHG2
coke
C3ART
62
CD4
LCK
CD4
CD3D
HLA-DRA
VAMP8
LCK
HEK
HLA-DRA
TNF
HLA-DRA
HLA-DRB5
TNF
Betweenness
BottleNeck
ITGB2
Degree
EPC
Closeness
TNF
RAC2
J
K
L
TACRI
CCL4
HPĢDS
APĻNR
RAB37
FOLR2
HLA-DOA
HỌK
HLA-ORA
CD4
LCK
LCK
GNG2
TEXAS1
P2RY12
SLO2A5
PRND
PIK3CO
TER
HLA-DRA
CD3D
IFITM2
CD4
APLNR
PLÁŽGIA
VAMP8
ClusteringCoefficient
Stress
P2RY12
EcCentricity
cell, macrophage, neutrophil, dendritic cell. In addition, the asso- ciation between tumor purity and hub genes was also analyzed.
2.7 Statistical analysis |
Survival package was applied in multivariate Cox regression analysis and K-M analysis. Limma pack- age was applied in differential analysis. Unpaired t test and ordinary one-way Analysis of Variance (ANOVA) were used to analysis data from two or more groups. Statistical analysis was implemented on R software (ver- sion 3.5.2). When we used FDR to evaluate the results,
FDR < 0.05 indicated that the results were statistically significant. In addition, when we used P to evaluate the results, P < . 05 indicated that the results were statistically significant.
3 RESULTS |
3.1 Lower immune scores significantly associated with malignant progression of ACC |
We downloaded gene expression profiles data of 79 samples from TCGA, including 31 male cases (39.2%) and 48 female
|
2
2.5
3
3.5
4
A
B
C
-LOG(FDR)
-LOG(FGR)
T cell receptor signaling pathway
-log(FDR)
MHC class Il protein complex
antigen processing and presentation of peptide or
polysaccharide antigen via MHC class II
integral component of lumenal side of
T cell costimulation
endoplasmic reticulum membrane
MHC class Il receptor activity
immune response
transport vesicle membrane
interferon-gamma-mediated signaling pathway
peptide antigen binding
antigen processing and presentation of exogenous peptide antigen via MHC class II
clathrin-coated endocytic vesicle
membrane
0
2
4
6
8
10
12
14
0
2
4
6
8
10
12
0
0.5
1
1.5
BP
CC
MF
D
KEGG pathway annotation
G
Environmental Information Processing
Enrichment plot: KEGG_CHEMOKINE_SIGNALING_PATHWAY
Enrichment plot: KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTIO
Signal transduction-
153
Signaling molecules and interaction -
148
0.45
N
Membrane transport -12
Organismal Systems-
Enrichment score (ES)
0.40
0.4
Immune system
226
0.30
Enrichment score (ES)
Endocrine system-
64
99
Development
52
1.15
9.2
Nervous system-
31
0.10
Digestive system-
26
0.1
Sensory system
19
0.00
Circulatory system-
16
Environmental adaptation -8
Excretory system -
8
Ranked list metric (Signal2Noise)
Ranked list metric (Signal2Noise)
Aging -
3
Human Diseases
lı
1.0
f (paneely corelated)
Infectious diseases-
166
1.0
Cancers
90
0.6
Zera cross at 906
.A.
Immune diseases
72
00
VA
Zeno cruce at 90d
Cardiovascular diseases-
47
o
Endocrine and metabolic diseases-
43
“(negatively comelated)
0.
” (negatively cenelatee)
Drug resistance-
15
0
100
200
100
400
800
doo
1.000 1,100
@
100
200
300
400
500
000
800
00 1.000 1.100
Neurodegenerative diseases-
13
Rank in Ordered Dataset
Rank in Ordered Dataset
Substance dependence-
11
Enrichment profile
Hits
Ranking metric scores
Enrichment profile
Hits
Ranking metric scores
Cellular Processes-
Transport and catabolism-
62
Enrichment plot: KEGG_JAK_STAT_SIGNALING_PATHWAY
Enrichment plot: KEGG_LEUKOCYTE_TRANSENDOTHELIAL_MIGRATION
Cell growth and death-
$6
Cellular community - eukaryotes-
$6
0.5
Cell motility
17
Metabolism-
Enrichment score (ES)
0.4
Enrichment score (ES)
0.4
Global and overview maps
55
9.3
Lipid metabolism-
24
0.3
Glycan biosynthesis and metabolism-
15
93
02
Carbohydrate metabolism
11
Xenobiotics biodegradation and metabolism-
10
2.1
Q.1
Metabolism of cofactors and vitamins-9
Amino acid metabolism-
18
9.0
0.0
Metabolism of other amino acids-
Nucleotide metabolism-
16
Energy metabolism-
Ranked list metric (Signal2Noise)
13
Ranked list metric (Signal2Noise)
Biosynthesis of other secondary metabolites -
11
1.0
(positively comelade d)
Metabolism of terpenoids and polyketides-
1
LA
Y (ponovely comelatee)
Genetic Information Processing Folding, sorting and degradation -
9.5
Zero cross at 906
0.6
Zara cross at 906
Translation -
11
2
50
100
150
200
5.5
”’ (negatively conelated)
” (negatively complained)
Number of Gene
a
100
200
200
400
800
800
10 1,000
1,100
900
Rank in Ordered Dataset
0
100
200
300
400
500
000
700
1.000
1.100
Rank in Ordered Dataset
Enrichment profile
Hits
Ranking metric scores
E
KEGG enrichment barplot
Enrichment profile Hits
Ranking metric scores
Cytokine-cytokine receptor interaction-
14.83
Hematopoietic cell lineage-
13.28
Enrichment plot: KEGG_HEMATOPOIETIC_CELL_LINEAGE
Enrichment plot: KEGG_T_CELL_RECEPTOR_SIGNALING_PATHWAY
Chemokine signaling pathway-
10.78
Cell adhesion molecules (CAMs)-
10.31
0.0
Malaria-
10.14
Staphylococcus aureus infection-
Enrichment score (ES)
0.5
Enrichment score (ES)
…
9.03
Tuberculosis-
8.32
B.4
0.4
Rheumatoid arthritis -
8.14
0.3
0.3
Toll-like receptor signaling pathway-
8.04
Leishmaniasis-
7.45
02
9.2
T cell receptor signaling pathway-
7.42
1.1
0.1
Chagas disease (American trypanosomiasis) -
6.63
Pertussis-
6.63
0.0
Complement and coagulation cascades-
6.23
ntestinal immune network for IgA production-
6.23
Ranked list metric (Signal2Noise)
Ranked list metric (Signa(2Noise)
Osteoclast differentiation-
6.19
NF-kappa B signaling pathway-
5.93
(positively & nielike d)
Th1 and Th2 cell differentiation-
1.0
1.0
5.13
Leukocyte transendothelial migration-
4.85
0.5
05
Inflammationy bowel disease (IBD)-
4.85
D.D
Zero croun at Pod
A
Zero croun at 906
Q
$
10
15
-Ig(Qvalue)
“L’ (negatively conelated)
“L’(negatively comelatee)
100
200
000
400
500
600
700
800
000 1,000 1,100
Q
1DO
200
100
400
500
800
700
000
1,000 1,100
F
Top 20 of Pathway Enrichment
Rank in Ordered Dataset
Rank in Ordered Dataset
Enrichment profite - Hits
Ranking metric scores
Enrichment profile
Hits
Ranking metric scores
Leukocyte transendothelial migration
Inflammationy bowel disease (IBD)
Enrichment plot: KEGG_NATURAL_KILLER_CELL_MEDIATED_CYTOTOXI
Th1 and Th2 cell differentiation
·
NF-kappa B signaling pathway
CITY
Osteoclast differentiation
Complement and coagulation cascades
D.A
Intestinal immune network for IgA production Chagas disease (American trypanosomiasis)
Enrichment score (ES)
D.5
GeneNumber
DJ
21
0.2
Pathway
Pertussis
54
T cell receptor signaling pathway
0.1
Leishmaniasis
Qvalue
Toll-like receptor signaling pathway
0.000010
Rheumatoid arthritis
Ranked list metric (Signal2Noise)
0.000005
Tuberculosis
Staphylococcus aureus infection
1.0
Malaria
OUS
OUR
Zeno cross at sea
Cell adhesion molecules (CAMs)
Chemokine signaling pathway
“L’(vegaQuely comelate )
Hematopoletic cell lineage
·
100
200
200
400
600
000
700
800
000
000 1_1
Cytokine-cytokine receptor interaction
Rank in Ordered Dataset
Enrichment profile - Hits
Ranking metrie scores
0.20
0.25
9.30
0.55
Rich Factor
|
Cancer Medicine
Open Access
A
B
100
ROC curve (AUC=0.887)
1.0
Percent survival
U
0.8
True positive rate
0.6
50
0.4
L-RS
0.2
H-RS
P <. 0001
0
0.0
0.0
0.2
0.4
0.6
0.8
1.0
0
1000
2000
3000
4000
5000
False positive rate
OS (d)
cases (60.8%). Selected clinical characteristics were shown in Table 1. Patients’ immune scores (range from -1181.25 to 1577.33), stromal scores (range from -1766.65 to 1519.45), and estimate scores (range from -2947.9 to 3096.78) were listed in Table S1. We explored differences of immune scores according to different clinical characteristics and presented the results in Figure 1A-E. Lower immune scores were observed in disease with distant metastasis (DM) and locoregional recurrence (LR) than other cases (P = . 0204). Besides, difference analysis of stromal scores and estimate scores in data groups of clinical characteristics were shown in Figure S1.
Furthermore, K-M analysis revealed that lower immune scores was significantly associated with poor overall survival (OS) (Figure 1F, P = . 0495). However, difference in survival data was not observed in stromal scores and estimate scores (Figure 1G-H, P = . 3064 and P = . 2210, respectively).
3.2 Comparison of gene expression profiles with immune scores and stromal scores |
We used limma package to analyze RNA-seq data. ACC pa- tients were divided into low-level group (n = 39) and high-level group (n = 40) in accordance with median immune scores. Clustering analysis was shown in Heatmap in Figure 2A. Volcano plot which exhibits significantly differentially ex- pressed genes was shown in Figure 2B. Based on immune scores, 1649 differentially expressed genes (DEGs) which contained 1209 upregulated DEGs and 440 downregulated DEGs were screened (|log2fold-change| > 1, FDR < 0.05, Table S2). Besides, 1521 DEGs which contained 1180 up- regulated DEGs and 341 downregulated DEGs were screened based on stromal sores (log2fold-changel > 1, FDR < 0.05, Table S3). The heatmap and volcano plot based on stromal scores were shown in Figure S2 A-B. Venn plots were ap- plied to screen 1122 intersect genes which including 907 upregulated genes in higher immune/stromal scores and 215 downregulated genes in lower immune/stromal scores (Figure 2C,D).
3.3 PPI, GO analysis, KEGG analysis, and GSEA |
About 1122 intersect genes were uploaded and analyzed by online SRING tool (https://string-db.org/, minimum re- quired interaction score: 0.9) and Cytoscape software. Then, 18 hub genes (CD4, HLA-DRA, HCK, CD53, RAC2, HLA- DRB5, RAB37, CD93, FOLR2, GRAP2, HLA-DOA, HLA- DPA1, HPGDS, LAIR1, PTPRB, TACR1, TBXAS1, WAS) were selected from the result analyzed by cytoHubba. The network generated by the 12 topological analysis methods was shown in Figure 3. The GO analysis revealed that hub genes might be associated with T-cell receptor signaling pathway, antigen processing and presentation of exogenous peptide antigen via MHC class II, T-cell costimulation, MHC class II protein complex, and MHC class II receptor activity (Figure 4A-C). In addition, KEGG pathway anno- tation and enrichment were analyzed and shown in Figure 4D-E. Top 20 pathway enrichment was revealed by bubble chart in Figure 4F.
GSEA analysis was performed to further screen the sig- nificant pathway between higher immune scores group and lower immune scores group. We selected seven immune-re- lated pathways: hematopoietic cell lineage, T-cell receptor signaling pathway, natural killer cell-mediated cytotoxicity, leukocyte transendothelial migration, cytokine receptor in- teraction, chemokine signaling pathway, and JAK-STAT stat signaling pathway (Figure 4G).
3.4 RS calculation and survival analysis |
We calculated the risk score (RS) as: RS = - 0.398*HLA- DPA1-0.656*RAC2-0.247*HLA-DRB5- 0.246*HPGDS + 0.985*PTPRB + 0.692*HCK + 0.092*HLA- DOA + 0.005*RAB37 + 0.016*FOLR2-0.538*GRAP2-0.55- 4*TACR1-1.981*TBXAS1 + 0.102*HLA-DRA + 1.802* CD53 + 1.545*LAIR1-0.035*WAS-0.645*CD93-0.33*CD4 for all 79 ACC patients, and divided those patients into low-RS group (n = 39) and high-RS group (n = 40) (Table S4). Survival
|
A
Overall Survival CD4
B
Overall Survival CD53
C
Overall Survival CD93
1.0
+Low (N=19)
1.0
-+-Low (N=19)
1.0
+Low (N=19)
Hi gh (N=20)
+High (N=20)
+Hỉ gh (N=20)
0.8-
Logrank=0. 032
0.8-
Lbgrank-01 045
0.8-
Logrank=0.049
Percent Survival
Percent Survival
Percent Survival
0.6-
0.6-
0.6-
0.4-
0.4-
0.4-
0.2-
0.2-
0.2-
0.0
0
730
1.460
2.920
3.650
0.0
0.0
2. 190
0
730
1.460
2.190
2.920
3.650
0
730
1.460
2.190
2.920
3.650
Days
Days
Days
D
Overall Survival FOLR2
E
Overall Survival GRAP2
F
Overall Survival HCK
1.0
-+-Low (N=19)
1.0
++Low (N=19)
1.0
+Low (N=19)
+Hỉ gh (N=20)
+High (N=20)
+Hì gh (N=20)
0.8-
Logrank=0.03
0.8-
0.8-
Lbgrank=0.031
Percent Survival
Percent Survival
Logrank=0.01
Percent Survival
0.6
0.6
0.6-
0.4-
0.4-
0.4-
0.2-
0.2-
0.2-
0.0
730
1.460
2. 190
2.920
3.650
0.0
0
730
1.460
2.190
2.920
3.650
0.0
0
0
730
1.460
2.190
2.920
3.650
Days
Days
Days
G
Overall Survival HLA-DOA
H
Overall Survival HLA-DPA1
Overall Survival HLA-DRA
1.0
+Low (N=19)
1.0
+Low (N=19)
|+ Hitch(1-26)
+Low (N=19)
1.0
Logrank=0.006
+High (N=20)
+High(N=20)
0.8
Percent Survival
0.8-
Lpgrank=0, 02|
0.8-
Logrank=0, 019
Percent Survival
Percent Survival
0.6-
0.6-
0.6-
0.4-
0.4-
0.4-
0.2-
0.2-
0.2-
0.0
0
730
1.460
2.190
2.920
3.650
0.0
0.0
Days
0
730
1.460
2.190
2.920
3.650
0
730
1.460
2. 190
2.920
3.650
J
Days
Overall Survival HLA-DRB5
K
Days
Overall Survival HPGDS
L
Overall Survival LAIRI
1.0
+Low (N=19)
1.0
+Low (N=19)
1.0
+Low (N=19)
+High (N=20)
+High (N=20)
++High (N=20)
0.8-
Logrank=0, 031
0.8-
0.8-
Lbgrank=0.044
Percent Survival
Percent Survival
Lpgrank=0. 014
Percent Survival
0.6-
0.6-
0.6-
0.4-
0.4-
0.4-
0.2-
0.2-
0.2-
0.0
0
730
1.460
2.190
2.920
3.650
0.0-
0.0
0
730
1.460
2.190
2.920
3.650
0
730
1.460
2.190
2.920
3.650
Days
Days
Days
M
Overall Survival PTPRB
N
Overall Survival RAB37
0
Overall Survival RAC2
1.0
-+-Low (N=19)
1.0
++Low (N=19)
1.0
+Low (N=19)
Hi gh (N=20)
++High(N=20)
+High (N=20)
0.8-
Percent Survival
Logrank=0.044
0. 8-
Logrank=0.002
0.8-
Logrank=0.028
Percent Survival
Percent Survival
0.6-
0.6-
0.6-
0. 4-
0. 4-
0. 4-
0.2-
0.2-
0.2-
0.0
0.0
0
730
1.460
2. 190
2.920
3.650
0.0-
0
730
1.460
2.190
2.920
3.650
0
730
1.460
2.190
2.920
3.650
Days
Days
Days
P
Overall Survival TACRI
Q
Overall Survival TEXAS1
R
Overall Survival WAS
1.0
++Low (N=19)
1.0
+Low (N=19)
1.0
+Low (N=19)
+High (N=20)
+Hỉ gh(N=20)
+-High (-20)
0.8-
Logrank=0.003
Percent Survival
0.8-
Logrank=0.018
0.8-
Logrank=0.049
Percent Survival
Percent Survival
0.6-
0.6-
0.6
0.4-
0. 4-
0.4-
0.2
0.2-
0.2
0.0-
0
730
1.460
2.190
2.920
3.650
0.0
0.0
Days
0
730
1.460
2.190
2.920
3.650
0
730
1.460
2. 190
2.920
3.650
Days
Days
analysis revealed that high RS is significantly related with poor survival outcomes (Figure 5A, P < . 0001). To further evaluate the prognostic value of RS, we performed the ROC curve and
the area under the curve (AUC) was 0.887, which indicated su- perior predictive accuracy in survival outcomes (Figure 5B). The survival curves of 18 hub genes were exhibited in Figure 6A-R.
|
Cancer Medicine
Open Access
A
HLA-DOA Expression Level [og& RSEM)
CD4+ T Cel
K
Party
Neutrophil
HPGOS Expression Level [log2 RSEM) HLA-DRBS Expression Level jog2 RiSEM
partial. cor . 0.224
P=4.620-02
partial cor - 0.5
P- 4.476-06
p - 1.080-90
p . 2.470-04
P= 3.014-01
P = 2.888-05
P = 6.460-04
P= 3.678-03
B
04 04 08 1.0
0.30
Infiltration Level
0.12
9.12
B2
0.4
0.12
0.30
8:00 0.91 0.13
Infiltration Level
0,00
0.12
014
0.15
0.49 0.50 0.31 0.52
B
COL Expression Level (Jog RSEM)
CD4+ T Cel
L
Purity
CDI- T Cell
Derchtc Cel
CD-T Col
CD4+ T Cel
Cendinc Del
p = 7.700-01
partial. por 4 0 2re
partial cor ; @ 500
partial por - 0 504
partial.cor . 0.481 p= 1.049-05
partialcor . 0.224 p = 5.100-03
p = 2 010-05
p. 7.31041
- 9.816-41
P= 1.979-04
0 - 1.05%-02
p = 3.300-01
5
02
04
0.20
0.4
6.8
6.20
0.25
6.35-800
0.15
4.14
Infiltrason Level
4.15
012
4.14
Infiltration Level
C
CDS3 Expression Level (og RSEM)
Punty
CD4+ T Cel
Nouroghi
M
LAIRE Expression Level [log2 RSEM)
CD8+ T Cel
Macrophage
Dentroc Gol
partial cor- 4
partial Go C 0-140 p = 2 214-01
partial cor - 0 Se a = 2.30%-07
partial por - 0.497
19
AOC
8.12
6.95
8.12
Na
.
0.42
0.30
1.35-0 47
Infiltration Lavel
0.00
0.12
0.12
Infiltration Level
D
CD90 Expression Level (log2 RSEM)
Puty
Denstrie Gel
N
PTPAB Expression Level (og? RSEM
GOT - - 0.575
partial cor . 028
parsal cor & D.129
p . 4.226-06
P . 3.800-0
P= 3.070-0
0
&
AR
64
0.25
9.30
Infiltration Level
0.06
4.12
9.12
9.49
8
G
10
0,11
0.00
Infiltration Level
4.15
F
FOUR2 Expression Level (log2 RSEM)
O
Purty
CD4+ T Cell
RAB37 Expression Level (og2 RSEM)
Macrophage
Dontrinc Col
5 Col
CD4+ T Cel
Denchức Del
partialloor & d.com
0= 1.879-03
partisi cor - 0.410 p = 2.2:26-04
partial cor . 0.122
partial cor . 0.4747
p = 1.999-05
+
+
4-
:
NO
02
0.4
M
0.12
0.25
0.30
9.12
0.4
84
8. 12
0.20
Infiltration Lavel
Infiltration Level
4.14
G
GRA/2 Expression Level [og2 RSEM)
P
Purity
0 Cal
CD4-T Cel
Macrophago
Neutrochi
Denchito Del
RAC2 Expression Level [og@ RSEN)
Party
Cp44 T Gel
Deretroin Cell
p= 1.77 — 00
50
partial cor 2 0 box P = 4.379-01
4
TE
4.25
4.14
0.30
Infiltration Level
4.12
8.18
9.12
Intitration Level
H
006- T Cell
CD4-T Col
Noutrophi
Dendri: Cel
Q
TAGRI Expression Level foge RISEM
CD4+ T Cel
Macrop/uage
Conditc Cal
HICK Expression Level (og2 RSEM
partial cor - 0,10 P = 1.270-01
parial por + 0.224 p = 5.630-40
partial por - 0.574
partial.cor . 0.400
p = 2.274-03
p= 6.670-0
3
02
0.20
8.25
030
0.35 637
0.09 |1
Intitration Level
414
150
0.82
12
04
€6
6
RLO
4.12
Infiltration Level
0.15
4.14
HLA-DPA1 Expression Level (log2 RSEM)
R
Purity
TEXAS: Expression Level jog2 RSEM)
Dendau Cell
partial cor - D.1.45
P = 1.986.03
partial por - 0.061
partial.cor - 0.184
p= 1.195-01
p = 2.260-00
p = 4.000-07
p - 7.060-08
p . 3.290-00
9
2
13
20
0.35 tờ
Institradon Level
115
2.14
0.18
0.4
Innanarion Level
@12
J
HLA-ORA Expression Level Jog2 RSEI
006- T Cell
CD4+ T Cel
Nextochl
Dendite Gel
S
WAS Expression Level (log2 RSEM)
Dendiac Cell
D - 2.160-01
P= 1.416-64
P= 1,97%-03
partial Por - 0 51
P = 4.908-01
P= 1.00%-02
:
02
84
0.20
0.25
Infiltration Level
8.19
0.56
8.14
4.15
0.50
0.52
6.4
0.12
Infiltration Level
0.94
A
3.5 Selected hub genes were associated with immune cell infiltration |
To further understand the relationship between hub genes and immune infiltration in ACC microenvironment, we explored the correlation ship in TIMER database. The results illumi- nated that HLA-DOA was significantly related with B cell, CD4+ T cell, CD8 + T cell, macrophage, neutrophils, and den- dritic cell infiltration (Figure 7A). Other genes were related to some of the above immune cells infiltration (Figure 7B-S).
4 DISCUSSION |
Recently, TME is getting increasing attention because of the profound research on the mechanism of immunotherapy and target therapy. Positive responses to immunotherapy usually depend on the interaction between tumor cells and immune regulation in TME. Previous studies had investigated the role
of TME in various cancers. For example, Wood et al32 dis- cussed the role of autocrine/paracrine signaling interactions, ECM remodeling, and cell-cell interactions between tumor cells and the surrounding stroma, and the mechanisms could provide an important theoretical basis for target therapy in patients with NSCLC. Florent Petitprez et al33 demonstrated that the interactions between tumor cells and TME compo- nents were dependent on the original organ, treatment type, and oncogenic process. Unlike other studies our study identi- fied specific signatures that associated with the infiltration of immune and stromal cells in ACC TME based on the tran- scriptional profiles in high-quality datasets rather than focus- ing on the immune molecule and nontumor cells in TME.
In 2013, ESTIMATE was firstly introduced by Yoshihara et al17 to assess the immune cells infiltration and the presence of stromal cells based on gene expression data. This tech- nique will help elucidate the role of TME to neoplastic cell and provide a new perspective on the occurrence of genomic alterations. In our study, immune scores, estimate scores, and
stromal scores were calculated using ESTIMATE algorithm. Results showed that lower immune scores were significantly associated with DM and LR (P = . 0204). Further K-M analy- sis reflected that a lower immune score was closely related to poor OS (P = . 0495), which was consistent with the previous research performed by Jia et al (GBM patients with lower immune scores had better OS than those with higher immune scores (442 d vs 394 d, P = . 0537)).
PPI analysis was performed by Cytoscape software and SRING tool, from which 18 related hub genes were identi- fied including CD4, CD53, CD93, FOLR2, GRAP2, HCK, HLA-DOA, HLA-DPA1, HLA-DRA, HLA-DRB5, HPGDS, LAIR1, PTPRB, RAB37, RAC2, TACR1, TBXAS1, WAS. Furthermore, potential pathways were obtained after GO- enriched analysis such as T-cell receptor signaling pathway, antigen processing and presentation of exogenous peptide antigen via MHC class II, T-cell costimulation, MHC class II protein complex, and MHC class II receptor activity. Based on 18 selected hub genes which were closely related to ACC TME, a predictive model was developed. ACC pa- tients with low RS had better survival outcomes than those with high RS (P < . 0001). Additionally, satisfactory predic- tive efficiency (AUC = 0.887) revealed that our predictive model may have its potential for clinical application.
Protein tyrosine phosphatase receptor type B (PTPRB), also known as VE-PTP, was regarded as an independent prognostic factor in patients with NSCLC. It was downreg- ulated and was significantly associated with OS in NSCLC patients. Besides, PTPRB regulated the tumorigenesis and Scr phosphorylation.34 Weng et al35 reported that PTPRB promoted tumor metastasis and invasion of colorectal cancer via inducing epithelial mesenchymal transition. Obviously, it would be a new therapeutic target for treat- ment. Increasing evidence revealed that dysregulated Rab proteins were related to cancer progression.36-38 RAB37 was a metastasis suppressor by regulating exocytosis to the extracellular compartment, resulting in inhibition of cancer metastasis and neo-angiogenesis. 39,40
In the analysis of the relationship between hub genes and im- mune infiltration, HLA-DOA was identified to have significant correlation with B cell, CD4 + T cell, CD8 + T cell, macrophage, neutrophils, and dendritic cell infiltration. Previous researches had reported that infiltration of macrophage, CD4+ T cell, and CD8 + T cell played an important role in tumor prognosis. HLA- DOA is a nonclassical class II MHC molecule. It exhibits very little sequence variation when compared with other classical HLA class II molecules. Studies on the role of HLA-DOA in tumor prognosis are rare. Yukinori Okada41 reported that HLA- DOA was an independent risk of anticitrullinated protein auto- antibody-positive rheumatoid arthritis.
Finally, some limitations should be taken into consider- ation in our research. First of all, lack of clinical data to val- idate the conclusions. We performed this study based on the
public database via biological algorithm approaches and fur- ther clinical trials are needed before its application. Moreover, further investigations of these selected 18 TME-related hub genes should be conducted, to have a better understanding of the regulatory mechanism in immune infiltrates, which may bring novel insights into the potential association of TME with ACC prognosis in a comprehensive manner.
5 CONCLUSIONS |
In our study, we screened a list of TME-related genes which predict poor survival outcomes in ACC patients from TCGA database. HLA-DOA was significantly related with immune infiltration in ACC microenvironment. Further studies are needed to verify our conclusion.
ACKNOWLEDGMENTS
This work was supported by grants from the Six Talents Peaks Program of Jiangsu Province (No. 2016-WSW- 021), National Natural Science Foundation of China (No. 81702520), Medical Research Project of Jiangsu Provincial Health and Family Planning Commission (No. H2018052), Research Project of Jiangsu Cancer Hospital (Nos. ZM201801, ZN201602), and the Young Talents Program of Jiangsu Cancer Hospital (No. 2017YQL-04).
CONFLICTS OF INTEREST
The authors have no conflicts of interest to declare.
ORCID Xiao Li ID https://orcid.org/0000-0002-9358-1136 Feng Qi ID https://orcid.org/0000-0003-1061-7058
REFERENCES
1. Mansmann G, Lau J, Balk E, Rothberg M, Miyachi Y, Bornstein SR. The clinically inapparent adrenal mass: update in diagnosis and management. Endocr Rev. 2004;25:309-340.
2. Else T, Kim AC, Sabolch A, et al. Adrenocortical carcinoma. Endocr Rev. 2014;35:282-326.
3. Fiorentini C, Grisanti S, Cosentini D, et al. Molecular drivers of potential immunotherapy failure in adrenocortical carcinoma. J Oncol. 2019;2019:6072863.
4. Kerkhofs TMA, Verhoeven RHA, Van der Zwan JM, et al. Adrenocortical carcinoma: a population-based study on inci- dence and survival in the Netherlands since 1993. Eur J Cancer. 2013;49:2579-2586.
5. Varghese J, Habra MA. Update on adrenocortical carcinoma man- agement and future directions. Curr Opin Endocrinol Diabetes Obes. 2017;24:208-214.
6. Ayala-Ramirez M, Jasim S, Feng L, et al. Adrenocortical carci- noma: clinical outcomes and prognosis of 330 patients at a tertiary care center. Eur J Endocrinol. 2013;169:891-899.
7. Miller BS, Ammori JB, Gauger PG, Broome JT, Hammer GD, Doherty GM. Laparoscopic resection is inappropriate in patients
with known or suspected adrenocortical carcinoma. World J Surg. 2010;34:1380-1385.
8. Leboulleux S, Deandreis D, Al Ghuzlan A, et al. Adrenocortical carcinoma: is the surgical approach a risk factor of peritoneal car- cinomatosis? Eur J Endocrinol. 2010;162:1147-1153.
9. Fassnacht M, Johanssen S, Quinkler M, et al. Limited prognos- tic value of the 2004 International Union Against Cancer staging classification for adrenocortical carcinoma: proposal for a Revised TNM Classification. Cancer. 2009;115:243-250.
10. Wang Q, Hu B, Hu X, et al. Tumor evolution of glioma-intrinsic gene expression subtypes associates with immunological changes in the microenvironment. Cancer Cell. 2017;32(42-56):e6.
11. Suvà M, Rheinbay E, Gillespie S, et al. Reconstructing and repro- gramming the tumor-propagating potential of glioblastoma stem- like cells. Cell. 2014;157:580-594.
12. Grivennikov SI, Greten FR, Karin M. Immunity, inflammation, and cancer. Cell. 2010;140:883-899.
13. Ren BO, Cui M, Yang G, et al. Tumor microenvironment partici- pates in metastasis of pancreatic cancer. Mol Cancer. 2018;17:108.
14. Huang J, Shen F, Huang H, et al. Th1high in tumor microenviron- ment is an indicator of poor prognosis for patients with NSCLC. Oncotarget. 2017;8:13116-13125.
15. Kitawaki Y, Nakamura Y, Kubota-Nakayama F, et al. Tumor mi- croenvironment in functional adrenocortical adenomas: immune cell infiltration in cortisol-producing adrenocortical adenoma. Hum Pathol. 2018;77:88-97.
16. Tomczak K, Czerwinska P, Wiznerowicz M. The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge. Contemp Oncol (Pozn). 2015;19:A68-77.
17. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tu- mour purity and stromal and immune cell admixture from expres- sion data. Nat Commun. 2013;4:2612.
18. Carter SL, Cibulskis K, Helman E, et al. Absolute quantification of somatic DNA alterations in human cancer. Nat Biotechnol. 2012;30:413-421.
19. Shah N, Wang P, Wongvipat J, et al. Regulation of the glucocorti- coid receptor via a BET-dependent enhancer drives antiandrogen resistance in prostate cancer. Elife. 2017;6:e27861. https://doi. org/10.7554/eLife.27861
20. Priedigkeit N, Watters RJ, Lucas PC, et al. Exome-capture RNA sequencing of decade-old breast cancers and matched decalcified bone metastases. JCI. Insight. 2017;2(17):e95703. https://doi. org/10.1172/jci.insight.95703
21. Alonso MH, Aussó S, Lopez-Doriga A, et al. Comprehensive analysis of copy number aberrations in microsatellite sta- ble colon cancer in view of stromal component. Br J Cancer. 2017;117:421-431.
22. Ritchie ME, Phipson B, Wu DI, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47.
23. Aalen OO. A linear regression model for the analysis of life times. Stat Med. 1989;8:907-925.
24. Aalen OO. Further results on the non-parametric linear regression model in survival analysis. Stat Med. 1993;12:1569-1588.
25. Chen H, Boutros PC. VennDiagram: a package for the genera- tion of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics. 2011;12:35.
26. Szklarczyk D, Franceschini A, Wyder S, et al. STRING v10: pro- tein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015;43:D447-D452.
27. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software en- vironment for integrated models of biomolecular interaction net- works. Genome Res. 2003;13:2498-2504.
28. Chin C-H, Chen S-H, Wu H-H, Ho C-W, Ko M-T, Lin C-Y. cyto- Hubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8(Suppl 4):S11.
29. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and ge- nomes. Nucleic Acids Res. 2000;28:27-30.
30. Tilford CA, Siemers NO. Gene set enrichment analysis. Methods Mol Biol. 2009;563:99-121.
31. Heagerty PJ, Lumley T, Pepe MS. Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics. 2000;56:337-344.
32. Wood SL, Pernemalm M, Crosbie PA, Whetton AD. The role of the tumor-microenvironment in lung cancer-metastasis and its relationship to potential therapeutic targets. Cancer Treat Rev. 2014;40:558-566.
33. Petitprez F, Vano YA, Becht E, et al. Transcriptomic analysis of the tumor microenvironment to guide prognosis and immunotherapies. Cancer Immunol Immunother. 2018;67:981-988.
34. Qi Y, Dai Y, Gui S. Protein tyrosine phosphatase PTPRB regulates Src phosphorylation and tumour progression in NSCLC. Clin Exp Pharmacol Physiol. 2016;43:1004-1012.
35. Weng X, Chen W, Hu W, et al. PTPRB promotes metastasis of col- orectal carcinoma via inducing epithelial-mesenchymal transition. Cell Death Dis. 2019;10:352.
36. Yang J, Liu W, Lu X, Fu Y, Li L, Luo Y. High expression of small GTPase Rab3D promotes cancer progression and metastasis. Oncotarget. 2015;6:11125-11138.
37. Pellinen T, Arjonen A, Vuoriluoto K, et al. Small GTPase Rab21 regulates cell adhesion and controls endosomal traffic of beta1-in- tegrins. J Cell Biol. 2006;173:767-780.
38. Hendrix AN, Maynard D, Pauwels P, et al. Effect of the secretory small GTPase Rab27B on breast cancer growth, invasion, and me- tastasis. J Natl Cancer Inst. 2010;102:866-880.
39. Tsai C-H, Cheng H-C, Wang Y-S, et al. Small GTPase Rab37 tar- gets tissue inhibitor of metalloproteinase 1 for exocytosis and thus suppresses tumour metastasis. Nat Commun. 2014;5:4804.
40. Tzeng H-T, Tsai C-H, Yen Y-T, et al. Dysregulation of Rab37- mediated cross-talk between cancer cells and endothelial cells via thrombospondin-1 promotes tumor neovasculature and metastasis. Clin Cancer Res. 2017;23:2335-2345.
41. Okada Y, Suzuki A, Ikari K, et al. Contribution of a non-classical HLA gene, HLA-DOA, to the risk of rheumatoid arthritis. Am J Hum Genet. 2016;99:366-374.
SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section.
How to cite this article: Li X, Gao Y, Xu Z, Zhang Z, Zheng Y, Qi F. Identification of prognostic genes in adrenocortical carcinoma microenvironment based on bioinformatic methods. Cancer Med. 2020;9:1161- 1172. https://doi.org/10.1002/cam4.2774