Abstract
Background
Clear cell renal cell carcinoma (ccRCC) is a common urinary cancer.
Although diagnostic and therapeutic approaches for ccRCC have been
improved, the survival outcomes of patients with advanced ccRCC remain
unsatisfactory. Fatty acid metabolism (FAM) has been increasingly
recognized as a critical modulator of cancer development. However, the
significance of the FAM in ccRCC remains unclear. Herein, we explored
the function of a FAM-related risk score in the stratification and
prediction of treatment responses in patients with ccRCC.
Methods
First, we applied an unsupervised clustering method to categorize
patients from The Cancer Genome Atlas (TCGA) and International Cancer
Genome Consortium (ICGC) datasets into subtypes and retrieved
FAM-related genes from the MSigDB database. We discern differentially
expressed genes (DEGs) among different subtypes. Then, we applied
univariate Cox regression analysis followed by least absolute shrinkage
and selection operator (LASSO) linear regression based on DEGs
expression to establish a FAM-related risk score for ccRCC.
Results
We stratified the three ccRCC subtypes based on FAM-related genes with
distinct overall survival (OS), clinical features, immune infiltration
patterns, and treatment sensitivities. We screened nine genes from the
FAM-related DEGs in the three subtypes to establish a risk prediction
model for ccRCC. Nine FAM-related genes were differentially expressed
in the ccRCC cell line ACHN compared to the normal kidney cell line
HK2. High-risk patients had worse OS, higher genomic heterogeneity, a
more complex tumor microenvironment (TME), and elevated expression of
immune checkpoints. This phenomenon was validated in the ICGC cohort.
Conclusion
We constructed a FAM-related risk score that predicts the prognosis and
therapeutic response of ccRCC. The close association between FAM and
ccRCC progression lays a foundation for further exploring FAM-related
functions in ccRCC.
Keywords: Fatty acid metabolism, ccRCC, Immune features, Prognostic
model, Treatment
1. Introduction
Clear cell renal cell carcinoma (ccRCC) accounts for approximately
70–80% of all renal malignancies and is one of the most prevalent
urinary tumors [[35][1], [36][2], [37][3], [38][4]]. More than 30% of
ccRCC patients have disease metastases at diagnosis, and the 5-year
survival rate of metastatic ccRCC is only 10–20% [[39][5], [40][6],
[41][7]]. Although targeted therapies focused on the vascular
epithelial growth factor (VEGF) pathway, such as sunitinib, axitinib,
and pazopanib, have improved the survival of patients with ccRCC,
acquired resistance and side effects remain a challenge during ccRCC
treatment [[42][8], [43][9], [44][10], [45][11]]. Recent studies have
shown that immune checkpoint inhibitor therapies, including
pembrolizumab, ipilimumab, and nivolumab, significantly prolong the
overall survival (OS) and disease-free survival (DFS) in patients with
localized ccRCC [[46][12], [47][13], [48][14]]. However, high
recurrence rates after immunotherapy and the lack of reliable
biomarkers to predict therapeutic effects still hinder the development
of ccRCC treatments [[49]12]. Therefore, building a powerful model to
precisely predict ccRCC characteristics, patient clinical stage, and
drug response is essential for formulating personalized treatment
regimens and improving overall patient survival.
Although cancer types and underlying etiologies greatly differ,
dysregulation of cellular metabolism is a common feature of
tumorigenesis [[50][15], [51][16], [52][17]]. Aerobic glycolysis is a
classic example of a metabolic perturbation, altering glucose
metabolism to aerobic glycolysis, called the Warburg effect [[53][18],
[54][19], [55][20]]. Abnormal glucose metabolism satisfies the
metabolic necessities for cell growth and supports high levels of
glycolytic intermediates for various biosynthetic pathways in cancer
progression [[56]21,[57]22]. Additionally, aerobic glycolysis
establishes a suppressive immune microenvironment that triggers
immunosuppression in the tumor microenvironment (TME) [[58][23],
[59][24], [60][25]]. Reprogramming fatty acid metabolism (FAM) is
another characteristic of cancer. Multiple critical enzymes involved in
FAM are dysregulated in several cancers [[61][26], [62][27], [63][28],
[64][29]]. Many studies have demonstrated that abnormally activated
fatty acid (FA) oxidation and de novo FA synthesis provide metabolic
energy, cell membranes, and signaling molecules for carcinogenesis
[[65][30], [66][31], [67][32], [68][33], [69][34]]. Most understanding
of FAM indicates that limiting the FA supply may be a therapeutic
strategy for cancer [[70][35], [71][36], [72][37]].
FAs in the TME influence the role and infiltration of immune cells
[[73][38], [74][39], [75][40], [76][41], [77][42]]. Metabolic
reprogramming is intricately linked to antitumor immunity [[78][43],
[79][44], [80][45], [81][46]]. ccRCC cells are characterized by
abundant glycogen and lipid accumulation [[82]47,[83]48] and exhibit
diverse remodeling of cellular metabolism, including glucose, FA, and
tricarboxylic acid cycle [[84][49], [85][50], [86][51], [87][52]].
Tumor-specific metabolic changes in ccRCC, including FAM, are
associated with poor patient outcomes [[88]53,[89]54]. Moreover,
enhanced lipid import and droplet formation benefit ccRCC tumorigenesis
[[90]55]. Multiple FAM enzymes, including anabolic enzymes such as
fatty acid synthase (FASN), and catabolic enzymes, such as carnitine
palmitoyltransferase I (CPT1A), hydroxyacyl-CoA dehydrogenase
trifunctional multienzyme complex subunit alpha (HADHA) and beta
(HADHB), and acetyl-CoA acetyltransferase 1 (ACAT1), have been
considered potential prognostic biomarkers for ccRCC [[91][56],
[92][57], [93][58]]. The von Hippel–Lindau (VHL) inactivation, commonly
present in ccRCC, stabilized the hypoxia-inducible factor (HIF)1α to
further repress FAM by inhibiting CPT1A and inducing lipid deposition
[[94]59]. Mechanistic investigations showed that the metabolic
vulnerabilities of key enzymes involved in FA metabolism are novel
therapeutic targets for ccRCC [[95]50]. Therefore, an in-depth
exploration of FAM may uncover the mechanisms underlying ccRCC
progression and new therapeutic strategies.
In this study, we explored the potential roles of key FAM-related genes
in ccRCC and provided a theoretical foundation for future ccRCC
research. We developed and validated a novel FAM-related risk score
using The Cancer Genome Atlas (TCGA) and International Cancer Genome
Consortium (ICGC) datasets. The FAM-related risk score predicted
clinical features, TME characteristics, and therapeutic efficacy in
patients with ccRCC. In summary, our results provide new insights into
FAM-related mechanisms in ccRCC. The established FAM-related risk score
is reliable for predicting OS and treatment response in patients with
ccRCC.
2. Materials and methods
2.1. Data acquisition and preprocessing
We retrieved RNA sequencing (RNA-Seq) profiles and clinical information
from TCGA (https://tcga‐data.nci.nih.gov/tcga/) and ICGC databases.
TCGA-KIRC cohort with 513 samples was used as the training set, and the
ICGC (RECA-EU) dataset with 91 patients was used as the validation set.
Samples with incomplete expression profiles, clinical follow-up data,
survival time, or statuses were excluded. OS is the period from the day
of diagnosis to death from any cause or the last follow-up and
expressed in years. We converted the ensemble gene symbols to gene
symbols for subsequent analyses and used the mean as the gene
expression value when several probes were matched to the same gene. We
retrieved 309 FAM-related genes from three FAM single-sample gene set,
including the hallmark, the Kyoto Encyclopedia of Genes and Genomes
(KEGG), and Reactome, from the molecular signature database (MSigDB)
([96]http://www.broadinstitute.org/gsea/msigdb/).
2.2. Enrichment analysis and correlations between FAM and clinicopathological
features
To explore whether FAM was involved in ccRCC pathogenesis, we performed
sample gene set enrichment analysis (ssGSEA) using the “GSVA” and
“clusterProfiler” R packages to compute FAM-related enrichment scores
and enrichment levels of FAM pathways, respectively, between tumor and
para-cancerous tissues in TCGA cohort. Furthermore, we analyzed the
association between FAM and clinicopathological features of patients
with ccRCC.
2.3. Molecular subtypes based on FAM-related genes
We conducted the consensus clustering of TCGA cohort based on FAM
levels using the “ConsensusClusterPlus” R package [[97]60]. We varied
the number of clusters k from 2 to 9 to determine the optimal k value
for stratifying the appropriate subtypes based on the delta area,
cumulative distribution function (CDF), and consensus matrix. We
applied the Spearman's rank correlation coefficient as the distance
measurement and completed 500 repetitions to ensure cluster stability.
We then performed principal component analysis (PCA) to demonstrate the
typing effect. We constructed Kaplan–Meier (K–M) survival curves to
assess the OS between molecular subtypes using the “survival” R package
[[98]61]. Finally, we evaluated the distribution of clinicopathological
features and differences in FAM between the different subtypes.
2.4. Gene mutational landscape of FAM-related subtypes
Tumor mutation characteristics are closely related to tumor progression
and immunotherapeutic efficacy. Therefore, we investigated mutational
characteristics, including genetic mutations and the mutant-allele
tumor heterogeneity (MATH) score of FAM-related subtypes using the
“maftools” R package. MATH is a tumor-specific score used to measure
intratumoral genetic heterogeneity based on variations in the variant
allele frequency of all mutations in the tumor. TCGA mutation data were
retrieved from Genomic Data Commons (GDC)
([99]https://portal.gdc.cancer.gov/). Furthermore, we estimated the
association of our subtypes with six previously published immune
subtypes (IFN-γ dominant, inflammatory, lymphocyte-depleted,
immunologically quiet, wound healing, and TGF-β dominant) using a
pan-cancer analysis [[100]62].
2.5. Functional enrichment analysis of ccRCC FAM-related subtypes
To identify the underlying characteristics of FAM-related subtypes in
ccRCC development, we performed GSEA using KEGG gene signatures to
explore the enriched pathways among different FAM-related subtypes. We
then conducted the enrichment score of FAM genes of the hallmarks gene
signatures to quantify the enrichment levels of different pathways
among FAM-related subtypes.
2.6. Immune features and pathway enrichment of FAM-related subtypes
To investigate the immune microenvironment of FAM-related subtypes in
ccRCC, we used the “ESTIMATE” R package to evaluate the immune status
of patients with ccRCC [[101]63]. The “CIBERSORT” R package was used to
measure the infiltration of immune cells [[102]64].
2.7. Therapeutic differences among FAM-related subtypes
Based on the immune characteristics among different FAM-related
subtypes, we explored the response of different molecular subtypes to
conventional ccRCC chemotherapy agents by the IC[50] value measurement
using “pRRophetic” R package, including cisplatin, bortezomib,
temsirolimus, axitinib, gefitinib, sunitinib, sorafenib, and
vinblastine. As immune checkpoint inhibitors (ICIs) have revolutionized
cancer immunotherapy, we used a heatmap to compare the distribution of
the published 45 immune coinhibitory/costimulatory genes among
FAM-related subtypes [[103]65].
2.8. Differentially expressed genes (DEGs) and relevant function analysis of
FAM-related subtypes
To classify key FAM genes, we conducted pairwise comparisons of DEGs
among FAM-related subtypes using the “Limma” R package and exhibited
the results in heatmaps [[104]66]. Significant DEGs were determined by
false discovery rate (FDR) < 0.05, and |log[2] [Fold Change (FC)]
| > 1. Subsequently, KEGG and Gene Ontology (GO) analyses were
performed using the “webgestalt” R package to evaluate pathway
differences for the key FAM genes identified [[105]67].
2.9. Establishment of a prognostic FAM-related risk score
We applied univariate Cox regression analysis to preliminarily screen
crucial prognostic genes from the DEGs identified preliminarily
[[106]68]. Then, we conducted LASSO analysis using the “glmnet” R
package to select the crucial genes in TCGA cohort [[107]69]. We used
10-time cross-validation to estimate the confidence interval for each
lambda and subsequently identified the optimal lambda with the lowest
average error. In addition, we conducted multivariate Cox regression in
a backward stepwise fashion to curtail the gene number and acquire the
regression coefficient for every gene. We assumed that p-selected genes
are input into a prognostic prediction model, denoted as (x[1],
…,x[p]). The risk score was a weighted sum of genes whose weights
reflect the degree of association and defined as
[MATH: Riskscore=β1×x1+…+βp×xp
:MATH]
[[108]70]. Then, patients were divided into high- and low-risk groups
based on the cutoff value determined by the “survminer” R package
([109]http://www.sthda.com/english/rpkgs/survminer). Survival curves
were constructed using the K–M to compare OS between the different risk
groups. The ICGC cohort was used for the external validation to verify
the robustness of the established prediction model. Finally, we
examined the relationships between risk scores and clinical
characteristics to elucidate their potential for clinical application.
We also evaluated the differences in risk scores among subtypes.
2.10. Immune status between different risk groups
We performed GSEA to compare the differences in enrichment scores of
FAM pathways between groups using genes selected from MSigDB, hallmark,
KEGG, and Reactome. We then applied the ESTIMATE R package to explore
the TME features between groups. We used CIBERSORT to analyze the
relative abundance of tumor-infiltrating immune cells between the two
groups. Thereafter, we collected 29 functional gene expression
signatures (FGES) previously published associated with the TME in
different risk groups [[110]71]. Additionally, we obtained previously
published gene signatures from 15 pathways characterized by immune,
stromal, DNA repair, and oncogenic signatures to investigate the
heterogeneity of the immune microenvironment between different risk
groups [[111]72]. We also assessed the correlation of the risk score
with the FAM pathways, 29 gene expression signatures, and gene
signatures from 15 pathways using Pearson’s correlation.
2.11. The predictive performance of the risk score on ccRCC therapy
First, we evaluated drug sensitivity using IC[50] for common
chemotherapy agents in TCGA risk groups. We then analyzed whether
immunotherapy sensitivity varied between risk groups. We analyzed the
differences in MATH scores using the Wilcoxon test to investigate
genomic heterogeneity between groups. Given the close relationship
between tumor immune regulation and immune checkpoint-related genes
(ICGs), we explored the relationship between the risk score and ICGs
and the expression profiles of ICGs between groups.
2.12. Expression validation of the selected genes in FAM-related prognostic
model
To further explore the FAM-related genes selected from the constructed
prognostic model at the mRNA level in ccRCC and normal renal cells, we
performed cell culture and real-time quantitative PCR (RT–qPCR)
[[112]73]. Human renal tubular epithelial HK2 and ccRCC ACHN cells were
cultured in Dulbecco’s modified eagle medium (DMEM, Gibco, Invitrogen,
Carlsbad, CA, USA). Total RNA was extracted using TRIzol Reagent
(Invitrogen, Carlsbad, CA, USA), following the standard protocol. cDNA
was synthesized using a PrimeScript™ RT reagent kit (Takara, Shiga,
Japan). RT–qPCR was conducted using an ABI7300 real-time PCR system
(Applied Biosystems) to validate the expression of the selected
FAM-related genes from the constructed prognostic model, following the
manufacturer’s instructions. The amplification conditions for qRT-PCR
were 95 °C for 30 s, 40 cycles of 95 °C for 5 s, and 60 °C for 30 s.
The primers used are listed in [113]Supplementary Table 1. Each mRNA
expression was normalized to GAPDH expression according to the 2^−ΔΔCt
method, and all reactions were performed thrice in triplicate to
calculate the average cycle threshold (Ct).
2.13. Statistical analysis
Statistical analyses were conducted using R (version 3.4.3) and
GraphPad Prism (version 8.4.2) software. To compare clinical features
between different groups, we applied the Wilcoxon test to evaluate two
independent non-parametric samples and the Kruskal–Wallis test to
compare three or more samples. Continuous variables with a normal
distribution are presented as means ± standard deviation. Pearson’s or
Spearman’s rank correlation coefficients were used to explore the
correlations between variables. Multivariate Cox regression analysis
was performed to explore the independence of the constructed
FAM-related model for prognostic performance. p-value < 0.05 indicated
statistical significance [[114]74].
3. Results
3.1. Influence of FAM in ccRCC
To elucidate the mechanisms of FAM in ccRCC progression, we found that
the enrichment scores of the genes in FAM pathways in TCGA cohort were
higher in para-cancer tissues than in ccRCC tissues (p < 0.05)
([115]Fig. S1A). We comprehensively assessed the association between
FAM and clinicopathological features. As expected, higher clinical
stages resulted in lower FAM scores, and women had more active FAM than
men (p < 0.05) ([116]Fig. S1B).
3.2. Identification of ccRCC subtypes associated with FAMGs
To explore the carcinogenic mechanisms of FAMGs in ccRCC, we performed
a consensus cluster analysis based on the expression of FAMGs in TCGA
cohort. Based on the relative CDF and delta area changes,
[MATH: k=3 :MATH]
(the number of clusters) achieved the flattest middle portion of the
CDF curve and a definitive boundary for the consensus matrix (Figs.
S2A–C). Therefore, we classified patients in TCGA cohort into subtypes:
subtype 1, subtype 2, and subtype 3, named S1, S2, and S3,
respectively. Furthermore, PCA analysis also proved the typing results
with distinct differences among the three subtypes (p < 0.05)
([117]Fig. 1A). K–M survival curves showed that ccRCC patients in S1
had shorter OS than those in S2 and S3 in both TCGA and ICGC cohorts
(p < 0.05) ([118]Fig. 1B). Additionally, we found clear differences in
the enrichment levels of FAM among the S1, S2, and S3 subtypes, with
the highest FAM scores in S3 and the lowest in S1 (p < 0.05) ([119]Fig.
1C). Next, we evaluated the differences in clinicopathological features
between the S1, S2, and S3 subtypes in TCGA cohort and found that S1
patients with poorer prognoses had higher TNM and clinical stages and
grades (p < 0.05) ([120]Fig. S3).
Fig. 1.
[121]Fig. 1
[122]Open in a new tab
Consensus clustering of FAMGs of ccRCC in TCGA cohort. (A) PCA plot of
the three FAM-related subtypes. (B) Kaplan–Meier survival curves for OS
of the S1, S2, and S3 subtypes. (C) Comparison of FAM of the S1, S2,
and S3 subtypes. Abbreviations: ccRCC, clear cell renal cell carcinoma;
CDF, cumulative distribution function; FAMGs, fatty acid metabolism
genes; KEGG, the Kyoto Encyclopedia of Genes and Genomes; OS, overall
survival; PCA, principal component analysis; ssGSEA, sample gene set
enrichment analysis; TCGA, The Cancer Genome Atlas.
3.3. Characteristics of gene mutations among FAMG-related subtypes
First, we found that VHL, PBRM1, TTN, SETD2, and BAP1 mutations were
commonly found in TCGA cohort and were associated with FAMG-related
subtypes ([123]Fig. S4A). To understand the influence of FAM-related
gene mutations in ccRCC, we observed that the S1, S2, and S3 subtypes
exhibited substantial differences in tumor heterogeneity, and S3 had
the lowest MATH scores (p < 0.05) ([124]Fig. S4B). Finally, we found
that S3 patients were mostly matched to the widely recognized immune
subtype S3 (p < 0.05) ([125]Fig. S4C).
3.4. Enriched pathways among FAMG-related subtypes
Furthermore, the GSEA results identified distinct enriched pathways in
the three ccRCC subtypes. For example, S1 was largely enriched in
pathways correlated to inflammation and immune response, including the
“complement and coagulation cascades,” “primary immunodeficiency,”
“systemic lupus erythematosus,” and “P53 signaling pathway” ([126]Fig.
2A). Besides “FAM pathways,” other metabolism-related pathways such as
“butanoate,” “propanoate,” and “pyruvate” were enriched in S2
([127]Fig. 2B). S3 was also enriched in several metabolism-related
pathways, including “FAM,” “pentose and glucuronate interconversions,”
and “ascorbate and aldarate metabolism” ([128]Fig. 2C). S1 was enriched
in several pathways, including immune-, cell cycle-, and cancer-related
pathways. S3 was enriched for some metabolism-related pathways such as
“FAM,” “oxidative phosphorylation,” “xenobiotic metabolism,” and “bile
acid metabolism” ([129]Fig. 2D).
Fig. 2.
[130]Fig. 2
[131]Open in a new tab
Pathway enrichment analysis among three subtypes in TCGA cohort. (A)
S1, (B) S2, and (C) S3. (D) Heatmap of differential pathway enrichment
among the three subtypes. Abbreviations: KEGG, the Kyoto Encyclopedia
of Genes and Genomes; TCGA, The Cancer Genome Atlas.
3.5. Immune features of FAMG-related subtypes
To evaluate the TME components among different FAMG-related subtypes,
we used ESTIMATE to assess the immune cell status. We found that S1 had
higher stromal, immune, and ESTIMATE scores than the other two subtypes
and a lower tumor purity, suggesting that S1 had a high immune
infiltration density and heterogeneity (p < 0.05) ([132]Fig. 3A). We
further used CIBERSORT to determine the abundance of 22 immune cells
among the FAMG-related subtypes. S1 showed an increased proportion of
M0 macrophages, regulatory T, and Treg cells (p < 0.05) ([133]Fig. 3B).
Fig. 3.
[134]Fig. 3
[135]Open in a new tab
Associations between the immune signatures of TME and FAM-related
subtypes in TCGA cohort. (A) Comparison of ESTIMATE scores in the S1,
S2, and S3 subtypes. (B) Comparison of the immune cell distribution in
the S1, S2, and S3 subtypes. Abbreviations: CD4, cluster of
differentiation 4; CD8, cluster of differentiation 8; FAM, fatty acid
metabolism; NK, natural killer; TCGA, The Cancer Genome Atlas; TME,
tumor microenvironment.
3.6. Differences in cancer therapy response among FAMG-related subtypes
Next, we examined whether the FAMG-related subtypes significantly
influenced cancer chemotherapy and immunotherapy. We analyzed the
response of the subtypes to conventional chemotherapy agents and found
that S1 was more sensitive to bortezomib and sunitinib, whereas S3 was
more sensitive to axitinib in TCGA cohort (p < 0.05) ([136]Fig. 4A).
Given the promising applications of immunotherapy with checkpoint
inhibitors, we also tested the levels of several representative immune
checkpoints in the three FAMG-related subtypes. The results revealed
that most immune checkpoint genes, characterized by “immune activation”
and “immune inhibition, were highly expressed in the S1 stage
([137]Fig. 4B). Furthermore, we identified 15 significantly
differentially expressed ICGs among FAMG-related subtypes: KIR3DL1,
ADORA2A, CD80, TNFSF18, TNFRSF18, CD28, TNFSF14, LAG3, CEACAM1, LAIR1,
IDO1, PVR, CD200, CD276, and HAVCR1 (p < 0.05) ([138]Fig. 4C).
Fig. 4.
[139]Fig. 4
[140]Open in a new tab
Estimation of chemotherapy and immunotherapy response of the S1, S2,
and S3 subtypes in TCGA cohort. (A) Estimated IC[50] for commonly used
chemotherapy drugs in the S1, S2, and S3 subtypes. (B) Expression of
ICGs in the S1, S2, and S3 subtypes. (C) Expression of 15
differentially expressed ICGs in the S1, S2, and S3 subtypes.
Abbreviations: ICGs, immune checkpoint-related genes; IC50, half
maximal inhibitory concentration; TCGA, The Cancer Genome Atlas.
3.7. Differential gene screening among FAMG-related subtypes
After establishing three distinct FAMG subtypes in ccRCC, we compared
their gene expression levels to determine the co-expressed differential
genes. We found 321 DEGs between the S1 and S2 subtypes, 146
upregulated and 175 downregulated genes in S1, 973 DEGs between the S1
and S3 subtypes, 305 upregulated and 668 downregulated genes in S1, 151
DEGs between the S2 and S3 subtypes, and 29 upregulated and 122
downregulated genes in S2 (p < 0.05) (Figs. S5A–B). We found that 86
DEGs overlapped among the three FAMG-related subtypes: SLC27A2,
CYP4A11, FMO1, CYP4A22, HAO2, CYP2J2, HMGCC2, and CA4 ([141]Fig. S5C).
Moreover, we performed a functional enrichment analysis of 305 DEGs
between S1 and S3 and found that these DEGs were mainly mapped to
stroma-related signaling pathways ([142]Figs. S5D–G).
3.8. Construction of a FAMG-related risk prediction model
To evaluate the relationship between FAMGs and the prognosis of ccRCC
patients, we used univariate Cox regression analysis and verified 852
DEGs with prognostic values; 214 were risk genes and 638 were
protective genes (p < 0.05) ([143]Fig. 5A). Next, we incorporated the
852 prognostic genes to establish a prognostic risk prediction model
using LASSO regression. First, regression coefficients approached zero
and gradually increased with an increase in lambda. The optimal lambda
was set to 0.0531 based on a 10-time validation, and 22 genes were
finally included in the risk score ([144]Fig. 5B). Subsequently,
multivariate stepwise linear regression analysis verified the final
nine prognostic genes and determined the coefficient for each targeted
gene ([145]Fig. 5C). The following formula was established using the
mRNA level of each risk gene and coefficients:
[MATH: Riskscore=(−0.151×ANK3)+(0.158×GYG2)+(0.157×SYCE
1L)+(0.378×SPAG
5)+(−0.108×PLG)+(−0.314×SLC2
A9)+(0.165×RNAS
E2)+(−0.154×CGN)+(−0.147×SEMA
3G) :MATH]
. Moreover, we validated the expression of nine genes in the HK2 and
ACHN cell lines. Compared to HK2 cells, the ccRCC cell line ACHN had
higher SEMA3G and SYCE1L mRNA levels and lower CGN, GYG2, RNASE2, and
SLC2A9 mRNA levels (p < 0.05) ([146]Fig. 5D). With the cutoff value of
−0.199807 set by the “survminer” R package, we allocated patients to
the high- and low-risk groups in TCGA cohort. We constructed K–M
survival curves and found that high-risk patients had significantly
poorer prognoses (p < 0.05) ([147]Fig. 5E). To verify the robustness of
risk scores, the same model was applied to the ICGC cohort for external
validation. High-risk patients showed shorter OS than low-risk ones in
the ICGC cohort (p < 0.05) ([148]Fig. 5F). These results demonstrated
that our risk score was predictive of OS in different cohorts,
suggesting robust performance of the constructed nine-gene risk score.
We also investigated the distribution of risk scores among clinical
characteristics (i.e., cancer grade, stage, and TNM) and found that the
more advanced the ccRCC clinical stage, the greater the risk score
(p < 0.05). In addition, the distribution of risk scores exhibited
significant differences among the three FAMG-related subtypes
(p < 0.05) ([149]Fig. S6).
Fig. 5.
[150]Fig. 5
[151]Open in a new tab
Construction and validation of FAM-related risk score for ccRCC in TCGA
and ICGC cohorts. (A) Univariate Cox regression analysis of DEGs
identifying OS-related FAMGs. (B) Partial likelihood deviance of LASSO
coefficient profiles. (C) Multivariate Cox regression analysis for nine
selected FAMGs of the risk score. (D) The mRNA expression levels of the
nine selected FAMGs by RT–qPCR. ROC curve of the FAM-related risk score
in the ICGC cohort. (E) Kaplan–Meier curves of the FAM-related risk
score in TCGA cohort. (F) Kaplan–Meier curves of the FAM-related risk
score in the ICGC cohort. Abbreviations: ccRCC, clear cell renal cell
carcinoma; DEGs, differentially expressed genes; FAM, fatty acid
metabolism; FAMGs, fatty acid metabolism genes; ICGC, International
Cancer Genome Consortium; LASSO, least absolute shrinkage and selection
operator; ROC, receiver operating characteristic; RT–qPCR, real-time
quantitative PCR; TCGA, The Cancer Genome Atlas.
3.9. Associations between the risk score and immune characteristics
We first evaluated the enrichment of FAM pathways between low- and
high-risk groups using the Hallmark, KEGG, and Reactome databases.
Low-risk patients had higher FAM enrichment scores than high-risk
patients (p < 0.05) ([152]Fig. 6A). To further illustrate the
differences in the TME between the different risk groups, we used
ESTIMATE to obtain tumor stroma and immune cell infiltration. The
immune, stromal, and ESTIMATE scores were higher in the high-risk group
than in the low-risk group. In contrast, the tumor purity was higher in
the low-risk group than in the high-risk one (p < 0.05) ([153]Fig. 6B).
We also compared the levels of 22 immune cell types between the risk
groups and found a remarkable difference in the proportion of 22 immune
cell types between the groups (p < 0.05) ([154]Fig. 6C). Moreover, the
enrichment scores between the high- and low-risk groups differed
significantly for 29 FGES. The FGES associated with tumor stroma and
promotion accounted for most of the enrichment in the high-risk group
([155]Fig. 6D). Compared to the low-risk group, several pathways
correlated with DNA damage repair were differentially enriched in the
high-risk group based on the 15 pathway-related enrichment score
distributions ([156]Fig. 6E). Additionally, we explored the
relationship between the risk score and the FAM pathways, 29 FGES, and
15 pathways. The risk score exhibited a negative correlation with FAM,
a positive correlation with FGES associated with tumor stroma and
promotion, and pathways correlated with DNA damage repair ([157]Fig.
6F).
Fig. 6.
[158]Fig. 6
[159]Open in a new tab
Relationship between immune infiltration features and the risk score in
TCGA cohort. (A) FAM enrichment scores in high- and low-risk groups.
(B) Immune, stromal, estimate, and tumor purity scores in high- and
low-risk groups. (C) Distribution of 22 immune cells in high- and
low-risk groups. (D) TME-related 29 gene signatures in high- and
low-risk groups. (E) Scores of 15 pathways in high- and low-risk
groups. (F) Correlation of risk score with TME-related gene signatures
and 15 pathways. Abbreviations: ssGSEA, sample gene set enrichment
analysis; TCGA, The Cancer Genome Atlas; TME, tumor microenvironment.
3.10. Comparisons of therapeutic responses between groups
For drug sensitivity analysis, we compared the IC[50] of several
commonly used chemotherapeutic agents between the low- and high-risk
groups, including bortezomib, axitinib, cisplatin, sorafenib,
temsirolimus, sunitinib, gefitinib, and vinblastine. The high-risk
group exhibited greater sensitivity to bortezomib, gefitinib, and
temsirolimus than the low-risk group (p < 0.05) ([160]Fig. 7A). Further
correlation analysis revealed that the risk score was significantly
correlated with the expression of several ICGs, including CTLA4, IDO1,
and ICOS (p < 0.05) ([161]Fig. 7B). Many ICGs were significantly
overexpressed in the high-risk group, including LAG3, IDO1, and CTLA4
(p < 0.05) ([162]Fig. 7C).
Fig. 7.
[163]Fig. 7
[164]Open in a new tab
Risk score predictions for ccRCC treatment sensitivity in TCGA cohort.
(A) Violin plots of the estimated IC[50] for chemotherapy agents in
high- and low-risk groups. (B) Association of risk score with ICGs. (C)
Heatmap of differentially expressed ICGs in high- and low-risk groups.
ccRCC, clear cell renal cell carcinoma; ICGs, International Cancer
Genome Consortium; TCGA, The Cancer Genome Atlas.
4. Discussion
Metabolic and epigenetic studies in the last decade have revealed a new
relationship between metabolic alterations and cancer progression
[[165][75], [166][76], [167][77], [168][78]]. Metabolic modifications
trigger oncogenic transformations and are hallmarks of cancer. For
example, upregulated glycolytic metabolism contributes to energy
production and essential metabolic intermediates in most solid tumors
[[169][79], [170][80], [171][81]]. The metabolic properties of tumors
induce immunosuppression in the TME and constrain cancer immunotherapy
[[172][82], [173][83], [174][84]]. Although most studies have focused
on the influence of FAM on cell physiological functions and biological
behaviors, the effects of FAM on ccRCC remain unclear [[175][85],
[176][86], [177][87]]. FAM is dysregulated in ccRCC and correlated with
a poor prognosis of patients. Various key enzymes related to FAM have
been found to participate in the malignance of ccRCC. However, the
pathogenic mechanism of FAM in ccRCC has not been elucidated [[178]53].
In addition, the prognostic performance of the previously published
FAM-related risk scores and gene signatures requires further clinical
verification [[179]88,[180]89]. The immune features and therapeutic
responses of FAM-related genes have also been inadequately illustrated
in the constructed models. Therefore, we constructed a novel
FAM-related risk score to explore the regulatory mechanisms of FAM in
ccRCC progression and identify high-risk patients with poor prognoses
to guide personalized treatment strategies.
FAM reprogramming is a ubiquitous metabolic process in cancer
[[181][90], [182][91], [183][92]]. Several key enzymes of the FAM play
a critical role in cancer cell survival and progression and are
currently attracting considerable attention as new therapeutic targets
[[184]93,[185]94]. Emerging omics techniques, including genomics,
epigenomics, transcriptomics, proteomics, and metabolomics, have shown
that ccRCC characteristically exhibits high lipid accumulation and
abnormal FAM [[186]53,[187]55,[188]95,[189]96]. In this study, we found
that the FAM scores differed between the primary tumor and
para-tumorous normal tissues. We also found a correlation between the
FAM and the clinicopathological features of patients with ccRCC. A
higher clinical stage was associated with lower FAM scores.
The three molecular subtypes S1, S2, and S3, categorized by FAM
features using consensus cluster analysis, exhibited significant
differences in prognosis, clinical grade, tumor genomic heterogeneity,
biological pathways, TME, and response to common treatments. Thus, to
better understand the functions of FAM-related genes in ccRCC, we
constructed a predictive risk score to discern different degrees of
risk in patients with ccRCC. We showed that the FAM-related risk score
predicted OS efficiently in both TCGA and ICGC cohorts. Consistent with
previous studies, high-risk patients had worse OS than low-risk
patients and were susceptible to posterior clinical staging, suggesting
that the FAM-related risk score can be used to identify high-risk
patients and provide personalized treatment to improve outcomes
[[190][97], [191][98], [192][99], [193][100]].
Previous studies demonstrated that ccRCC is highly infiltrated by
different immune cells and has high immune heterogeneity
[[194]98,[195][101], [196][102], [197][103]]. We found that high-risk
patients had more plentiful immune and stromal components, lower tumor
purity in the TME, and more abundant immunosuppressive cells. As
expected, the high-risk group had elevated enrichment scores for
cancer-associated fibroblasts (CAFs), tumor- and matrix-associated
macrophages, and tumor proliferation rates. Additionally, pathways
involved in oncogenesis, such as the DNA damage repair pathway, were
identified in high-risk groups.
Chemotherapy has remained the most effective treatment for advanced
ccRCC for decades [[198]104,[199]105]. Herein, we found that high-risk
patients were more sensitive to bortezomib, gefitinib, and temsirolimus
than other treatments, indicating that they might benefit from
chemotherapy. The significant genomic heterogeneity differences between
the groups reflect features of tumor promotion and treatment resistance
in the high-risk group, highlighting the importance of immunotherapy in
ccRCC [[200]106]. Therefore, we explored immunotherapy responses
between the low- and high-risk groups and found that the
characteristics reflecting the immunotherapy response were increased in
the high-risk group compared to the low-risk group. Thus, high-risk
patients may show an improved immunotherapy response. In summary, these
results demonstrate that the constructed FAM-related risk score
exhibits a promising prognostic performance for patients with ccRCC and
provides insights into the function of FAMGs in ccRCC.
This study has some limitations. First, we used multiple online
databases and bioinformatics methods without prospective clinical
validation. Thus, we will conduct an additional validation cohort study
at our center to confirm the predictive values of the constructed
nine-gene risk score for ccRCC. Moreover, we validated the mRNA levels
of the nine genes in the risk prediction model. Therefore, future
studies should investigate the underlying mechanisms of the nine FAMGs
in ccRCC progression.
In conclusion, we successfully developed and verified a robust risk
score associated with the FAM in ccRCC that can be applied to predict
the prognosis of patients with ccRCC. The FAM-related risk score
described the clinicopathological features of the patients and
predicted their susceptibility to chemotherapy and immunotherapy. The
FAM-related risk score has contributed to research on promising ccRCC
prognostic predictors and future individualized treatments for ccRCC.
Author contribution statement
Conceived and designed the experiments: Qinfan Yao, Dajin Chen, and
Jianghua Chen.
Performed the experiments: Xiuyuan Zhang, Chunchun Wei, Qiannan Xu, and
Hongjun Chen.
Analyzed and interpreted the data: Qinfan Yao, Chunchun Wei, and
Xiuyuan Zhang.
Contributed reagents, materials, analysis tools or data: Qiannan Xu,
and Jianghua Chen.
Wrote the paper: Qinfan Yao, Dajin Chen, and Xiuyuan Zhang.
Funding statement
This work was funded by the National Natural Science Foundation of
China (81802085), and the Natural Science Foundation of Zhejiang
(LY21H160033).
Data availability statement
The datasets generated in the current study are available from the
corresponding author.
Declaration of competing interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to
influence the work reported in this paper
Footnotes
^Appendix A
Supplementary data to this article can be found online at
[201]https://doi.org/10.1016/j.heliyon.2023.e17224.
Contributor Information
Jianghua Chen, Email: chenjianghua@zju.edu.cn.
Dajin Chen, Email: zju2001@zju.edu.cn.
Appendix A. Supplementary data
The following are the Supplementary data to this article.
Multimedia component 1
[202]mmc1.docx^ (13.8KB, docx)
Multimedia component 2
[203]mmc2.pdf^ (1,012KB, pdf)
References