Identification of multi-omics landscape and prognostic LMAGs in RPLS
The whole study process, including consensus clustering, immune landscape visualization, DEGs analysis, nomogram construction and WGCNA analysis were systematically evaluated and depicted in the workflow chart, as shown in Table S1.
To systematic appraise the lipid metabolism dysregulation in RPLS, 741 LMAGs were obtained from the MSigDB database. Initially, we encapsulated the incidence of copy number variations and somatic mutations in cohort-TCGA, where 14862 amplified genes were screened to identify potent antigens (Fig. 1A), while missense mutation was the most common variant classification (Fig. 1B). Additionally, we investigated the incidence of CNV among these LMAGs. By evaluating the frequency of CNV status, existing CNV alterations of LMAGs were also clarified, and the top 35 genes in amplified or deleted CNV status were summarized (Fig. S1A).
Subsequently, according to the transcriptomic data, consensus clustering was performed to cluster RPLS patients into two lipid metabolism subgroups (LMSs) (K = 2) (Fig. 1C) based on 135 overall survival (OS) related genes by univariable Cox analysis (Fig. S1B). 31 and 27 patients were clustered into LMS1 and LMS2, respectively (Fig. S1C). Heatmap visualization also indicated that prognostic LMAGs profiles differed significantly between LMS1 and LMS2, while LMS2 was enriched with LMAGs (Fig. 1D), but predicted poor prognosis (Fig. 1E). Interestingly, two subtypes harbored heterogenous somatic mutations profiles, demonstrated ATRX and B4GALNT1 as the genes with the highest mutation frequency in LMS1, but MUC16 in LMS2 (Fig. S1D, E). However, not significant difference was found in genome altered fraction and mutation counts between LMS1 and LMS2 (Fig. S1F, G).
Moreover, we also investigated 26 disease free survival (DFS) related LMAGs through univariate Cox analysis (Fig. S1H). Intersecting the results of OS-related and DFS-related genes, 8 overlapping LMAGs (GSTM4, GRHL1, PI4K2B, GK3P, ARSK, ELOVL2, NR1H4, and ABHD4) were excavated and eligible for further screen prognostic relevant antigens (Fig. S1I). The location on chromosomes and expression levels were visualized in Fig. 1F. The regulatory network described the comprehensive landscapes of the 8 LMAGs indicating their interactions, correlation feature (Fig. 1G) and prognostic values (Fig. 1H). These findings indicated that the LMAGs classified RPLS patients into two subtypes with different molecular features and prognosis.
Heterogeneous Functional Enrichment And Immune Landscape In Lmags Subtypes
To better understand the innate difference of survival and underlying signaling mechanisms between LMS1 and LMS2, DEGs and functional enrichment analyses were performed, respectively. A total of 4144 DEGs were identified, of which 3586 genes were downregulated and 558 genes were upregulated in LMS2, as compared with LMS1 (Fig. 2A). GO enrichment analysis indicated that these up-regulated DEGs were involved in positive regulation of immune system process, immune response and regulation of immune system process (Fig. 2B). Similarly, KEGG enrichment analysis also validated these pathways associated with cytokine-cytokine receptor interaction, chemokine signaling pathway and B cell receptor signaling pathway, of which are part of the immune response (Fig. 2C). Meanwhile, PPI analysis further confirmed 15 submodels, all of which were closely associated with immune and metabolism (Fig. S2A). These results depict that the expression of LMAGs is closely related to immune-related biological processes, confirming that lipid metabolic reprogramming is significantly associated with tumor immune microenvironment (TIME) in RPLS.
To further evaluate the dysregulation of immune and metabolism remodeling involved in RPLS, a series of TIME profiles was conducted. Firstly, we illustrated the distribution feature of previously reported six pan-cancer immune subtypes (C1-C6)[18], of which LMS1 and LMS2 were mostly clustered into C1 (Wound Healing), C2 (IFN-γ Dominant), C3 (Inflammatory), C4 (Lymphocyte Depleted) and C6 (TGF-β Dominant) (Fig. S2B). Intriguingly, C3 and C6 presented a higher proportion in LMS1, while they were associated with better outcomes. Accordingly, immune score, stromal score and ESTIMATE score were also calculated using the ESTIMATE algorithm (Fig. 2D). Remarkably, LMS1 demonstrated significantly higher TIME scores and better prognosis than those in LMS2, implying that favorable immune components and immune-related molecules were abundant in LMS1.
Next, the relationship between LMSs and 34 infiltrated immune cells was further explored. In concordance with previous TIME scores, LMS1 harbored more infiltrated DC, B cells, CD4 + Tem, macrophages, monocytes and NKT than those in LMS2 (Fig. 2E), indicating the significant impacts of these immune cells in the progression of RPLS. However, the complete immune response involves a close combination of multiple events, not only the infiltrated immune cells[19]. Thereafter, we further calculated and compared the immune activity score of each step through TIP analysis. Similarly, the abundance of anti-tumor immune cells was significant higher in LMS1 than those in LMS2, as well as in Step1, Step3 and Step5 (Fig. S2C). In addition, we explored substantial differences in TLS-associated 12-chemokine signature between LMS1 and LMS2 (Fig. 2F), while the expression of lymphoid-structures-associated B-cell-specific chemokine CXCL13 was notably higher in LMS1 (Fig. 2G).
Accumulating evidence indicated that tumor with high TMB level predicted better efficacy of immunotherapy[20]. We then estimated the value of TMB in both LMS1 and LMS2, but not significant difference was found (Fig. S2D). Intriguingly, patients in LMS1 with low TMB demonstrated a satisfactory survival benefit. Considering the importance of immune checkpoint inhibitors in the treatment of solid tumor[21], we further examined the differences in immune checkpoint profiles and found notably substantial differences in CD28, CD40, CD86, HAVCR2 and PD-1, between these two subtypes (Fig. S2E). Immunogenic cell death (ICD) has been classified as a form of regulated cell death (RCD) that is sufficient to activate an adaptive immune response[22]. We next identified ICD-related genes and analyzed the expression patterns. Importantly, we discovered that significant higher expression of FPR1, TLR4 and CXCL10 were enriched in LMS1 (Fig. S2F). Taken together, these findings demonstrated the unique characteristics of TIME within two LMSs, offering a conducive complement to previous studies.
Identification Of Immune Gene Co-expression Modules And Immune Hub Genes Of Rpls
The immune gene co-expression module was designed and applied to classify immune-related genes, whose expression may significantly influenced the effectiveness of potential mRNA vaccine[23]. Therefore, we re-analyzed and enrolled immune-related genes to establish gene modules through WGCNA (Fig. 3A). The soft threshold was set at four in the scale-free network (Fig. S3A). The representation matrix was converted to adjacency and next to a topological matrix. The average-linkage hierarchy clustering approach was applied with a minimum of 30 genes for each network according to the standard of a hybrid dynamic shear tree. Eigengenes of each module were calculated and the close modules were integrated (height = 0.25, deep split = 3 and min module size = 30). Notably, six gene modules were identified (Fig. 3B, C), and correlation feature was also visualized (Fig. S3B). In addition, the module eigengenes in LMS1 were significantly higher in yellow and blue modules (Fig. 3D). Moreover, the prognostic analysis indicated that eigengenes in the brown module was significantly associated with OS in RPLS (Fig. 3E).
Further functional enrichment analysis illustrated that genes involved in cytokine-cytokine receptor interaction and T cell receptor signaling were enriched in the brown module (Fig. 3F). The brown module (361 immune-related genes) was further selected, and three of them (PLCG1, ZC3HAV1L and NFAT5) were filtered to build the risk score through LASSO algorithm (Fig. S3C, D). Patients were classified into the high-risk and low-risk groups, while high-risk group predicted unfavorable OS (P = 0.02) (Fig. 3G). The area under the receiver operating curve (AUC) was 0.75, indicating a good accuracy of the model (Fig. S3E). Taken together, this risk model may serve as a novel tool for prognostic predicting in RPLS, based on the immunogenic genes co-expression network.
Development Of Survival And Relapse Risk Models And Nomograms Based On Lmags
Given the significant biological roles of LMAGs in lipid metabolic reprogramming, the association between LMAGs-related risk score and the prognosis needed thoroughly study. Thus, two prognostic models were conducted for OS and DFS, respectively.
24 LMAGs were found to be considerably linked to the OS of patients through LASSO regression analysis (Fig. 4A, B), 13 of which were tested and selected for the risk score model from multivariate Cox analysis (Fig. S4A). We also investigated the relationship between risk score and survival status, and the low-risk subgroup harbored significant more alive statuses (Fig. 4C) and better OS than those in high-risk subgroup (Fig. 4D). Specifically, this OS related model indicated a great accuracy with AUC values of 0.94 in 1 year, 0.97 in 3 years and 0.97 in 5 years (Fig. 4E).
To evaluate and validate the universality of this LMAGs-related prognostic model from cohort-TCGA, an independent dataset (cohort-FD) was performed as a validation cohort. The risk score of each case in cohort-FD was calculated using the same formula as that for the cohort-TCGA. Similarly, patients in high-risk group suffered unfavorable OS than those in low-risk group (Fig. S4B). In addition, the AUC values of this model according to ROC analysis were 0.7 in 1 year, 0.8 in 3 years, and 0.85 in 5 years (Fig. S4C).
Since this model predicted great potency in clinical prognosis based on LMAGs, we further investigated the correlation between risk score and TIME in cohort-TCGA. In concordance with previous findings, the high-risk subgroup was characterized with a significantly lower immune score (Fig. 4F), and infiltrated less CD8 + T cells and plasma cells, but more resting memory CD4 + T cells and Tregs than those in low-risk subgroup (Fig. S4D), implying that high-risk subgroup was considered with the characteristics of immunosuppressive status.
Meanwhile, we also explored the value of model between the two subgroups stratified by different clinical features. Univariate Cox analysis indicated that patients with dismal OS were authenticated with larger tumor size, chemotherapeutic efficacy and high risk score (Table 1). Multivariate Cox analysis further confirmed that all of them were independent risk factors (Table 1). Subsequently, we developed a nomogram for OS prediction using these two clinical parameters and the LMAGs-based risk scores. A calibration plot for internal validation of the nomogram presented excellent consistency between the nomogram-predicted probability and actual observations of the 1-, 3-, and 5-year OS (Fig. 4G).
Table 1
Univariate and multivariate analysis of OS in cohort-TCGA (n = 58).
Variables | Univariate analysis | Multivariate analysis |
HR (95%CI) | P value | HR (95%CI) | P value |
Age, years (> 60 vs. ≤60) | 3.778 (0.888–16.075) | 0.072 | | |
Sex (female vs. male) | 0.531 (0.240–1.176) | 0.119 | | |
Tumor size, cm (> 20 vs. ≤20) | 2.579 (1.089–6.106) | 0.031 | 5.299 (1.302–21.567) | 0.02 |
Multifocal (single vs. multiple) | 2.059 (0.899–4.715) | 0.088 | | |
Residual tumor (yes vs. no) | 2.194 (0.909–5.299) | 0.081 | | |
Treatment outcome (CR vs. PD) | 0.170 (0.047–0.611) | 0.007 | 0.220 (0.052–0.935) | 0.04 |
Recurrence (yes vs. no) | 2.388 (0.953–5.983) | 0.063 | | |
Riskscore (high vs. low) | 26.477 (5.915–118.520) | < 0.001 | 10.296 (1.960-54.073) | 0.006 |
Bold values identify statistical significance (p < 0.05). |
Abbreviations: OS, Overall Survival; HR, hazard ratio; CI, confidential interval; CR, Complete Response; PD, Progressive Disease. |
Variables with P < 0.05 in the univariate analysis were included in the multivariate analysis. |
Given the significant high local recurrence rate in clinical treatment of RPLS[24], 21 LMAGs were also found to be considerably linked to the DFS of patients through LASSO regression analysis (Fig. S4E), 2 of which were tested and selected for the prediction model in the multivariate Cox analysis (Fig. S4F). The association between risk score and recurrence status was next evaluated, and the low-risk subgroup harbored significant more alive statuses (Fig. S4G) and better DFS than those in high-risk subgroup (Fig. S4H). Similarly, this DFS related model also presented a great accuracy with AUC values of 0.72 in 1 year, 0.89 in 3 years and 0.79 in 5 years (Fig. S4I). In addition, univariate Cox analysis indicated that patients with worse DFS were authenticated with tumor residue, poor chemotherapeutic efficacy and high risk score (Table 2). Multivariate Cox analysis further confirmed that dismal chemotherapeutic efficacy and high risk score were independent risk factors (Table 2). Furthermore, a nomogram for DFS prediction was conducted with excellent consistency (Fig. S4J).
Table 2
Univariate and multivariate analysis of DFS in cohort-TCGA (n = 58).
Variables | Univariate analysis | Multivariate analysis |
HR (95%CI) | P value | HR (95%CI) | P value |
Age, years (> 60 vs. ≤60) | 2.686 (0.941–7.666) | 0.065 | | |
Sex (female vs. male) | 1.217 (0.561–2.642) | 0.619 | | |
Tumor size, cm (> 20 vs. ≤20) | 1.619 (0.747–3.5008) | 0.222 | | |
Multifocal (single vs. multiple) | 1.594 (0.766–3.319) | 0.213 | | |
Residual tumor (yes vs. no) | 2.169 (1.017–4.626) | 0.045 | 1.880 (0.740–4.779) | 0.185 |
Treatment outcome (CR vs. PD) | 0.109 (0.045–0.265) | < 0.001 | 0.144 (0.056–0.372) | < 0.001 |
Riskscore (High vs. Low) | 3.662 (1.748–7.670) | < 0.001 | 2.854 (1.088–7.491) | 0.033 |
Bold values identify statistical significance (p < 0.05). |
Abbreviations: DFS, disease free survival; HR, hazard ratio; CI, confidential interval; CR, Complete Response; PD, Progressive Disease. |
Variables with P < 0.05 in the univariate analysis were included in the multivariate analysis. |
Identification Of Lipid Metabolism-associated Tumor Antigens
To explore key genes that functioned as potential candidates for mRNA vaccine targets, we further systematically screened and identified two candidates (NR1H4 and ELOVL2) with both gene amplification and mutation, which were also associated both OS and DFS from the 135 LMAGs (Fig. 5A). Given the essential role of antigen-presenting cells (APCs) in the function of mRNA vaccines, we also evaluated the relationship of these two genes with APCs using TIMER analysis[25]. Intriguingly, ELOVL2, but not NR1H4, was identified with closely related to APCs (Pearson correlation coefficient > 0.3; Fig. 5B, S5A), which could serve as a potential tumor antigen and triggered strong immune response. In addition, similar results were also found in TCGA-SARC (Fig. S5B). Notably, survival analysis demonstrated that high mRNA expression of ELOVL2 was associated with unfavorable OS and DFS (Fig. 5C, D), suggesting ELOVL2 was of importance in RPLS development and progression. In concordance with cohort-TCGA, the mRNA expression of ELOVL2 validated similar prognostic efficiency in both cohort-FD and cohort-GSE30929 (Fig. S5C, D).
Dedifferentiated liposarcoma (DDLPS) was often progressed from primary or recurrent well-differentiated liposarcomas (WDLPS), which constituted the most common pathological type of RPLS[26]. Thus, we also investigated the heterogeneous expression of ELOVL2 in WDLPS and DDLPS. Interestingly, the expression of ELOVL2 was significantly higher in LMS2 than that in LMS1 (Fig. 5E). Specifically, ELOVL2 exhibited a potential enrichment in DDLPS, as compared with WDLPS in both cohort-FD and cohort-GSE30929 (Fig. S5E, F), indicating the significant impacts of ELOVL2 in the dedifferentiation evolution of RPLS.
To further identify the association between immune status and ELOVL2 expression, a series of TIME profiles was also conducted. Accordingly, immune score, stromal score, ESTIMATE score (Fig. 5F), infiltrating immune cells and TLS signature (Fig. S6A, B) were all calculated. As expected, we discovered that significant difference was found in the immune status between ELOVL2high and ELOVL2low subgroups, suggesting that RPLS patients with high expression of ELOVL2 might towards a status of immune desert. In addition, GO enrichment analysis indicated that high expression of ELOVL2 was involved in immune system process and cell surface reportor signaling. Similarly, KEGG enrichment analysis also validated cytokine-cytokine receptor interaction and cytokine signaling pathway (Fig. 5G, H).
Accumulating evidence indicated that the increase of DNA methylation of ELOVL2 leaded to the decrease of its protein expression and polyunsaturated fatty acid synthesis, but the accumulation of short chain fatty acids, which is closely related to aging[27]. Thus, we investigated the genome and epigenome landscape of ELOVL2 in cohort-TCGA (Fig. S6C), and 13 ELOVL2 related CpG sites were identified. However, only the cg20462795 exhibited significantly survival patterns (Fig. S6D), depicting that ELOVL2 associated epigenetic metabolic axis could be a novel therapeutic target in RPLS.
Transcription factor (TF) is one of the most common tool that involving in regulating gene expression[28]. Thereafter, we further systematically screened and identified three ELOVL2 related TFs (TFDP1, TP73 and MYBL2) in TCGA-SARC cohort (Fig. S6E). Interestingly, high expression of them were significantly associated with decreased OS (Fig. S6F). Consistent with survival prediction occurring in our cohort, TFDP1 and MYBL2 exhibited similar survival patterns in both TCGA-SARC and GTEx databases (data not shown). Taken together, these results demonstrated that gene expression was determined by the synergistic regulation of both transcription and epigenetic factors.
Elovl2 Dominated Lipid Metabolism Reprogramming And Executive Time Affect Prognosis In Rpls
Considering the dysregulation of immune status and metabolism remodeling involved in RPLS, we discovered that ELOVL2 played an essential role in tumor progression and even dedifferentiation, and served as the only one lipid metabolism related tumor antigen for potential mRNA vaccine. Moreover, PLCG1, ZC3HAV1L and NFAT5 were selected as immune hub genes through immune gene co-expression modules analysis and predicted great prognostic efficacy. Therefore, we re-analyzed the association between ELOVL2 and three hub genes.
Firstly, we validated these hub genes in cohort-FD and cohort-GSE30929 through survival analysis. Interestingly, the expression of PLCG1 was positive associated with ELOVL2 in all cohorts (Fig. 6A). In addition, PLCG1 was significantly higher in LMS2 than that in LMS1, and exhibited similar survival patterns in all cohorts (Fig. 6B). Next, the combined prognostic analysis of ELOVL2 and PLCG1 was also performed. Patient stratification based on these three groups presented that the Group I was associated with favorable prognosis, whereas the Group III was associated with dismal prognosis, and Group II was associated with medium prognosis (Fig. S7A). Subsequently, we developed another prognostic model for OS prediction based on ELOVL2 and PLCG1, which suggested a great prognostic efficacy (Fig. 6C).
In order to systematically comprehend the extent of relevancy between ELOVL2 and PLCG1 within TIME, we re-analyzed the expression profile in cohort-TCGA. Based on comparative analysis, we observed remarkable positive correlation between ELOVL2 and PLCG1. In addition, the infiltration of T cells, monocytics, dendritic cells, MDSC and M2-TAM were significantly lower with high expression of ELOVL2 and PLCG1. Moreover, the chemokine enrichment of CCL3, CCL4, CCL5 and CCL18, as well as PDCD1LG2 and HAVCR2 were significantly negative associated with high expression of ELOVL2 and PLCG1 (Fig. 6D). These results demonstrated that highly expression of ELOVL2 and PLCG1 in TIME of RPLS was involved with a immune-exculded phenotype.
To fully characterize the specific cellular localization ELOVL2 and PLCG1 in RPLS, we first evaluated them in the human protein atlas and further validated them in situ single cell spatial phenotype analysis through single cell RNA-sequencing and immunohistochemistry (IHC) in cohort-FD. Interestingly, ELOVL2 was seemed enriching in fibroblasts (no images in database), while PLCG1 in T cells (Fig. S7B). At the single-cell level, ELOVL2 was expressed in tumor cells, cancer associated fibroblasts (CAFs) and smooth muscle cells (SMCs), while PLCG1 was presented in CD4 + and CD8 + T cells, as expected (Fig. 6E; S7C). Moreover, the densities of CD3 + and cytotoxic CD8 + T cells in the area of tumor and invasive margin were quantified by IHC (Fig. S7D), and the immunoscore was evaluated for each patient (Fig. S7E). We observed that the density of CD3+ and cytotoxic CD8+ T cells were significantly but negatively associated with ELOVL2+ cells (Fig. 6F). In concordance with previous findings, the mRNA expression level of CD3 and CD8 was negatively associated with ELOVL2 (Fig. S7F). Taken together, the dysregulation of lipid metabolism might remodel an immune desert of TIME further supported their crucial roles in the evolution of RPLS.