Lugang Deng1,2
Peixi Wang- 1School of Public Health, Guangzhou Medical University, Guangzhou, China
- 2South China Hospital, Health Science Center, Shenzhen University, Shenzhen, China
- 3Institute of Chronic Disease Risks Assessment, School of Nursing and Health, Henan University, Kaifeng, China
Background: Kidney renal clear cell carcinoma (KIRC) has the highest invasion, mortality and metastasis of the renal cell carcinomas and seriously affects patient’s quality of life. However, the composition of the immune microenvironment and regulatory mechanisms at transcriptomic level such as ceRNA of KIRC are still unclear.
Methods: We constructed a ceRNA network associated with KIRC by analyzing the long non-coding RNA (lncRNA), miRNA and mRNA expression data of 506 tumor tissue samples and 71 normal adjacent tissue samples downloaded from The Cancer Genome Atlas (TCGA) database. In addition, we estimated the proportion of 22 immune cell types in these samples through “The Cell Type Identification by Estimating Relative Subsets of RNA Transcripts.” Based on the ceRNA network and immune cells screened by univariate Cox analysis and Lasso regression, two nomograms were constructed to predict the prognosis of patients with KIRC. Receiver operating characteristic curves (ROC) and calibration curves were employed to assess the discrimination and accuracy of the nomograms. Consequently, co-expression analysis was carried out to explore the relationship between each prognostic gene in a Cox proportional hazards regression model of ceRNA and each survival-related immune cell in a Cox proportional hazards regression model of immune cell types to reveal the potential regulatory mechanism.
Results: We established a ceRNA network consisting of 12 lncRNAs, 25 miRNAs and 136 mRNAs. Two nomograms containing seven prognostic genes and two immune cells, respectively, were successfully constructed. Both ROC [area under curves (AUCs) of 1, 3, and 5-year survival in the nomogram based on ceRNA network: 0.779, 0.747, and 0.772; AUCs of 1, 3, and 5-year survivals in nomogram based on immune cells: 0.603, 0.642, and 0.607] and calibration curves indicated good accuracy and clinical application value of both models. Through co-correlation analysis between ceRNA and immune cells, we found both LINC00894 and KIAA1324 were positively correlated with follicular helper T (Tfh) cells and negatively correlated with resting mast cells.
Conclusion: Based on the ceRNA network and tumor-infiltrating immune cells, we constructed two nomograms to predict the survival of KIRC patients and demonstrated their value in improving the personalized management of KIRC.
Introduction
Renal cell carcinoma is a malignant kidney tumor that the American Cancer Society reports was diagnosed in 65,340 patients, and caused 14,970 deaths, in the United States during 2018 (Siegel et al., 2018). It is a serious health problem worldwide, and the economic burden of the disease has steadily increased during the last century (Hsieh et al., 2017). Kidney renal clear cell carcinoma (KIRC) is the subtype of renal cell carcinoma with the highest invasiveness, mortality and metastasis rates (Xiao et al., 2020); it is typically screened and diagnosed through computed tomography (CT), magnetic resonance imaging (MRI), and pathological section testing in the clinic (Morrissey et al., 2015). Several difficulties continue to make KIRC a challenge to treat, including its stealth in early stages, a lack of effective biomarkers, radiation risks and the high cost of diagnostic imaging required. Underlying biomarkers and molecular mechanisms for the prevention, diagnosis and prognosis of KIRC require further studies.
Long non-coding RNA (lncRNA) is a type of RNA that contains more than 200 nucleotides but does not have the function of encoding proteins. Due to their potential for exploring novel biomarkers in diseases and elucidating mechanisms of biological processes, they have gradually achieved attention recently (Nagano and Fraser, 2011; Quinn and Chang, 2016). Current research reveals that lncRNA plays an important role in the genesis and development of tumors through interactions with competing endogenous RNA (ceRNA) or by acting as microRNA sponges (Tay et al., 2014). Tumor-infiltrating immune cells are a major component of the tumor microenvironment and affect the clinical outcomes and pathological staging of multiple tumor types (Gao et al., 2020; Hong et al., 2020). Some studies suggest lncRNA can affect tumor development and tumor immune cell microenvironment through the ceRNA network (Huang et al., 2018; Cao et al., 2019; Xu et al., 2019). Therefore, exploring interactions between the ceRNA network and various immune cells in the development of KIRC is essential.
In this study, we construct a ceRNA network related to the occurrence and development of KIRC using transcriptome and clinical data from The Cancer Genome Atlas (TCGA)1 to explore potential molecular mechanisms. In addition, we used ‘‘The Cell Type Identification by Estimating Relative Subsets of RNA Transcripts algorithm’’ (CIBERSORT)2 (Newman et al., 2015) algorithm to assess differences in the composition of immune cells within tumor tissues and normal tissues. Furthermore, predictive nomograms were constructed based on the ceRNA network and on significant immune cells. Finally, we evaluated the relationship between tumor-infiltrating immune cells and the ceRNA network to identify the underlying regulatory mechanisms.
Materials and Methods
Data Source and Selection
A flowchart illustrating the research process of our study is provided in Figure 1. Data on lncRNA, miRNA, and mRNA expression of KIRC patients were downloaded from TCGA database (see text foot note 1) through the GDC data transfer tool. Relevant clinical information including age, gender, survival data (vital status, days to death, and days to last follow-up), grade and TNM stage of tumor tissues was also obtained for further analysis. We excluded samples with incomplete survival or pathological staging data, as well as duplicate samples. Only samples with both mRNA and miRNA expression profiles were included in this study. Correspondingly, the adjacent and paired normal tissues collected from the included patients was set as the control group.
Figure 1. Flowchart of the analysis process.
Differential Expression Analysis
Based on the HTSeq-Counts data of included samples, we used the “DESeq2” package (Love et al., 2014)3 of Bioconductor to screen for differentially expressed lncRNAs (DElncRNAs), miRNAs (DEmiRNAs) and mRNAs (DEmRNAs) between KIRC tissue and normal adjacent tissue. The log fold change (FC) criterion for differential expression was set as |log2(FC)| > 1 and the false discovery rate (FDR) at < 0.05. For the differential expression results of each type of RNA, volcano plots and heatmaps were generated through R packages “ggpolt2” and “pheatmap.”
Construction of lncRNA-miRNA-mRNA ceRNA Network and Survival Analysis
Using the differential expression results, we predicted the miRNAs-lncRNAs interactions and the miRNAs-mRNAs interactions based on the starBase v2.0 database (Yang et al., 2011; Li et al., 2014),4 which contained the miRNA-target interactions overlapped with CLIP-Seq data. We conducted hypergeometric tests and correlation analysis to filter the predicted interactions with FDR < 0.05 as the screening criteria. In addition, the lncRNA-miRNA-mRNA regulatory network of KIRC was visualized using Cytoscape v.3.8.15 (Shannon et al., 2003). Besides, we conducted a Kaplan–Meier (K–M) survival curve analysis to assess the prognostic value of genes in the KIRC ceRNA network.
ceRNA Network: Cox Proportional Hazards Regression Model
Firstly, significant variables were screened for integration into the initial Cox model using univariate Cox analysis. Next, a least absolute shrinkage and selection operator (Lasso) regression was employed to further evaluate the fitness of multifaceted models and filter variables. Finally, we generated a nomogram of from multivariate Cox proportional hazards regression models to predict the prognosis of KIRC patients. To confirm the discrimination and precision of the nomogram, receiver operating characteristic curves (ROC) and calibration curves were employed. In addition, we calculated the risk score of each patient and then divided patients into high and low risk groups based on the median to conduct a K–M survival curve analysis.
Tumor-Infiltrating Immune Cells: CIBERSORT Estimation and Survival Analysis
To analyze differences between the two groups at the level of immune cells, and explore the connection between key biomarkers in the ceRNA network and immune cells, we estimated the proportion of 22 kinds of immune cell types in 506 tumor tissues and 71 normal adjacent tissues by CIBERSORT (Newman et al., 2015), an algorithm for characterizing cell composition of complex tissues from their gene expression profiles. Only samples with CIBERSORT output of P < 0.05 were included in further analysis. After screening samples, we searched for immune cells with significant differences between cancer tissues and normal adjacent tissues, as indicated by a Wilcoxon rank-sum test. A K–M survival curve analysis was performed to evaluate the relationship between the content of different immune cells and the overall survival of KIRC patients.
Tumor-Infiltrating Immune Cells: Cox Proportional Hazards Regression Model
Immune cells with significant effects in univariate Cox analysis were integrated into the initial Cox model and the Lasso regression again used to evaluate the fitness of multifaceted models and filter variables further. As above, we constructed a nomogram based on the multivariate model of immune cells, and used ROC and calibration curves to confirm the discrimination and precision of the nomogram. Similarly, we conducted a K–M survival curve analysis by dividing patients into high and low risk groups according to the median of risk score. Finally, a Pearson correlation analysis was used to explore the relationship between genes and immune cells.
Statistical Analysis
Only two-sided P-values < 0.05 was considered statistically significant. All statistical analyses were conducted in R version 4.0.2 software (Institute for Statistics and Mathematics, Vienna, Austria)6 (packages: corrplot, DESeq2, GDCRNATools, ggplot2, glmnet, pheatmap, rms, survival, survminer, and timeROC).
Results
Differential Expression in lncRNAs, miRNAs, and mRNAs
After screening the samples of TCGA-KIRC, we constructed an experimental group and a control group containing of 506 tumor tissues and 71 normal adjacent tissues, respectively. Compared to the control group, we identified 357 DElncRNAs (287 up-regulated and 70 down-regulated), 132 DEmiRNAs (61 up-regulated and 71 down-regulated) and 3092 DEmRNAs (2,032 up-regulated and 1,060 down-regulated) using the cutoffs of the |log2(FC)| > 1 and FDR < 0.05 (Figures 2A–G).
Figure 2. Genes differentially expressed in 506 KIRC tissues versus 71 normal adjacent tissues. Provided are the heatmap (A) and the volcano plot (B) of 357 differentially expressed lncRNAs between the two groups; the heatmap (C) and the volcano plot (D) of 132 differentially expressed miRNAs between the two groups; and the heatmap (E) and the volcano plot (F) of 3092 differentially expressed protein-coding genes between the two groups; and the composition of differentially expressed genes (G). All differentially expressed genes had a log (fold-change) > 1.0 or < −1.0 and the overall of FDR < 0.05.
Construction of lncRNA-miRNA-mRNA ceRNA Network and Survival Analysis
Using the starBase v2.0 database and ceRNA theory, a lncRNA-miRNA-mRNA ceRNA network consisting of 173 genes (12 lncRNAs, 25 miRNAs and 136 mRNAs) was established (Figure 3A and Supplementary Table 1). A K–M survival curve analysis to assess the prognostic value of nodes in the ceRNA network revealed 83 genes (10 lncRNAs, 7 miRNAs, and 66 mRNAs) related to the overall survival of KIRC patients. Survival curves of the top three (according to P-value) lncRNAs (PVT1, AC005154.1, AC015813.1), miRNAs (hsa-miR-21-5p, hsa-miR-130b-3p, and hsa-miR-204-5p), and mRNAs (MXD3, VPS13D, and KCNN4) are shown in Figures 3B–J, respectively.
Figure 3. (A) The KIRC related ceRNA network consisting of 12 lncRNAs, 25 miRNAs and 136 mRNAs. (B–J) The K–M survival curves of the top three lncRNAs, miRNAs, and mRNAs: PVT1 (B), AC005154.1 (C), AC015813.1 (D), hsa-miR-21-5p (E), hsa-miR-130b-3p (F), hsa-miR-204-5p (G), MXD3 (H), VPS13D (I), and KCNN4 (J).
ceRNA Network: Cox Proportional Hazards Regression Model
Univariate Cox analysis applied to filter genes in the ceRNA network retained 97 genes for incorporation into the initial model (Supplementary Table 2). Lasso regression analysis indicated that 15 genes were statistically significant in this model. Subsequently, a Cox proportional hazards regression model with seven genes was established through further screening of key genes based on the Akaike Information Criterion (AIC). Finally, we constructed a nomogram to predict the 1, 3, and 5-year overall survival probability of KIRC patients. ROC [area under curves (AUCs)] for the model gave 1, 3, and 5-year survivals of 0.779, 0.747, and 0.772, respectively, which reflected the discrimination and precision of the nomogram with calibration curves (Figure 4). Overall survival (based on K–M curves) of the high-risk group is lower than that of the low-risk group (Supplementary Figure 1). Survival curves of the genes in the model are also provided in Supplementary Figure 1.
Figure 4. The Cox proportional hazards regression model (A) based on RNAs screened by Lasso regression analysis (B,C) in the ceRNA network. Seven potential prognosis related RNAs were integrated into the Cox proportional hazards regression model. The nomogram (E) was constructed based on the model. The ROC (D) and calibration curves (F) indicate the acceptable accuracy (AUCs of 1, 3, and 5-year survivals: 0.779, 0.747, and 0.772) and discrimination of the nomogram. *P < 0.05; **P < 0.010; ***P < 0.001.
Tumor-Infiltrating Immune Cells: CIBERSORT Estimation and Survival Analysis
The analysis by CIBERSORT identified five normal tissue samples and 206 tumor tissue samples for inclusion in subsequent analysis (P < 0.05) and the proportions of 22 immune cells in each sample are displayed in the histogram and heat map of Figure 5. A Wilcoxon rank-sum test showed that naive B cells (P < 0.001), plasma cells (P < 0.001), CD8 T cells (P = 0.014), CD4 naive T cells (P < 0.001), CD4 memory resting T cells (P = 0.013), regulatory (Tregs) T cells (P = 0.015), gamma delta T cells (P = 0.021), resting dendritic cells (P = 0.016), and resting mast cells (P = 0.048) all varied between the two tissue sample types (Figure 5C). Based on the results of K–M survival curve analysis, five immune cells, including memory B cells, plasma cells, follicular helper T (Tfh) cells, T regulatory (Tregs) cells, and resting mast cells were related to overall survival of KIRC patients (Supplementary Figure 2). The clinical correlation of 22 immune cells is displayed in Supplementary Figure 3.
Figure 5. The composition (A) and heatmap (B) of immune cells estimated by the CIBERSORT algorithm in five normal adjacent tissues and 206 KIRC tissues. The violin plot of immune cells (C) compares cells’ proportion between the two groups.
Tumor-Infiltrating Immune Cells: Cox Proportional Hazards Regression Model
After filtering based on univariate Cox analyses (Supplementary Table 3), two of 22 immune cells (Tfh and resting mast cells) were incorporated into the initial regression model. Lasso regression analysis and AIC were then employed for further optimization of the model, and both cells were kept in the final Cox proportional hazards regression model. A nomogram was developed to predict the 1, 3, and 5-year overall survival probability of KIRC patients based on the model. Both the ROC (AUCs of 1, 3, and 5-year survivals: 0.603, 0.642, and 0.607) and the calibration curves suggested that the nomogram had good accuracy (Figure 6). The K–M survival curves suggested that significant differences existed in the overall survival of the high versus low-risk group (Supplementary Figure 4).
Figure 6. The Cox proportional hazards regression model (A) based on immune cells screened by Lasso regression analysis (B,C). Tfh cells and resting mast cells were integrated into the Cox proportional hazards regression model. The nomogram (E) was constructed based on the model. The ROC (D) and the calibration curves (F) indicate the acceptable accuracy (AUCs of 1, 3, and 5-year survivals: 0.603, 0.642, and 0.607) and discrimination of the nomogram. *P < 0.05.
Co-expression Analysis of Key Genes and Immune Cells
Pearson correlation analysis was used to explore the co-expression pattern between 22 types of immune cells (Figure 7A). Similarly, the co-expression relationship between the key biomarkers and immune cells of the two models is illustrated in Figure 7B. Tfh cells are positively correlated with LINC00894 (Figure 7C, R = 0.36, P < 0.001) and KIAA1324 (Figure 7D, R = 0.39, P < 0.001), while resting mast cells are negatively correlated with LINC00894 (Figure 7E, R = −0.32, P < 0.001) and KIAA1324 (Figure 7F, R = −0.30, P < 0.001).
Figure 7. Co-expression heatmap of all immune cells in KIRC (A); the co-expression heatmap of the RNAs and immune cells in the two Cox proportional hazards regression models (B); the significant results of Pearson correlation coefficients between the RNAs and immune cells: Tfh cells and LINC00894 (C); Tfh cells and KIAA1324 (D); resting mast cells and LINC00894 (E); and resting mast cells and KIAA1324 (F).
Validation of the Correlation Between Key Genes and Immune Cells
Multiple databases were searched to verify the correlation between the key genes (LINC00894 and KIAA1324) and immune cells (Tfh cells and resting mast cells) identified. B cell lymphoma 6 (BCL6), CXCL13, CD10 (MME), ICOS, and PD-1 (PDCD1) were the most common surface markers of Tfh cells in the CellMarker database (Supplementary Figure 5). In the OncomineTM platform, we explored the expression of KIAA1324 and Tfh markers among various cancers. KIAA1324, BCL6, CXCL13, ICOS, and PDCD1 were highly expressed and MME was downregulated in KIRC compared to normal kidney tissue (Supplementary Figure 6); median rank and P-values are listed in Supplementary Table 4.
Using the gene expression profiling interactive analysis (GEPIA), KIAA1324 is positively correlated with CXCL13, ICOS, and PDCD1 in KIRC (Supplementary Figure 7), and LINC00894 is positively correlated with BCL6, ICOS, and PDCD1 in KIRC (Supplementary Figure 8). Similarly, the results of the LinkedOmics database showed that KIAA1324 was positively correlated with BCL6, CXCL13, ICOS, and PDCD1 in KIRC. In addition, the expression of KIAA1324 was significantly different among tumor stages and related to tumor purity in KIRC (Supplementary Figure 9).
Discussion
Kidney renal clear cell carcinoma is one of the most common and malignant subtypes of kidney tumors (Xiao et al., 2020). However, because of limitations in early detection and a lack of sensitive biomarkers, patients with KIRC typically have a poor prognosis. Recently, ceRNA was identified as playing a potentially important role in the molecular regulatory function of tumorigenesis and growth, and in affecting differences in immune-infiltrating cells within multiple tumors. In this study, we investigated the role and connection of the ceRNA network and immune infiltration in KIRC to explore key prognostic biomarkers and potential regulatory mechanisms.
We established a ceRNA network consisting of 12 lncRNAs, 25 miRNAs, and 136 mRNAs and compared the composition differences of immune cells related to KIRC. Based on the results, two prediction nomograms consisting of seven genes and two immune cells, respectively, were constructed to evaluate overall survival. A correlation analysis of the key factors in both nomograms indicated that both LINC00894 (lncRNA) and KIAA1324 (protein-coding RNA) were positively correlated with Tfh cells and negatively correlated with resting mast cells. Consequently, we infer that ceRNA regulatory function, in which both LINC00894 and KIAA1324 participate, and Tfh cells and resting mast cells may play a crucial role in the development and treatment of KIRC.
LINC00894 is lncRNA derived from a locus on the X chromosome that is differentially expressed in various cancers (Zhang et al., 2018). Similar to our work, Liang (Liu et al., 2018) reported that LINC00894 was a tumor-special lncRNA in KIRC that was upregulated in tumor tissues and correlated with overall survival time. Epithelial-to-mesenchymal transition (EMT) regulating miRNAs regulate drug resistance and cancer progression and metastasis (Park et al., 2008; Ward et al., 2013). LINC00894 can inhibit the TGF-β2/ZEB1 signaling pathway by acting as the sponge of EMT-regulating miR-200, or can be upregulated by ERα activation and positively regulate the expression of miR-200a-3p and miR-200b-3p, which results in induction of the TGF-β2/ZEB1 signaling pathway to reduce the occurrence of drug resistance in breast cancers (Zhang et al., 2018; Farhan et al., 2019; Huang et al., 2020). A large number of miRNAs are dysregulated in the development and treatment of KIRC, and affect tumor proliferation, metastasis and treatment resistance, etc. through degrading target mRNA (Qu et al., 2016; Zhai et al., 2018; Incorvaia et al., 2020). In our study, hsa-miR-130b-3p is the only miRNA in the model we have constructed related to patient prognosis. It showed high abundance in the serum of patients with Wilms tumor, and reflected valuable diagnostic potentials with an AUC of 0.94 (Ludwig et al., 2015). In addition, its up-regulation in hepatocellular carcinoma leads to a decrease in the expression of HOXA5, which in turn promotes tumor growth and angiogenesis, is associated with a poor prognosis (Liao et al., 2020). In the ceRNA network we constructed related to KIRC, a large number of lncRNA/miRNA/mRNA axes were revealed, and further experiments are needed to verify their function and mechanism.
KIAA1324, also known as EIG121, encodes a 1013 amino acid transmembrane protein that shows high sequence conservation among different species. In type I endometrial cancer (estrogen-related), the expression of KIAA1324 is upregulated, but it is downregulated in type II endometrial cancer (not estrogen-related), which is more malignant than type I (Deng et al., 2005). KIAA1324 may be a novel suppressor, being epigenetically downregulated in gastric cancer (Kang et al., 2015); in addition, its high expression in endometrial and pancreatic carcinoma predicts improved survival (Oh et al., 2006; Westin et al., 2009). However, increased expression of KIAA1324 in ovarian cancer is associated with a poor prognosis (Schlumbrecht et al., 2011). These studies indicate that KIAA1324 may exhibit different biological functions in various cancers. Furthermore, recent mechanistic studies reveal that KIAA1324 is localized to endosome-lysosome compartments and associated with autophagy, which helps protect cells from cell death by upregulating the autophagy pathway under unfavorable conditions, such as starvation or chemotherapy (Deng et al., 2010). Therefore, we infer that overexpressed KIAA1324 helps cancer cells proliferate and resist chemotherapy in KIRC.
In this study, we explored five among 22 types of immune cells related to patient prognosis, and screened out the Tfh cells and resting mast cells to construct a Cox proportional hazards regression model, which evaluation of its AUC suggests has clear clinical value. Gou et al. (Zhu et al., 2019) also explored the relationship between immune cells and renal cell carcinoma, and similar to our research, found Tfh cells were associated with poor prognosis and resting mast cells were positively associated with long term survival.
Under the action of myeloid professional antigen presenting cells (APCs), naive CD4 + T cells could differentiate into Tfh cells and non-Tfh effector cells (such as Th1, Th2, or Th17 cells) (Crotty, 2011). Therefore, Tfh cells are a specialized subset of CD4 + T cells that help B cells respond to antigens (Crotty, 2019), which is, in turn, depends on the expression of transcriptional factor B cell lymphoma 6 (BCL6) (Johnston et al., 2009). In the past decade, several studies have revealed a potential connection between Tfh cells and infiltrated tumors. Amé-Thomas et al. (2012) demonstrated that Tfh cells protect follicular lymphoma malignant B cells from spontaneous and rituximab-induced apoptosis by upregulating CD40L and IL-4, which might lead to a worse prognosis. In non-small cell lung cancer (NSCLC) tumor tissues, the proportion of Tfh cells is increased, suggesting that Tfh cells might play an important role in estimating the poor survival of NSCLC patients (Guo et al., 2017). In addition, Tfh cell abundance is known to improve prognosis of patients with breast cancer or colorectal cancer (Bindea et al., 2013; Gu-Trantien et al., 2013), which indicates the value of Tfh cells against tumors through immune responses.
Mast cells, both resting and activated, are large granulated innate immune cells found predominantly in sites between the host and its external environment, such as skin, respiratory mucosa, and the gastrointestinal tract (Frossi et al., 2017). The activation of mast cells is related to various stimuli received by numerous receptors on the cell surface, such as pathogens, neuropeptides, cytokines, growth factors, toxins, basic compounds, complement, immune complexes, etc. (da Silva et al., 2014). Recent studies have shown that the association between mast cells and cancer is two-sided, including both involvement in, and protection against, tumor progression. Xu et al. (Wang et al., 2020) constructed a prognosis-associated microRNA–mast cell network in lung adenocarcinoma, and found resting mast cell numbers were closely related to better overall survival and disease-free survival (DFS), while activated mast cells were related to poor survival. miR-30a and miR-550a are also involved in the regulation of immune response by acting as the promoter and suppressor of resting mast cells, respectively. In contrast, mast cells in hepatocellular carcinoma are mostly resting, inactivation of mast cells can promote immune escape, which is beneficial to tumor growth (Rohr-Udilova et al., 2018).
Our research inevitably has some limitations that should be recognized. First, our research data was obtained from public databases. The limited data means that the clinical pathological parameters are incomplete, which can lead to potential errors and biases. Secondly, we have not taken the heterogeneity of the immune microenvironment associated with the location of immune infiltration into consideration. In other words, the accuracy and applicability of the prediction models will be affected by the heterogeneity of tissue subtypes. In addition, the data series used in this study are all from Western countries, so applying the research conclusions to Asian countries has certain bias and may be inappropriate. As well, the lack of cell markers for resting mast cells leads to incomplete validation in our study. Lastly, this study is not a biological mechanism research, as it lacks experiments on the interaction mechanisms between ceRNA and immune cells. Despite limitations of the study, we did combine tumor infiltrating immune cells with a ceRNA network for analysis in KIRC, which not only constructed two normograms with clinical predictive value, but also identified a correlation between two key prognostic molecules and two prognosis-related immune cells, providing novel directions for future research.
Conclusion
Based on tumor-infiltrating immune cells and ceRNA networks, we constructed two nomograms to predict survival of KIRC patients. The high AUC values we obtained reflect the utility of this approach, which provides clinicians with more comprehensive clinical information through the models to improve individual management of KIRC patients. In addition, our study suggests that LINC00894, KIAA1324, Tfh cells and resting mast cells are significantly related to the prognosis of KIRC patients and we infer that LINC00894 and KIAA1324 might interact with Tfh cells and resting mast cells to affect the progression of KIRC. This study also provides new ideas for the pathogenesis and clinical treatment of KIRC.
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.
Author Contributions
ZQ and NL conceived the idea, designed the study, and contributed to the revision of the manuscript draft. LD conducted the data analysis and visualization. PW performed the data interpretation and literature collection. LD and NL wrote the manuscript. All authors read and approved the final manuscript.
Funding
This work was supported by National Natural Science Foundation of China (Nos. 81872584 and 81472941), the National 863 Young Scientist Program (No. 2015AA020940), the Key Scientific Research Project Plan of Henan Province (No. 21A330001), the Natural Science Foundation of Guangdong Province (No. 2016A030313138), the Key Projects of Guangzhou Science and Technology Program (No. 201704020056), the Interdisciplinary Research for First-class Discipline Construction Project of Henan University (No. 2019YLXKJC04), and the Scientific Research Project for University of Education Bureau of Guangzhou (No. 201831841).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.667610/full#supplementary-material
Footnotes
- ^ https://portal.gdc.cancer.gov/
- ^ http://cibersort.stanford.edu/
- ^ http://bioconductor.org/packages/DESeq2/
- ^ http://starbase.sysu.edu.cn/starbase2/index.php
- ^ http://cytoscape.org/
- ^ www.r-project.org
References
Amé-Thomas, P., Le Priol, J., Yssel, H., Caron, G., Pangault, C., Jean, R., et al. (2012). Characterization of intratumoral follicular helper T cells in follicular lymphoma: role in the survival of malignant B cells. Leukemia 26, 1053–1063. doi: 10.1038/leu.2011.301
Bindea, G., Mlecnik, B., Tosolini, M., Kirilovsky, A., Waldner, M., Obenauf, A. C., et al. (2013). Spatiotemporal Dynamics of Intratumoral Immune Cells Reveal the Immune Landscape in Human Cancer. Immunity 39, 782–795. doi: 10.1016/j.immuni.2013.10.003
Cao, J., Dong, R., Jiang, L., Gong, Y., Yuan, M., You, J., et al. (2019). LncRNA-MM2P Identified as a Modulator of Macrophage M2 Polarization. Cancer Immunol. Res. 7, 292–305. doi: 10.1158/2326-6066.Cir-18-0145
Crotty, S. (2011). Follicular Helper CD4 T Cells (TFH). Annu. Rev. Immunol. 29, 621–663. doi: 10.1146/annurev-immunol-031210-101400
Crotty, S. (2019). T Follicular Helper Cell Biology: A Decade of Discovery and Diseases. Immunity 50, 1132–1148. doi: 10.1016/j.immuni.2019.04.011
da Silva, E. Z. M., Jamur, M. C., and Oliver, C. (2014). Mast cell function: a new vision of an old cell. J. Histochem. Cytochem. 62, 698–738. doi: 10.1369/0022155414545334
Deng, L., Broaddus, R. R., McCampbell, A., Shipley, G. L., Loose, D. S., Stancel, G. M., et al. (2005). Identification of a Novel Estrogen-Regulated Gene, <em>EIG121</em>, Induced by Hormone Replacement Therapy and Differentially Expressed in Type I and Type II Endometrial Cancer. Clin. Cancer Res. 11:8258. doi: 10.1158/1078-0432.CCR-05-1189
Deng, L., Feng, J., and Broaddus, R. R. (2010). The novel estrogen-induced gene EIG121 regulates autophagy and promotes cell survival under stress. Cell Death Dis. 1, e32–e32. doi: 10.1038/cddis.2010.9
Farhan, M., Aatif, M., Dandawate, P., and Ahmad, A. (2019). “Non-coding RNAs as Mediators of Tamoxifen Resistance in Breast Cancers,” in Breast Cancer Metastasis and Drug Resistance: Challenges and Progress, ed. A. Ahmad (Cham: Springer International Publishing), doi: 10.1007/978-3-030-20301-6_11
Frossi, B., Mion, F., Tripodo, C., Colombo, M. P., and Pucillo, C. E. (2017). Rheostatic Functions of Mast Cells in the Control of Innate and Adaptive Immune Responses. Trends Immunol. 38, 648–656. doi: 10.1016/j.it.2017.04.001
Gao, Y., Chen, S., Vafaei, S., and Zhong, X. (2020). Tumor-Infiltrating Immune Cell Signature Predicts the Prognosis and Chemosensitivity of Patients With Pancreatic Ductal Adenocarcinoma. Front. Oncol. 10:557638–557638. doi: 10.3389/fonc.2020.557638
Guo, Z., Liang, H., Xu, Y., Liu, L., Ren, X., Zhang, S., et al. (2017). The Role of Circulating T Follicular Helper Cells and Regulatory Cells in Non-Small Cell Lung Cancer Patients. Scand. J. Immunol. 86, 107–112. doi: 10.1111/sji.12566
Gu-Trantien, C., Loi, S., Garaud, S., Equeter, C., Libin, M., de Wind, A., et al. (2013). CD4+ follicular helper T cell infiltration predicts breast cancer survival. J. Clin. Investigat. 123, 2873–2892. doi: 10.1172/JCI67428
Hong, W., Gu, Y., Guan, R., Xie, D., Zhou, H., and Yu, M. (2020). Pan-cancer analysis of the CASP gene family in relation to survival, tumor-infiltrating immune cells and therapeutic targets. Genomics 112, 4304–4315. doi: 10.1016/j.ygeno.2020.07.026
Hsieh, J. J., Purdue, M. P., Signoretti, S., Swanton, C., Albiges, L., Schmidinger, M., et al. (2017). Renal cell carcinoma. Nat. Rev. Dis. Prim. 3:17009. doi: 10.1038/nrdp.2017.9
Huang, D., Chen, J., Yang, L., Ouyang, Q., Li, J., Lao, L., et al. (2018). NKILA lncRNA promotes tumor immune evasion by sensitizing T cells to activation-induced cell death. Nat. Immunol. 19, 1112–1125. doi: 10.1038/s41590-018-0207-y
Huang, L., Liang, G., Zhang, Q., and Zhao, W. (2020). The Role of Long Noncoding RNAs in Antiestrogen Resistance in Breast Cancer: An Overview and Update. J. Breast Cancer 23, 129–140. doi: 10.4048/jbc.2020.23.e10
Incorvaia, L., Fanale, D., Badalamenti, G., Brando, C., Bono, M., De Luca, I., et al. (2020). A “Lymphocyte MicroRNA Signature” as Predictive Biomarker of Immunotherapy Response and Plasma PD-1/PD-L1 Expression Levels in Patients with Metastatic Renal Cell Carcinoma: Pointing towards Epigenetic Reprogramming. Cancers 12:3396.
Johnston, R. J., Poholek, A. C., DiToro, D., Yusuf, I., Eto, D., Barnett, B., et al. (2009). Bcl6 and Blimp-1 are reciprocal and antagonistic regulators of T follicular helper cell differentiation. Science 325, 1006–1010. doi: 10.1126/science.1175870
Kang, J. M., Park, S., Kim, S. J., Kim, H., Lee, B., Kim, J., et al. (2015). KIAA1324 Suppresses Gastric Cancer Progression by Inhibiting the Oncoprotein GRP78. Cancer Res. 75:3087. doi: 10.1158/0008-5472.CAN-14-3751
Li, J.-H., Liu, S., Zhou, H., Qu, L.-H., and Yang, J.-H. (2014). starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein–RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 42, D92–D97. doi: 10.1093/nar/gkt1248
Liao, Y., Wang, C., Yang, Z., Liu, W., Yuan, Y., Li, K., et al. (2020). Dysregulated Sp1/miR-130b-3p/HOXA5 axis contributes to tumor angiogenesis and progression of hepatocellular carcinoma. Theranostics 10, 5209–5224. doi: 10.7150/thno.43640
Liu, T., Sui, J., Zhang, Y., Zhang, X. M., Wu, W. J., Yang, S., et al. (2018). Comprehensive analysis of a novel lncRNA profile reveals potential prognostic biomarkers in clear cell renal cell carcinoma. Oncol. Rep. 40, 1503–1514. doi: 10.3892/or.2018.6540
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550. doi: 10.1186/s13059-014-0550-8
Ludwig, N., Nourkami-Tutdibi, N., Backes, C., Lenhof, H.-P., Graf, N., Keller, A., et al. (2015). Circulating serum miRNAs as potential biomarkers for nephroblastoma. Pediatr. Blood Cancer 62, 1360–1367. doi: 10.1002/pbc.25481
Morrissey, J. J., Mellnick, V. M., Luo, J., Siegel, M. J., Figenshau, R. S., Bhayani, S., et al. (2015). Evaluation of Urine Aquaporin-1 and Perilipin-2 Concentrations as Biomarkers to Screen for Renal Cell Carcinoma: A Prospective Cohort Study. JAMA Oncol. 1, 204–212. doi: 10.1001/jamaoncol.2015.0213
Nagano, T., and Fraser, P. (2011). No-Nonsense Functions for Long Noncoding RNAs. Cell 145, 178–181. doi: 10.1016/j.cell.2011.03.014
Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457. doi: 10.1038/nmeth.3337
Oh, D. S., Troester, M. A., Usary, J., Hu, Z., He, X., Fan, C., et al. (2006). Estrogen-Regulated Genes Predict Survival in Hormone Receptor–Positive Breast Cancers. J. Clin. Oncol. 24, 1656–1664. doi: 10.1200/JCO.2005.03.2755
Park, S.-M., Gaur, A. B., Lengyel, E., and Peter, M. E. (2008). The miR-200 family determines the epithelial phenotype of cancer cells by targeting the E-cadherin repressors ZEB1 and ZEB2. Genes Dev. 22, 894–907. doi: 10.1101/gad.1640608
Qu, L., Ding, J., Chen, C., Wu, Z.-J., Liu, B., Gao, Y., et al. (2016). Exosome-Transmitted lncARSR Promotes Sunitinib Resistance in Renal Cancer by Acting as a Competing Endogenous RNA. Cancer Cell 29, 653–668. doi: 10.1016/j.ccell.2016.03.004
Quinn, J. J., and Chang, H. Y. (2016). Unique features of long non-coding RNA biogenesis and function. Nat. Rev. Genet. 17, 47–62. doi: 10.1038/nrg.2015.10
Rohr-Udilova, N., Klinglmüller, F., Schulte-Hermann, R., Stift, J., Herac, M., Salzmann, M., et al. (2018). Deviations of the immune cell landscape between healthy liver and hepatocellular carcinoma. Sci. Rep. 8, 6220–6220. doi: 10.1038/s41598-018-24437-5
Schlumbrecht, M. P., Xie, S.-S., Shipley, G. L., Urbauer, D. L., and Broaddus, R. R. (2011). Molecular clustering based on ERα and EIG121 predicts survival in high-grade serous carcinoma of the ovary/peritoneum. Modern Pathol. 24, 453–462. doi: 10.1038/modpathol.2010.211
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303
Siegel, R. L., Miller, K. D., and Jemal, A. (2018). Cancer statistics, 2018. CA Cancer J. Clin. 68, 7–30. doi: 10.3322/caac.21442
Tay, Y., Rinn, J., and Pandolfi, P. P. (2014). The multilayered complexity of ceRNA crosstalk and competition. Nature 505, 344–352. doi: 10.1038/nature12986
Wang, C., Tang, X., Wang, J., and Xu, Y. (2020). Patterns of immune infiltration in lung adenocarcinoma revealed a prognosis-associated microRNA–mast cells network. Human Cell 33, 205–219. doi: 10.1007/s13577-019-00300-1
Ward, A., Balwierz, A., Zhang, J. D., Küblbeck, M., Pawitan, Y., Hielscher, T., et al. (2013). Re-expression of microRNA-375 reverses both tamoxifen resistance and accompanying EMT-like properties in breast cancer. Oncogene 32, 1173–1182. doi: 10.1038/onc.2012.128
Westin, S. N., Broaddus, R. R., Deng, L., McCampbell, A., Lu, K. H., Lacour, R. A., et al. (2009). Molecular clustering of endometrial carcinoma based on estrogen-induced gene expression. Cancer Biol. Ther. 8, 2126–2135. doi: 10.4161/cbt.8.22.9740
Xiao, W., Wang, X., Wang, T., and Xing, J. (2020). Overexpression of BMP1 reflects poor prognosis in clear cell renal cell carcinoma. Cancer Gene Ther. 27, 330–340. doi: 10.1038/s41417-019-0107-9
Xu, M., Xu, X., Pan, B., Chen, X., Lin, K., Zeng, K., et al. (2019). LncRNA SATB2-AS1 inhibits tumor metastasis and affects the tumor immune cell microenvironment in colorectal cancer by regulating SATB2. Mol. Cancer 18:135. doi: 10.1186/s12943-019-1063-6
Yang, J.-H., Li, J.-H., Shao, P., Zhou, H., Chen, Y.-Q., and Qu, L.-H. (2011). starBase: a database for exploring microRNA–mRNA interaction maps from Argonaute CLIP-Seq and Degradome-Seq data. Nucleic Acids Res. 39(Suppl._1), D202–D209. doi: 10.1093/nar/gkq1056
Zhai, W., Li, S., Zhang, J., Chen, Y., Ma, J., Kong, W., et al. (2018). Sunitinib-suppressed miR-452-5p facilitates renal cancer cell invasion and metastasis through modulating SMAD4/SMAD7 signals. Mol. Cancer 17:157. doi: 10.1186/s12943-018-0906-x
Zhang, X., Wang, M., Sun, H., Zhu, T., and Wang, X. (2018). Downregulation of LINC00894-002 Contributes to Tamoxifen Resistance by Enhancing the TGF-β Signaling Pathway. Biochemistry 83, 603–611. doi: 10.1134/S0006297918050139
Keywords: kidney renal clear cell carcinoma, TCGA, competing endogenous RNA network, immune cell, prognosis
Citation: Deng L, Wang P, Qu Z and Liu N (2021) The Construction and Analysis of ceRNA Network and Immune Infiltration in Kidney Renal Clear Cell Carcinoma. Front. Genet. 12:667610. doi: 10.3389/fgene.2021.667610
Received: 29 May 2021; Accepted: 09 August 2021;
Published: 08 September 2021.
Edited by:
Yadong Zheng, Zhejiang Agriculture and Forestry University, ChinaReviewed by:
Gopal Pandi, Madurai Kamaraj University, IndiaDaniele Fanale, Azienda Ospedaliera Universitaria Policlinico Paolo Giaccone, Italy
Copyright © 2021 Deng, Wang, Qu and Liu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Nan Liu, MTM2ODg4Njk4NzVAMTYzLmNvbQ==; orcid.org/0000-0002-8895-3169