Tumor Immune Estimation Resource (TIMER, http://timer.cistrome.org/) is a comprehensive resource for systematic analysis of immune infiltration across diverse cancer types, including gene expression, clinical outcomes, somatic mutations, somatic copy number alterations, and the characterization of the tumor-immune system. The database includes 10,897 samples across 32 cancer types from TCGA to estimate the abundance of immune infiltration (Li et al., 2017; Li et al., 2020). We input PPP1R14B in the “Gene_DE” module of TIMER2.0 web and analyzed the expression difference of PPP1R14B between tumor and adjacent normal tissues for the 32 cancer types derived from the TCGA project.

GEPIA (http://gepia.cancer-pku.cn/) is a web-based tool to deliver fast and customizable functions, including differential expression analysis, correlation analysis, patient survival analysis, similar gene detection, and dimensionality reduction analysis (Tang et al., 2017; Tang et al., 2019). We analyzed the expression of PPP1R14B across TCGA tumors, and the matched TCGA normal and GTEx data were included as controls.

UALCAN (http://ualcan.path.uab.edu), an interactive web-portal, uses TCGA RNA-seq and clinical data from 31 cancer types and can be easily used for in-depth analyses of TCGA gene expression data (Chandrashekar et al., 2017). Herein, we analyzed the relative expression of PPP1R14B across tumor and normal samples, as well as in various tumor models based on individual tumor grade. The protein expression was analyzed using the data from the Clinical Proteomic Tumor Analysis Consortium (CPTAC, http://ualcan.path.uab.edu/analysis-prot.html).

The GEPIA2.0 web-based tool was used flor patient survival analysis. We input the PPP1R14B gene in the “survival analysis” module, selected a custom cancer type, and used the log-rank test for overall survival analysis. High-expression and low-expression cohorts were splitted by the high (50%) and low (50%) cutoff values.

The Kaplan–Meier plotter (http://kmplot.com/analysis/) is capable of assessing the effect of 54k genes (mRNA, miRNA, and protein) on survival in 21 cancer types. Sources of the databases include GEO, EGA, and TCGA (Nagy et al., 2021). We analyzed the correlation between PPP1R14B expression and survival in different cancer types in the “pan-cancer RNA-seq” module.

The cBioPortal web (https://www.cbioportal.org/) was used for genetic alteration analysis (Cerami et al., 2012). We chose the “curated set of non-redundant studies” (184 studies, 48,035 samples) in the “query” section and entered “PPP1R14B” for queries about the genetic alteration characteristics of PPP1R14B. The results of the alteration frequency, mutation type, and CNA (copy number alteration) across all tumors were observed in the “cancer types summary” module. Then, the data of the PP1R14B mRNA expression level vs. putative copy number alterations (mRNA VS CNA) were obtained in the “plots” module and then analyzed by using Graphpad software. The mutated site information of PPP1R14B can be displayed in the “mutations” module. We also used the “comparison/survival” module to obtain the data on the overall and relapse-free survival differences for the cancer cases with or without PPP1R14B genetic alteration.

We used the TIMER2.0 (http://timer.cistrome.org/) resource for systematic analysis of immune infiltrates across diverse cancer types (Li et al., 2020). We input the PPP1R14B gene in the “gene” module in the “immune” section and then selected the immune infiltrate type to obtain the data about the correlation of its expression with the immune infiltration level in diverse cancer types. We analyzed the important immune infiltrating cells, including B cells, CD4⁺ T cells, CD8⁺ T cells, Tregs, Tfh cells, γδ T cells, monocytes, macrophages, neutrophils, DC, and MDSCs.

We used the STRING website (https://string-db.org/) to obtain the PPP1R14B-binding proteins for protein–protein interaction network functional enrichment analysis (von Mering et al., 2003). We searched PPP1R14B in the query of “Protein by name” and chose “Homo sapiens”. Subsequently, we set the following main parameters: meaning of network edges (“evidence”), active interaction sources (“textmining”, “experiments”, “databases”, “co expression”, “neighborhood”, “gene fusion”, and “co-occurrence”), network type (“full network”), minimum required interaction score [“medium confidence (0.400)”], and maximum number of interactors to show (“no more than 20 interactors” in 1st shell).

We used TIMER2.0 to explore the correlation between PPP1R14B with a list of genes of the 20 interactors of PPP1R14B-binding in various cancer types. We input the PPP1R14B gene in the “Gene_Corr” section in the “exploration” module and then input a list of genes of the 19 interactors (ENSG00000108825 is not retrieved in the database) to obtain the data about the correlation of its expression.

Jvenn (http://bioinfo.genotoul.fr/jvenn) is an open-source component for web environments to process lists and produce Venn diagrams (Bardou et al., 2014). We pasted up our lists with one element per row to conduct an intersection analysis to compare the PPP1R14B-binding and interactive genes.

We used ShinyGO (http://bioinformatics.sdstate.edu/go/), a graphical gene set enrichment tool (Ge et al., 2020) for enrichment analysis by combining the PPP1R14B-correlated and interacted genes. Then we used the Gene Ontology chord plot tool in bioinformatics analyses (http://www.bioinformatics.com.cn/) to visualize the enrichment results.

Data about the association of PPP1R14B copy number alteration with its mRNA expression was obtained from the cBioPortal web (https://www.cbioportal.org/) and then analyzed using unpaired t-tests by GraphPad Prism (Version 8.0.1) for Windows. P-values less than 0.05 were considered statistically significant. The following annotations were used to show statistical significance: *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001.

We used the TCGA project and GEO databases to perform pan-cancer analysis of PPP1R14B, including expression differences, correlations between expression levels and survival, genetic alteration, immune infiltration, and relevant cellular pathways, to investigate the functions and potential mechanisms of PPP1R14B in the pathogenesis or clinical prognosis of different cancers. The setup of this study is shown in Figure 1.

To determine the differences of PPP1R14B expression between tumor and adjacent normal tissues, we used the TIMER2.0 database to analyze PPP1R14B mRNA expression levels across all TCGA tumors. The results showed that PPP1R14B was highly expressed in BLCA (bladder urothelial carcinoma), BRCA (breast invasive carcinoma), CESC (cervical squamous cell carcinoma and endo-cervical adenocarcinoma), CHOL (cholangio carcinoma), COAD (colon adenocarcinoma), ESCA (esophageal carcinoma), GBM (glioblastoma multiforme), HNSC (head and neck squamous cell carcinoma), KICH (kidney chromophobe), KIRC (kidney renal clear cell carcinoma), KIRP (kidney renal papillary cell carcinoma), LIHC (liver hepatocellular carcinoma), LUAD (lung adenocarcinoma), LUSC (lung squamous cell), PRAD (prostate adenocarcinoma), READ (rectum adenocarcinoma), STAD (stomach adenocarcinoma), THCA (thyroid carcinoma), and UCEC (uterine corpus endometrial carcinoma) compared with their adjacent normal tissues (Figure 2A).

In addition, considering that the TIEMR2.0 database lacks normal controls for some tumors, we chose the normal tissues of the GTEx dataset as controls to evaluate the expression difference of PPP1R14B between the normal and tumor tissues of BLCA (bladder urothelial carcinoma), DLBCL (lymphoid neoplasm diffuse large B-cell lymphoma), GBM (glioblastoma multiforme), LGG (brain lower grade glioma), OV (ovarian serous cystadenocarcinoma), and UCS (uterine carcinosarcoma). As shown in Figure 2B, PPP1R14B was highly expressed in the abovementioned tumor tissues compared with normal tissues.

Proteins are undoubtedly the main molecule related to disease, and the change in the protein expression level is directly related to disease, drug action, or toxin action. Therefore, we further analyzed the expression difference of the PPP1R14B protein between tumor and normal tissues using data from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) (Ellis et al., 2013). The results showed that the PPP1R14B protein was highly expressed in ccRCC (Clear cell renal cell carcinoma), COAD, LUAD, OV, and UCEC (Uterine corpus endometrial carcinoma) compared with the normal tissues (Figure 2C).

Next, we investigated the expression of PPP1R14B according to the pathologic stage of the patients in the TCGA cancer type. We found that in BLCA, BRCA, CESC, COAD, ESCA, HNSC, KICH, KIRC, LIHC, LUAD, LUSC, READ, STAD, THCA, and UCEC, PPP1R14B expression levels were significantly higher in early stages. But there were no significant differences between the late stages and early stages. This indicates a possible involvement of PPP1R14B in the initiation but not the progression of cancer (Figure 3).

We divided cancer cases into high-expression and low-expression groups according to the mRNA expression level of PPP1R14B and investigated the correlation between PPP1R14B expression and the prognosis of patients with different tumors in the TCGA and GEO datasets. These results showed that highly expressed PPP1R14B was linked to poor prognosis of overall survival (OS) for cancers of ACC (p = 8.1e-07), MESO (P = 8e-04), SKCM (p = 0.02), and UVM (p = 0.042) within the TCGA project in GEPIA2 (Figures 4A–D). In addition, the analysis of survival data using the Kaplan–Meier plotter tool showed a correlation between high expression of PPP1R14B and poor OS for BRCA (p = 2.7e-07), CESC (p = 0.0081), HNSC (p = 0.0024), KIRC (P = 9e-08), LIHC (2e-06), LUAD (p = 0.0041), PAAD (p = 0.00073), SARC (p = 0.013), STAD (p = 0.018), and UCEC (p = 0.013) within the TCGA project (Figures 4E–N). All of the abovementioned data showed that high expression of PPP1R14B was closely related to poor prognosis of most tumors, which may be a promising pan-cancer prognostic marker.

In order to more systematically clarify the mutation feature and biological functions of PPP1R14B in tumor progression, we used the cBioPortal web to investigate the genetic alteration status of PPP1R14B in human pan-cancer. We observed that the highest alteration frequency of PPP1R14B (>4%) appears for patients with bladder/urinary tract cancer with “Amplification” as the primary type. In addition, the “amplification” type of CNA was also the primary type of genetic alteration in endometrial carcinoma, pancreatic cancer, cholangio carcinoma, and prostate cancer, of which the alteration frequency exceeded 2%. We also observed that the “mutation” type of CNA was the primary alteration type in cervical adenocarcinoma and lung cancer with an alteration frequency of 1–2% (Figure 5A). Then, we retrieved the PPP1R14B expression data from TCGA pan-cancer dataset and found that PPP1R14B was amplified along with the significantly high mRNA expression (Figure 5B). The types, sites, and case numbers of the PPP1R14B genetic alteration are further presented in Figure 5C. We found that the missense mutation of PPP1R14B was the main type of genetic alteration, which was detected in 18 cases. In addition, five cases contained the splice mutation and two cases carried “In_Frame_Del” (IF del) mutation. The fusion, FS ins, and nonsense mutations were detected in one case alone. Additionally, we explored the association between genetic alteration of PPP1R14B and clinical survival in different cancer. As shown in Figures 5D,E, the BRCA cases with altered PPP1R14B showed poor OS (p = 4.031e-04) and relapse-free survival (p = 9.337e-06). The prostate cancer cases with altered PPP1R14B also showed poor OS (p = 2.59e-14) (Figure 5F). These results suggested that genetic alteration of PPP1R14B in tumors is involved in PPP1R14B mRNA expression alteration and poor clinical survival, which deserved more in-depth research.

The level of immune infiltration in the TME is closely correlated with the initiation, progression, or metastasis of cancer. We investigated the various immune cell infiltration levels in diverse cancer types in the TIMER database. As shown in Supplementary Table S1, the PPP1R14B expression level was significantly correlated with MDSC infiltration in BLCline 277-280A, BRCA, CESC, LIHC, LUAD, PAAD, PRAD, STAD, and THCA. However, there was no universally significant correlation between PPP1R14B expression and the other subgroup immune cell infiltrations, which included B cells, CD4⁺ T cells, CD8⁺ T cells, Tregs, Tfh cells, γδ T cells, monocytes, macrophages, neutrophils, and DC. Therefore, we focused on the MDSC in the TME. As shown in Figure 6A, the heatmap and scatter plot presented the relationship between the infiltrate estimation value and PPP1R14B gene expression. The results showed that in most cancer types, the expression of PPP1R14B is positively correlated with the MDSC infiltration level (Figure 6B). Furthermore, we investigated the correlation of the MDSC infiltration level and the prognosis of patients with different tumors in TCGA datasets in TIMER2.0. These results suggested that high levels of MDSC infiltration predict poor clinical outcomes in pan-cancer (Supplementary Figure S1). Our studies suggested that PPP1R14B is involved in the tumor immunology process.

To further study the molecular mechanism of PPP1R14B in tumorigenesis, we tried to screen out PPP1R14B-binding proteins for protein–protein interaction network analysis. Using the STRING online tool, we obtained a total of 20 PPP1R14B binding proteins, including members which were supported by experimental evidence or were predicted (Figure 7A). Then, we analyzed the correlation of PPP1R14B with the 20 interactors of PPP1R14B-binding in various cancer types. As shown in Figure 7B, the expression of NUBP2, PLCB3, PPP1CA, TBCB, and UBE2C were positively correlated with PPP1R14B in the majority of cancer types. Then, we used the GEPIA2 tool to combine all tumor expression data of TCGA and obtained the top 100 genes related to the expression of PPP1R14B. Through intersection analysis of the abovementioned two groups, we obtained two common members, namely, PPP1CA and UBE2C (Figure 7C). PPP1CA and UBE2C have strong positive correlations with PPP1R14B expression (Figure 7D).

In addition, we used the ShinyGO tool for enrichment analysis based on the two sets of data. Here, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Gene Ontology (GO) enrichment analysis. Then, the Gene Ontology chord plot tool in was used in bioinformatics to visualize the enrichment results. The KEGG pathway analysis data further indicated that PPP1R14B may be involved in oncogenesis through the “pyrimidine metabolism”, “ribosome biogenesis in eukaryotes”, “spliceosome” and “mRNA surveillance pathway” (Figure 7E). The GO enrichment analysis suggested that these genes were mainly related to “cellular component biogenesis”, “RNA processing”, “cell cycle”, “cell division”, and “ribonucleoprotein complex biogenesis” biological processes (Figure 7F).