Immune-Related RNA-Binding Protein-Based Signature With Predictive and Prognostic Implications in Patients With Lung Adenocarcinoma.

Background: Dysregulation of RNA-binding proteins (RBPs) in cancers is associated with immune and cancer development. Here, we aimed to profile immune-related RBPs in lung adenocarcinoma (LUAD) and construct an immune-related RBP signature (IRBPS) to predict the survival and response to immunotherapy.
Methods: A correlation analysis was performed to establish a co-expression network of RBPs and immune-related genes (IRGs) to characterize immune-related RBPs in the TCGA-LUAD cohort (n = 497 cases). Then, a combination of the Random survival forest (RSF) and Cox regression analysis was performed to screen the RBPs and establish IRBPS. This was followed by independent validation of IRBPS in GSE72094 (n = 398 cases), GSE31210, (n = 226 cases), and GSE26939 (n = 114 cases). Differences between the low- and high-risk groups were compared in terms of gene mutations, tumor mutation burden, tumor-infiltrating lymphocytes, and biomarkers responsive to immunotherapy.
Results: DDX56, CTSL, ZC3H12D, and PSMC5 were selected and used to construct IRBPS. The high-risk scores of patients had a significantly worse prognosis in both training and testing cohorts (p < 0.0001 and p < 0.05, respectively), and they tended to be older and have an advanced TNM stage. Furthermore, IRBPS was a prognostic factor independent of age, gender, smoking history, TNM stage, and EGFR mutation status (p = 0.002). In addition, high-risk scores of IRBPS were significantly correlated with tumor-infiltrating lymphocytes (p < 0.05). They also had a high level of PD-L1 protein expression (p < 0.01), number of neoantigens (p < 0.001), and TMB (p < 0.001), implying the possible prediction of IRBPS in the immunotherapy of LUAD.
Conclusion: The currently established IRBPS encompassing immune-related RBPs might serve as a promising tool to predict survival, reflect the immune microenvironment, and predict the efficacy of immunotherapy among LUAD patients.

Introduction

Worldwide, 2.2 million new lung cancer cases were diagnosed, and 1.8 million deaths occurred in 2020 (Sung et al., 2021). Although the incidence of lung cancer in the United States has decreased in recent years, it remains the leading cause of cancer deaths, with approximately 350 deaths per day projected by 2022 (Siegel et al., 2022). Non-small-cell lung cancer (NSCLC) accounts for about 85% of all lung cancer cases (Bade and Dela Cruz, 2020). Lung adenocarcinoma (LUAD) is the most common pathological subtype of NSCLC, accounting for about 63% of NSCLC (Gridelli et al., 2015). Although the multidisciplinary cooperative treatment strategy based on traditional treatment methods (surgery, chemotherapy, radiotherapy, and targeted therapy) reduces mortality, the LUAD 5-year overall survival rate remains 15.9% (Chen et al., 2014). Compared with traditional treatment methods, immunotherapy, especially immune checkpoint inhibitors, has changed the treatment pattern of lung cancer. Some patients can obtain a long-lasting treatment response, but this only accounts for 20–40% of all patients (Garon et al., 2019; Brahmer et al., 2018; Article, 2019; Garon et al., 2019; Brahmer et al., 2018; Article, 2019). Based on a high-throughput and micro-assay sequencing development, identifying and recognizing single markers or multifactor markers has become a hotspot in our current research (Ahluwalia et al., 2021; Yi et al., 2021; Zhang et al., 2020a). RNA-binding proteins (RBPs) are critical regulators of posttranscriptional processes (RNA splicing, modification, transport, cell localization, stability, translation, and degradation) (Wang et al., 2018; Martin and Ephrussi, 2009; Moore and Proudfoot, 2009; Fu and Ares, 2014; Martin and Ephrussi, 2009; Moore and Proudfoot, 2009; Fu and Ares, 2014; and Wang et al., 2018) and play an essential role in the occurrence and development of diseases, including cancer (Pereira et al., 2017). Identifying RBPs with potential, predictive, or prognostic functions may further broaden and open up new molecular markers.

Recently, the role of RBPs in the field of cancer has been recognized. The dysregulated expression of certain RBPs can lead to cancer. For example, the mutation of SF3B1 is observed in about 10–15% of patients with chronic lymphocytic leukemia (Wang et al., 2011); overexpression of mRNA 5′cap-binding protein eIF4E promotes the malignant transformation of human and mouse cells (Avdulov et al., 2004). Moreover, there are dozens of RBPs in the currently known cancer driver genes (Zhang et al., 2020b). In addition, some studies have proven that the dysregulation expression of RBPs in cancer relative to adjacent normal tissues is related to the efficacy and prognosis (Wang et al., 2018). For example, overexpression of MEX3B in tumors can reduce the lethality of tumor-infiltrating lymphocytes, thereby mediating tumor immunotherapy resistance (Huang et al., 2018); NONO promotes cancer proliferation by regulating SKP2 and E2F8, significantly affecting patient prognosis (Iino et al., 2020). More importantly, RBPs are considered essential regulators in the immune system (Jeltsch and Heissmeyer, 2016; Turner and Díaz-Muñoz, 2018; Newman et al., 2016; Jeltsch and Heissmeyer, 2016; Newman et al., 2016; and Turner and Díaz-Muñoz, 2018) and may be involved in immune modulation in the tumor microenvironment. However, RBPs associated with immune-related genes (IRGs) are rarely reported.

In this study, by establishing a co-expression network between RBPs and IRGs, we undertook a systematic and comprehensive biomarker discovery and validation effort to identify and develop an IRBPS for robust, prognostic, and efficacious prediction of LUAD patients. Herein, we first report a novel 4-gene immune-infiltrating and proliferation-associated IRBPS that not only identifies patients with poor prognosis but is also unique because of the close association with immune cell infiltration and immunotherapeutic biomarkers. In summary, the IRBPS gene signature provides an attractive platform for risk stratification of prognosis and efficacy in patients with LUAD, which has important implications for the clinical management of patients with this fatal malignancy.

Methods

Transcript Data Acquisition and Preprocessing

Genome-wide mRNA expression and clinicopathological information (age, gender, stage, smoking, and follow-up) of the TCGA–LUAD were downloaded as training cohorts from UCSC Xena (https://xenabrowser.net/); 497 LUAD primary tumor tissues matched with the complete follow-up recording were collected. Before further analysis, we converted the quantitative gene expression of fragment per kilobase million (FPKM) into transcripts per million (TPM) values processed by log2(value + 1) for all samples. Gene copy number data of the TCGA–LUAD quantified by the GISTIC2 approach were acquired from the cbioportal database (https://www.cbioportal.org/). Gene-level mutation (VarScan) of the TCGA–LUAD was obtained from the TCGA portal (https://www.cancer.gov/) by the R package “TCGAbiolinks.” The proteomic data of the TCGA-LUAD were downloaded from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) (https://cptac-data-portal.georgetown.edu/). Invalidation sets and normalized transcriptomic data recorded by series matrix files in gene expression series (GSE) including GSE72094, GSE31210, and GSE26939 were downloaded from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/), respectively. The probe IDs in the microarray datasets were converted into gene symbols. Only cases with complete follow-up recordings were collected; 6,094 unique candidate RBPs in Homo sapiens were offered from RBP2GO (https://rbp2go.dkfz.de/); 1,793 unique IRGs were extracted from the Immunology Database and Analysis Portal (ImmPort) database (https://www.immport.org/home).

Identification of Immune-Related RBPs With Prognosis Significance

Significant prognostic RBPs and IRGs in 497 LUAD samples were identified using the univariate Cox proportional-hazards regression analysis, respectively; genes were considered significant at p < 0.01. Immune-related RBPs were identified using the Pearson correlation analysis between RBP and IRG expression value. The threshold was set as a p-value < 0.01 and |correlation coefficient| > 0.7.

Construction and Validation of IRBPS at the Transcription Level

The random survival forest (RSF) algorithm is suitable for survival analysis of high-dimensional genomic data by reducing the dimensionality of features. The selection of genes was based on the variable importance (VIMP) and the minimum depth, where both were evaluation indexes to measure the predictive ability of variables for survival prediction. Here, integration of immune-related RBPs into an RFS approach was used to select genes and only choose the intersection of the top 15 genes in both VIMP and minimum depth. Before constructing the Cox regression model, RBPs with a high mutual covariance in the presence of 15 genes were removed by correlation analysis (the cutoff value of p < 0.01 and |correlation| > 0.7) to prevent overfitting of the constructed model, and the remaining genes were included in the multivariate analysis. The best model for predicting the prognosis of LUAD patients was determined based on the Akaike information criterion (AIC). The following formula calculated the risk score of each LUAD patient: risk score = expression of gene a * coefficient a + expression of gene b * coefficient ß + expression of gene c * coefficient ? + …… + expression of gene n * coefficient n. The median value of the risk scores stratified high- and low-risk groups. The Kaplan–Meier survival curves were applied for survival comparison between low- and high-risk groups, and the log-rank test was used to test statistical significance (p-value < 0.05). The performance of the IRBPS was evaluated and validated in the TCGA–LUAD cohorts and GEO–LUAD cohorts, respectively.

Expression Profile of Genes in IRBPS at the Transcription Level

In the TCGA–LUAD cohort, the Wilcoxon test was used to calculate the differential expression levels of genes in IRBPS between LUAD and normal lung tissues. Copy number variation (CNV) is a necessary part of genome structure variation, leading to variation in gene expression levels. Then, to explore whether the differences in gene expression in IRBPS are affected by CNV, we extracted CNV data defined by the GISTIC2 method and performed a correlation analysis between the gene expression values and CNV values.

Validation of Expression and Predictive Value of IRBPS at the Protein Level

Protein is the ultimate biological and functional unit for gene function. We separated proteomic data of genes in IRBPS from the Clinical Proteomic Tumor Analysis Consortium (CPTAC, https://proteomics.cancer.gov/programs/cptac) LUAD cohort to verify these gene expression features between LUAD and normal lung tissues. Immunohistochemical staining images were isolated from the Human Protein Atlas project (https://www.proteinatlas.org/) to verify the actual expression of these genes. Furthermore, 111 LUAD cases were obtained with complete follow-up data at an unshared log-ratio-quantified protein expression level. IRBPS at the protein level was also calculated using the aforementioned formula, and multivariate Cox regression analysis identified the independent role of the signature.

Clinical and Molecular Characteristics Between Low- and High-Risk Groups

The different prognostic characteristics between low- and high-risk groups might demonstrate that these patients have other clinical features and molecular profiles. Therefore, we performed a differential analysis incorporating clinical information such as age, gender, TNM stage, and survival status, as well as carving out the mutational characteristics of patients in different subgroups using the R package “maftools,” all of which were used to assess clinical and genomic differences between groups.

GO and KEGG Analyses

To further explore the differences in molecular functions and pathways between low- and high-risk groups, the Limma algorithm was used to perform differentially expressed genes (DEGs), with a false discovery rate (FDR) < 0.01 and | log2 (fold change) | > 0.5 as thresholds. Then, the potential biological processes and pathways involved in DEGs were explored through WebGestalt (http://www.webgestalt.org/).

Profile of Tumor-Infiltrating Immune Cells

In order to explore the relationship between IRBPS and tumor-infiltrating immune cells, the tumor immune score of the TCGA–LUAD cohort was extracted from the Tumor Immune Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/). The abundance of 24 immune cells in the TCGA–LUAD cohort was obtained from the immune cell abundance identifier (ImmuCellAI, https://cibersort.stanford.edu/). Thus, the immune score and immune cell components between the low- and high-risk groups were compared.

Exploring the Potential of IRBPS in Immunotherapy Response Prediction

To explore relationships between IRBPS and immunotherapy biomarkers, PD-L1 protein expression, neoantigen, and tumor mutation burden (TMB) data were collected for TCGA–LUAD. PD-L1 protein expression values were extracted from cBioPortal (https://www.cbioportal.org/) calculated by the reverse-phase protein array (RPPA) analysis. TMB and neoantigen datasets were attained from the Cancer Immunome Atlas (TCIA) (https://tcia.at/home). The aforementioned immunotherapy biomarker data were pooled and compared between the low- and high-risk groups.

Statistical Analysis

The Wilcoxon test and the chi-square test were used in comparisons of the two groups for continuous variables and categorical variables, respectively. All statistical analyses were performed using R software (version 4.0.2). p < 0.05 was considered statistically significant.

Results

Construction of a Co-Expressed Network Between Selected RBPs and IRGs

Overall, 1,367 RBPs and 194 IRGs were significantly associated with prognosis (p < 0.01) in 497 LUAD patients, respectively (Supplementary Tables S1, S2). A two-by-two correlation line analysis of RBPs and IRGs was performed to establish co-expression networks (Supplementary Table S3). Finally, 309 immune-related RBPs were identified based on screening criteria (|correlation| > 0.7 and p < 0.01).

Establishment of IRBPS for Predicting Prognosis in the TCGA–LUAD Cohort

309 RBPs were fed into the RFS algorithm. The forest prediction error stabilized as the number of trees increased (Figure 1A). The top 15 RBPs ranked by minimum depth are listed in Figure 1B, where the yellow markers (S100A16, DDX56, S100A10, CTSL, ZC3H12D, and PSMC5) were genes that were also in the top 15 of the VIMP ranking. A correlation analysis between the genes was performed for S100A16, DDX56, S100A10, CTSL, ZC3H12D, and PSMC5, and all genes satisfied mutual independence (Supplementary Figure S1). The stepwise Cox analysis illustrated that the model combining DDX56, CTSL, ZC3H12D, and PSMC5 had the smallest AIC as the best model. Therefore, the multivariate Cox regression model combining four genes was developed as IRBPS. The risk score for each patient was calculated according to the formula: 0.2127 * DDX56 + 0.2453 * CTSL - 0.2487 * ZC3H12D + 0.1431 * PSMC5. Patients were classified into the low- and high-risk groups based on the median risk score as the cutoff point (cutoff value = −0.01). The Kaplan–Meier survival curves demonstrated that patients in the high-risk group had a significantly worse overall survival (OS) in patients with higher risk scores (p < 0.0001) (Figure 1C). The area under the ROC curves (AUC) for 1-, 3-, and 5-years was 0.684, 0.625, and 0.662, respectively (Figure 1D). Furthermore, the results of the univariate and multivariate Cox model analyses confirmed that the risk score-based IRBPS was a significant prognostic factor independent of age, gender, smoking history, TNM stage, and EGFR mutation status (HR = 1.63, 95% CI 1.19-2.22, p = 0.002) (Table 1).

FIGURE 1

Gene selection and model construction. (A) Error plots using the out-of-bag prediction of error estimator based on trees. (B) Overlap genes (annotated yellow) ranked by minimum depth and VIPM. (C) Kaplan–Meier curves and log-rank test of the overall survival based on the low- and high-risk groups in the TCGA–LUAD. (D) ROC curves for evaluating the prediction performance of the signature in the TCGA–LUAD.

TABLE 1

Univariate and multivariate Cox analysis for prognosis in the TCGA-LUAD cohort.

	Univariate model			Multivariate model
Variables	HR	CI 95%	P Value	HR	CI 95%	P Value
Age (years)	1.01	0.99–1.02	0.305	-	-	-
Gender	-	-	-	-	-	-
female	1	-	-	-	-	-
male	0.75	0.78–1.41	0.747	-	-	-
Stage	-	-	-	-	-	-
I	1	-	-	1	-	-
II	2.47	1.71–3.57	0.000	2.31	1.60–3.34	0.000
III	3.55	2.42–5.22	0.000	3.19	2.16–4.71	0.000
IV	3.88	2.24–6.75	0.000	3.61	2.07–6.30	0.000
Smoking	-	-	-	-	-	-
no	1	-	-	-	-	-
yes	1.05	0.77–1.43	0.747	-	-	-
EGFR	-	-	-	-	-	-
wild	1	-	-	-	-	-
mutation	0.69	0.26–1.87	0.472	-	-	-
Risk Score	-	-	-	-	-	-
low	1	-	-	1	-	-
high	1.89	1.39–2.55	0.000	1.63	1.19–2.22	0.002

Validating the Predictive Robustness of IRBPS in Independent Cohorts

We assessed the robustness of IRBPS derived from the TCGA–LUAD cohort in three independent GEO cohorts, including GSE72094, GSE31210, and GSE26939. Patients in these cohorts were classified into low- and high-risk groups using the median risk score cutoffs calculated with the same formula. The results of the Kaplan–Meier curve showed that IRBPS significantly stratified patients in the GSE72094 cohort (cutoff value = −0.03, p = 0.026) (Figure 2A), GSE31210 cohort (cutoff value = −0.02, p = 0.029) (Figure 2B), and GSE26939 cohort (cutoff value = 0.01, p = 0.005) (Figure 2C). AUC values illustrated the good predictive performance of IRBPS in the GSE72094 cohort (0.609, 0.525, and 0.617 for 1-, 3-, and 5-years), GSE31210 cohort (0.584, 0.623, and 0.693 for 1-, 3-, and 5-years), and GSE26939 cohort (0.730, 0.659, and 0.661 for 1-, 3-, and 5-years).

FIGURE 2

Validation of the signature in independent cohorts. (A) GSE72094, (B) GSE31210, and (C) GSE26939.

We collected 497 LUAD samples and 58 normal lung samples with the transcription data. The expression levels of DDX56 (Figure 3A), ZC3H12D (Figure 3C), and PSMC5 (Figure 3D) were all significantly increased in LUAD tissues compared to normal lung tissues. The CTSL expression had a higher level in normal lung tissues than LUAD tissues (Figure 3B). A diploid implied the absence of CNV at this gene location. Compared with the diploid, other types (deep deletion and shallow deletion gain amplification) of CNV were significantly correlated with gene expressions of DDX56 (Figure 3E), ZC3H12D (Figure 3G), and PSMC5 (Figure 3H) (p < 0.05), except for CTSL (Figure 3F) (p > 0.05), indicating that CNV might regulate the expression of these genes in cancer.

FIGURE 3

Expression profile of RBPs in IRBPS at the transcription level. (A–D) Expression differences between LUAD tissues and normal lung tissues. (E–H) The correlation of CNV with the four genes’ expression levels. Diploid implied the absence of CNV at this gene location. Compared with diploid, investigating the effect of other types (deep deletion and shallow deletion gain amplification) of CNV on their gene expression.

Expression Profile of Genes in IRBPS at the Protein Level

At the protein level, proteomic data of 111 LUAD tissues and 106 normal tissues were collected. As shown in Figure 4A, LUAD tissues had a significantly higher expression level of these four genes (DDX56, ZC3H12D, and PSMC5) than normal lung tissues. The CTSL expression had a higher level in normal lung tissues than in LUAD tissues. Immunohistochemical images of the LUAD tissues and normal lung tissues for DDX56, ZC3H12D, and PSMC5 are shown in Figure 4B, respectively. CTSL, due to lack of documentation, was not available for immunohistochemical images. As illustrated in Figure 4C, the high protein level of CTSL and PSMC5 had a significant adverse impact on the survival of LUAD patients. The effect of DDX56 and ZC3H12D on survival did not show statistical significance and might be limited by the sample size. In addition, the high-risk score of IRBPS still adversely affected the survival of LUAD patients and was not affected by age, gender, and TNM stage at the protein level. This further demonstrated the ability of IRBPS for predicting the patient prognosis at the protein level.

FIGURE 4

Expression profile of genes in IRBPS at the protein level. (A) Gene expression differences between LUAD tissues and normal lung tissues at the protein level. (B) Immunohistochemical staining for DDX56, ZC3H12D, and PSMC5 in LUAD tissues and normal lung tissues. (C) Kaplan–Meier curves, log-rank test, and multivariate Cox analysis of the overall survival based on the expression level of these genes and IRBPS at the protein level.

Clinical and Genomic Mutational Features Between IRBPS Subgroups

As shown in Figures 5A–G, patients in the high-risk group tended to be older and had advanced TNM stage than the low-risk group. The 20 most frequently mutated genes in the low- and high-risk groups were calculated, respectively (Figures 5H,I). By comparing between subgroups, 19 genes including TP53, KEAP1, PAPPA2, LPA, TRPC7, CAD, SCN1A, GCC2, LRRC7, SMARCA4, CSMD3, BMPER, DYNC2H1, GTF3C1, SORCS1, TIE1, PPAPDC1A, AGO2, and WDFY3 were significantly mutated in the high-risk group (Figure 5J).

FIGURE 5

Clinical and mutational features between IRBPS subgroups. (A) Distribution of survival status and the four-gene expressions trend between IRBPS subgroups. Clinical characteristic differences between IRBPS subgroups in terms of age (B), gender (C), stage (D), T-stage (E), N-stage (F), and M-stage (G). (H,I) Ranked the top 20 of gene mutations in low- and high-risk groups. (J) The significant mutation genes in the high-risk group compared to the low-risk group.

Biological Processes Between IRBPS Subgroups

According to the criteria (FDR <0.01 and | log 2 (Fold Change) | > 0.7), there were 389 DEGs between the low- and high-risk groups (Supplementary Table S3), and they were applied to GO and KEGG analyses. The enrichment results (Figure 6A) indicated that these DEGs were more involved in the biological functions of the regulation of cell proliferation and cell cycle processes. KEGG analysis (Figure 6B) illustrated that these genes were closely associated with cell cycle, metabolism, p53 signaling pathways, and the B-cell reporter signaling pathway.

FIGURE 6

Enrichment results of GO (A) and KEGG (B) analyses. The color of the bar indicated the significance level; the length of the bar indicated the enrichment ratio.

Correlation of Tumor-Infiltrating Immune Cells and Immunotherapy Biomarkers With IRBPS Subgroups

Immune cell panorama was significantly different between low- and high-risk patients (p < 0.05) (Figure 7A), and the immune score calculated by the TIMER algorithm was significantly higher in the low-risk group (p < 0.001) (Figure 7B). As shown in Figure 7A, the abundance of nTreg (p = 0.000), Th17 (p = 0.018), monocytes (p = 0.000), macrophages (p = 0.013), and neutrophils (p = 0.000) were significantly higher in the high-risk group than those in the low-risk group. On the contrary, the proportion of CD8_naive (p = 0.047), cytotoxic (p = 0.000), Tfh (p = 0.000), central_memory (p = 0.000), B cells (p = 0.000), NK (p < 0.001), CD4_T (p = 0.000), and CD8_T (p = 0.000) in the low-risk group was significantly higher. In addition, the predictive value of IRBPS for widely recognized immunotherapy biomarkers was further evaluated. The results demonstrated that these biomarkers had a higher level in the high-risk group than in the low-risk group in terms of the PD-L1 protein expression level (Figure 7C, p < 0.01), TMB (Figure 7D, p < 0.001), and number of neoantigens (Figure 7E, p < 0.001).

FIGURE 7

Correlation of tumor-infiltrating immune cells and immunotherapy biomarkers with IRBPS subgroups. (A) Abundances of immune cells between subgroups. (B–E) Correlation of immunotherapy response biomarkers with subgroups.

Discussion

RBPs are a group of conserved proteins in eukaryotes that play an essential role in co-transcriptional and post-transcriptional gene regulation. RBPs can interact with RNA to form protein–RNA complexes and participate in biological processes such as cell differentiation, proliferation, and cell fate transition. The imbalance of RBPs has been implicated in various diseases. In particular, several studies have demonstrated that RBPs are closely associated with cancer development, and immune and therapeutic responses. However, studies on RBPs, especially immune-related RBPs in LUAD, are still rarely reported. The present study is the first comprehensive analysis of the prognostic impact of 6,094 RBP candidates in LUAD based on a bioinformatics analysis. In addition, by screening immune-related RBP genes, we established a risk score module (IRBPS) based on these genes. We validated their prognostic value at the transcriptional level with four datasets, respectively. Notably, at the protein level, IRBPS still played an essential role in predicting the survival of LUAD patients. Furthermore, IRBPS was closely associated with tumor-infiltrating immune cells and immunotherapeutic biomarkers, which implied that this signature has an essential role in predicting the efficacy of immunotherapy and deserves further investigation.

In this study, RFS was used for variable selection, followed by the traditional Cox regression analysis to construct the model instead of the most commonly used least absolute shrinkage and selection operator. First, RSF is a combinatorial survival tree method that inherits the advantages of random forest, such as noise resistance, over-fitting prevention, and nonlinear correlation processing. It can be used for variable selection for the survival prediction analysis of high-dimensional data. In addition, the RFS algorithm combines two forms of randomization methods: case resampling and variable sub-setting, to prevent the over-fitting problem, make the results more robust, and predict accurately. Then, the combination of the Cox regression analysis provides an exact formula to calculate risk scores, which is convenient for clinical application. In this study, 309 immune-related RBPs were screened out by establishing a co-expression network of RBPs and IRGs, and four RBPs (DDX56, CTSL, ZC3H12D, and PSMC5) were finally identified to construct IRBPS by RSF and Cox regression analysis. The robustness of IRBPS in predicting the prognosis was validated in three independent GEO datasets, respectively. DDX56, ZC3H12D, and PSMC5 were up-regulated in LUAD tissues compared to normal lung tissues. In contrast, CTSL is down-regulated in LUAD tissues. CNV might be involved in the process of regulating the mRNA expression of DDX56, ZC3H12D, and PSMC5. Here, the number of tumor tissues was greater than the number of normal tissues, which had limited utility for assessing differences in gene expression between the two. Therefore, a comparable and sufficient sample size between the two was required. The protein expression features of the four genes between LUAD tissues and normal lung tissues followed the same trend at the transcription level. Although the CPTAC–LUAD data recorded the protein expression values of four RBPs in 111 patients, there were many missing values of the ZC3H12D gene. We here filled in the missing values by the KNN method to initially assess the prognostic impact of each gene and IRBPS at the protein level. The results showed that high expression of CTSL and PSMC5 still negatively affected the prognosis of LUAD patients, and that IRBPS was an adverse prognostic factor independent of age, gender, and the TNM stage. These results directly or indirectly provide evidence that RBPs are dysregulated in LUAD and affect the prognosis of LUAD patients. Notably, patients in the high-risk group of IRBPS exhibited older age, more advanced tumor stages, had more genetic mutations, lower immune scores, and higher levels of immune efficacy marker indices. This also explains, in one way, why the prognosis is worse in the high-risk group of IRBPS, and it also suggests that patients with high-risk scores of IRBPS may have a higher possibility of immunotherapeutic efficacy.

Previously, many studies have shown that RBPs are associated with immune responses and immune regulation. In the innate immune response, some RBPs (such as TTP, Roquin, and Regnase-1) promote inflammation by eliminating the mRNA of pro-inflammatory cytokines (Mino and Takeuchi, 2018), controlling the adaptive immune responses, and maintaining immune homeostasis during T-cell activation (Jeltsch and Heissmeyer, 2016). This study found that IRBPS was significantly associated with the tumor immune score (TIMER) and immune cell infiltration (ImmuCellAI). Based on these results, we explored the potential of IRBPS in predicting the immune response. The results showed higher TMB, PD-L1 protein expression, and the number of tumor neoantigens in the high-risk group than in the low-risk group. These indicators are well-established biological biomarkers for predicting immunotherapeutic efficacy in LUAD. DDX56 encodes a member of the DEAD-box protein family. The pathway analysis showed that DDX56 regulated the p53 signaling pathway, affected cell cycle progression, and promoted cancer cell proliferation (Zhu et al., 2020). The high expression level of the DDX56 gene was correlated with a high level of TP53 mutation and was positively correlated with the MYC signal (Zhu et al., 2020; Cui et al., 2021). MYC could up-regulate PD-L1 protein by interacting with PD-L1 enhancer genes (Casey et al., 2016). TP53 mutation also promoted the PD-L1 protein expression (Cha et al., 2016; Dong et al., 2017). CTSL promoted tumor angiogenesis by regulating the CDP/Cux/VEGF-D pathway (Pan et al., 2020a). VEGF was the most important angiogenesis factor, which inhibited antigen presentation, promoted regulatory T cell infiltration, and induced PD-L1 expression on tumor-infiltrating T cells (Voron et al., 2015). ZC3H12D was a tumor suppressor gene associated with memory T lymphocytes and macrophages, involved in the regulation of inflammation (Wang et al., 2007; Huang et al., 2012). Previous studies showed that ZC3H12D was associated with a variety of immune-related genes, including CD274 (PD-L1) and CTLA4 (Chen et al., 2022). ZC3H12D could regulate IL-6 expression, and IL-6 was positively correlated with the PD-L1 expression (Wawro et al., 2017; Chan et al., 2019). IL-6 enhanced the glycosylation initiation and stability of PD-L1, which was necessary to maintain the stability of PD-L1 and its interaction with PD-1 (Chan et al., 2019). PSMC5 not only has the proteasome function but also participates in transcriptional regulation. PSMC5 regulated the epithelial–mesenchymal transition, hypoxia, immune response, and other pathways, and was positively correlated with gene expression such as CD274 and CTLA4 (He et al., 2021). Obviously, high expression of the four genes (DDX56, CTSL, ZC3H12D, and PSMC5) could promote the expression of PD-L1 as they have been proven to be promising biomarkers and immunotherapy targets in a variety of tumors, including colorectal cancer (Kouyama et al., 2019; He et al., 2021), NSCLC (Wu et al., 2021; Chen et al., 2022), and gastric cancer (Pan et al., 2020b).

Although sufficient verification was achieved in the TCGA–LUAD database and the GEO–LUAD database, the performance of IRBPS in predicting the survival status was still limited (AUC<0.7). IRBPS combined with clinical features to jointly build the model might further enhance the predictive power of the model. Second, this study’s incorporation of multi-omics data was limited to the TCGA portal, which prevented us from extensively validating the robustness of IRBPS. Furthermore, the absence of transcriptome data from patients receiving immunotherapy prevented further evaluation of the potential of IRBPS to predict its immune response in the real world. In addition, future cellular and animal experiments were also warranted to prove this concept in this study. Despite such limitations, it is undeniable that our study provides important clues to elucidate the prognostic value of RBPs in LUAD, to facilitate the selection of beneficiary groups for immunotherapy, and even to develop new therapeutic strategies for patients with LUAD.

In summary, this study first comprehensively characterizes immune-related RBPs with a prognostic value in LUAD and has constructed an IRBPS which offers an easy tool for identifying LUAD patients with worse prognosis and reflecting the efficacy of immunotherapy patients with this lethal malignancy.