A novel model based on liquid-liquid phase separation-Related genes correlates immune microenvironment profiles and predicts prognosis of lung squamous cell carcinoma.

OBJECTIVE: The aim of the study was to construct and validate a robust prognostic model based on liquid-liquid phase separation (LLPS)-related genes in lung squamous cell carcinoma (LUSC).
METHODS: The Cancer Genome Atlas dataset was used as the discovery set to identify the LLPS-related differentially expressed genes (DEGs) between LUSC and normal tissue. These DEGs were screened by the LASSO Cox regression analysis to identify the genes with nonzero coefficient, which were next included in the multivariate Cox regression analysis to construct the prediction model. The dataset GSE41271 was adopted as the validation set to verify the efficacy of the model. Enrichment analysis and the CIBERSORT were performed to illustrate potential immune mechanisms underlying the prediction model.
RESULTS: A total of 48 LLPS-related genes were aberrantly expressed in LUSC. Among them, 7 genes were selected by the LASSO Cox regression analysis to construct the prediction model. Risk index (RI) was calculated according to the model for each patient. The prognosis was significantly different between the patients with high and low RI in the discovery set and the validation set (p < 0.001 and p = 0.028, respectively). The multivariate survival analysis confirmed RI as an independent prognostic factor in LUSC (in the discovery set: p < 0.001, HR = 2.643, 95% CI = 1.986-3.518; in the validation set: p = 0.042, HR = 2.144, 95% CI = 1.026-4.480). A series of pathways involving immune cells were found to be related to RI. The distribution pattern of immune cells and chemokines varied according to the value of RI.
CONCLUSION: The prediction model based on LLPS-related genes was constructed and validated as a robust prognostic tool for LUSC using multiple datasets. LLPS might have an impact on LUSC through immune pathways.

INTRODUCTION

Non–small‐cell lung cancer (NSCLC), one of the most common malignancies worldwide, is characterized by high morbidity and mortality. Lung squamous cell carcinoma (LUSC) is known as a main pathological subtype of NSCLC. , Although molecular targeted therapy has been constantly applied in clinical practice, its effect against LUSC remained limited. , The prognosis of LUSC is still dismal, especially in patients with advanced‐stage disease. Recently, the promising outcome of immune checkpoint inhibitors (ICIs) such as anti‐PD‐1 agents has shed light on the anti‐LUSC strategies. , Its treatment response is mainly predicted by expression level of PD‐L1. , However, only a small group of patients with NSCLC can benefit from ICIs. More predictors for the efficacy of immunotherapy and more profound understanding of the underlying mechanisms become an urgent requirement for clinical practice.

The important roles and comprehensive mechanisms of liquid‐liquid phase separation (LLPS) of proteins and nucleic acids in cellular activities have been continuously uncovered by recent studies. , , Its correlation with the development of some human diseases has been indicated by a growing number of studies and has become a new paradigm in relevant researches. , Emerging evidence also revealed the correlation of LLPS‐related activities with anti‐tumoral immune response, providing a new perspective for understanding tumor immunology. , All these findings indicated the potential impact of LLPS on cancer development by its direct and indirect effects.

The aim of this study was to construct and validate a model based on LLPS‐related genes for predicting the prognosis of LUSC, and to preliminarily uncover the relationship between LLPS‐ and tumor‐related immune activities in LUSC.

MATERIALS AND METHODS

Study cohorts and LLPS‐related genes

The list of LLPS‐related genes was downloaded from PhaSepDB (a database provides a collection of manually curated phase separation proteins and membraneless organelles–related proteins; http://db.phasep.pro/). The gene expression data of LUSC patients were downloaded from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/). Among them, 550 cases with sufficient clinical and survival data were enrolled as the discovery set to identify the LLPS‐related differentially expressed genes (DEGs) between LUSC and normal tissue, and to construct the prediction model. Another dataset (GSE41271) obtained from Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) was used to validate the prediction model. The workflow of this study is shown in Figure 1.

FIGURE 1

Workflow of the present study

Construction of a LLPS‐related prediction model

TCGA dataset was adopted to calculate the LLPS‐related DEGs between LUSC and normal tissue using “limma” package in R. The genes with fold changes >1.5 and p value <0.05 were considered to be statistically significant.

All the LLPS‐related DEGs in the discovery set were subjected to least absolute shrinkage and selection operator (LASSO) Cox regression. Then, the genes with nonzero coefficients found in the LASSO Cox analysis by using the “glmnet” package were put into the multivariate Cox regression analysis to construct the prediction model. The model was expressed as:

In the formula, “Exp” represented the expression value of a gene, and “Coef” represented the coefficient of the gene.

Prognostic value of the prediction model

RI was calculated for each individual in the discovery set according to our prediction model. The prognosis of 2 patient groups divided by the level of RI was compared to illustrate the prognostic value of the model. Then, RI and other clinical characteristics were subjected to the multivariate Cox regression analysis to confirm the independence of RI (cases with missing data for these variables were excluded).

All the aforementioned analyses were additionally performed in the validation set to confirm the efficacy of the model.

Enrichment analysis and the CIBERSORT method

Enrichment analysis was conducted in order to identify the related biological pathways through the Gene Ontology (GO) database and Kyoto Encyclopedia of Genes and Genomes (KEGG) database by using “clusterProfiler” and “enrichplot” packages in R. The significant cutoff threshold was defined as adjusted p value <0.05 and |NES|>1.5.

The CIBERSORT method, known as a useful tool to convert gene expression profile to relative proportion of different types of immune cells, was performed to assess the distribution of tumor‐infiltrating immune cells (TICs) in LUSC using TCGA data.

Statistical analysis

All the statistical analysis was performed using R (version 4.1.1) (https://www.r‐project.org). Statistical significance was defined as a two‐sided p value <0.05 for all results in this study, unless otherwise stated.

RESULTS

Construction of the prediction model based on LLPS‐related genes

LLPS‐related genes were derived from PhaSepDB and matched in the discovery set. A total of 43 LLPS‐related genes were identified to be up‐regulated, and 5 of them were down‐regulated in LUSC compared with normal tissues (Figure 2A). Next, all the LLPS‐related DEGs of 493 LUSC patients were included in the LASSO Cox regression model. After the screening, 7 genes remained with nonzero regression coefficient, including HOXA13, SSX1, MAGEA4, DPPA2, PLA2G1B, CALML5, and H2BC9 (Figure 2B,C). These selected genes were used to construct a prognostic model for LUSC patients using multivariate Cox regression analysis (Figure 2D). The formula of the model was expressed as follows: Risk index (RI) = −0.03359×HOXA13+0.11907×SSX1+0.01246×MAGEA4+0.06373×DPPA2+0.09617×PLA2G1B +0.02794×CALML5‐0.13094×H2BC9.

FIGURE 2

Construction of the prognostic model based on liquid‐liquid phase separation (LLPS)‐related genes. (A) The volcano plot showing the different expression of LLPS‐related genes between lung squamous cell carcinoma (LUSC) and normal lung tissue. (B) The coefficient profiles of the 7 LLPS‐related genes selected by the least absolute shrinkage and selection operator (LASSO) Cox regression analysis. (C) Selection of the optimal parameter in the LASSO Cox regression analysis with 10‐fold cross‐validation. (D) The 7 LLPS‐related genes constituting the prediction model in the discovery set

The area under the curve (AUC) of the model was calculated to reveal its efficacy. The AUC at 1 year, 3 years, and 5 years equaled to 62.76, 64.24, and 65.40, respectively (Figure 3A). The result of the log‐rank test indicated that RI was significantly related to OS in LUSC (p < 0.001). Patients with high RI had a poorer prognosis in comparison with those with low RI (Figure 3B). Additionally, RI and the other clinical characteristics including age, gender, and cancer stage were included in the multivariate Cox analysis. As a result, the prognostic value of RI remained significant (p < 0.001, hazard ratio [HR] = 2.643, 95% confidence interval [CI] = 1.986–3.518; Figure 3C).

FIGURE 3

Predictive value of the prognostic model. (A) The area under the curve (AUC) of the prediction model in the discovery set. (B) Significantly different overall survival (OS) between patients with high and low risk index (RI). (C) RI as an independent prognostic factor in lung squamous cell carcinoma (LUSC) in comparison with routine clinical characteristics

Validation of the prediction model

Patients were divided into two groups according to the level of RI in the validation set. As shown in Figure 4A, high RI was significantly related to shorter OS (p = 0.028). After multivariate analysis, RI remained as an independent predictor for the prognosis of LUSC (p = 0.042, HR = 2.144, 95% CI = 1.026–4.480; Figure 4B).

FIGURE 4

Validation of the liquid‐liquid phase separation (LLPS)‐related prediction model. (A) The worse prognosis of patients with high risk index (RI) compared to those with low RI in the validation set. (B) The independent predictive value of RI in the validation set

Correlation between LLPS and immune in LUSC

In order to identify the plausible mechanisms to illustrate the prognostic effect of RI, enrichment analysis was performed between patients with low and high RI in TCGA dataset. Interestingly, a series of immune‐related pathways were found to be associated with RI in LUSC, such as “cytokine production” (NES = 1.818, adjusted p = 0.003), “activation of immune response” (NES = 1.704, adjusted p = 0.003), and “leukocyte‐mediated immunity” (NES = 1.781, adjusted p = 0.003) (Figure 5A,B). Next, the distribution pattern of tumor‐infiltrating immune cells was compared according to the level of RI. The results indicated that natural killer cells and CD8+ T cells were more abundant in the low‐RI group (p < 0.05 and p < 0.01, respectively), and more resting CD4+ memory T cells were found in the high‐RI group (p < 0.05; Figure 6A). The phenotypes of CD8+ cells were further analyzed between low‐ and high‐RI groups. Exhaustion phenotype markers such as PDCD1 (PD‐1) and HAVCR2 (TIM‐3) were highly expressed in the patients with high RI (p < 0.05 and p < 0.001, respectively; Figure 6B). The expression of chemokines from CC, CXC, XC, and CX3C subfamily was also analyzed according to the level of RI. Among them, CXCL12, CCL23, and CCL14 were found to be correlated with RI (|r|>0.3, p < 0.05; Figure 6C–E).

FIGURE 5

Immune‐related pathways involving in risk index (RI) identified by the Gene Ontology (GO) term enrichment analysis (A) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (B)

FIGURE 6

Distribution patterns and phenotypes of tumor‐infiltrating immune cells (TICs) and immune molecules in patients with high and low risk index (RI). (A) Distribution differences in TICs according to the value of RI. (B) Higher density of exhausted CD8+ T cells in the high‐RI group than in the low‐RI group. (C–E) The correlation between the level of RI and CXCL2, CCL23, and CCL14

DISCUSSION

Aberrant LLPS of molecules in living cells might trigger abnormal cellular activities, leading to human diseases. Previous studies presented robust evidence regarding the various LLPS‐related biological processes in the progression of neurodegenerative diseases. , As for the impact of LLPS on the development of cancer, relevant studies were limited. In the present study, the expression of LLPS‐related genes in LUSC was significantly different from normal tissues, which indicated the potential roles of LLPS processes in the cancer development.

The LASSO Cox regression analysis was known as an efficient method dealing with a high‐dimensional data and avoiding overfitting problems. , Thus, LASSO Cox was adopted in this study to provide the optimal solution for the prediction model by screening the suitable genes from a large number of candidates with low mutual correlation. In addition, the efficacy of the prediction model constructed by the selected genes was verified by the validation set.

The indispensable role of TICs in different types of cancer including LUSC has been reported by numerous studies. , , Complex immune mechanisms illustrated the various effects of immune cells and immune molecules on cancer development, providing potential therapeutic targets and anti‐tumor strategies. , Recently, some studies tried to interpret immune activities in tumor microenvironment from the perspective of LLPS. , However, this field is still in its infancy. In the present study, our analysis identified a series of immune‐related processes were correlated with LLPS in LUSC, including pathways involving different immune cells, adaptive and innate immune response, and immune checkpoint. Thus, we further investigated the distribution pattern of TICs and immune molecules with different LLPS‐related gene signatures. As a result, the abundance of some critical cells in anti‐cancer immune activities, including CD8+ T cells, NK cells and CD4+ T cells, varied according to the value of RI. Additionally, the functions of CD8+ T cells were altered in patients with high RI in comparison with low RI. It was noteworthy that these biomarkers indicating the function of CD8+ T cells were also known as immune checkpoint, such as PD‐1 and TIM‐3. The differential expression of immune checkpoint was correlated with the level of RI, which implied the potential predictive value of RI for the therapeutic effectiveness of ICIs. The significant correlation between the expression of some chemokines and RI further indicated the immune effect of LLPS‐related processes in LUSC.

There were several limitations of this study. First, all the datasets were obtained from retrospective studies. A prospective study should be performed in the future to confirm these findings. Second, our conclusion was based on bioinformatic analysis without verification of molecular biological experiments.

In conclusion, the LLPS‐related gene–based prediction model was a robust prognostic tool for LUSC patients. The correlation between LLPS processes and immune activities might have an impact on the development of LUSC.

CONFLICT OF INTEREST

All authors have completed the ICMJE uniform disclosure form. The authors have no conflicts of interest to declare.

AUTHOR CONTRIBUTIONS

LZ and KZ conceptualized and designed the study. QB provided administrative support. ZZ and YL contributed to the provision of study materials or patients. WG collected and assembled the data. LZ analyzed and interpreted the data. All authors wrote the article and gave final approval of the article.