Development of multivariable risk signature based on four immune-related RNA-binding proteins to predict survival and immune status in lung adenocarcinoma.

Background: This study aimed to construct a risk signature with predictive power based on immune-related RNA-binding proteins (RBPs) for lung adenocarcinoma (LUAD) patients.
Methods: The Cancer Genome Atlas (TCGA) database was used as the data source. Immune genes (IGs) were obtained from the Immunology Database and Analysis Portal (immPort) database. Differentially expressed RBPs and IGs between tumor and normal tissues were screened. For external validation, an independent cohort from the Gene-Expression Omnibus (GEO) database was used. The accuracy of the risk signature prediction was evaluated using Cox regression analysis and the receiver operating characteristic (ROC) curve.
Results: The risk signature was constructed from four immune-related and prognostic RBPs (OAS3, PCF11, TLR7, and EXO1). The patients were divided into the low- and high-risk groups, with the low-risk group having a higher survival rate than the high-risk group. The risk signature outperformed other clinical parameters, with a multivariable hazard ratio of 1.862 (95% confidence interval: 1.292-2.683). The tumor immune microenvironment, stemness index, immune checkpoint, immune infiltration, and proportion of immune cells were significantly different between the low- and high-risk groups (all P<0.05). Conclusions: The risk signature of immune-related RBPs can provide the basis for clinical decisions regarding diagnosis, prognosis, and immunotherapy in LUAD patients. 2022 Translational Cancer Research. All rights reserved.

Introduction

Lung cancer is the first leading cause of cancer-related deaths worldwide, accounting for almost a quarter of all cancer deaths (1). Lung cancer can be classified into two types based on histology: small cell lung cancer (SCLC) and non-small cell lung cancer (NSCLC); they are responsible for 15% and 85% of lung cancers, respectively. Lung adenocarcinoma (LUAD) accounts for >40% of lung cancers and is the most common type of NSCLC (2). In the last decade, immunotherapy has become an effective therapy for several cancers, and LUAD is an immune-sensitive cancer (3,4). Indeed, the immune status of cancer predicts the effectiveness of immunotherapy and is closely related to patients’ prognoses (5).

Studies have confirmed that several biomarkers are prognostic factors of LUAD, such as the epidermal growth factor receptor (EGFR), V-Ki-Ras2 Kirsten rat sarcoma 2 viral oncogene homolog (KRAS), and ALK receptor tyrosine kinase (ALK) (6-8). Many drugs have subsequently been developed to target driver gene mutations in these genes. Unfortunately, most patients will eventually develop resistance to targeted therapies after receiving them. For example, patients’ development of resistance to EGFR-tyrosine kinase inhibitors is part due to a secondary mutation in tumors (9). Moreover, changes at the genetic level of the tumor occur before the clinical characteristics, and patients with the same stage of cancer may have varying prognoses. Most patients are at an advanced stage when diagnosed and have thus missed the opportunity to undergo standard treatment. For these reasons, the 5-year survival rate of patients with LUAD is only 23% (10,11). Therefore, it is imperative to identify novel diagnostic and prognostic biomarkers of LUAD. Then, management schemes can be designed for distinct subsets of patients to allow for precise treatment.

RNA-binding protein (RBP) is a class of proteins that control RNA stabilization, degradation, and modification at the post-transcriptional level by selectively binding to RNA (12). So far, 1,542 RBP genes have been identified by genome-wide screening in the human genome, accounting for 7.5% of all protein-coding genes (13). Post-transcriptional regulation plays a critical role in cellular processes. Therefore, RBP dysregulation might lead to abnormal gene expression in cells and eventually result in various diseases, including cancer and immune disorders (14). For instance, by regulating the translation efficiency of Cyclin D1, CDK2, and CDK4, m6A modified RBP YTHDF1 deficiency inhibits the proliferation of NSCLC cells and the formation of xenograft tumors, and YTHDF1 decrease inhibits de novo LUAD (15). RBPs play important roles in various cancers (16-19). No study has tried to develop a risk signature for LUAD based on immune-related RBPs.

Lately, a study constructed a prognostic signature of liver cancer using immune-related RBPs and demonstrated good predictive power (20). To the best of our knowledge, no research has been reported to construct a risk signature for LUAD based on immune-related RBPs. Therefore, this study aimed to construct and validate a risk signature with predictive power based on immune-related RBPs for LUAD patients. Such a risk signature might help provide information for accurate prognosis and individualized immunotherapy of patients with LUAD. We present the following article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-22-698/rc).

Methods

Data sources

We used public data from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov) and Gene-Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo) databases to conduct a comprehensive analysis. As a training cohort, the expression data and clinical information of LUAD patients (486 tumor samples and 54 normal samples) were downloaded from TCGA, after excluding those with unknown survival data, overall survival (OS) <30 days, or died from non-cancer causes (e.g., myocardial infarction, hemorrhage, infection, etc.). GSE68465 from GEO (including 443 LUAD patients) was used as a validation cohort (Table S1). Gerstberger et al. (13) reviewed the literature to create a census of RBPs and identified 1,542 RBPs, which were all examined in the present study. Based on the Immunology Database and Analysis Portal (immPort) database (https://immport.niaid.nih.gov), we acquired 2,498 immune genes (IGs) (21). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Identification of differentially expressed immune-related RBPs

Using Wilcoxon’s test, the R package “limma” was used to identify the differentially expressed RBPs and IGs between the LUAD and normal samples. The adjusted P value <0.05 and |log2 fold-change| >0.5 were considered as the thresholds. Based on their expression in LUAD samples, Pearson correlation analysis was used to determine the RBPs related to IGs. The threshold was set to P value <0.01 and |correlation coefficient| >0.5. Genes that met the criteria were displayed in a regulatory network using the Cytoscape software. Gene Ontology (GO) (http://geneontology.org) is an important bioinformatics tool to annotate and illustrate genes and their biological process (BP), cellular component (CC), and molecular function (MF) (22). The Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg) is a comprehensive database resource for the biological interpretation of genome sequences (23). The analysis of GO and KEGG pathway enrichment was performed in the R package “clusterProfiler” (24).

Screening and analysis of prognostic immune-related RBPs

Univariable Cox regression analyses were used to determine the immune-related RBPs that were associated with LUAD prognosis. Using the R package “glmnet”, Lasso regression analysis was used to exclude genes with high correlation, thus reducing the complexity of the risk signature (25). Multivariable Cox regression analysis was ultimately performed to evaluate the immune-related RBPs as independent predictive indicators for LUAD (P<0.05). The final immune-related RBPs were defined as risk genes for the construction of the risk signature. A series of online databases were selected for systematic analysis of the above risk genes. First, genetic alteration conditions of risk genes in LUAD were downloaded from cBioPortal (http://www.cbioportal.org) (26). Subsequently, mRNA expression levels of risk genes in LUAD and normal tissues were compared in the UALCAN (http://ualcan.path.uab.edu/) database (27). Finally, we conducted OS analysis using The Gene Expression Profiling Interactive Analysis 2 (GEPIA2) (http://gepia2.cancer-pku.cn) database (28).

Construction and verification of prognostic risk model

The risk score was derived from the following formula: Risk score = Exp1 * B1 + Exp2 * B2 + Expi * Bi. Expi was the expression level of risk genes, and Bi was the multivariable Cox regression coefficient of risk genes. According to the median risk score, LUAD patients were divided into the low- and high-risk groups. The R package “survival” was used to compare the survival differences between the low- and the high-risk groups. External validation of the prognostic model was performed using the GEO data set GSE68465 to test stability. Group-based survival analyses were performed based on age, sex, stage, and TNM to determine the predictive power of the risk signature for different patients. In order to investigate whether the risk signature can be used as an independent prognostic factor, we integrated the risk score with clinical parameters for univariable and multivariable Cox regression analyses. The R package “survivalROC” was used to calculate the area under the curve (AUC) to verify the prognostic performance. A nomogram was drawn using the R package “rms”.

Assessment of tumor immune microenvironment, stemness index, and immune checkpoint

Understanding the tumor immune microenvironment is a prerequisite for effective tumor immunotherapy. Using the R package “estimate”, the tumor stromal cell and immune cell scores were calculated based on the ESTIMATE algorithm to analyze the purity of the tumor. A correlation analysis between the risk score and DNA stemness score (DNAss) and RNA stemness score (RNAss) was performed through the Spearman method. Tumor immune evasion is involved in immune checkpoint pathways, and anticancer immunity can be enhanced by immune checkpoint inhibitors (29). Thus, we analyzed the expression differences of some common immune checkpoints between low- and high-risk groups.

Immune infiltrates analysis

The tumor-infiltrated immune cells are an important component of the tumor immune microenvironment. The tumor immune estimation resource (TIMER) database (https://cistrome.shinyapps.io/timer/) can comprehensively analyze the infiltrating immune cells in tumors (30). According to the data of immune infiltration in LUAD patients in TIMER, we analyzed the correlation between risk score and infiltrating immune cells. Based on CIBERSORT (http://cibersort.stanford.edu/), the deconvolution algorithm was employed to evaluate the proportion of 22 immune cell subtypes in LUAD patients (31) to compare the differences in immune cell subtypes between the low- and high-risk groups, as well as the association between immune cells.

Principal component analysis (PCA) and gene set enrichment analysis (GSEA)

PCA was used to show the distribution patterns of different groups. The 3-dimensional graph was plotted by the R package “pca3d”. GSEA was used to assess whether there were differences in immune-related pathways between low- and high-risk patients. GSEA was used to analyze the KEGG enrichment pathways in the low- and high-risk groups. Each analyzed genome was repeated 1,000 times (32).

Drug sensitivity

Finally, we wanted to explore which drugs were sensitive to the immune-related RBPs constructed in the risk signature. Drug sensitivity processed data was obtained from the CellMiner database (https://discover.nci.nih.gov/cellminer/home.do). According to the analysis results, the top 16 sensitive drugs were selected for display using the R package “ggplot2”.

Statistical analysis

Unless otherwise specified, P<0.05 was considered statistically significant. All analyses were performed using R software (version 4.0.4), as described above.

Results

Differentially expressed genes

First, the differentially expressed RBPs and IGs between LUAG and normal samples were identified. Based on the threshold, 245 up-regulated and 139 down-regulated differentially expressed RBPs were identified (). In addition, 484 up-regulated and 337 down-regulated IGs were identified ().

Figure 1

Differentially expressed genes in LUAD patients. (A) Heatmap and volcano plot of differentially expressed RBPs in LUAD patients. (B) Heatmap and volcano plot of differentially expressed IGs in LUAD patients. (C) Regulatory network and functional enrichment analysis. The regulatory network of RBPs and IGs. GO (D) and KEGG (E) enrichment analyses of the immune-related RBPs. (F) Univariable Cox regression analyses. (G) Multivariable Cox regression analysis. LUAD, lung adenocarcinoma; RBP, RNA-binding protein; IG, immune gene; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; FC, fold change; ncRNA, non-coding RNA; UTR, untranslated region; BP, biological process; CC, cellular component; MF, molecular function; FDR, false discovery rate.

Co-expression regulatory network and enrichment analysis

Then, the role of the RBPs in the immune system was examined. A total of 278 pairs of interacting genes were identified through Pearson correlation analysis, among which the number of immune-related RBPs was 91 (90 significantly up-regulated genes, and one significantly down-regulated gene). A co-expression regulatory network was constructed to show the regulatory relationship between differentially expressed RBPs and IGs by Cytoscape (). GO enrichment () and KEGG pathway () analysis of immune-related RBPs were performed to understand their roles better. GO enrichment analysis indicated that the BPs were mainly associated with the non-coding RNA (ncRNA) metabolic process and RNA splicing. The CCs were abundant in the external side of the plasma membrane and spliceosomal complex. The MFs were most significantly enriched in catalytic activity, acting on RNA and nuclease activity (Table S2). Moreover, the KEGG pathways analysis showed that the RBPs were mainly enriched in cytokine-cytokine receptor interaction and viral protein interaction with cytokine and cytokine receptor (Table S3).

Prognostic immune-related RBPs

Univariable Cox regression was used to analyze genes associated with OS. From 91 immune-related RBPs, we identified 15 RBPs with prognostic value ( and Table S4). After lasso Cox regression and multivariable Cox regression analysis, we finally got four independent prognostic RBPs ( and Table S5). Then the risk genes were validated from multiple angles. First, gene mutations of the risk genes were detected in the cBioPortal database. Risk genes expression patterns were altered in 23% of 586 LUAD patients. EXO1 had the highest mutation rate of 10%. The most common forms of genetic alterations were amplification and missense mutations (Figure S1A). Next, we verified the differences in the expression of risk genes between tumor and normal tissues. Based on UALCAN database, the mRNA expressions of OAS3 and EXO1 in tumor tissues were increased while the mRNA expressions of PCF11 and TLR7 were opposite (Figure S1B). Based on the GEPIA2 database, we compared the prognosis of high and low expression of risk genes in LUAD patients (Figure S1C). We can see that the expression of risk genes and the results of their survival analysis were consistent with the effect of the signs of these coefficients. This validated our previous analysis.

Validation of the immune-RBP-based prognostic risk model for LUAD

According to the expression of the risk genes and regression coefficients, the formula for constructing a risk signature was: risk score = (OAS3*0.185096492873965) + (PCF11*−0.387198865656384) + (TLR7*−0.15166316298057) + (EXO1*0.214581605857755). The training group consisted of 486 LUAD patients from TCGA, then divided into low- and high-risk groups based on the median risk score. The median value of TCGA was between 0.987 and 0.935. The median value of GEO was between 0.986 and 0.988. Then, the risk score range of TCGA patients was 0.303 to 3.045, and the risk score range of GEO patients was 0.248 to 3.306. Compared with the high-risk group, OS was significantly higher in the low-risk group (). The AUC of the training group model was 0.731 (). We analyzed the survival distribution among patients in the low- and high-risk groups and observed that mortality increased with the increase in risk score (). Finally, a heat map showed the expression changes of risk genes in different groups ().

Figure 2

Construction of the prognostic model in the TCGA cohort based on risk scores. (A) Overall survival of the low- and high-risk groups. (B) ROC curve of the TCGA cohort. (C) The risk score distribution and survival status of LUAD patients. (D) Heat map of risk genes expression. AUC, area under the curve; TCGA, The Cancer Genome Atlas; ROC, receiver operating characteristic; LUAD, lung adenocarcinoma.

Independent external verification was performed using GSE68465 from the GEO database, which contains a total of 443 LUAD patients. The above formula was used again to obtain a risk score for each patient. The low-risk group also showed a better prognosis (). The AUC of the test group model was 0.707 (). Besides, the survival overview and gene expression heat map were shown in . Therefore, the above validation results proved the stability of the model prediction and its good predictive ability of LUAD patients’ OS. Finally, we further analyzed the relationship between the risk signature and the grouping of patients with different clinical parameters. Except for metastasis, all the others showed significant differences ().

Figure 3

Validation of the prognostic model in the GSE68465 cohort based on risk scores. (A) Overall survival of the low- and high-risk groups. (B) ROC curve of the GSE68465 cohort. (C) The risk score distribution and survival status of LUAD patients. (D) Heat map of risk genes expression. (E) The association between risk score and clinical parameters (age, gender, stage, T, N, and M). AUC, area under the curve; ROC, receiver operating characteristic; LUAD, lung adenocarcinoma.

Construction of a nomogram for clinical practice

The prognostic value of the risk model was further evaluated using hazard ratios (HRs) of the risk scores in univariable and multivariable Cox regression analyses (, and Table S6). The HRs were 2.11 and 1.86 in the univariable and multivariable Cox regression analyses, respectively (P<0.001). This means that the risk signature can be used as an independent prognostic factor, and the higher the risk score, the worse the prognosis. ROC curves were used to evaluate the specificity and sensitivity of our risk signature, and the AUC was 0.733, suggesting the risk signature for LUAD was highly reliable (). We used four immune-related RBPs to construct a nomogram for predicting prognosis, which can easily calculate the OS of LUAD patients at 1, 3, and 5 years ().

Figure 4

Independent prognostic analysis and nomogram for prediction. Univariable (A) and multivariable (B) cox regression analysis. (C) Comparison of risk signature and clinical parameters in predicting the prognosis of LUAD. (D) Development of a nomogram for predicting 1-, 3-, and 5-year OS of LUAD patients. AUC, area under the curve; LUAD, lung adenocarcinoma; OS, overall survival.

Tumor immune microenvironment, stemness index, and immune checkpoint

We performed risk signature analysis on tumor immune microenvironment, Stemness index, and immune checkpoints. The result showed that stromal cell score (R=−0.29, P<0.001) and immune cell score (R=−0.35, P<0.001) were negatively correlated with a risk score. Moreover, the risk score was positively correlated with DNAss (R=0.21, P<0.001) and RNAss (R=0.55, P<0.001) (). Cancer stemness represents the similarity of cancer cells with stem cells. Immune checkpoints are proteins involved in tumor immune escape. Among traditional immune checkpoints such as PD-L1 (Figure S2A), PD-1 (Figure S2B), and CTLA4, only CTLA4 (P=0.046, ) showed significant differences. Still, other potential emerging immune checkpoints under investigation all showed significant differences. For instance, LAG3 (P=0.0044), TIGIT (P=0.04), and TIM3 (P=0.0015) ().

Figure 5

Analysis of risk score with the tumor immune microenvironment, stemness index, and immune checkpoint in LUAD. (A) Based on the ESTIMATE algorithm, correlation analysis risk score with stromal cell scores, immune cell scores, DNAss, and RNAss. (B) R means correlation value. The immune checkpoint gene expression, including CTLA4, LAG3, TIGIT, and TIM3, in the low- and high-risk groups. (C) Immune cell infiltration analysis. Correlations between risk signature and immune cell infiltration levels, including B cell, CD4+ T cell, CD8+ T cell, dendritic cells, neutrophil, macrophage. LUAD, lung adenocarcinoma; DNAss, DNA stemness score; RNAss, RNA stemness score.

Immune infiltration and immune cell subtypes

Based on the data downloaded from the TIMER database, the Spearman correlation test was performed with the R software to obtain the correlation between the risk signature and the level of immune cell infiltration in LUAD. All immune cells were negatively correlated with risk score: B cell (Cor =−0.357, P<0.001), CD4+ T cell (Cor =−0.392, P<0.001), CD8+ T cell (Cor =−0.195, P=0.004), dendritic cells (Cor =−0.380, P<0.001), neutrophil (Cor =−0.284, P<0.001), and macrophage (Cor =−0.303, P<0.001) (). Using CIBERSORT, we further assessed whether there were differences in the abundance of 22 different immune cell infiltrates based on the risk signature. In the relatively high-risk group, B cell memory macrophages M2, resting dendritic cells, resting mast cells, and eosinophils were highly expressed in the low-risk group, while natural killer cells resting, Macrophages M0 and dendritic cells activated were the opposite (P<0.05; ). Heatmap showed the interaction between different cells based on Pearson correlation analysis ().

Figure 6

Distribution pattern and gene enrichment of low- and high-risk groups. (A) Composition of 22 immune cells in low- and high-risk groups. (B) Correlation heat map of 22 immune cells in LUAD. (C) PCA for low- and high-risk groups based on risk genes, whole RBPs, whole IGs, and whole genome. (D) GSEA suggested differences in immune response and immune system process. (E) GSEA suggested differences in the KEGG pathway between the low- and high-risk groups. LUAD, lung adenocarcinoma; PCA, Principal components analysis; RBP, RNA-binding protein; IG, immune gene; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Comparison of patients in the high- and low-risk groups

We used PCA to show the different distribution patterns of the high- and low-risk groups. Based on our immune-related RBP, the low- and high-risk groups were clearly distinguished from each other (). According to other standards, such as whole RBPs, whole IGs, and whole-genome sets (), all fail to distinguish between low- and high-risk groups of patients better. Normalized enrichment score (NES) is the standardized result, with positive values indicating that the pathway is associated with the first phenotype group and negative values indicating that the pathway is associated with the second phenotype. As shown in , immune response [NES =−1.69, false discovery rate (FDR) =0.049] and immune system process (NES =−1.71, FDR =0.044) were more annotated in the low-risk group. These results are consistent with our previous immune-related analyses, indicating that the low-risk scores were correlated with the enhanced immunophenotype. We further performed GSEA on the gene sets of patients in the low- and high-risk groups to identify the KEGG pathways involved. The top five enrichment pathways were selected to display, respectively. The high-risk group was mainly enriched in cell cycle, DNA replication, pathways in cancer, p53 signaling pathway, and SCLC. The low-risk group was mainly enriched in allograft rejection, autoimmune thyroid, B cell receptor signaling, jak stat signaling, and natural killer cell-mediated cytotoxicity (). Different enrichment pathways of patients in low- and high-risk groups may indicate different prognoses and appropriate treatment methods.

In order to explore the potential correlation between immune-related RBPs in the signature and drug sensitivity, we used human cancer cell lines from the CellMiner database for correlation analysis. TLR7 was positively correlated with alectinib, denileukin diftitox (ontak), fluphenazine, isotretinoin, LDK-378, imiquimod, megestrol acetate, nelfinavir, estramustine, celecoxib, dromostanolone propionate, and negatively correlated with irofulven. EXO1 was positively correlated with vorinostat and 6-thioguanine. PCF11 was positively correlated with nelarabine and vorinostat (). All these drugs modulate the tumor immune microenvironment to induce cancer cell death (33-46).

Figure 7

Drug sensitive to risk genes. (A) Alectinib. (B) Denileukin diftitox (ontak). (C) Irofulven. (D) Fluphenazine. (E) Isotretinoin. (F) LDK-378. (G) Imiquimod. (H) Vorinostat. (I) Megestrol acetate. (J) Nelarabine. (K) Nelfinavir. (L) Estramustine. (M) 6-thioguanine. (N) Vorinostat. (O) Celecoxib. (P) Dromostanolone propionate.

Discussion

RBPs play critical roles in cancer. Immunotherapy is effective in many cancer types, LUAD. This study aimed to construct and validate a risk signature with predictive power based on immune-related RBPs for LUAD patients. The results suggest that the risk signature of immune-related RBPs established in this study can provide the basis for clinical decisions regarding diagnosis, prognosis, and immunotherapy in LUAD patients.

The expression level of oncogenes and tumor suppressor genes can be affected by changes in RBPs expression and location. In addition, various cancer-associated phenotypes such as metastasis, invasion, aging, angiogenesis, apoptosis, and proliferation can be regulated by RBP-mediated genes. Therefore, there must be some relationship between cancer prognosis and RBP (47). Researchers are committed to developing anti-cancer drugs that target RBP, which is expected to become a new direction for cancer treatment in the future (48). Still, the biological roles and mechanisms of most RBPs have not been fully elucidated in LUAD. RBPs are also important regulators of the immune response. More and more scholars have focused on the study of RBP in tumor immunotherapy (49). For instance, it has been reported recently that cold-inducible RBP can enhance immunotherapeutic responses against hepatocellular carcinoma (50). The specific deletion of RBP PCBP1 in T cells facilitates the differentiation of Treg cells and generates a variety of inhibitory immune checkpoints, including VISTA, TIGIT, and PD-1 on tumor-infiltrating lymphocytes, thereby inactivating anti-tumor immunity (51). Nevertheless, many immunomodulatory RBPs have not been identified in LUAD. Consequently, we aimed to investigate the immune-related RBPs in LUAD diagnosis, prognosis, and immunity.

Based on the TCGA and immPort databases, we identified 384 RBPs and 488 IGs that were differentially expressed in LUAD and normal tissues. Through the correlation analysis between IGs and RBPs, 91 immune-related RBPs in LUAD were identified. GO enrichment analysis revealed the immune-related RBPs were primarily enriched in ncRNA metabolic process, external side of the plasma membrane, catalytic activity, and acting on RNA. KEGG enrichment analysis showed that the immune-related RBPs were mainly enriched in cytokine-cytokine receptor interaction.

We eventually identified four independent prognostic immune-related RBPs (OAS3, PCF11, TLR7, and EXO1) as risk genes to construct the risk signature. The expression levels of OAS3 and EXO1 were up-regulated in LUAD, while PCF11 and TLR7 were down-regulated in LUAD. Simultaneously, higher expression of OAS3 and EXO1 and lower expression of PCF11 and TLR7 in LUAD patients were associated with worse OS. These results suggest that they may become potential biomarkers for the diagnosis and prognosis of LUAD. OAS3 (2'-5'-oligoadenylate synthetase 3) has an antiviral effect on alphavirus (52). Recent study has identified several molecules as inhibitors of the OAS family, especially interacting with Tyr230, Asp75, Asp77, and Asp75 in OAS3, indicating that OAS3 could be expected to become a new drug target (53). PCF11 (PCF11 cleavage and polyadenylation factor subunit) directs neural differentiation, apoptosis, proliferation, and cell cycle by shaping the inputs that converge on WNT signaling. A recent study reported that down-regulation of postnatal PCF11 promotes neurodifferentiation, and low levels of PCF11 were associated with spontaneous tumor regression and good prognosis (54). Toll-like receptor (TLR) ligand’s effect on immune activation makes it a drug target for treating cancer and infectious diseases (55). Because TLRs bind to cell membranes and present a wide range of expression patterns. Even local application may cause the release of systemic cytokines, resulting in inhibitory side effects (56). In contrast, TLR7 has a more limited range of action, including only certain types of immune cells, and can serve as a treatment target (57). Preclinical studies have shown that TLR7 agonists have antitumor and immunomodulatory effects. For example, studies have demonstrated that the TLR7 agonist R-848 can enhance anti-tumor effects in cetuximab-coated colorectal cancer (CRC) cell lines through antibody-dependent cellular cytotoxicity (ADCC) (58). EXO1 is an evolutionarily highly conserved exonuclease (59). Compared with the above genes, EXO1 had the highest mutation rate. As previously reported, abnormal expression of EXO1 may cause genomic instability, including homologous recombination, telomere maintenance, and DNA mismatch repair (60-62). Some studies claimed that increased expression of EXO1 promotes tumor invasion and metastasis, leading to a worse prognosis for cancer patients (63,64). In terms of immune infiltrate, overexpression of EXO1 might lead to decreased antibody diversity and immunoglobulin maturation malfunction, lowering the number of infiltrating B cells even further (62).

Through gradual screening, a risk signature was finally established based on the above four immune-related RBPs. According to the median risk score, LUAD patients were divided into the high- and low-risk groups. Low-risk patients had better OS than high-risk patients. It was validated in another independent cohort from the GEO database. In order to explore the feasibility of the prognostic model in clinical application, we compared this prognostic model with the clinical parameters of LUAD patients, such as age, gender, and pathological stage. After univariable and multivariable Cox and ROC analyses, the risk signature was an independent prognostic factor of LUAD. To enable our prognostic model to gain more reliable and valuable predictive power in clinical applications, we developed a nomogram that included OAS3, PCF11, TLR7, and EXO1 to calculate the scores of different patients to predict the prognosis.

Recent studies have associated the immune microenvironment with the immunotherapeutic response of several cancers, including LUAD (65-67). By calculating immune scores, we could accurately determine the tumor purity and immune cell infiltration in the tumor microenvironment. DNAss and RNAss reflect epigenetic characteristics and gene expression, respectively (68). They were found to be positively correlated with the risk score in our analysis. An increase in stemness index indicates that the BPs of cancer stem cells were active, which might cause the dedifferentiation of differentiated tumors. Interestingly, a study found that elevated DNAss in some tumors led to decreased PD-L1 expression and leukocyte fraction, including BLCA, HNSC, LUSC, and GBM. Because of inherent down-regulation of the PD-L1 pathway or insufficient immune cell infiltration of tumors, such tumors were expected to be less susceptible to immune checkpoint blockade (68), supporting the present study.

With the successful application of immunotherapy, immune checkpoints have become the new focus of molecular oncologists. Although anti-cytotoxic T lymphocyte-associated protein 4 (anti-CTLA4) and anti-PD-1 therapies have achieved initial success, limitations were also exposed, such as the emergence of drug resistance in most patients. Other immune checkpoints that have the potential to become targets for immunotherapy have aroused the interest of researchers and were under study, such as TIM3, LAG3, and TIGIT (69). Coincidentally, our results suggested significant differences in immune checkpoints between the high- and low-risk groups, especially emerging immune checkpoints.

Based on the CellMiner database, we explored drugs with higher sensitivity to immune-related RBPs. Drug sensitivity was proportional to the amount of gene expression. Among them, TLR7 showed a close correlation with various drugs, and recent study has shown a link. For example, studies indicated that imiquimod, as a TLR7 agonist, promotes the immunogenicity of mesenchymal stem cells (70). In addition, EXO1 was positively correlated with vorinostat and 6-thioguanine. PCF11 was positively correlated with nelarabine and vorinostat. Of course, we only selected the top 16 sensitive drugs, and there are still many sensitive drugs worth exploring. Therefore, detecting the expression of TLR7, EXO1, and PCF11 has a certain guiding significance for selecting clinical drugs.

As far as we know, this is the first study to establish a risk signature based on immune-related RBPs in LUAD. This risk signature should be examined in the future for guiding immunotherapy in patients with LUAD. Nevertheless, several limitations of this study should also be noted. First, patient information was taken from the TCGA and GEO databases, and more sufficient clinical samples were needed for verification. Second, transcriptome analysis can only reflect certain aspects of the immune status and lacks overall changes. Third, this study was a retrospective study, and prospective studies were needed to confirm our conclusions.

In conclusion, this study identified differentially expressed and prognostic immune-related RBPs and used them to construct a risk signature of LUAD. Low-risk LUAD patients showed a better prognosis and better immunotherapy responses than high-risk patients. This study provided several biomarkers and a potential model for the prognosis and individualized immunotherapy of LUAD patients.

The article’s supplementary files as