FMR1 is identified as an immune-related novel prognostic biomarker for renal clear cell carcinoma: A bioinformatics analysis of TAZ/YAP

Abbreviations: AUC: Area under the curve; ANOVA: One-way analysis of variance; AGO2: Argonaute RISC catalytic component 2; BP: Biological processes; CC: Cellular components; CI: Confidence interval; DEGs: Differential expression genes; DFS: Disease-free survival; DFELMDA: Deep forest ensemble learning based on autoencoder; FC: Fold change; FDR: False discovery rate; FMR1: Fragile X mental retardation 1; GEO: Gene expression omnibus; GTEx: Genotype-tissue expression; GEPIA: Gene expression profiling interactive analysis; GSEA: Gene set enrichment analysis; GO: Gene ontology; HPA: Human protein atlas; KIRC: Kidney renal clear cell carcinoma; KIRP: Kidney renal papillary cell carcinoma; KICH: Kidney chromophobe; KEGG: Kyoto encyclopedia of genes and genomes; Limma: Linear models for microarray data; LMFNRLMI: Logistic matrix factorization with neighborhood regularized; MCODE: Molecular complex detection; MSigDB: Molecular signatures database; MF: Molecular functions; NES: Normalized enrichment score; NDALMA: LncRNA–miRNA association prediction; OS: Overall survival; PPI: Protein-protein interaction; PDL1: Programmed cell death ligand 1; RCC: Renal cell carcinoma; ROS: Reactive oxygen species; ROC: Receiver operating character; STRING: Search tool for the retrieval of interacting genes; STAU2: Staufen double-stranded RNA binding protein 2; TCGA: The cancer genome atlas; TIMER: Tumor immune estimation resource; TAZ: WW domain-containing transcription regulator 1; VEGFA: Vascular endothelial growth factors A; YAP: Yes-associated protein 1

1.   Introduction

Renal cell carcinoma (RCC) is the most prevalent and malignant solid tumor in the urinary system ^[1]. According to the statistics data of the 2020 global cancer report, the cases of birth and death are 431,288 and 179,368, respectively ^[2]; while the China statistic data show that there are approximately 73,587 birth cases and 43,196 death cases ^[3]. Generally, there are mainly three types of histologic features: kidney renal clear cell carcinoma (KIRC, 70%), kidney renal papillary cell carcinoma (KIRP, 10–15%), and kidney chromophobe (KICH, 5%) ^[4]. Currently, the first treatment choice is surgery in most KIRC patients; however, there is still a 30–40% recurrence rate in patients after surgery ^[5], with a 5-year overall survival (OS) rate of only approximately 12% ^[3]. Recently, emerging scholars have found that the changes in the tumor immune microenvironment were closely associated with the process of tumor ^[6,7]. Some researchers suggested that KIRC was an immune-infiltrated tumor ^[8]. The expression of immune-related genes was correlated with patients' prognosis ^[9,10]. The immune infiltration mediated approximately 1% of the spontaneous regression of patients ^[11]. Even though KIRC is one of the most sensitive tumors for immunotherapy, we are still unsatisfied with the clinical outcomes of patients ^[12]. Therefore, it is important to understand the regulatory mechanisms of the KIRC tumor immune microenvironment and identify immune-related prognostic biomarkers.

The Hippo pathway is a classic tumor suppressor signaling pathway to maintain the homeostasis of cells and organizations via regulating cell proliferation and differentiation, programmed cell death, immunity, metastasis, and chemotherapy suppression ^[13,14,15]. WW domain-containing transcription regulator 1 (TAZ, or WWTR1)/Yes-associated protein 1 (YAP) (homologous paralogues, partially common structures, and overlapping regulatory functions) are a type of transcriptional activators and both are essential Hippo pathway effectors ^[16]. When the Hippo pathway is activated, they are phosphorylated and degraded in the cytoplasm to inhibit the interaction with the TEA domain (TEAD) ^[17,18] or other transcription factors (paired box gene 3 (PAX3) ^[19,20], T-box transcription factor 5 (TBX5) ^[21,22], AP-1 ^[23] and SMADs ^[24]) in the nucleus, thereby influencing the homeostasis process of cells and organs through regulating the expression of common target genes. Recently, emerging evidence indicated that TAZ and YAP were the most common oncogenes. For example, a report stated that TAZ/YAP was inclined to interact with conserved enhancers and exhibited high chromatin profile activity in various cancers, which possibly was the fundamental reason to perturb the gene expression of cells ^[25]. Another report indicated that the TAZ/YAP mutation was tightly associated with six tumor development processes and immune cell infiltration ^[26]. However, it is intriguing that TAZ/YAP also possesses anti-cancer capacity. For example, TAZ/YAP contribute to ferroptosis and autophagy by regulating epithelial membrane protein 1 (EMP1)/angiopoietin-like 4 (ANGPTL4)-nicotinamide adenine dinucleotide phosphate (NADPH) oxidase 2 (NOX2) axis ^[27,28] and the expression of S-phase kinase-associated protein 2 (SKP2) ^[29] and LC3 ^[30], respectively. Presently, some studies, the combined analysis of TAZ/YAP, had revealed their common regulatory mechanisms. In cholangiocarcinoma, their low expressions inhibit the genes' expression of the cell cycle, DNA replication, and the Hippo pathway ^[31]. In RCC with neurofibromatosis type 2 (NF2) mutations, their depletions weaken the glycolysis and increase mitochondrial respiration and reactive oxygen species (ROS) to induce oxidative stress cell death ^[32]. In conclusion, it is complicated that TAZ/YAP was involved in the regulatory mechanisms. There are growing reports trying to understand their mechanisms and signaling pathways. However, due to some objective factors (cells' tolerance and sensitivity), it is not enough that we conduct the combined research of TAZ/YAP, so that we are incompletely clear about their specific regulatory mechanisms, especially in KIRC.

Bioinformatics is a novel interdisciplinary subject that involves the various processes of biological information by combining life sciences with computer technologies ^[33]. It is a great method to analyze complex biological data and understand the biological networks and regulatory mechanisms in theory. To date, numerous studies have been performed to uncover the roles of differential expression genes (DEGs) and identify many new biomarkers by using the aforementioned technologies in the corresponding diseases, including cancer. In KIRC, we screen many DEGs, identify various biomarkers, and construct the correlated prognostic risk models applying the methods to evaluate the status of patients from different aspects and elucidate the corresponding pathological mechanisms. Based on the research direction, the associated DEGs, markers, and models could be divided into the following types: tumor stage-related ^[34], immune-related ^[35,36,37], epigenetic-related ^[38,39,40], programmed cell death-related ^[41,42,43], drug-related ^[44,45], pathway-related ^[46,47], gene mutation-related ^[48], miRNA-related ^[49,50], lncRNA-related ^[51,52], etc. Therefore, it is significant that we found the diagnostic and prognostic markers and constructed the biological effect models to evaluate the KIRC patients using bioinformatics methods.

In this study, we carried out the comprehensive coalition analysis of TAZ/YAP using the bioinformatic methods to explore their co-regulatory mechanisms and screen new prognostic biomarkers in KIRC. Firstly, we combined GSE121689 with the GSE146354 database and identified the co-expressed DEGs. Secondly, we explored their shared regulatory network, constructed a protein-protein interaction (PPI) network, screened important modules, and identified the hub genes. Thirdly, we discussed the signaling pathways of the hub genes. Fourthly, the hub genes were verified in the samples of the human protein atlas (HPA), UALCAN, and gene expression profiling interactive analysis (GEPIA). Finally, we explored the diagnostic value of the candidate gene using receiver operating character (ROC) analysis and the relationship between the candidate gene and physical immunity using the tumor immune estimation resource (TIMER) database. Based on the results of the aforementioned analysis, we found that the DEGs, influenced by TAZ/YAP, took part in the common regulatory interaction networks (TGF-β signaling pathway and endocytosis) and identified an immune-related novel biomarker (FMR1) that was positively correlated with immune infiltration and patients' prognosis. This study may provide new insight into the understanding of molecular mechanisms and the clinical prognosis of KIRC.

2.   Materials and methods

The detailed process of the research was exhibited in Figure 1. The dataset, software, and website of the study were shown in Table 1.

Figure 1.  The analyzed flowchart of the research.

DownLoad: Full-Size Img PowerPoint

Table 1.  The information of databases, software, and website.

Name	Source
Cytoscape v3.8.2	http://cytoscape.org/
GEO	https://www.ncbi.nlm.nih.gov/geo/
GEPIA	http://gepia.cancer-pku.cn/
GDC	https://portal.gdc.cancer.gov/
GSEA v4.1.0	https://www.gsea-msigdb.org/gsea/index.jsp
HPA	https://www.proteinatlas.org/
MSigDB	https://www.gsea-msigdb.org/gsea/msigdb/index.jsp
RStudio v4.1.0	https://www.r-project.org/
STRING	https://www.string-db.org/
SangerBox	http://www.sangerbox.com/
TIMER	https://cistrome.shinyapps.io/timer/
UCSC Xena	https://xenabrowser.net/
UALCAN	http://ualcan.path.uab.edu/

---

 | Show Table

DownLoad: CSV

2.1.   Selecting database and recombining microarray data

The gene expression omnibus (GEO) database is an international public repository that is easy to flexibly and freely obtain genomics data ^[53]. In this study, we downloaded the microarray datasets of GSE121689 (annotation information: Affymetrix human genome GPL571 platform) and GSE146354 (annotation information: Affymetrix human genome GPL570 platform) from GEO. GSE121689 was associated with TAZ and had 12 samples that were treated with DMSO, siRNA, erastin, and siRNA-erastin, successively. GSE146354 was correlated with YAP and also had 12 samples that the processing method was the same as GSE121689. Here, we want to explore the co-regulatory mechanism of TAZ/YAP and identify novel biomarkers. According to our intention, we selected 12 samples that were handled with DMSO (as control) and siRNA in GSE121689 and GSE146354, respectively. Then, we recombined gene expression data for subsequent analysis.

2.2.   Processing data and identifying co-expressed DEGs

The aforementioned microarray data was converted probe IDs to the corresponding gene symbols and was analyzed using RStudio v4.1.0. The Linear Models for Microarray data (limma) package was utilized to identify DEGs ^[54]. The filtered conditions were referenced in the related literature ^[55]. The results were saved as CSV text and were applied to draw volcano plots and heatmaps using packages of ggplot2 and pheatmap. Subsequently, the Venn diagram was drawn with SangerBox which was a simple visualization webpage tool to identify the overlapping genes.

2.3.   Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis of DEGs

The clusterProfiler, an R package, was applied to execute GO and KEGG enrichment analyses to predict the functional roles of the co-expressed DEGs ^[56]. Among them, GO enrichment analyses consist of three aspects, which were described in the reference ^[57].

2.4.   Constructing a PPI network and identifying the hub genes

To construct networks, the search tool for the retrieval of interacting genes (STRING) database was employed to analyze and predict the interactions of DEGs ^[58]. Cytoscape software v3.8.2 was applied to visualize the prediction network of DEGs proteins ^[59]. The threshold of the combined score was set as 0.4. The molecular complex detection (MCODE) and cytoHubba, plugins of Cytoscape, were used to filter the significant module and screen the nodes of the highest score and strongest connectivity to serve as the hub genes. The parameters of MCODE have referenced the related research ^[60].

2.5.   Performing gene set enrichment analysis (GSEA) to identify the involved pathways of the hub genes

GSEA software v4.1.0 is applied to analyze the potential biological pathways of the different phenotypes based on statistically significant and concordant differences in RNA-seq profiles ^[61]. The KIRC expression profile (535 cancer samples) was downloaded from the UCSC Xena functional genomics explorer (UCSC Xena) and was separated into two groups according to the median level of the hub genes. Single-gene GSEA was executed to obtain their biological pathways by ranking the Spearman correlation coefficient of the hub gene and genes in the gene list. The selection of reference was the annotated genes set (c2.cp.kegg.v7.4.symbols.gmt) in Molecular Signatures Database (MSigDB). The threshold was set based on the related study ^[62].

2.6.   Validation and OS and disease-free survival (DFS) analysis of the hub genes

The HPA is an opening database that is intended to map proteins in organs, tissues, and cells via various omics technologies ^[63]; UALCAN is a simply and conveniently online analyses tool to carry out the cancer genome atlas (TCGA) data mining ^[64]; while GEPIA is a formidable analysis website to perform some functions that are the analysis of gene expression, correlation, and survival based on the combination of TCGA database and genotype-tissue expression (GTEx) database ^[65]. In this research, immunohistochemistry images of the hub genes were obtained from HPA. The corresponding antibodies were CAB020836 (STAU2), HPA050118 (FMR1), and CAB019309 (AGO2). Their expressed data were obtained from the "expression" module and "KIRC" dataset of UALCAN. The "single gene analysis" module of GEPIA was used to execute the analyses of the pathological stage, OS, and DFS.

2.7.   The ROC curves of the candidate genes

TCGA is an open-access database and collects a total of 33 cancer omics data, which is usually used by cancer-related researchers. Here, we downloaded the clinical data of KIRC from the genomic data commons (GDC). The collation of the data was executed by RStudio v4.1.0 to obtain a corresponding CSV text, which was analyzed and was visualized by the R packages of pROC and ggplot2, respectively. The area under the curve (AUC) of ROC represents the predicted accuracy of the candidate gene. The 0.5–1 of AUC is considered to predict sensitivity, which is divided into three degrees; including the low accuracy (0.5–0.7), the moderate accuracy (0.7–0.9), and the high accuracy (0.9–1) ^[66].

2.8.   Immune cell infiltrates analysis of the candidate gene

TIMER is a webpage tool to systematically evaluate the relationships between gene and immune cell infiltration in diverse cancers ^[67]. In this research, the "Gene module" was applied to evaluate the relation between FMR1 and immune cells infiltration; the "SCNA module" was employed to assess the correlation between FMR1 mutation and immune cells infiltration; the "survival module" was utilized to explore the relationships between clinical outcomes and immune infiltrates or FMR1 expression; the "correlation module" was performed to identify the relevance between immune biomarkers and FMR1 expression.

2.9.   Statistical analysis

All data was conducted in RStudio v4.1.0 and corresponding bioinformatic webpage tools. The statistical differences were calculated by Student's t-test, one-way analysis of variance (ANOVA), Spearman correlation analysis, cox regression analysis, and Kaplan-Meier analysis, respectively. The threshold was P-value < 0.05; *P < 0.05, **P < 0.01 and ***P < 0.001.

3.   Results

3.1.   Choosing the DEGs and identifying the co-expressed DEGs based on GSE121689 and GSE146354

The data of GSE121689 and GSE146354 were normalized and analyzed using the limma package. We obtained a total of 165 DEGs (99 upregulated and 66 downregulated) in the GSE121689 dataset (Table S1); similarly, a total of 341 DEGs (117 upregulated and 224 downregulated) were selected in GSE146354 (Table S2). The volcano plots and heatmaps of DEGs were shown in Figure 2A–D. The Venn diagram was utilized to identify and visualize 51 overlapping DEGs (46 upregulated and 5 downregulated) that were shown in Figure 2E, F and Table S3, respectively.

Figure 2.  Identification of overlapping DEGs in the GSE121689 and GSE146354 database. (A, B) The volcano plots of DEGs; the upregulated genes and downregulated genes were shown in red and blue, respectively. (C, D) The heatmap of DEGs. (E, F) The Venn plots of overlapping upregulated and downregulated genes in two datasets, respectively. Log |fold change (FC)| > 1 and ^*P < 0.05.

DownLoad: Full-Size Img PowerPoint

3.2.   GO and KEGG functional enrichment analysis of co-expressed DEGs

To comprehensively understand the biological functions and pathways of co-expressed DEGs, we executed the GO and KEGG analysis using clusterProfiler. The KEGG results indicated that the DEGs were significantly enriched in endocytosis and TGF-β signaling pathway (Figure 3A and Table S4). Meanwhile, the GO results mainly displayed the following three aspects: in the biological processes (BP) cluster, the enrichment terms were negative regulation of transport, signal release, negative regulation of secretion, endocrine process, regulation of receptor-mediated endocytosis, and negative regulation of exocytosis (Figure 3B and Table S4); while in the molecular functions (MF) cluster, the enrichment terms were receptor ligand activity, signaling receptor activator activity, hormone activity, growth factor activity, and structural molecule activity conferring elasticity (these functions were associated with the characteristics of TGF-β signaling pathway ^[68]) (Figure 3C and Table S4); in the cellular components (CC) cluster, the DEGs diversifications mainly focused on collagen-containing extracellular matrix and synapse components (asymmetric synapse, neuron to neuron synapse, presynapse, axon terminus, neuron projection terminus, distal axon, and postsynaptic density) (Figure 3D and Table S4). Subsequently, we also explored the interacting strength of signaling pathways of BP, MF, and CC, which found that there were different degrees of connection of the aforementioned pathways (Figure 3E–G). The aforementioned results concluded that the DEGs had enormous potential to regulate the activity of endocytosis and TGF-β signaling pathway to control the cellular procession.

Figure 3.  GO and KEGG analysis results of 51 co-expressed DEGs. (A) The KEGG results of the DEGs. (B) The GO BP results of the DEGs. (C) The GO MF results of the DEGs. (D) The GO CC results of the DEGs. (E–G) The interaction of pathways of BP, MF, and CC; the degrees of the line-thickness represented the interacting strength of pathways. ^*P < 0.05.

DownLoad: Full-Size Img PowerPoint

3.3.   Constructing a PPI network and identifying the hub genes

To intuitively understand and explore the interactive relationships of co-expressed DEGs, we performed the PPI network analyses using the STRING database and Cytoscape software. In this network, there were 14 nodes and 11 edges; 14 DEGs exhibited the interaction, including 10 upregulated genes (SLC4A8, STAU2, FMR1, CPEB3, DESI2, DENND5B, MFSD6, CSGALNACT2, ACT, and CFH) and 4 downregulated genes (NUPL1, AGO2, INHBE, and ADM) (Figure 4A). Then, the CytoHubba and MCODE of Cytoscape plugins were employed to confirm the hub genes; in the CytoHubba module, there were six genes which were AGO2, FMR1, STAU2, CPEB3, ADM, and MFSD6 (Figure 4B); in the MCODE module, there were three genes which were STAU2, FMR1, and AGO2 (Figure 4C). Finally, we found that the hub genes were STAU2, FMR1, and AGO2 according to the exhibition of the Venn diagram (Figure 4D and Table 2).

Figure 4.  The PPI network of co-expressed DEGs and the confirmation of hub genes. (A) The PPI network of 14 DEGs; red nodes are upregulated genes, whereas blue nodes are downregulated genes. (B) The top six genes were filtered by CytoHubba; the variational colors of nodes represent the different statistical significance. Red is the highest significance, orange is the middle significance, and yellow is the lowest significance. (C) The significant module was identified by MCODE, including three genes. (D) The Venn diagram was applied to screen the hub genes in MCODE and CytoHubba.

DownLoad: Full-Size Img PowerPoint

Table 2.  The Venn analysis results of the MCODE and CytoHubba.

Names	Total	Elements	Common Elements
MCODE	3	STAU2, FMR1, AGO2	STAU2, AGO2, FMR1
CytoHubba	6	CPEB3, FMR1, AGO2, STAU2, ADM, MFSD6

---

 | Show Table

DownLoad: CSV

3.4.   The biological signaling pathways of the hub genes

To further elucidate the biological networks and pathways of the aforementioned hub genes in KIRC, we executed a KEGG analysis using GSEA. The results indicated that the hub genes were associated with metabolic pathways (ubiquitin-mediated proteolysis, lysosome, endocytosis, valine leucine and isoleucine degradation, and peroxisome) and physical immune system (allograft rejection, graft versus host disease, systemic lupus erythematosus, primary immunodeficiency, intestinal immune network for IGA production, B cell receptor signaling pathway, T cell receptor signaling pathway, and Toll-like receptor signaling pathway) (Figure 5A–C). Based on enriched similar functional terms, we found that the overexpression of STAU2 and AGO2 made the active states of ubiquitin-mediated proteolysis and endocytosis (Figure 5A, B), while AGO2 and FMR1 enhanced the immune system activity to probably cause immune function disorder in the high expression groups (Figure 5B, C). In conclusion, the abovementioned results demonstrated that the three genes and their expression-related genes had affected physical functions by regulating the procession of metabolism and immune; especially, FMR1 paid high attention to the regulation of the immune system.

Figure 5.  The KEGG-related pathways of three hub genes. (A) The GSEA results of STAU2. (B) The GSEA results of AGO2. (C) The GSEA results of FMR1. |NES| > 1, NOM ^*P-value < 0.05, and FDR q-value < 0.25.

DownLoad: Full-Size Img PowerPoint

3.5.   Validation and survival analysis of the hub genes to determine FMR1 as a candidate gene

To confirm the reliable roles of the hub genes in KIRC development, the datasets of HPA, UALCAN, and GEPIA were utilized to explore the relationships between expressional changes of the hub genes and the prognosis and clinical outcomes of the patients. The immunohistochemical results of the hub genes in the HPA database showed that STAU2 and FMR1 were high-expressed in the cytoplasm and membrane of normal kidney cells, while STAU2 and FMR1 were weakly expressed in the nucleus and cytoplasm/membrane of tumor cells, respectively (Figure 6A, B). However, the expression level of AGO2 was unchanged in the cytoplasm and membrane of normal and KIRC tissues (Figure 6C). Similarly, the UALCAN boxplot results indicated that STAU2 and FMR1 were significantly downregulated in the tumor, whereas AGO2 variation was not statistically significant (Figure 7A). Meanwhile, we found that STAU2 and FMR1 were negatively associated with KIRC stages (Figure 7B). Subsequently, we also performed the Kaplan-Meier curve and log-rank test analyses to exhibit the correlation between the hub genes (STAU2, FMR1, and AGO2) and OS and DFS of patients. The results showed that only FMR1 was meaningful for patients' OS and DFS. The deregulated FMR1 was correlated with poor OS and DFS (Figure 7C, D). FMR1 may be served as an anti-oncogene.

Figure 6.  The expressions of the hub genes were verified by the HPA database in KIRC. (A) The immunohistochemistry images of STAU2. (B) The immunohistochemistry images of FMR1. (C) The immunohistochemistry images of AGO2.

DownLoad: Full-Size Img PowerPoint

Figure 7.  The relationships between three hub genes and clinical characteristics and prognostic outcomes of KIRC patients. (A) The differential expression of STAU2, FMR1, and AGO2 between adjacent normal tissues and tumor tissues. (B) The clinical stages of STAU2, FMR1, and AGO2. (C) The OS curves of STAU2, FMR1, and AGO2. (D) The DFS curves of STAU2, FMR1, and AGO2. ^*P < 0.05, ^**P < 0.01, and ^***P < 0.001.

DownLoad: Full-Size Img PowerPoint

3.6.   The ROC curves of FMR1

To identify the diagnostic value of FMR1 in distinguishing KIRC patients from the normal, we executed the ROC analysis using pROC. The result showed that the AUC of FMR1 was 0.645 and the confidence interval (CI) was 0.59–0.7 (Figure 8), which implied that FMR1 had a certain predictive ability.

Figure 8.  The ROC analysis of FMR1 (normal vs tumor).

DownLoad: Full-Size Img PowerPoint

3.7.   FMR1 serves as an immune-related gene in KIRC

Based on the above-mentioned GSEA and prognosis analysis results, FMR1 was probably a potential prognostic biomarker and had a strong influence on immunomodulation in KIRC. Therefore, we performed the co-expressed analysis to explore the relationship between FMR1 and immune cell infiltration using TIMER. The results demonstrated that FMR1 was positively correlated with the B cells (cor = 0.241, P = 1.77e-07), CD8⁺ T cells (cor = 0.135, P = 4.69e-03), CD4⁺ T cells (cor = 0.356, P = 3.59e-15), macrophages (cor = 0.437, P = 2.41e-22), neutrophils (cor = 0.395, P = 1.64e-18) and dendritic cells (DCs, cor = 0.256, P = 3.18e-08) (Figure 9A). In 57 immune cell biomarkers, FMR1 had co-expressed relationships with 33 immune markers, including 27 positive correlations and 6 negative relations (Figure 9D and Table S5). Moreover, we found that FMR1 mutation negatively regulated immune cell infiltration (Figure 9B). The low expression of CD8⁺ T cells, macrophages, and FMR1 were correlated with poor prognosis through survival analysis (Figure 9C). In addition to the above analysis, we also selected more than 40 common immune checkpoint genes to perform the co-expressed analysis. The results indicated that there were 31 co-expressed connections, including 23 positive correlations and 8 negative correlations (Figure 9E). Some research revealed that the chemokine-chemokine receptor interactions promoted the immune cells' enrichment of the tumor microenvironment to influence tumor development ^[69]. Subsequently, we discussed the co-expressed relationships between FMR1 and chemokines or chemokine receptors. The results exhibited that FMR1 was associated with the expression of the chemokines (21, 11 positive correlations and 10 negative correlations) (Figure 9F) and commonly positively correlated with the chemokine receptors' expression (14, 12 positive correlations and 2 negatively correlations) (Figure 9G). In summary, FMR1 plays a vital role in immune infiltration to affect the KIRC progression and is a potential immune-related prognostic gene.

Figure 9.  Immunological analysis of FMR1. (A) The association between FMR1 and immune cell infiltration. (B) The relation between FMR1 mutation and immune cell infiltration. (C) The survival analysis of FMR1 and immune cells. (D) The heatmap of co-expressed analysis of FMR1 and immune cells biomarkers. (E) The heatmap of co-expressed analysis of FMR1 and immunological checkpoint genes. (F) The heatmap of co-expressed analysis of FMR1 and chemokines-related genes. (G) The heatmap of co-expressed analysis of FMR1 and chemokine receptors related genes. ^*P < 0.05, ^**P < 0.01, and ^***P < 0.001.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

TAZ/YAP is the most common oncogenes in various cancers; while there are double functions of TAZ/YAP in RCC, including the promotor and inhibitor of cancer. Mechanically, the overwhelming research indicated that TAZ/YAP promoted cell proliferation, migration, and invasion to influence the process of cancer ^[70,71]. Furthermore, TAZ and YAP could regulate the expression of the related immune molecule (programmed cell death ligand 1 (PDL1) ^[72]) or related angiogenic growth factor (vascular endothelial growth factors A (VEGFA) ^[73]), respectively; thereby creating a suitable microenvironment for the aforementioned pathological process. It is indubitable that the nuclear accumulation of TAZ/YAP did contribute to the formation of the aggressive RCC ^[74]. Interestingly, some studies demonstrated that TAZ/YAP prevented the above-mentioned functions by regulating F-actin, P53, apoptosis, and ferroptosis ^{[27,29,75,76]}. While a combined TAZ/YAP analysis revealed that they influenced the cell metabolism process to promote the RCC development of NF2 mutation ^[32]. However, there isn't a combined analysis of TAZ/YAP to explore their common regulatory mechanisms in mutation-free RCC, especially in KIRC.

Bioinformatics essentially is an analyzed method that uses related computer technology to process genomics data. It theoretically can screen the diagnostic and prognostic biomarkers and explain the corresponding mechanism of diseases, which is convenient for developing effective diagnosis and treatment measures. Recently, emerging scholars have paid attention to the development of the related analytical tools with more effective and accurate, so that the analytical methods have become more diverse and simpler. Presently, the R packages of limma ^[54], edgeR ^[77], and DESeq2 ^[78,79] were the common methods to screen DEGs; the webpage tools of DAVID ^[80], KOBAS ^[81], Metascape ^[82], and WebGestalt ^[83] could be applied to enrichment analysis to obtain the corresponding signaling pathways; the tools of STRING ^[58], LncBase ^[84], miRTarBase ^[85], Starbase ^[86], and Cytoscape ^[59] were used to predict and construct the interaction of gene-gene; and the associated computational models (IMBDANET ^[87], lncRNA–miRNA association prediction (NDALMA) ^[88], logistic matrix factorization with neighborhood regularized (LMFNRLMI) ^[89], deep forest ensemble learning based on autoencoder (DFELMDA) ^[90], and the computational models of predicted circRNA-diseases ^[91]) were proposed to predict the correlation between gene and gene or diseases. These methods revealed the related regulatory networks and the corresponding molecular mechanisms form different aspects, which made us have comprehensive understanded to the pathological process of diseases. Some above-analyzed results have confirmed by the relevant experiments ^[92,93,94], which verified the accuracy of the bioinformatic results.

In our study, we performed the comprehensive conjoint analysis of TAZ/YAP utilizing bioinformatic methods to reveal their co-regulatory mechanisms and identify new prognostic biomarkers in KIRC. Here, we identified 51 co-expressed DEGs and three hub genes (STAU2, FMR1, and AGO2), but only FMR1 was consistent with the clinical characteristics and prognosis outcomes of patients. The KEGG terms of the co-expressed DEGs were TGF-β signaling pathway and endocytosis, which was supported by the GO results. Research indicated that the TGF-β signaling pathway took part in the regulation of complicated mechanical signal networks and exhibited double functions (tumor-suppressing and tumor-promoting) to maintain physical homeostasis or promote tumor development by influencing cell proliferation, migration, invasion, drug resistance, and tumor escape ^[95]; while endocytosis was an internalized process that involved in absorbing nutrients, signaling proteins, and other cargoes, which were correlated with cellular metabolism. This study showed that TAZ/YAP commonly regulated TGF-β signaling pathway and endocytosis to control cellular procession in KIRC.

FMR1 is a multifunctional polyribosome-associated RNA-binding protein that is involved in RNA maturity and stability ^[96]. Its epigenetic transformation is tightly associated with the progression of cancers ^{[97,98,99,100]}, neurodegenerative disorders, and immunomodulation dysfunctions ^{[101,102,103]}. In GO results, FMR1 was involved in some GO pathways, including the negative regulation of transport, signal release, negative regulation of secretion, regulation of receptor−mediated endocytosis, negative regulation of exocytosis, and synapse components (Figure 10A, B and Table S4); which indicated that this gene was closely associated with the regulation of neuroendocrine functions and the development of the nervous system. In addition, some research showed that FMR1 created an abnormal immune microenvironment by affecting the production of cytokines (IL-β, IL-6, TNF-α, and IL-10) ^[104,105] and the surface biomarkers activity of immune cells (T cells) ^[103], which was consistent with its GSEA results. Therefore, we performed the co-expressed analysis of this gene and immune cells, immune checkpoints, chemokines, and chemokine receptors. We found that the gene was closely correlated with the aforementioned biomarkers. Meanwhile, we also found that FMR1 was significantly deregulated in KIRC tissue and exhibited a negative relationship with the stage of patients; the poor prognosis was associated with its weak expression; the ROC result indicated that FMR1 owned a certain predictive sensibility to distinguish KIRC patients from the normal. These results implied that FMR1 could be served as a tumor suppressor gene and a potential novel immune-related prognostic molecule in KIRC, which provided a new target and theoretical basis for the KIRC treatment strategy. Although the above conclusions were confirmed by the public databases, more experiments are still needed for further analysis and validation.

Figure 10.  The connection between the pathways and genes in GO enrichment pathways. (A) The relation between the BP pathways and genes. (B) The association between the CC pathways and genes.

DownLoad: Full-Size Img PowerPoint

5.   Conclusions

In conclusion, TAZ/YAP had influenced 51 common DEGs to take part in the regulation of TGF-β signaling pathway and endocytosis, thereby controlling the KIRC development procession. Even though we identified three hub genes (STAU2, FMR1, and AGO2), only FMR1 was consistent with the clinical characteristics and prognosis outcomes of the patients. We found that FMR1 was downregulated in KIRC tissues, exhibited an unfavorable prognosis, and had a certain predictive sensibility to distinguish KIRC patients from the normal. Meanwhile, FMR1 was tightly associated with immunomodulation. The prementioned results implied that FMR1 could be served as an antitumor gene and was a potential novel immune-related prognostic molecule in KIRC.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (81802785 [Y. Jiang]), Hunan Provincial Natural Science Foundation of China (2020JJ5382 [Y. Jiang], 2020JJ8051 [J. Huang]), Scientific Research Program of Hunan Provincial Health Commission (202202065445 [Y. Jiang]). We sincerely thank Ph.D. Chanjuan Zheng for her suggestions in this research. We thank all participants.

Conflict of interest

The authors declare that there are no potential conflicts of interest.