Fatty acid metabolism-related signature predicts survival in patients with clear cell renal carcinoma

Data source

Two clear cell renal carcinoma patients’ datasets with survival results were used in this study: TCGA-KIRC dataset, including RNA sequences and clinical phenotypes, downloaded from the UCSC Xena [8] and GSE29609 dataset, including expression matrix and survival results, downloaded from the GEO DataSets [9]. These two datasets were publicly available from the databases. We used TCGA-KIRC dataset as a training set to perform a series of survival analyses and identify a survival-related signature. In contrast, GSE29609 was designed as a test set for testing the signature’s accuracy and stability.

Survival-related fatty acid genes

Based on a review of previous literature, 90 fatty acid-related genes were selected [10]. To find the survival value of each of the 90 fatty acid genes, we performed Lasso cox regression analysis on the TCGA-KIRC dataset. Subsequently, univariate and multivariate Cox regressions were performed on genes with survival value to find survival-related fatty acid genes.

Signature of survival-related fatty acid genes

We predicted a new signature containing survival-related fatty acid genes. The following formula was used to calculate each sample risk score:

$$\text{Risk Score}={\sum }^{\text{i}}(\text{Coefi}\times \text{ExpGenei})$$

“Coef” represents non-zero regression coefficients calculated with univariate Cox regression analysis. “ExpGene” represents expression values of genes from the prognostic risk score model. We calculated risk score level for TCGA-KIRC dataset’ patients and verified the survival value of the signature in both the training set and test set using KM analysis.

Relationship between risk score and clinical features

First, using test set, we performed a time-dependent ROC analysis to predict risk score on 1-, 3- and 5-year survival results. Second, important clinical prognostic features, such as age and TMN stage, were compared in the training set using Cox regression analysis. R package “survival”, and “glmnt” were used for ROC and Cox regression analyses. We constructed a forest plot showing these results. In addition, prognostic characteristics with significantly different risk scores selected for further analysis. A nomogram of overall survival for clear cell renal carcinoma patients in the training set with risk scores and other important clinical prognostic indexes was drawn.

Enrichment analysis for low- and high-risk score patient groups

Patients included in the training set were divided into two groups: high risk-score group and low risk-score group. R package ‘limma’ was used to identify differentially expressed genes between two patient groups. Genes with FDR <0.05 were considered statistically significant. Three functional enrichment analyses were performed on differentially expressed genes for the two groups: GO analysis [11], KEGG analysis [12], and Gene set variation analysis (GSVA) [13]. R packages “cluster Profile” and “GSVA” were used for enrichment analysis. Based on enrichment analysis results, we identified candidate different functions between the high-score group and low-score group patients.

Immune-related characteristics in low- and high-risk score groups

We performed multi-dimension immune-related analysis on the two groups of patients. First, we used cibersort website (https://cibersort.stanford.edu/) to analyze the immune component of the two patients groups and divided patients into immune subtypes [14]. Second, we used R package “GSVA” to analyze immune-related score. Third, we explored immune status of the two groups using R package “estimate”. Fourth, the TAM score and Cytotoxic score, important prognostic indexes generated using R package “GSVA”, were compared between the two patient groups. Finally, we compared the expression levels of two important immune checkpoint genes: CTLA4 and PDL1.

Chemotherapy drug sensitivity analysis

Ridge regression model was constructed with R package “pRRophetic” and “oncoPredict” using the remaining data. Drug sensitivity (IC50) was chosen as an outcome index. Chemotherapeutic response was predicted using tumor genes expression and drug sensitivity data of cell lines from the cancer genome project (CGP) and drug sensitivity in cancer (GDSC) [15].

Verification of survival-related fatty acid genes in cell lines

Increasing the expression level of fatty acid synthase (FASN) was positively correlated with aggressive cell proliferation, migration, apoptosis, and lipid droplet formation. It was also shown to regulate metabolic disorders associated with clear cell renal carcinoma. Researchers also proved that pharmacological inhibitor of FANS could suppress the growth and invasiveness of clear cell renal carcinoma. Therefore, we used two human clear cell renal carcinoma cell line (caki-1 and 786-0 cell lines) from authenticated cell cultures of the Chinese national collection. We down-regulated the expression of FANS in caik-1 cells by the way of siRNA virus and exposed FASN inhibitorC75 (HY-12364, MedChemExpress, China) with two cell line (caki-1 and 786-0). We verified the expression level of survival-related fatty acid genes with quantitative real-time polymerase chain reaction (Q-PCR) [16].

Data availability statement

The original data presented in this study are included in the article. All data were retrieved from public databases.

Identity of survival-related fatty acid genes

We downloaded 90 fatty acid genes, 525 patient clinical records and RNA sequences from the TCGA-KIRC dataset (the training set) and survival results of 39 renal cancer patients and the matrix of genes’ expression from the GSE29609 (the test set). We performed Lasso regression analysis on the training set for all fatty acid genes (Figure 1A, 1B). The results showed that only 14 genes had survival value for renal cancer patients. After univariate and multivariate Cox regression analyses (Figure 1D), only five survival-related genes (AGR2, HAO2, IGF2BP1, MCCD1 and OLFM4) were significantly different to have the statistical value (p < 0.05). The correlation among the five genes was also established in the Figure 1C.

Figure 1. (A, B) The lasso regression analysis; (C) The correlation of 5 key genes; (D) The univariate and multivariate COX regression analysis; (E) The risk score distribution in the training set; (F) The KM analysis in the TCGA cohort; (G) The KM analysis in the GSE cohort.

Risk score construction and survival value verification

The risk score was calculated using the following formula:

Risk score = (0.0726303) × AGR2 + (−0.17735540) × HAO2 + (0.18400723) × IGF2BP1 + (−0.08972776) × MCCD1 + (−0.08952658) × OLFM4.

The patients in the training set were divided into low- and high-risk score groups based on their risk score level. We also presented the distribution of risk score in the training set (Figure 1E). To explore the survival value of each risk score, we tested it on both training and test sets using KM analysis. The results showed that renal cancer patients with high risk score had shorter survival time for both the training set (p < 0.001, Figure 1F) and test set (p = 0.008, Figure 1G). We also used time-dependent ROC curve to analyze the risk score in the test set. The AUC value of risk score was 0.813, with satisfactory results for ROC curves (Figure 2A, 2B).

Figure 2. (A, B) The ROC time dependent analysis; (C) The forest plot for clinical features in the training set; (D–I) The comparative results for age, gender, T stage, N stage, M stage, and pathology stage.

Relationship between risk score and clinical features

We compared clinical features (like TMN stages) in the training set using Cox regression analysis. The results showed that six clinical features (age, T stage, N stage, M stage, gender and pathology stage) had significant survival value (p < 0.05, Figure 2C) and were presented in a forest plot. A further separate comparison of these six indexes between the two groups showed that patients with high-risk score had worse T stage, N stage, M stage, and pathology stage (P < 0.001) than low-risk score patients (Figure 2D–2I). A nomogram containing several clinical features and risk score was drawn predicting 1-, 3- and 5-year overall survival (Figure 3A–3D).

Figure 3. (A) The nomogram for all the clinical features; (B–D) The nomogram to predict probability 1-/3-/5- year OS.

Enrichment analysis for low- and high-risk score patient groups

Using R package “limma”, we identified differentially expressed genes, which were represent in a volcano plot (Figure 4A). Three forms of functional enrichment analysis were used to determine potential function of risk score in two groups of renal cancer patients. First, GO analysis showed that functions of extracellular matrix organization, extracellular structure organization, and urogenital system development would be activated in high-risk score patients (Figure 4C). In contrast, the functions of negative thymic T cell selection and actin cytoskeleton reorganization were inhibited in high-risk score patients. Second, KEGG analysis showed that pathways of PIK3, papillomavirus infection, and MAPK signaling would be activated in high-risk score patients (Figure 4D). Third, GSVA analysis ultimately found that a high-risk score was strongly correlated with inflammatory response and epithelial mesenchymal response in patients (Figure 4B).

Figure 4. (A) The volcano plot for the different expression level genes; (B) The GSVA analysis; (C) The GO analysis; (D) The KEGG analysis.

Immune-related characteristics in low- and high-risk score groups and chemotherapy drug sensitivity analysis

Using cibersort website, we identified immune components of two patient groups. However, differences between main immune cells could not be identified by scanning whole immune landscapes of the two groups (Figure 5A). However, high-risk score patient group were more likely to exhibit immune insensitivity subtypes, such as immune C2 and C5, than low-risk score patients. Subsequently, we used R package “estimate” to analyze the estimate score and immune score (Figure 5B, 5C). The results showed that immune scores were significantly different between the two patient groups (p = 0.043). Scanning the immune scores with R package “GSVA” revealed a strong correlation of TAM score and cytotoxic score were with risk score. TAM score and cytotoxic score have been reported to have a prognostic value for many cancers. Here, we analyzed the two scores for renal cancer patients. The results showed a strong positive survival value of TAM score and cytotoxic score in renal cancer patients (Figure 5D). KM analysis further revealed that patients with high-risk score and TAM score or cytotoxic score had poorer survival results than low-risk score patients (p < 0.001, Figure 5E, 5F). We further compared the correlation of risk score with two important immune therapy genes: PDL1 and CTLA4. The results showed a strong positive correlation between risk score and expression level of PDL1 or CTLA4 (Figure 6B, 6C). Finally, we performed chemotherapy drug sensitivity analysis on CGP and GDSC platforms. The landscape of drug sensitivity analysis was presented in Figure 6A.

Figure 5. (A) The landscape of immunity response by the cibersort and GSVA; (B, C) The immune score and estimate score by the R package “estimate”; (D) The correlation between the TAM score and cytotoxic score; (E, F) The KM analysis of TAM score and cytotoxic score.

Figure 6. (A) The landscape of drug response by the CGP and GDSC; (B, C) The correlation between the risk score and the expression level of CTLA4 and PDL1; (D, E) The results of 5 key genes by the test of Q-PCR among three groups: control group; fatty acid inhibitor group, and the FASN siRNA group, caki-1 cell line (D); 786-0 cell line (E).

Verification of survival-related fatty acid genes in caki-1 cell line

We verified survival related fatty acid genes in caki-1 cell and 786-0 cell (two human renal clear cell carcinoma cell lines) in two ways: down-regulating FASN expression and exposure to pharmacological inhibitor of FASN. We measured expression levels of five survival-related fatty acid genes. The results showed that three genes (AGR2, IGF2BP1 and OLFM4) had been inhibited whereas two genes (HAO2 and MCCD1) had been activated in both two ways (Figure 6D, 6E). These changes suggested that the five genes regulated fatty acid synthase function in renal cancer patients. Therefore, the five-gene signature could be used to predict survival outcomes of renal cancer patients.