Prognostic Value and Biological Functions of RNA Binding Proteins in Stomach Adenocarcinoma.

PURPOSE: To investigate the prognostic value and biological function of RNA binding proteins (RBPs) in stomach adenocarcinoma (STAD).
MATERIALS AND METHODS: Datasets of the differentially expressed genes (DEGs) were downloaded from the TCGA-based (The Cancer Genome Atlas) GEPIA database, from which the differentially expressed RBPs were determined. Functions and prognostic values of these determined RBPs were systematically investigated by a series of methods in bioinformatics analysis. In addition, transwell assays were performed to explore the effect of PTBP1 in STAD cells.
RESULTS: Three hundred and sixty-two differentially expressed RBPs were determined, with 331 up-regulated and 31 down-regulated. Seven RBPs (PTBP1, PPIH, SMAD5, MSI2, RBM15, MRPS17, and ADAT3) were identified to be prognosis-related and adopted to construct a prognostic model. Compared with low-risk patients in TCGA training cohort, TCGA testing cohort and GEO cohort, high-risk patients had poorer overall survival (OS). The area under the ROC curves of this prognostic model were 0.804, 0.644 and 0.581 for TCGA training cohort, TCGA testing cohort and GEO cohort, respectively, justifying itself as a good prognostic model with reliable predictive ability. Using the seven identified RBPs, we then constructed a nomogram to generate a clinical utility model. The regulatory networks and functions of the seven RBPs were then investigated, the results of which demonstrated that MRPS17 and PTBP1 reduced the number of infiltrated immune cells. In-vitro experiments showed that the downregulation of PTBP1 weakened the migration and invasion capability of AGS and HGC27 cells.
CONCLUSION: The seven-gene signature can be used as a reliable STAD prognostic biomarker, and these findings help us better understand the prognostic roles and functions of RBPs in STAD.

Introduction

Globally, stomach adenocarcinoma (STAD) is a common malignant tumor with high mortality.1 Survival of patients with early STAD has been remarkably improved due to progress made in diagnosis and treatment over the last few decades.2 Unfortunately, the prognosis of advanced STAD patients remains dismal due to the lack of effective treatment. Little progress has been made in improving patients, survival despite the application of chemotherapy, targeted therapy and immunotherapy.3 The role of immunotherapy still needs further verification.3 Thus, it is still an urgent task for scientists to elucidate molecular mechanisms of STAD that may help us design new therapeutic strategies and improve the prognosis of patients with advanced STAD.

RNA binding proteins (RBPs) are a class of proteins that interact with various types of RNA, such as mRNAs, tRNAs, snRNAs, snoRNAs and others4,5. In human cells, more than 1500 RBPs have been identified through genome-wide screening, accounting for more than 7.5% of protein-encoding genes.4,5 By regulating processing, translation, export, and localization of target RNAs, RBPs play crucial roles in post-transcriptional regulation.6 It has demonstrated that the abnormal expression of RBPs is associated with cancer and many other diseases.7 For example, as one of the most RBP implicated in tumorigenesis, HuR has been reported to be associated with the prognosis of many malignant tumors, such as lung cancer, gallbladder cancer, ovarian cancer, breast cancer and colon cancer.8 Nucleolin promotes the occurrence of breast cancer through activating ERBB2.9 FOXK2 promotes colorectal cancer metastasis by upregulating mRNA expression of ZEB1 and EGFR.10 ESPR1 promotes the mesenchymal-epithelial transition of ovarian cells.11 AGO2 facilitates cancer progression by increasing oncogenic miR-19b biogenesis.12 Taken together, these studies indicate that RBPs are closely related to the occurrence and development of human tumors. However, till today, the roles of RBPs in STAD have only been superficially investigated.6 Thus, a systematic investigation of RBPs, roles in STAD is warranted since it can help us better understand novel mechanisms of SATD development and improve the prognosis of STAD patients.

In the present study, data of STAD were downloaded from TCGA, from which differentially expressed RBPs were identified and the potential functions and mechanisms of these RBPs were explored. An RBPs-based prognostic model was also developed and validated as some RBPs could be used as potential prognostic biomarkers.

Materials and Methods

Data Processing

STAD RNA-sequencing and clinical dataset were downloaded from TCGA () and GEO datasets (GSE84437). GSE84437 dataset included gene expression data from 433 samples, which was obtained through the Illumina HumanHT-12 V3.0 Expression BeadChip platform. We used the ANOVA method to obtain the differently expressed genes (DEGs) from the GEPIA online database () in view of |log2fold change (FC)|>1 and p<0.01. We extracted 1542 RBPs gene names from Gerstberger’s article.5 Volcanic maps of DEGs and differently expressed RBPs were produced via an online tool ().

GO Enrichment and KEGG Pathway Analysis

Gene ontology (GO) enrichment and kyoto encyclopedia of genes and genomes (KEGG) pathway analysis were performed by utilizing an online WEB-based Gene Set Analysis Toolkit (WebGestalt, ).13

Construction and Validation of Prognostic Models

STAD database in TCGA was randomly divided into training and testing cohorts. Univariate Cox regression analysis was performed on all differently expressed RBPs in the TCGA training dataset using the survival R package. The significance of candidate genes was verified by Log rank test. Subsequently, a multivariate Cox proportional hazards regression model was constructed based on the above significant candidate genes, and the risk score (RS) was calculated to assess patient prognosis outcomes. The RS formula for each sample was as follows: RS= β1*Exp1 + β2*Exp2 + βi*Expi, where β represented the coefficient value, and Exp represented the gene expression level. According to the median RS survival analysis, STAD patients were divided into low-risk and high-risk groups. OS of the two subgroups were compared by Log rank test. The predictive capability of the aforementioned prognostic model was evaluated by SurvivalROC package.14 TCGA testing cohort and GEO cohort were used as the validation groups to confirm the predictive capability of this prognostic model. The prognostic value of hub RBPs in STAD was also verified by the Kaplan-Meier plotter ()15 online tool. Finally, the nomogram was conducted using R package to forecast the likelihood of OS. P<0.05 was considered as a significant difference.

Biology Network

To understand the mechanisms and functions of hub RBPs in STAD, we created tissue-specific protein–protein interaction (PPI), transcription factor (TF)–gene interaction and miRNA–gene interaction network of prognosis-related hub RBPs under the condition that betweenness was greater than 1 using Cytoscape 3.7.2, based on the data from NetworkAnalyst 3.0 ().16 Function enrichment analysis of prognosis-related hub RBPs was performed by the ClueGO plug to Cytoscape 3.7.2.

Tumor Immune Infiltrating Analysis

To fully understand the function of hub RBPs, we assessed the relationship between hub RBPs expression and tumor-infiltrating lymphocytes (TILs) abundance via TISIDB ().17 TISIDB is a web portal for tumor and immune system interaction, which integrates multiple heterogeneous data types.

Cell Culture and siRNA Transfection

AGS and HGC27 cell lines were purchased from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). HGC27 cells were cultured in RPMI-1640 medium (Biological Industries, Shanghai, China), while AGS cells were cultured in DMEM/F-12 medium (Biological Industries, Shanghai, China). All media were supplemented with 10% fetal bovine serum (FBS; Gibco, Logan, UT, USA). Cells were grown at 37°C in an atmosphere of 5% CO2. The genOFF PTBP1 siRNA and negative control siRNA oligonucleotides were designed and synthesized by Ribobio (Guangzhou, China). Sequences of si-1 and si-2 were as follows: GCCTCAACGTCAAGTACAA, and GCACAGTGTTGAAGATCAT. The siRNA transfections were performed with Ribo FECT™CP Transrection Kit (Ribobio, Guangzhou, China) in accordance with the manufacturer’s instructions. In-vitro experiments were approved by the Ethics Committee of The Seventh Affiliated Hospital of Sun Yat-sen University.

RNA Extraction and Quantitative Real-Time PCR

Under instructions provided by the manufacture, total RNA was retrieved with AG RNAex Pro Reagent (Accurate Biology, Changsha, China) and was reverse‐transcribed into cDNA by Evo M-MLV RT Premix Kit (Accurate Biology). Quantitative real-time PCR (RT-PCR) assays were performed using SYBR® Green Premix Pro Taq HS qPCR Kit (Accurate Biology) by CFX96 Touch Deep Well Real-Time PCR System (Bio-Rad, Guangzhou, China) according to the manufacturers’ protocols. The PCR primer sequences were: GAPDH(F):ACAACTTTGGTATCGTGGAAGG; GAPDH(R):GCCATCACGCCACAGTTTC; PTBP1(F):AGCGCGTGAAGATCCTGTTC; PTBP1(R):CAGGGGTGAGTTGCCGTAG.

Western Blotting

Cells were washed with PBS and lysed with RIPA lysis buffer containing phosphatase and protease inhibitor reagents (Beyotime Biotechnology, Shanghai, China), which were then separated with the method of SDS-PAGE. Then, the proteins separated by SDS-PAGE were transferred onto PVDF membranes (Merck Millipore, Billerica, MA, USA). After being loaded with transferred proteins, these PVDF membranes were soaked in 5% BSA (Beyotime Biotechnology) solution for at least 1 hour at room temperature, which was immediately followed by incubation with different primary antibodies at 4°C overnight. Then on the second day, after being washed three times, the PVDF membranes were soaked into a solution containing HRP-conjugated goat anti-rabbit IgG (SA00001-2, Proteintech, Wuhan, China) for at least 1 hour at room temperature. After being washed with TBST solution for 20 minutes, these PVDF membranes containing antibody-conjugated protein bands were visualized with a BeyoECL Plus Kit (Beyotime Biotechnology). Anti-GADPH (10,494-1-AP) and anti-PTBP1 antibody (12,582-1-AP) were purchased from Proteintech.

Cell Migration and Invasion Assays

In order to assess the effects of PTBP1 expression on the migration and invasion capability of STAD cells, we subsequently performed transwell assays. For migration assay, 5 × 10^4 cells were seeded into the upper chamber containing serum-free medium (8-μm pore size, Corning, NY, USA), while the lower chamber was filled with 10% fetal bovine serum-supplemented medium. For invasion assay, 5 × 10^4 cells were seeded into the upper chamber covered with 10% Matrigel (Corning). AGS and HGC27 cells were incubated for 24 h. Cells left in the upper chamber after the incubation period were wiped off by swabs. Then, cells that have transpassed the membrane were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet. These cells migrating or invading the lower chamber were counted using five different fields under a microscope.

Statistical Analysis

For in-vitro experiments, all quantitative data are presented as the mean ± standard deviation of three independent experiments. Differences between the three groups were analyzed with one-way ANOVA using GraphPad Prism 8.0 (GraphPad, La Jolla, CA, USA). A P-value <0.05 was considered statistically significant.

Results

Differentially Expressed RBPs in STAD

The study design is illustrated in Figure 1. It was revealed by GEPIA analysis that DEGs in STAD included 896 downregulated genes and 3475 upregulated genes (Figure 2A). Among the 1542 RBPs, 362 were differentially expressed with 331 upregulated and 31 downregulated (Figure 2B, ).

Figure 1

Framework for analyzing RBPs in STAD.

Figure 2

Volcano plot of related DEGs in STAD. (A) Volcano plot of all DEGs in STAD. (B) Volcano plot of differentially expressed RBPs in STAD.

GO and KEGG Pathway Enrichment Analysis of the Differentially Expressed RBPs

To explore the functions and mechanisms of these differentially expressed RBPs, we performed the functional analysis for downregulated and upregulated RBPs via WebGestalt.

As shown in Table 1, significant differences were observed in functional enrichment of downregulated and upregulated RBPs. As for localization within the cell, downregulated RBPs were enriched in U2-type spliceosomal complex, spliceosomal snRNP complex, U2-type prespliceosome, prespliceosome and U1 snRNP, and upregulated RBPs in U1 snRNP, U4 snRNP, U12-type spliceosomal complex, box H/ACA snoRNP complex and histone pre-mRNA 3ʹend processing complex. Differences in localization within the cell meant different molecular function, which was further confirmed by molecular functional analysis. Molecular functional analysis demonstrated that downregulated RBPs participated in RNA binding, mRNA binding and snRNA binding and upregulated RBPs in RNA 7-methylguanosine cap binding, exoribonuclease activity, transforming growth factor-beta receptor, box H/ACA snoRNA binding and snRNP binding (Table 1). Biological process analysis showed that downregulated RBPs were related to positive regulation of mRNA processing, mRNA splice site selection, positive regulation of RNA splicing, positive regulation of mRNA splicing and regulation of chaperone-mediated autophagy while upregulated RBPs to ribosomal subunit export from the nucleus, termination of RNA polymerase II transcription, positive regulation of cytoplasmic mRNA processing body assembly, regulation of ribonuclease activity and mRNA cleavage involved in gene silencing (Table 1).

Table 1

GO Enrichment and KEGG Pathway Analysis of Differently Expressed RBPs

	GO Term	Name	P value	FDR
Down-regulated RBPs	Cellular Component	U2-type spliceosomal complex	2.07E-08	6.09E-06
		Spliceosomal snRNP complex	4.95E-07	9.51E-05
		U2-type prespliceosome	5.67E-07	9.51E-05
		Prespliceosome	5.67E-07	9.51E-05
		U1 snRNP	3.68E-09	1.44E-06
	Molecular Function	RNA binding	0	0
		mRNA binding	5.81E-9	5.45E-6
		snRNA binding	4.60E-5	2.88E-2
	Biological Process	Positive regulation of mRNA processing	2.27E-05	1.15E-02
		mRNA splice site selection	2.06E-05	1.10E-02
		Positive regulation of RNA splicing	4.19E-07	4.23E-04
		Positive regulation of mRNA splicing, via spliceosome	7.14E-06	4.06E-03
		Regulation of chaperone-mediated autophagy	8.14E-05	3.89E-02
	KEGG pathway	Spliceosome	3.19E-5	1.04E-2
Up-regulated RBPs	Cellular Component	U1 snRNP	2.17E-12	1.06E-10
		U4 snRNP	2.81E-9	8.48E-8
		U12-type spliceosomal complex	1.33E-15	1.04E-13
		Box H/ACA snoRNP complex	3.73E-5	6.09E-4
		Histone pre-mRNA 3ʹend processing complex	5.43E-9	1.59E-7
	Molecular Function	RNA 7-methylguanosine cap binding	1.50E-5	6.40E-4
		Exoribonuclease activity, producing 5ʹ-phosphomonoesters	6.66E-16	1.25E-13
		Transforming growth factor beta receptor, pathway-specific cytoplasmic mediator	6.62E-5	2.44E-3
		Box H/ACA snoRNA binding	6.62E-5	2.44E-3
		snRNP binding	4.30E-6	2.18E-4
	Biological Process	Ribosomal subunit export from nucleus	1.59E-07	1.24E-05
		Termination of RNA polymerase II transcription	1.38E-14	1.90E-12
		Positive regulation of cytoplasmic mRNA processing body assembly	1.18E-04	4.84E-03
		Regulation of ribonuclease activity	2.40E-07	1.80E-05
		mRNA cleavage involved in gene silencing	5.97E-05	2.66E-03
	KEGG pathway	Aminoacyl-tRNA biosynthesis	3.13E-8	1.70E-6
		Ribosome biogenesis in eukaryotes	0	0
		RNA transport	2.22E-16	2.41E-14
		mRNA surveillance pathway	5.33E-15	4.34E-13
		RNA degradation	6.75E-14	4.40E-12
		Spliceosome	0	0

Moreover, downregulated RBPs were only enriched in the spliceosome pathway, while upregulated RBPs significantly in aminoacyl-tRNA biosynthesis, ribosome biogenesis in eukaryotes, RNA transport, mRNA surveillance pathway, RNA degradation and spliceosome (Table 1).

Prognosis-Related Hub RBPs

To further analyze the effects of RBPs on the prognosis of STAD patients, we assessed the relationship between the differentially expressed RBPs and OS through the univariate Cox regression analysis and Kaplan–Meier method in the TCGA training cohort, the results of which suggested that 25 candidate hub RBPs were significantly associated with OS (Table 2). Afterwards, the impacts of these 25 candidate hub RBPs on OS were evaluated by multivariate analysis, the results of which demonstrated that seven hub RBPs were independent prognostic predictors for STAD patients (Table 2).

Table 2

Cox Regression Analysis for Identification of Hub RBPs in the Training Dataset

Gene	Univariate Analysis				Multivariate Analysis				Coef
	P value	HR	95.0% CI		P value	HR	95.0% CI
			Lower	Upper			Lower	Upper
EZH2	0.040*	0.943	0.892	0.997
RAVER1	0.006**	0.959	0.931	0.988
IFIT1	0.042*	1.020	1.001	1.039
PTBP1	0.037*	0.990	0.981	0.999	0.000***	1.041	1.020	1.062	0.040
TTF2	0.038*	0.878	0.777	0.993
CSTF2	0.025*	0.956	0.919	0.994
PATL1	0.046*	0.980	0.961	0.999
RNASE6	0.036*	1.011	1.001	1.022
PPIH	0.040*	0.975	0.951	0.999	0.011*	0.949	0.912	0.988	−0.052
SMAD5	0.009**	1.077	1.018	1.140	0.000***	1.187	1.078	1.307	0.172
MSI2	0.004**	0.881	0.808	0.962	0.003**	0.812	0.708	0.931	−0.208
RBM15	0.003**	0.848	0.760	0.946	0.000***	0.623	0.509	0.762	−0.474
REPIN1	0.010*	0.989	0.981	0.997
METTL2B	0.002**	0.802	0.697	0.922
MRPS17	0.000***	1.062	1.030	1.096	0.000***	1.075	1.042	1.109	0.072
YBX1	0.019*	0.998	0.996	0.999
DAZAP1	0.009**	0.944	0.904	0.986
FTO	0.013*	1.160	1.032	1.304
SF3A3	0.044*	0.953	0.909	0.999
ADAT3	0.001**	0.854	0.780	0.934	0.001**	0.781	0.675	0.903	−0.248
BICC1	0.005**	1.072	1.021	1.125
RPS4Y1	0.015*	1.006	1.001	1.010
LARP6	0.029*	1.099	1.010	1.197
DDX3Y	0.014*	1.047	1.010	1.087
POLR2J2	0.045*	202,696.903	1.120	3.67E+10

Notes: Statistical significance was described as follows: *p <0.05, **p<0.01, ***p<0.001.

Construction of Prognostic Model

A predictive model based on the aforementioned seven hub RBPs was then established. According to the formula: RS= (0.040* Exp PTBP1) + (−0.052* Exp PPIH) + (0.172* Exp SMAD5) + (−0.208* Exp MSI2) + (−0.474* Exp RBM15) + (0.072* Exp MRPS17) + (−0.248* Exp ADAT3), RS of each individual patient was assessed. The predictive capability of RS was evaluated by survival analysis. According to the median RS, the 187 patients from TCGA training cohort were assigned into low-risk group and high-risk group. Results of survival analysis demonstrated that compared with patients of the low-risk group, those of high-risk group had significantly poorer OS (p<0.001, Figure 3A). In order to further evaluate the prognostic capability of the seven identified hub RBPs, we subsequently performed a time-dependent ROC analysis, the results of which demonstrated that the area under the ROC curve (AUC) of this RBPs RS model was 0.804 (Figure 3B), indicating the moderate diagnostic performance of this model. The survival status of patients, RS and expression heat map of the signature consisting of seven hub RBPs in the low- and high-risk subgroups are displayed in Figure 3C–E.

Figure 3

Risk score analysis of seven-genes prognostic model in the TCGA training cohort. (A) Survival curve for low- and high-risk subgroups. (B) ROC curves for forecasting OS based on risk score. (C) Survival status. (D) Risk score. (E) Expression heat map.

Validation of Prognostic Model

To further verify the validity of the seven RBPs-based predictive model, we analyzed the TCGA testing cohort included 184 STAD patients and GEO cohort (GSE84437) included 433 STAD patients, the results of which showed that compared with patients with low-risk score, those with high-risk score had significantly worse OS (p<0.05, Figure 4A; p<0.05, ). AUC of the TCGA testing cohort and GEO cohort were 0.644 (Figure 4B) and 0.581 (), which suggested good sensitivity and specificity of the predictive model. The survival status of patients, RS and expression heat map of seven hub RBPs in the TCGA testing cohorts and GEO cohort are shown in Figure 4C–E and . Additionally, the prognostic significance of different variables was assessed among patients of TCGA training cohort, TCGA testing cohort and GEO cohort by Cox regression analysis, the results of which demonstrated that for three cohorts, RS were independent prognostic factors of OS (p<0.01, Figure 5A; p<0.05, Figure 5B; p<0.01, ). The prognostic values of the seven hub RBPs were also further investigated by Kaplan-Meier plotter online tool, the results of which revealed that all the seven hub RBPs were not only significantly associated with OS but also with relapse free survival (RFS) (). In summary, all the aforementioned results suggested that the seven RBPs-based prognostic model was reliable in predicting outcomes of STAD patients.

Figure 4

Risk score analysis of seven-genes prognostic model in the TCGA testing cohort. (A) Survival curve for low- and high-risk subgroups. (B) ROC curves for forecasting OS based on risk score. (C) Survival status. (D) Risk score. (E) Expression heat map.

Figure 5

The prognostic value of different clinical parameters. (A, B) Univariate and multivariate COX regression analysis of different clinical parameters in TCGA training and TCGA testing cohort.

Building a Predictive Nomogram

A nomogram was constructed to generate a clinically practical model that would enable physicians to evaluate the prognosis of STAD patients using the seven hub RBPs (Figure 6). Based on the results of multivariate Cox analysis, the corresponding points were assigned to each individual variable according to the point scale in the nomogram. A horizontal line was drawn to determine the point of each variable. The total point for each patient was calculated by adding up the points of all variables, based on which we estimated the survival rate of each patient at 1 year, 2 years, and 3 years.

Figure 6

Nomogram for predicting 1, 2, and 3 year OS of STAD patients in the TCGA training cohort.

Biology Network and Functions of the Seven Hub RBPs

To further investigate the functions of the seven hub RBPs in STAD, we created TF-gene, micro-RNA and tissue-specific PPI network of the seven hub RBPs. TF-gene-specific network included 68 nodes, 152 edges and 7 seeds (PTBP1, PPIH, SMAD5, MSI2, RBM15, MRPS17, and ADAT3) (Figure 7A), miRNA-gene-specific network 36 nodes, 64 edges and 4 seeds (MRPS17, MSI2, SMAD5 and PTBP1) (Figure 7B), and tissue-specific PPI network 14 nodes, 19 edges and 5 seeds (PPIH, RBM15, PTBP1, SMAD5 and MSI2) (Figure 7C). Finally, function enrichment analysis revealed that they were enriched in tRNA wobble base modification, thrombopoietin-mediated signaling pathway, Mullerian duct regression, poly-pyrimidine tract binding and U4/U6 snRNP (Figure 7D). All these results suggested that the seven hub RBPs were widely involved in many biological processes.

Figure 7

Biology network and functions of the seven hub RBPs. (A) TF-gene network. (B) miRNA-gene network. (C) Tissue-specific PPI network. (D) Functional enrichment network.

Assessment of Tumor Immune Infiltration

The effects of the seven RBPs on tumor immune infiltration were also explored, the results of which demonstrated that expression of MRPS17 and PTBP1 was negatively significantly correlated with the abundance of TILs. It was also revealed that expression of MRPS17 was negatively correlated with Th1 cell, Tem-CD8 cell, NKT cell and NK cell abundance (Figure 8A). The expression of PTBP1 was also revealed to be negatively correlated with Th1 cell, Tem-CD4 cell, NKT cell and eosinophil abundance (Figure 8B). All these results indicated that MRPS17 and PTBP1 may reduce the infiltration of immune cells.

Figure 8

The relationship between the hub RBPs expression and TILs abundance. (A) The relationship between MRPS17 expression and TILs abundance. (B) The relationship between PTBP1 expression and TILs abundance.

The Downregulation of PTBP1 Weakened the Migration and Invasion Capability of STAD Cells

Since the role of PTBP1 in STAD is still unclear, we further explored the effect of PTBP1 on the migration and invasion capability of STAD cells. PTBP1 mRNA expression in multiple STAD cell lines was detected by RT-PCR (Figure 9A). Subsequently, siRNAs (si-NC, si-1, and si-2) were transfected into AGS and HGC27 cells. RT-PCR and Western Blotting assays showed that PTBP1 expression was significantly down-regulated in STAD cells (AGS and HGC27) transfected with si-1 or si-2 compared with those transfected with si-NC (Figure 9B). Then, it was proven by transwell assays that both migration and invasion of AGS (p<0.01, Figure 9C) and HGC27 (p<0.01, Figure 9D) were significantly reduced by knockdown of PTBP1. Therefore, we could infer that PTBP1 played important roles in promoting metastasis of STAD.

Figure 9

The downregulation of PTBP1 weakened the migration and invasion capability of STAD cells. (A) PTBP1 mRNA expression level in STAD cells. (B) PTBP1 knockdown AGS and HGC27 cells were constructed and then confirmed by RT-PCR and Western blotting. (C, D) Migration and invasion capability of AGS and HGC27 cells were significantly weakened by downregulation of PTBP1. (**p<0 0.01, ***p <0.001).

Discussion

In the present study, 362 differentially expressed RBPs were identified from GEPIA based on STAD TCGA data. We comprehensively analyzed the relevant biological functions and pathways. Additionally, seven hub RBPs were found to be prognosis-related by univariate and multivariate Cox regression analysis, and ROC analysis. A risk model based on the seven hub RBPs was established in TCGA training cohort and validated in TCGA testing cohort and GEO cohort. The risk model was proved reliable in predicting OS. To improve its clinical practicality, we then constructed a nomogram to evaluate the survival of STAD patients at 1 year, 2 years and 3 years. Furthermore, to understand the mechanisms and functions through which they contributed to STAD, we investigated TF gene, miRNA gene, and tissue-specific PPI network and assessed their effects on tumor immune infiltration. In addition, we assessed the effects of PTBP1 expression on the migration and invasion capability of STAD cells. These findings may enable us to better understand the mechanisms involved in the occurrence and progression of STAD and develop novel biomarkers for diagnosis and predicting prognoses of STAD patients.

Although 362 differentially expressed RBPs were identified in STAD, functional enrichment analysis demonstrated that they were significantly enriched in post-transcriptional regulation of RNA, including mRNA processing, RNA splicing, termination of RNA polymerase II transcription, cytoplasmic mRNA processing body assembly, ribonuclease activity, and mRNA cleavage involved in gene silencing. As known to all, it has been reported that post-transcriptional regulation of RNA plays an important role in many types of cancers.6 A recently published study revealed that RBPs significantly affected the mRNA alternative splicing process and led to altered expression of tumor-associated genes.18 As one kind of RBP, SRSF1 promoted the occurrence and progression of non-small cell lung cancer through stabilizing LIG1 mRNA.19 Additionally, it was revealed by KEGG pathway analysis that it was by affecting spliceosome, aminoacyl-tRNA biosynthesis, RNA transport, mRNA surveillance pathway and RNA degradation that RBPs contributed to the occurrence and progression of STAD.

To explore the prognostic significance of the hub RBPs, we performed Cox survival analysis for these differentially expressed RBPs in TCGA training cohort. A total of seven RBPs that included PTBP1, PPIH, SMAD5, MSI2, RBM15, MRPS17, and ADAT3 were proven to be prognosis-related. Subsequently, based on the seven identified hub RBPs, we established a risk model for predicting prognoses of STAD patients in TCGA training cohort. The ROC analysis demonstrated that the seven RBPs-based risk model was reliable in selecting out STAD patients with poor prognosis, which was further validated by the TCGA testing cohort and GEO cohort. The aforementioned reliability of the risk model was also revealed by Kaplan-Meier plotter online tool. Consistent with our findings, Su et al20 demonstrated that high expression of RBM15 was associated with better survival of STAD patients. Additionally, revealed by multiple Cox regression analysis of the three cohorts, RS was an independent prognostic factor for STAD patients. All the aforementioned findings of our study suggested the clinical practicality of this seven RBPs-based prognostic model. To enable the physicians to predict patients, survival at 1 year, 2 years, and 3 years more intuitively, we constructed a nomogram.

Nevertheless, the molecular mechanisms through which these seven RBPs contribute to the pathogenesis of STAD remains poorly understood. To solve this puzzle, we further investigated the potential mechanisms at different levels. As shown in the network of TF-gene, these RBPs can be regulated by various TFs. As we all have known TF, MYC promotes the progression of STAD through regulating downstream genes and different pathways.21 FOXM1 has also been testified to play a crucial role in the occurrence and chemoresistance of gastric cancer.22 As the potential targets of TFs, RBPs and their roles in STAD deserve our further efforts. Network analysis of miRNA-gene revealed that many miRNAs could target PTBP1, MSI2, MRPS17, and SMAD5. The roles miRNAs play in cancer need no further illumination as it has been revealed by many studies that miRNA acts as oncogene or tumor suppressor gene in many kinds of cancers.23 It has been reported that miRNA-133b inhibited the growth of gastric cancer cells through silencing PTBP1.24 Despite the crucial roles played by miRNA155-5p,25 miRNA-12426 and miRNA-17,27 it remains unclear whether they regulate RBPs. Thus, the roles of miRNA-RBPs axis in STAD deserve to be further investigated. In our study, tissue-specific PPI was also constructed to confirm the proteins that interacted with the seven RBPs, the results of which demonstrated that PTBP1, MSI2, PPIH, RBM15, and SMAD5 interacted with nine other proteins that included HMGA2, FN1, SOX2, CUL3, CUL7, RPA2, RNF2, CALCOCO2, and MAP1LC3A. HMGA2 promoted the metastasis and epithelial-to-mesenchymal transition of chemoresistant gastric cancer.28 Down-regulation of FN1 suppressed the proliferation, migration and invasion of gastric cancer cells.29 The role of SOX2 in gastric cancer is still controversial.30,31 Thus, it is necessary to confirm the roles of these hub RBPs and other proteins to help us understand the mechanisms contributing to STAD. For the last part, the effects of the seven RBPs on tumor immune infiltration were explored, the results of which demonstrated that both MRPS17 and PTBP1 were negatively associated with the abundance of TILs, such as Th1 cell, NKT cell, Tem CD8 cell, and Tem CD4 cell. As known to all, these infiltrated cells are indispensable in killing cells,32 which may explain the phenomenon that MRPS17 and PTBP1 affect the prognosis of STAD. Multiple studies have reported that PTBP1 is closely associated with immune response. PTBP1 is crucial to B-cell receptor-mediated activation33 and knockdown of PTBP1 inhibits tumor-promoting functions of SASP and influences immune surveillance.34 All these studies remind us that RBPs play multiple roles at different levels.

To further verify the role of these hub RBPs in STAD, we performed functional experiments in STAD cells. PTBP1 was selected for in-vitro experiments because its role in STAD is still unclear. Transwell assays showed that the downregulation of PTBP1 significantly weakened the migration and invasion capability of STAD cells. These results suggested that PTBP1 played important roles in the metastasis of STAD, which provided more powerful support for us to further explore the roles of RBPs in STAD. Although the mechanisms through which RBPs contribute to STAD remain to be solved, comprehensive investigation of abnormally expressed RBPs that are closely related to the occurrence and progression of STAD may bring us the potential targets to design novel strategies and the chances to improve outcomes of STAD patients.

Conclusion

In conclusion, we comprehensively investigated the prognostic values and potential functions of differentially expressed RBPs in STAD. A risk model that can reliably predict the prognosis of STAD patients was constructed based on the seven identified RBPs and validated by the TCGA testing cohort and GEO cohort. Our study may provide important evidences for future studies on the roles of RBPs in STAD.