A novel immune-related RNA-binding proteins signature to predict survival and therapeutic responses in prostate cancer.

BACKGROUND: The present study aims to identify immune-related RBPs signature to predict prognosis and therapy response in prostate cancer.
METHODS: Differentially expressed RBPs were compared and visualized using R packages. Immune-related RBPs were selected by Pearson correlation analysis. The prognostic immune-related RBPs were identified using the Kaplan-Meier method and LASSO regression. A multivariable Cox regression model was used to construct immune-related RBPs signature.
RESULTS: We constructed a prognostic predictive risk model of prostate cancer containing ten immune-related RBP genes. We found that high-risk prostate cancer patients presented poorer prognosis, higher tumor immune cell infiltration, higher rates of genomic alterations, and were more sensitive to targeted and immunotherapy than the low-risk group.
CONCLUSIONS: The immune-related RBPs' signature is an independent prognostic marker that could help screen patients with advanced prostate cancer who are better suited for targeted and immunotherapy.

Background

Prostate cancer is the second most common cancer and the main leading cause of cancer-related death in humans, with almost 1.3 million new cases and 359,000 deaths in 2018[1]. Due to the widespread use of prostate-specific antigen (PSA) testing and biopsy in the past 30 years, prostate cancer incidence rates have increased rapidly. PSA-based screening indeed decreased mortality. When newly diagnosed with prostate cancer, patients first want to know the risk of progression to fatal disease and the benefits of surgical treatment. Urologists perform risk stratification based on stage, Gleason score (GS), serum PSA, imaging, and sometimes genomic profiling. Therefore, it is vital to explore new biomarkers for both the early detection and prognosis of prostate cancer.f

RNA-binding proteins (RBPs) play a crucial role in posttranscriptional events. From the conventional view, RBPs control gene expression in the form of ribonucleoproteins (RNPs). They interact with other proteins or RNAs, regulating RNA polyadenylation, stability, translation, splicing, and degradation [2]. Many RBPs bind to their target RNA through specific RNA-binding domains (RBDs). Different RBDs assemble like puzzle pieces and give RBPs the conformational flexibility and adaptability properties. Scientists have presented a census of 1542 human RBPs and detailed their essential roles in developmental processes [3]. The loss-of-function of RBPs causes specific disease phenotypes. For example, ribosomal protein defects may cause Shwachman–Diamond syndrome [4]. Many studies have shown that RBP dysregulation may be related to cancer genesis or progression by impacting the expression of tumor suppressor proteins or oncoproteins. Some studies have compared RBP expression between cancer tissue and adjacent normal tissue and found some clinically relevant RBPs. RBP regulation mechanisms such as genomic alterations [5], posttranscriptional regulation [6], and posttranslational modifications (PTMs) have been revealed [7]. However, there are still many questions concerning RBPs in tumor biology that need to be explained. Exploring RBPs and their target interaction networks will undoubtedly help understand cancer biology better and provide new targets for tumor treatment.

The relative genes and molecular mechanisms of prostate cancer initiation, progression, and metastasis are still being studied globally. Identifying innovative functional genes or pathways can help us identify new biomarkers for early diagnosis and therapy targets. High-throughput microarray technology and various online datasets facilitate the screening of potentially differentially expressed genes and key players or hub genes controlling the deep interaction network. In the present study, we downloaded the RNA sequencing (RNA-Seq) data of all the RBPs and clinical information of patients in The Cancer Genome Atlas (TCGA)-prostate adenocarcinoma (PRAD) dataset. Differentially expressed RBPs were compared. Gene Ontology (GO), gene set enrichment analysis (GSEA), and prognostic analysis were then performed for downstream analysis. Our study reported the key RBPs of prostate cancer using integrated bioinformatics analysis and may help discover potential diagnostic markers or prognostic indicators for prostate cancer.

Methods

Data collection

We downloaded the RNA-Seq data of all the RBPs and clinical information of 495 patients in the TCGA-PRAD dataset (https://portal.gdc.cancer.gov/) and MSKCC dataset (https://www.mskcc.org/) as validation prostate cancer cohort. Gene alteration data were visualized using R “maftoosl” package. A list of RBP genes was obtained based on previous literature [3, 7, 8]. A list of immune-related genes was downloaded from the Immunology Database and Analysis Portal (ImmPort) database (https://www.immport.org) [9]. The Human Protein Atlas (HPA) (https://www.proteinatlas.org/) database was used to validate the protein expression of ten hub RBPs in normal prostate tissues and prostate cancer tissues.

Identification of differentially expressed RBPs and immune-related RBPs

The “edgeR” R package was used to compare mRNA expression differences between 55 normal prostate tissues and 495 tumor tissues. All the differentially expressed RBPs were visualized by the “pheatmap” R package. All immune-related RBPs genes were screened according to the Pearson correlation between differentially expressed RBP genes and immune-related genes in 495 samples with the p-value less than 0.01 and |correlatin coefficient| more than 0.5.

Visualization and enrichment analyses

To explore the mechanisms that may exist downstream of mining, we performed a functional enrichment analysis of all differentially expressed RBPs. We applied the “clusterProfiler” R package to implement and visualize the downstream analyses of the differentially expressed RBPs, including GO and GSEA.

Kaplan–Meier estimator and Cox regression analyses

Before performing cox survival analysis, we screened RBPs with Lasso's method to filter out those with prognostic significance. The Kaplan–Meier method was used to evaluate the effects of RBPs on progression-free survival (PFS) in prostate cancer. Then, univariable and multivariable Cox survival regression analyses were carried out to identify significantly valuable RBPs. We set a p-value of less than 0.0005 to indicate significance in the Kaplan–Meier estimator analysis.

Tumor microenvironment and infiltrated immune cell evaluation

We applied the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) tool to calculate the presence of stromal/immune cells infiltration and prostate cancer purity based on transcriptome data [10]. Three principal scores were generated by ESTIMATE algorithm based on the single sample Gene Set Enrichment Analysis (ssGSEA), including stromal score that captures the presence of stromal cells in tumor tissue, immune score that represents the infiltration of immune cells in tumor tissue, and estimate score that infers tumor purity. Then, we calculate the scores of 28 different immune cell types by ssGSEA method [11]. Panels of gene markers for each immune cell type were acquired from the previous publication [12]. Each ssGSEA enrichment score represents the degree to which the genes in a particular immune cell type are coordinately up- or down-regulated within a sample.

Target therapy and immunotherapy response prediction

To further elucidate the efficacy of target therapy drugs in high- and low-risk patients in prostate cancer, the “pRRophetic” R package was used to calculate IC50 of all 138 types of target therapy drug in prostate cancer patients according to the Genomics of Drug Sensitivity in Cancer database. The Cancer Immunome Atlas (TCIA, https://tcia.at/) was applied to evaluate the benefits of response to anti-PD1 and anti-CTLA4 immunotherapy in prostate cancer.

Statistics analyses

All data analyses were performed in R v3.6 platform. Wilcoxon rank-sum tests were applied to compare the mRNA expression level of RBPs between normal and tumor tissue. The hub genes of RBPs were filtered by Lasso and Cox regression method. The risk score for prostate cancer patients was calculated by multiplying each RBP expression in the Cox regression model by the corresponding coefficient and then summing cumulatively. Differences with p-value less than 0.05 were considered statistically significant.

Results

Differentially expressed RBPs between tumor and normal prostate tissues

We compared all the mRNA expression levels of the RBPs in normal prostate and tumor tissues in the TCGA-PRAD dataset and found a total of 110 differentially expressed RBPs with a log fold change greater than 1 (Fig. 1A, C). At the same time, 13 RBPs were identified as differentially expressed with a log fold change greater than 2 (Supplementary Fig. 1A). Most of these differentially expressed RBPs were down-regulated in prostate cancer compared with normal tissue. However, in prostate cancer patients with different GSs, the RBPs did not show different expression levels (Supplementary Fig. 1B).

Fig. 1

Differentially expressed RBPs and enrichment analysis in Prostate cancer. A Heatmap of differentially expressed RBPs genes. The rows include 110 differentially expressed RBPs with a log fold change greater than 1. The columns include normal and cancer tissue samples. The colors indicate RBPs’ expression, which increased from a relatively low level (blue) to a relatively high level (red). B GO analysis of the differentially expressed RBPs. C Volcano plot of differentially expressed RBPs genes. RBPs, RNA-binding proteins. D GSEA analysis of the 110 proteins. GO gene ontology, GSEA gene set enrichment analysis

To elucidate the biological functions of the differentially expressed RBPs between normal and tumor tissues, enrichment analyses were carried out. The results showed that the differentially expressed RBPs could participate in some vital biological processes, such as mRNA processing, RNA splicing, and DNA modification (Fig. 1B). The GSEA results revealed that MYC target genes, Apical Junction, and Epithelial–Mesenchymal transition pathway played essential roles in the regulatory mechanisms of RBPs in prostate cancer (Fig. 1D).

Identification and construction immune-related RBPs’ signature in prostate cancer

From 110 RBP genes that were differentially expressed in prostate cancer tumor tissues and normal tissues, we screened 71 immune-related RBP genes based on gene correlation analysis. Then, LASSO regression analysis was applied to screen 15 RPB genes (Fig. 2 A, B). Then, using univariable and multivariable Cox survival analyses, we identified ten immune-related RBPs to construct the model (Fig. 2C, Supplementary Table 1). The Concordance index of ten immune-related RBP genes model is 0.75, which is relatively high. Concretely, the risk scores for the ten immune-related RBP gene models were calculated as follows: risk score = (− 0.159 × EXPCSRP1) + (0.403 × EXPTFRC) + (0.417 × EXPZC3HAV1L) + (− 0.712 * EXPOGN) + (− 1.06 × EXPTPT1) + (0.055 × EXPRNASE2) + (0.252 × EXPDHX58) + (0.591 × EXPDDX17) + (0.1 × EXPTLR8) + (0.452 × EXPTMSB4X). To further assess the discrimination of the model and its ability to predict the prognosis of prostate cancer patients, we compared the prognostic differences between high- and low-risk patients. Survival analysis shows shorter PFS in high-risk prostate cancer patients compared to low-risk group (p < 0.001, Fig. 2D). The AUC of PFS at one year, two years, and three years was 0.846, 0.753, and 0.744, respectively (Fig. 2E). Moreover, compared with conventional biomarkers of prostate cancer, such as prostate-specific antigen (PSA) and androgen receptor (AR), the performance of RBPs risk model shows better predictive accuracy in MSKCC prostate cancer cohort (c-index 0.8, p < 0.001). In addition, we compared the expression differences of these ten immune-related RBP genes in high- and low-risk groups (Fig. 3). CSRP1, OGN, and TPT1 were significantly high- expressed in the low-risk group, while other RBPs were remarkably upregulated in the high-risk group (p < 0.001). To validate the differential expression of ten immune-related RBPs in protein level, we obtained and compared the immunohistochemistry figures of these hub RBPs from HPA database (Fig. 4). Furthermore, we found that these differential expressed RBPs between normal and tumor prostate tissues were consistent with our former RNA expression analysis.

Fig. 2

Construction and evaluation of an immune-related RBPs signature in prostate cancer. A The results of feature selection in LASSO regression. B The relationship between coefficients and log lambda in LASSO regression. C Multivariable Cox regression analysis of ten immune-related RBPs genes. D Kaplan–Meier curves of recurrence-free survival between high-risk and low-risk patients by risk score in TCGA. E Time-dependent ROC curves in TCGA. F Concordance index curves compared RBPs model with AR and PSA model in MSKCC cohort. RBPs RNA-binding proteins, LASSO least absolute shrinkage, and selection operator, TCGA The Cancer Genome Atlas, MSKCC Memorial Sloan Kettering Cancer Center

Fig. 3

The expression levels of ten immune-related RBPs in the high- and low-risk groups in the TCGA cohort. A CSRP1. B TFRC. C ZC3HAV1L. D OGN. E TPT1. F RNASE2. G DHX58. H DDX17. I TLR8. J TMSB4X. RBPs RNA-binding proteins, TCGA The Cancer Genome Atlas

Fig. 4

The immunohistochemistry figures of ten hub RBPs from HPA database

Establishment and validation of the predictive nomogram

To improve the prognostic power of patients with prostate cancer by combining their risk scores with other important clinicopathological features, we first compared the AUC of risk scores with other clinicopathological features and found that the AUC values of risk scores for the ten immune-related RBP gene prediction models (AUC = 0.846) were much higher than those of other clinicopathological markers (Fig. 5A). Then, we constructed a nomogram of PFS including risk scores and confirmed that the model has good stability by consistency analysis (Fig. 5B, C). To further validate the stability of our constructed model, we compared the ten immune-related RBP gene prediction models with four other published prediction models based on the TCGA prostate cancer database and found that our model had the best 1-year AUC to performance, which fully illustrated the better predictive ability of our signature (Supplementary Table 2).

Fig. 5

Construction of predictive nomogram combined the immune-related RBP signature with clinicopathologic features of prostate cancer. A ROC curves of the immune-related RBP signature and clinicopathological features. B Nomogram validation. C Predictive nomogram of RBP model risk score with critical clinicopathologic characteristics in prostate cancer

TME characteristics in different immune-related RBPs signatures

We found that patients in the high-risk group had slightly higher immune scores (p < 0.05) compared to the low-risk group, and then, stromal scores and tumor purity did not differ significantly between the two groups (Fig. 6A). Furthermore, immune cell infiltrations were accessed in high- and low-risk groups to demonstrate the immune landscape (Fig. 6B). The high-risk group presented higher abundance of activated CD4 T cell, activated dendritic cell, central memory CD8 T cell, gamma delta T cell, immature B cell, memory B cell, regulatory T cell, and type 2 T helper cell. In contrast, the low-risk group had a more infiltration of CD56bright natural killer cell, effector memory CD4 T cell, immature dendritic cell, mast cell, and type 17 T helper cell (Fig. 6C).

Fig. 6

Tumor microenvironment in high- and low-risk groups. A Comparison of the immune score, stromal score, and ESTIMATE score between high-risk and low-risk groups. B The immune landscapes C immune cell infiltration

Enrichment analysis and mutation landscape in different immune-related RBPs signatures

The GSEA results showed that the top six pathways of androgen response, cholesterol homeostasis, fatty acid metabolism, IL2-STAT5 signaling, peroxisome, and α-TNF signaling were highly enriched in the high-risk group (Fig. 7A). Moreover, we describe the distribution of the top 30 genes with mutation frequency in the high- and low-risk groups (Fig. 7B, C). The frequency of mutations in prostate cancer was significantly higher in the high-risk group than in the low-risk group (80.66% vs. 46.98%). In the high- and low-risk groups, the top three genes with the highest mutation frequencies were SPOP (15% vs. 5%), TP53 (15% vs. 6%), and TTN (13% vs. 4%). Notably, the frequency of mutations in FOXA1, a critical transcription factor in prostate cancer, was significantly higher in the high-risk group (7% vs. 2%).

Fig. 7

Gene set enrichment analysis and genomic mutations analysis. A Gene set enrichment analysis results in high-risk and low-risk groups. The top 30 most frequently mutated genes in high-risk (B) and low-risk (C) groups

Drug sensitivities and immunotherapy responses of different immune-related RBPs signatures

To better explore the implications of our constructed risk model for guiding prostate cancer treatment, we evaluated the differences in the efficacy of various target therapy agents and immunotherapy in high- and low-risk groups. The results showed that 21 targeted therapies were more effective in treating patients with high-risk prostate cancer (Fig. 8, Supplementary Table 3). Interestingly, we found that conventional chemotherapy regimens based on cisplatin or docetaxel are more suitable for low-risk patients (Fig. 9A, B). Furthermore, the immunophenoscore (IPS) of PD1- and CTLA4-based immunotherapy was higher in the low-risk group, suggesting that the risk model helps screen prostate cancer patients more suitable for immunotherapy and improves the efficacy of comprehensive treatment for advanced prostate cancer (Fig. 9C, D).

Fig. 8

Drug sensitivity comparison in high- and low- risk groups of 21 kinds target therapy drug in prostate cancer

Fig. 9

Chemotherapy and immunotherapy responses in high- and low-risk groups of prostate cancer. A Cisplatin. B Docetaxel. C. IPS of CTLA4 and PD1 negative. D IPS of CTLA4 negative and PD1 positive. IPS immunophenoscore

Discussion

The exact role of RBPs in the tumorigenesis and progression of prostate cancer is still unclear and needs further investigation. However, some RBPs are abnormally expressed in cancer tissues compared with normal prostate tissues, implying their potential role in tumorigenesis and progression. In the present study, we downloaded RNA-Seq profiling data from TCGA and found 110 differentially expressed RBPs that served as candidates for downstream data analysis. We performed GO analysis to acquire the biological functions of these differentially expressed RBPs.

First, regarding biological processes, these differential RBPs mainly participate in mRNA processing, RNA splicing, DNA modification, and DNA methylation or demethylation. Some of them are also related to base conversion or substitution editing and cytidine metabolic processes. Next, these RBPs are involved in cellular components, such as RNP granules or cytoplasmic RNP granules. RNP granules are responsible for the storage, transportation, and degradation of their target transcripts by physically separating or connecting mRNAs with proteins. Thus, they have a far-reaching influence on the regulation of oncogene or tumor suppressor gene expression [13]. Finally, regarding molecular functions, most of the RBPs are involved in mRNA binding, including 3´-UTR binding and AU-rich region binding. The GSEA results show that together, these prostate cancer-related differential RBPs play a vital part in the ribosome pathway. Some ribosomal proteins are related to prostate cancer, and those upregulated proteins may serve as potential biomarkers [14].

In this study, we constructed risk prediction models for ten RBP genes by machine learning methods and univariate and multifactor Cox regression. CSRP1, cysteine, and glycine-rich protein 1 are a specific biomarker for various tumors including prostate cancer, breast cancer, and colon cancer [15-17]. TFRC induces iron death in tumor cells in neuroblastoma through upregulation of iron ion transport by transferrin receptors on cell membranes, which regulates tumor cell growth [18]. YTHDF1 as an m6A reader can recognize the m6A methylation site in the mRNA of TFRC, which in turn stabilizes the mRNA level of TFRC and promotes the protein translation of TFRC, thus upregulating the ferritin content and promoting iron metabolism in tumor cells [19]. Meanwhile, TFRC has a role in predicting tumor prognosis and as a therapeutic target in bladder and lung cancers [20, 21]. ZC3HAV1L is also named Zinc finger CCCH-type antiviral protein 1-like. The function and role of ZC3HAV1L in tumors are unknown, and this study confirms for the first time the prognostic value of this gene in prostate cancer. OGN could induce bone formation by interacting with TGF and is highly associated with prognosis in multitype cancer [22-25]. The transcription factor NFIC can bind to the promoter region of OGN and promote the transcriptional activity of OGN, which in turn inhibits the NF-κB signaling pathway and suppresses the ability of proliferation in bladder cancer [26]. OGN expression was reduced in paired breast cancer samples compared to normal tissue. Downregulated OGN expression was associated with higher pathological grade, more aggressive tumor subtypes, and poor overall survival in breast cancer. Previous study has also found that OGN expression reverses EMT through inhibition of the PI3K/Akt/mTOR pathway [27]. The ability of OGN highly expressed in tumor tissues to increase the infiltration of T lymphocytes in the tumor microenvironment by inhibiting VEGF also suggests that a similar role may exist in prostate cancer [28].

In prostate cancer, although TPT1 can be used as a biomarker to distinguish CRPC from HSPC, the prediction model including TPT1 constructed by Jun et al., is generally effective in predicting prostate cancer prognosis and needs to be improved [29]. Moreover, the immune-related RBP prediction model constructed in our study including TPT1 is better. Like prostate cancer, RNASE2 was able to participate in the construction of genetic risk prediction models in various tumors, indicating a significant correlation between RNASE2 and tumor development and prognosis [30, 31]. DHX58, also named Probable ATP-dependent RNA helicase DHX58, regulates DDX58/RIG-I and IFIH1/MDA5-mediated antiviral signaling. DHX58 was mainly involved in innate immune responses to some RNA viruses and some DNA viruses such as SARS-CoV-2 coronavirus and the bacterial pathogen Listeria monocytogenes [32]. Meanwhile, DHX58 in dendritic cells is crucial for activating the body's immune anti-tumor response after tumor radiotherapy [33]. As an RNA helicase, DDX17 can deconvolute RNA by binding and hydrolyzing ATP to change the tertiary structure of RNA and is involved in various cellular processes such as immature mRNA splicing, ribosomal RNA and microRNA synthesis, and transcriptional regulation [34-37]. TLR8, Toll-like receptor 8, plays an important role in innate and acquired immunity as an intranuclear receptor. TLR8 induces the body's immune response by recognizing specific RNA degradation products of pathogenic microbial origin [38-40]. The above literature review shows that the ten immune-related RBP genes included in our constructed prediction model are indeed significant in tumor or prostate cancer.

We found that the high- and low-risk groups of prostate cancer patients distinguished by this risk prediction model had significantly different genomic characteristics and tumor immune microenvironment. The frequency of genetic mutations was significantly higher in the high-risk group than in the low-risk group. The type and proportion of immune cell infiltration were also different in the high-risk group than in the low-risk group. To better evaluate the predictive model for the treatment of advanced prostate cancer, we evaluated 138 different targeted therapies and ultimately identified 21 drugs that were more suitable for patients in the high-risk group and 72 drugs that were more suitable for patients in the low-risk group.

The limitation of our study is obvious since our results were based on pure bioinformatics analysis. Future molecular biological experiments need to be performed to provide more evidence for our conclusion.

Conclusions

The risk prediction model we constructed for ten immune-related RBP genes was able to significantly distinguish between the high-risk and low-risk groups of prostate cancer. We also found that patients in the high-risk group had a poor prognosis and were more sensitive to targeted therapy and immunotherapy.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 3062 KB)