Yu Wang1,2
Ke Hao- 1Laboratory Medicine Center, Allergy Center, Department of Transfusion Medicine, Zhejiang Provincial People’s Hospital (Affiliated People’s Hospital, Hangzhou Medical College), Hangzhou, Zhejiang, China
- 2School of Laboratory Medicine and Life Science, Wenzhou Medical University, Wenzhou, Zhejiang, China
- 3Department of Scientific Research, Zhejiang Provincial People’s Hospital, People’s Hospital of Hangzhou Medical College, Hangzhou, Zhejiang, China
Pancreatic adenocarcinoma (PAAD) is one of the most common malignant tumors with poor prognosis worldwide. Mounting evidence suggests that the expression of lncRNAs and the infiltration of immune cells have prognostic value for patients with PAAD. We used Gene Expression Omnibus (GEO) database and identified six genes (COL1A2, ITGA2, ITGB6, LAMA3, LAMB3, and LAMC2) that could affect the survival rate of pancreatic adenocarcinoma patients. Based on a series of in silico analyses for reverse prediction of target genes associated with the prognosis of PAAD, a ceRNA network of mRNA (COL1A2, ITGA2, LAMA3, LAMB3, and LAMC2)–microRNA (miR-15a-5p)–long non-coding RNA (LINC00511, LINC01578, PVT1, and TNFRSF14-AS1) was constructed. We used the algorithm “CIBERSORT” to assess the proportion of immune cells and found three overall survival (OS)–associated immune cells (monocytes, M1 macrophages, and resting mast cell). Moreover, the OS-associated gene level was significantly positively associated with immune checkpoint expression and biomarkers of immune cells. In summary, our results clarified that ncRNA-mediated upregulation of OS-associated genes and tumor-infiltration immune cells (monocytes, M1 macrophages M1, and resting mast cell resting) correlated with poor prognosis in PAAD.
Introduction
Pancreatic adenocarcinoma (PAAD) is one of the malignant tumors of the digestive system, and it is difficult to diagnose and treat (Siegel et al., 2018). Due to insufficient early symptoms, more than 50% of PAAD patients have metastatic disease at diagnosis (Balsano et al., 2019). Despite the significant improvement in diagnostic imaging and surgical mortality in the past 20 years, the 5-year survival rate is still less than 5% (Philip et al., 2009). The treatment options for pancreatic cancer are limited. About 80% of patients with PAAD lose the chance of surgical resection (Wang et al., 2020). In addition, patients undergoing complete tumor resection usually have a local or distant recurrence within 2 years after the surgery (Hessmann et al., 2020). Other treatment strategies such as radiation, combination chemotherapy, and biological therapy also have limited efficacy (Iacobuzio-Donahue et al., 2009; Pipas et al., 2012). The most fundamental reason is that the pathogenesis and prognostic markers of PAAD are not clear. For better performing early diagnosis, treatment, and prognostic evaluation of PAAD, it is particularly important to find effective biomarkers.
With advances in high-throughput sequencing technology, more biomarkers can be discovered to determine prognosis and improve treatment strategies (Zhang et al., 2018; Ping et al., 2020; Qiu et al., 2020). Non-coding RNAs (ncRNAs) have an important role in the occurrence and development of cancer and are promising biomarkers and therapeutic targets (Cai et al., 2020; Mu et al., 2020; Tseng et al., 2020). lncRNAs are broadly defined as non-coding RNA with a length greater than 200 nucleotides, which previously were thought to have limited protein-coding ability (Mi et al., 2020). However, accumulating evidence demonstrated that specific lncRNAs are related to cancer recurrence, progression, and metastasis (Lin et al., 2020a; Zhou et al., 2020a; Peng et al., 2020; Teng et al., 2020). Salmena et al. proposed the competitive endogenous RNA (ceRNA) hypothesis that lncRNA can affect mRNA through competitively combining different miRNAs (Salmena et al., 2011). miRNAs are a type of small single-stranded ncRNA with regulatory functions, with a length of 21–24 nucleotides (Zhang et al., 2020a; Lin et al., 2020b; Zhou et al., 2020b). One hypothesis of ceRNA considers that miRNAs can regulate multiple target genes, and the same target genes can be regulated by different miRNAs. Also, the ceRNA hypothesis also illustrates that the expression levels of lncRNA and miRNA are negatively correlated and positively correlated with mRNA expression (Wang and Liu, 2017; Li et al., 2020). Previous studies have focused on elucidating the involvement of the ceRNA network in pancreatic adenocarcinoma progression, metastasis, and prognosis, such as LINC00958 and TP73-AS1 (Chen et al., 2019; Miao et al., 2021). In addition, tumor cells and aggressive immune cells participate in the occurrence and development of cancer (Hanahan and Weinberg, 2011). It has become clear that evaluating the type and degree of tumor-infiltrating immune cells is very important for predicting progression and prognosis. However, only a few studies focus on the interaction mechanism between ceRNA networks and immune infiltrating cells. Therefore, a deeper understanding of ceRNA networks and pancreatic adenocarcinoma immune infiltrating cells is needed.
In this study, we established a significant ceRNA network that involved four lncRNAs (LINC00511, LINC01578, PVT1, and TNFRSF14-AS1), one miRNA (hsa-miR-15a-5p), and six mRNAs (COL1A2, ITGA2, ITGB6, LAMA3, LAMB3, and LAMC2). Moreover, we used CIBERSORT to infer the proportion of immune cells in pancreatic adenocarcinoma samples (Chen et al., 2018). In addition, we determined the relationship of OS-associated gene expression with biomarkers of immune cells and immune checkpoints in PAAD.
Materials and Methods
Data Collection and Differential Gene Expression Analysis
The four pancreatic cancer gene chip data used in this study (GSE63111, GSE23397, GSE62165, and GSE43795) were derived from the NCBI-GEO database (https://www.ncbi.nlm.nih.gov/geo). Among them, GSE63111 was included in 35 cases, including in seven tumor samples and 28 non-tumor samples; GSE23397 was included in 21 cases, including in 15 tumor samples and six non-tumor samples; GSE62165 was included in 131 cases, including in 118 tumor samples and 13 non-tumor samples; GSE43795 was included in 31 cases, including in 26 tumor samples and five non-tumor samples. We used the limma package in the R software to screen differentially expressed genes (DEGs) (Ritchie et al., 2015). The screening threshold was set to p < 0.05, and the absolute value of log-fold change | log2FC| ≥ 2. We used the online tool “Venny” to construct the Venn diagram of the DEGs and identified common DEGs (Oliveros, 2007-2015).
Functional Enrichment Analysis, Interaction Network Analysis, and Hub Gene Identification
To further elucidate the potential functional annotation and pathway enrichment related to DEGs, we used the clusterProfiler package to perform the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses (Yu et al., 2012). Significance was defined as p < 0.05. We used an online database Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/) to construct the protein–protein interaction (PPI) network of DEGs, and the confidence score ≥0.9 was set as the threshold. We removed protein nodes that did not interact with other proteins. In addition, we used Cytoscape v.3.8.2 to analyze the PPI network for screening the hub genes. The CytoHubba (version: 0.1) plug-in was used to identify hub genes in the PAAD.
Survival Analysis and Verification
Kaplan–Meier survival analysis curves (p < 0.05) and a univariate Cox proportional hazards model were used to assess the prognostic value of all biomarkers. GEPIA (http://gepia.cancer-pku.cn/) is a web server for gene expression analyses based on TCGA and GTEx databases (Tang et al., 2017). Using this database to perform differential expression analysis between tumor samples and non-tumor samples. Kaplan–Meier Plotter (http://kmplot.com/analysis/) is an online server that can access the effects of genes on survival in multiple cancer types (Nagy et al., 2021). The statistical significance of expression analyses and survival analysis is defined as p < 0.05. Differential expressions of these hub genes at the transcription level were inspected by matching pancreatic tumors and adjacent normal pancreatic tissues from the microarray data we selected in the GEO database. To verify our results, We used the Human Protein Atlas (HPA) online database (www.proteinatlas.org/) to determine the expression of these genes at the translational level (Thul et al., 2017).
Construction of the ceRNA Network
According to the ceRNA hypothesis, we constructed ceRNA networks by backward prediction. First, we used multiple target gene prediction programs to predict the upstream binding miRNAs of OS-associated genes, including miRbase, miRNet, starBase, Tarbase, PITA, and RNA22. Only the predicted miRNAs that appeared in more than two programs previously were regarded as potential miRNAs for subsequent analyses. Then, the mRNA–miRNA network was established by Cytoscape, and the CytoHubba plug-in was used to determine the top potential miRNAs. Next, Starbase (http://starbase.sysu.edu.cn/) was used to perform survival analysis and differential expression on miRNAs. The upstream lncRNAs of miRNAs identified earlier were predicted in the Starbase. Finally, expression analysis and survival analysis were performed on the identified lncRNA. p < 0.05 was considered statistically significant.
Cibesort Estimation and Survival Analysis of Key Members of the Immune Cells
CIBERSORT (http://cibersort.stanford.edu/) was utilized to evaluate the 22 kinds of immune cell types in PAAD. The p-value < 0.05 was chosen as the criterion to identify the immune-infiltrating cells possibly affected by the expression of OS-associated genes. To identify the interrelation between 22 types of immune cells, we performed a correlation-based heat map analysis between every two immune cells. To identify OS-associated immune cells, we used Kaplan–Meier survival analysis to determine the prognostic value.
Immune Infiltrate Analysis and Immune Checkpoint Expression Analysis
We used the GEPIA database to identify the interrelation between OS-associated genes and biomarkers of immune cells. To investigate the association between the expression of OS-associated genes and immune checkpoint expression, we applied TIMER (https://cistrome.shinyapps.io/timer/), which is an online server that can analyze immune infiltrates across multiple cancer types (Li et al., 2017).
Result
Identification of Differential Genes
Through the screening of gene expression microarrays associated with PAAD from the GEO database, four datasets (GSE63111, GSE23397, GSE62165, and GSE43795) were chosen. The DEGs in four datasets are shown in Figures 1A–D. The co-expressed DEGs were confirmed by the Venn Diagram online tool, including 72 upregulated and 58 downregulated genes (Figures 1D–F).
FIGURE 1. Screening of differentially expressed genes. Volcano plot of GSE63111 (A), GSE23397 (B), GSE62165 (C), and GSE43795 (D). (E,F) Upregulated and downregulated DEGs in the GSE63111, GSE23397, GSE62165, and GSE43795. DEGs, differentially expressed genes.
Functional Enrichment Analysis of DEGs, and Protein-Protein Interaction Network
To further investigate the biological functions of DEGs, we conducted GO and KEGG functional enrichment analyses. In the biological process (BP) group, DEGs were mainly enriched in extracellular matrix organization, extracellular structure organization, and external encapsulating structure organization (Figure 2A). In the cellular component (CC) group, DEGs were mainly enriched in the endoplasmic reticulum lumen and collagen-containing extracellular matrix (Figure 2A). In the molecular function (MF) group, DEGs were mainly enriched in endopeptidase activity, extracellular matrix structural constituent, and glycosaminoglycan binding (Figure 2A). Moreover, KEGG pathway analysis results demonstrated that DEGs were significantly enriched in pancreatic secretion and protein digestion and absorption (Figure 2B).
FIGURE 2. GO enrichment analysis, KEGG pathway analysis, and PPI network of DEGs. (A) Bubble chart shows GO enrichment significance items of DEGs. (B) Bubble chart shows enrichment of DEGs in signaling pathways. (C) PPI network of DEGs. PPI, protein–protein interaction. GO, Gene Ontology. KEGG, Kyoto Encyclopedia of Genes and Genomes.
The PPI network of DEGs was analyzed by the STRING online database. A total of 68 nodes and 122 edges were obtained in the PPI network (Figure 2C), and this network was then analyzed using Cytoscape software. Fourteen DEGs were determined as hub genes through the degree algorithm in the CytoHubba plug-in.
Survival Analysis and Verification
Cox regression, Kaplan–Meier, and the log-rank test (p < 0.05) were used to assess the correlation between DEGs and prognosis. Twenty-five genes were found to be significant in the Kaplan–Meier analysis. However, only six promising prognosis-related biomarkers were hub genes of the PPI network (Table 1). We identified these six genes (COL1A2, ITGA2, ITGB6, LAMA3, LAMB3, and LAMC2) that could obviously affect the OS rate of PAAD patients. We analyzed the expression of OS-related genes with tumor stage for PAAD. The results showed that the expression levels of ITGA2, ITGB6, LAMA3, LAMB3, and LAMC2 are significantly correlated with the tumor stage of PAAD patients, while the expression levels of COL1A2 are not correlated with the tumor stage of PAAD patients (Figure 3).
TABLE 1. Cox proportional hazard regression model for OS in patients with PAAD.
FIGURE 3. Correlation between OS-related gene expression and tumor stage in PAAD patients (GEPIA).
The difference between the tumor samples and the non-tumor samples was statistically significant (Figures 4A–F). We used R software to analyze the expression of six hub genes at the transcriptional level. The results showed that the expression of COL1A2, ITGA2, ITGB6, LAMA3, LAMB3, and LAMC2 in PAAD tissue was higher than that in adjacent normal pancreatic tissues (Figures 5A–F). Subsequently, we used the HPA database to explore the expression of the hub gene in protein at the translational level. The results illustrated that the protein level of COL1A2, ITGA2, ITGB6, LAMA3, LAMB3, and LAMC2 was higher in PAAD tissues than in normal pancreatic tissues, matching their mRNA expression levels (Figures 6A–F).
FIGURE 4. Survival analysis results of six DEGs. (A–F) Differential expression and survival analysis of the six selected genes via the GEPIA database. From the abovementioned figures, it was found that these six genes were differentially expressed in the normal group and cancer group and obviously reduce the OS of patients. OS, overall survival.
FIGURE 5. Comparison of the hub gene mRNA levels in paired adjacent normal tissues and PAAD tissues from the GEO database. (A) COL1A2, (B) ITGA2, (C) ITGB6, (D) LAMA3, (E) LAMB3, and (F) LAMC2. PAAD, pancreatic adenocarcinoma.
FIGURE 6. Verification of hub DEG expression in PAAD and normal tissue using the HPA database. (A) COL1A2, (B) ITGA2, (C) ITGB6, (D) LAMA3, (E) LAMB3, and (F) LAMC2.
Construction of the ceRNA Network
To construct a prognostic ceRNA in PAAD, we first predicted upstream miRNAs that could potentially bind to six mRNAs with prognostic values and finally found 160 miRNAs (Figure 7). Cytohubba identified the top fourteen highly connected miRNAs (Figure 8A). The starBase database identified one miRNA, hsa-miR-15a-5p, which was obviously downregulated in PAAD, and its up-regulation was positively correlated with patient prognosis (Figures 8B,C). Next, we use the starBase database to predict the upstream lncRNAs of hsa-miR-15a-5p. 286 possible lncRNAs were obtained. Then, the expression levels of these lncRNAs in PAAD were identified by GEPIA. As shown in Figure 9, LINC00511, LINC01578, PVT1, and TNFRSF14-AS1were significantly upregulated in PAAD compared with normal controls (Figures 9A,C,E,G). Overexpressed LINC00511, LINC01578, PVT1, and TNFRSF14-AS1 indicated poor OS in patients with PAAD (Figures 9B,D,F,H).
FIGURE 7. Result of the use of mRNA to reverse-predict miRNA. Red represents highly expressed mRNA, and blue represents miRNA.
FIGURE 8. Identification of hsa-miR-15a-5p as a potential miRNA in PAAD. Cytoscape software was used to screen the miRNAs that were closely related to hub-mRNAs. Through the differential expression and survival analysis of the starbase database, one miRNA was selected. (A) miRNA–mRNA regulatory network was set up by Cytoscape software. (B) Expression of hsa-miR-15a-5p in PAAD and normal samples was examined by the starBase database. (C) Prognostic value of hsa-miR-15a-5p in PAAD was evaluated by the Kaplan–Meier plotter.
FIGURE 9. Continued.
Composition and Co-Expression Analysis of the Immune Cells in Pancreatic Adenocarcinoma
We use the CIBERSORT algorithm to evaluate the composition of tumor-infiltrating immune cells in PAAD tissue. The histogram (Figure 10A) and the heat map (Figure 10B) illustrate the proportions of the 22 immune cells. The violin plot (Figure 10C) reveals the results of the Wilcoxon rank-sum test. The data showed that the fraction of follicular helper T cells (p < 0.001), monocytes (p = 0.017), M1 macrophages (p = 0.022), resting mast cells (p = 0.001), and activated mast cells (p = 0.008) was obviously different between the pancreatic adenocarcinoma samples and normal pancreatic samples. Then, we used Kaplan–Meier analysis to assess the correlation between the different immune cell subtypes and prognosis. The fraction of monocytes (p = 0.00921), M1 macrophages (p = 0.0315), and resting mast cell (p = 0.0499) was found to be obviously associated with OS (Figures 11A–C). We perform co-expression analysis on prognostic tumor-infiltrating immune cells. (Figure 11D). Finally, Figure 11E showed the relationship of co-expression between hub genes and immune cells. The fraction of M1 macrophages was positively correlated with COL1A2 expression (R = 0.41, p < 0.001) (Supplementary Figure S1A). The M1 macrophages were positively correlated with LAMA3 expression (R = 0.28, p < 0.001) (Supplementary Figure S1B). The monocytes were negatively correlated with LAMA3 expression (R = −0.38) (Supplementary Figure S1C). The monocytes were negatively correlated with ITGA2 expression (R = −0.3, p < 0.001) (Supplementary Figure S1D).
FIGURE 10. Proportion of the 22 immune cells detected by the CIBERSORT algorithm (A,B). Results of the Wilcoxon rank-sum test (C). Proportions of follicular helper T cells (p < 0.001), monocytes (p = 0.017), M1 macrophages (p = 0.022), resting mast cells (p = 0.001), and activated mast cells (p = 0.008) were significantly different between pancreatic cancer samples and normal pancreatic samples.
FIGURE 11. Correlation between immune cell fraction and OS (A–C). Significant co-expression patterns between proportions of immune cells (D). Significant co-expression patterns between hub mRNAs and immune cells (E).
Association of Biomarkers of Immune Cells and Immune Checkpoints With OS-Associated Gene Expression in PAAD
We researched the relationship between the expression of OS-associated genes and biomarkers of immune cells. Monocytes, M1 macrophages, and resting mast cells are the significant tumor-infiltrating immune cell. The results were strongly correlated with OS-associated immune cells (Table 2). PD1/PD-L1 and CTLA-4 are significant immune checkpoints for immune escape. We evaluated the relationship of OS-associated genes with PD1, PD-L1, and CTLA-4 by the TIMER platform. As suggested in Figure 12, OS-associated gene expression was clearly positively associated with PD1, PD-L1, and CTLA-4 in PAAD, which was adjusted by purity. These data indicate that immune escape might be concerned with the key gene–mediated poor prognosis of PAAD.
TABLE 2. Relationship between hub genes and biomarkers of immune cells in PAAD detected by GEPIA.
FIGURE 12. Correlation of hub mRNA expression with PD-1, PD-L1, and CTLA-4 expression in PAAD. Spearman relationship of COL1A2 with an expression of PD-1 (A), PD-L1 (B), and CTLA-4 (C) in PAAD adjusted by purity using TIMER. Spearman relationship of ITGA2 with an expression of PD-1 (D), PD-L1 (E), and CTLA-4 (F) in PAAD adjusted by purity using TIMER. Spearman relationship of LAMA3 with an expression of PD-1 (G), PD-L1 (H), and CTLA-4 (I) in PAAD adjusted by purity using TIMER. Spearman relationship of LAMB3 with an expression of PD-1 (J), PD-L1 (K), and CTLA-4 (L) in PAAD adjusted by purity using TIMER. Spearman relationship of LAMC2 with an expression of PD-1 (M), PD-L1 (N), and CTLA-4 (O) in PAAD adjusted by purity using TIMER.
Discussion
Pancreatic cancer is the deadliest form of cancer. Despite different types of treatment options that have been adopted for pancreatic cancer patients, the prognosis is still poor (Siegel et al., 2020). Molecularly targeted therapies targeting oncogenic genes associated with pancreatic cancer genesis and development are increasingly considered a promising way to cure pancreatic cancer (Paulson et al., 2013). Therefore, illustrating the molecular mechanisms of the pathogenesis of pancreatic cancer and determining potential biomarkers are crucial to improve patient outcomes.
Our study first performed bioinformatics analysis based on GEO datasets, and 130 genes of differentially expressed genes were selected by analyzing the microarray (GSE63111, GSE23397, GSE62165, and GSE43795) related to pancreatic cancer, of which 72 were upregulated and 58 were downregulated. Enrichment analysis showed that DEGs were mainly concentrated in external encapsulating structure organization, extracellular structure organization, and extracellular matrix organization. Next, six key genes in the PPI network were screened by Cox regression and Kaplan–Meier analysis. According to the ceRNA hypothesis, we gradually determined upstream miRNAs and lncRNAs from mRNAs and finally successfully established the ceRNA network. All genes of the ceRNA network are associated with poor prognosis of PAAD. In addition, we found that hsa-miR-15a-5p and hub genes in ceRNA and immune cell subtypes (monocytes, M1 macrophages, and resting mast cell) could predict prognosis effectively. Finally, our study found that high expression of OS-associated genes was closely related to PD1, PD-L1, or CTLA-4 in PAAD, suggesting that targeting mRNAs might improve the efficacy of immunotherapy in PAAD.
COL1A2 belongs to the family of ECM proteins. The ECM provides structural support, mediates cell-matrix adhesion, and integrates complex, multivalent signals for cells (Egusa et al., 2007; Hynes, 2009). Wu et al. Wu et al. (2019a) reduced the expression of COL1A2 to reduce the invasion and migration ability of pancreatic cancer cells and found that COL1A2 was highly expressed in pancreatic cancer patients (p < 0.05), which was related to poor prognosis in these patients. ITGA2 belongs to the integrin alpha chain family. Integrin is a membrane protein with cell adhesion and signal transduction functions (Desgrosellier and Cheresh, 2010). Studies have found that the expression of ITGA2 is an independent prognostic marker of PAAD and is related to the grading and aggressiveness of PAAD (Shang et al., 2019; Yan et al., 2020). They also found that the expression of ITGA2 could significantly reduce the survival rate of PAAD patients, which is in accordance with our research. ITGB6 belongs to a family of β2 integrins. Some studies have illustrated that the expression of ITGB6 is negatively related to prognosis in PAAD patients (Wu et al., 2019b), which is consistent with our survival analysis. LAMA3, LAMB3, and LAMC2 belong to the laminin gene family. Laminin is known to stimulate cell migration in cancers (Tripathi et al., 2008). Yang et al. Yang et al. (2019) found that four subunits of the laminin gene family (LAMA3, LAMA4, LAMB3, and LAMC2) were connected with the outcomes of PAAD patients. The findings in our study are consistent with previous reports as we found that the LAMA3, LAMB3, and LAMC2 genes were upregulated in patients with PAAD, with hazard ratios of 3.86, 2.18, and 3.06, respectively. Combined with patient survival curves, these genes were found to be poor prognostic markers in patients with PAAD.
Hsa-miR-15a-5p is a microRNA whose family has already become promising cancer-specific biomarkers (Shao et al., 2020). Compared with PAAD patients with high miR-15a-5p expression, those patients with low miR-15a-5p levels have a shorter overall survival (Guo et al., 2020), which is consistent with the results of our Figure 8C. Functionally, exogenous miR-15a-5p expression decreases cell growth, epithelial-to-mesenchymal transition, and metastasis and increases cell cycle arrest in PAAD (Guo et al., 2014; Feng et al., 2019). Mechanistically, It has been confirmed that miR-15a-5p exerts its anti-pancreatic cancer activity by directly targeting FGFR1, which is highly expressed and executes cancer-promoting actions in PAAD (Turner et al., 2008; Lehnen et al., 2013; Coleman et al., 2014). A study focusing on the mechanism of hsa-miR-15a-5p in pancreatic cancer independent from the CERS6-AS1 (Shen et al., 2021). In addition, research shows a relationship between patient prognosis and tumor-associated macrophage infiltration in pancreatic cancer (Di Caro et al., 2016). In our present study, M1 macrophages which were chosen from CIEBERSORT and the co-expression correlation of COL1A2, LAMA3, and M1 macrophages. The correlation analysis results were consistent with the previous research that indicated that there is a regulatory connection between M1 macrophages and COL1A2 (Zhang et al., 2020b). Cells of the monocyte lineage contribute to tumor angiogenesis (Dalton et al., 2014). One study suggested monocytes are recruited into Pancreatic Ductal Adenocarcinoma tumors via IL35 and increased expression of genes that promote angiogenesis in IL35 products (Huang et al., 2018). Nevertheless, there are few research studies before about the connections of monocytes and Hsa-miR-15a-5p. Thus, we deduced that Hsa-miR-15a-5p, as a miRNA, might connect with monocytes and took active part in the progression of PAAD.
Mishra et al.s’ study represents the first TCGA-based PDAC methylome data analysis and first reported that several genes (B3GNT3, DMBT1, and DEPDC1B) and lncRNAs (PVT1, and GATA6-AS) are strongly correlated with pancreatic ductal adenocarcinoma survival (Kumar Mishra et al., 2019). The conclusions regarding PVT1 are consistent with our results in Figure 9E, F. Raman et al. identified a 5-gene panel (ADM, ASPM, DCBLD2, E2F7, and KRT6A) that performed well against previous signatures across multiple datasets and captures subtype-specific differences in patient prognosis (Raman et al., 2018). This study provides a framework for similar assessments in other cancers.
It should be acknowledged that our study has some unavoidable limitations. First, the quantity of samples we gather from GEO databases is indeed limited, which suggests that clinico pathological information was not comprehensive and might bring about some potential errors or biases. In addition, only a limited number of links were predicted. The ceRNA network is highly complicated and diverse, and its complex regulatory network needs to be further elucidated. Finally, we did not consider the heterogeneity of the immune microenvironment which is allied to the immune infiltration location. The heterogeneity and complexity of the immune microenvironment may affect the generalization of results. However, to reduce the bias, we used numerous databases to reveal the expression of crucial biomarkers at the gene level, protein level, and tissue levels, and it showed that crucial biomarkers were upregulated in PAAD. This study provides valuable material for further research, in which we would further explore the expression levels of key members of the ceRNA network in clinical samples and the interaction mechanism between ceRNA and cell communication through experiments.
In summary, through the integration of several bioinformatics analyses, we constructed a reverse mRNA model based on the LINC00511, LINC01578, PVT1, TNFRSF14-AS/hsa-miR-15a-5p/COL1A2, ITGA2, ITGB6, LAMA3, LAMB3, and LAMC2 ceRNA network, which can enhance our understanding of the development of PAAD. HR values of all genes in ceRNA are greater than 1, indicating that they are risk factors for PAAD, so they might be utilized as potential prognostic biomarkers and therapeutic targets for PAAD. Moreover, our study reveals the potential mechanism that hsa-miR-15a-5p regulates COL1A2, ITGA2, and LAMA3 so that monocytes and M1 macrophages are activated, which may have great effects on pancreatic cancer prognosis.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author Contributions
KH and LJ were responsible for study design; YW and ZW were involved in data collection; KH and LJ analyzed the data; YW wrote the manuscript. KH, LJ, ZW, KL, WX, and BC revised the manuscript.
Funding
This work was supported by grants from the National Natural Science Foundation of China, NO. 81600595, and NO. 82072366; the Medicine and Health Research Foundation of Zhejiang Province, Nos. 2019RC113, 2019KY017, 2019RC124, and 2020KY012; Key projects jointly constructed by the Ministry and the province of Zhejiang Medical and Health Science and Technology Project, Nos. WKJ-ZJ-2019 and WKJ-ZJ-2101.
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.2022.874667/full#supplementary-material
Supplementary Figure S1 | Fraction of M1 macrophage cells was positively associated with COL1A2 expression (A) and LAMA3 expression (B). Monocyte cells were positively associated with LAMA3 (C) and ITGA2 (D).
References
Balsano, R., Tommasi, C., and Garajova, I. (2019). State of the Art for Metastatic Pancreatic Cancer Treatment: Where Are We Now? Anticancer Res. 39 (7), 3405–3412. doi:10.21873/anticanres.13484
Cai, L., Ye, Y., Jiang, Q., Chen, Y., Lyu, X., Li, J., et al. (2020). Author Correction: Epstein-Barr Virus-Encoded microRNA BART1 Induces Tumour Metastasis by Regulating PTEN-dependent Pathways in Nasopharyngeal Carcinoma. Nat. Commun. 11, 3437. doi:10.1038/s41467-020-17048-0
Chen, B., Khodadoust, M. S., Liu, C. L., Newman, A. M., and Alizadeh, A. A. (2018). Profiling Tumor Infiltrating Immune Cells with Cibersort. Methods Mol. Biol. 1711, 243–259. doi:10.1007/978-1-4939-7493-1_12
Chen, S., Chen, J. Z., Zhang, J. Q., Chen, H. X., Qiu, F. N., Yan, M. L., et al. (2019). Silencing of Long Noncoding RNA LINC00958 Prevents Tumor Initiation of Pancreatic Cancer by Acting as a Sponge of microRNA-330-5p to Down-Regulate PAX8. Cancer Lett. 446, 49–61. doi:10.1016/j.canlet.2018.12.017
Coleman, S. J., Chioni, A. M., Ghallab, M., Anderson, R. K., Lemoine, N. R., Kocher, H. M., et al. (2014). Nuclear Translocation ofFGFR1 andFGF2 in Pancreatic Stellate Cells Facilitates Pancreatic Cancer Cell Invasion. EMBO Mol. Med. 6, 467–481. doi:10.1002/emmm.201302698
Dalton, H. J., Armaiz-Pena, G. N., Gonzalez-Villasana, V., Lopez-Berestein, G., Bar-Eli, M., and Sood, A. K. (2014). Monocyte Subpopulations in Angiogenesis. Cancer Res. 74, 1287–1293. doi:10.1158/0008-5472.CAN-13-2825
Desgrosellier, J. S., and Cheresh, D. A. (2010). Integrins in Cancer: Biological Implications and Therapeutic Opportunities. Nat. Rev. Cancer 10, 9–22. doi:10.1038/nrc2748
Di Caro, G., Cortese, N., Castino, G. F., Grizzi, F., Gavazzi, F., Ridolfi, C., et al. (2016). Dual Prognostic Significance of Tumour-Associated Macrophages in Human Pancreatic Adenocarcinoma Treated or Untreated with Chemotherapy. Gut 65, 1710–1720. doi:10.1136/gutjnl-2015-309193
Egusa, H., Iida, K., Kobayashi, M., Lin, T. Y., Zhu, M., Zuk, P. A., et al. (2007). Downregulation of Extracellular Matrix-Related Gene Clusters during Osteogenic Differentiation of Human Bone Marrow- and Adipose Tissue-Derived Stromal Cells. Tissue Eng. 13, 2589–2600. doi:10.1089/ten.2007.0080
Feng, H., Wei, B., and Zhang, Y. (2019). RETRACTED: Long Non-coding RNA HULC Promotes Proliferation, Migration and Invasion of Pancreatic Cancer Cells by Down-Regulating microRNA-15a. Int. J. Biol. Macromol. 126, 891–898. doi:10.1016/j.ijbiomac.2018.12.238
Guo, S., Fesler, A., Huang, W., Wang, Y., Yang, J., Wang, X., et al. (2020). Functional Significance and Therapeutic Potential of miR-15a Mimic in Pancreatic Ductal Adenocarcinoma. Mol. Ther. - Nucleic Acids 19, 228–239. doi:10.1016/j.omtn.2019.11.010
Guo, S., Xu, X., Tang, Y., Zhang, C., Li, J., Ouyang, Y., et al. (2014). miR-15a Inhibits Cell Proliferation and Epithelial to Mesenchymal Transition in Pancreatic Ductal Adenocarcinoma by Down-Regulating Bmi-1 Expression. Cancer Lett. 344, 40–46. doi:10.1016/j.canlet.2013.10.009
Hanahan, D., and Weinberg, R. A. (2011). Hallmarks of Cancer: The Next Generation. Cell. 144, 646–674. doi:10.1016/j.cell.2011.02.013
Hessmann, E., Buchholz, S. M., Demir, I. E., Singh, S. K., Gress, T. M., Ellenrieder, V., et al. (2020). Microenvironmental Determinants of Pancreatic Cancer. Physiol. Rev. Rev. 100, 1707–1751. doi:10.1152/physrev.00042.2019
Huang, C., Li, Z., Li, N., Li, Y., Chang, A., Zhao, T., et al. (2018). Interleukin 35 Expression Correlates with Microvessel Density in Pancreatic Ductal Adenocarcinoma, Recruits Monocytes, and Promotes Growth and Angiogenesis of Xenograft Tumors in Mice. Gastroenterology 154, 675–688. doi:10.1158/0008-5472.CAN-17-161210.1053/j.gastro.2017.09.039
Hynes, R. O. (2009). The Extracellular Matrix: Not Just Pretty Fibrils. Science 326, 1216–1219. doi:10.1126/science.1176009
Iacobuzio-Donahue, C. A., Fu, B., Yachida, S., Luo, M., Abe, H., Henderson, C. M., et al. (2009). DPC4 Gene Status of the Primary Carcinoma Correlates with Patterns of Failure in Patients with Pancreatic Cancer. Jco 27, 1806–1813. doi:10.1200/jco.2008.17.7188
Kumar Mishra, N., Southekal, S., and Guda, C. (2019). Survival Analysis of Multi-Omics Data Identifies Potential Prognostic Markers of Pancreatic Ductal Adenocarcinoma. Front. Genet. 10, 624. doi:10.3389/fgene.2019.00624
Lehnen, N. C., von Mässenhausen, A., Kalthoff, H., Zhou, H., Glowka, T., Schütte, U., et al. (2013). Fibroblast Growth Factor Receptor 1 Gene Amplification in Pancreatic Ductal Adenocarcinoma. Histopathology 63, 157–166. doi:10.1111/his.12115
Li, J., Wang, X., Wang, Y., and Yang, Q. (2020). H19 Promotes the Gastric Carcinogenesis by Sponging miR-29a-3p: Evidence from lncRNA-miRNA-mRNA Network Analysis. Epigenomics 12, 989–1002. doi:10.2217/epi-2020-0114
Li, T., Fan, J., Wang, B., Traugh, N., Chen, Q., Liu, J. S., et al. (2017). TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res. 77, e108–e110. doi:10.1158/0008-5472.can-17-0307
Lin, L., Que, Y., Lu, P., Li, H., Xiao, M., Zhu, X., et al. (2020). Chidamide Inhibits Acute Myeloid Leukemia Cell Proliferation by lncRNA VPS9D1-AS1 Downregulation via MEK/ERK Signaling Pathway. Front. Pharmacol. 11, 569651. doi:10.3389/fphar.2020.569651
Lin, L., Que, Y., Lu, P., Li, H., Xiao, M., Zhu, X., et al. (2020). Chidamide Inhibits Acute Myeloid Leukemia Cell Proliferation by lncRNA VPS9D1-AS1 Downregulation via MEK/ERK Signaling Pathway. Front. Pharmacol. 11, 569651. doi:10.3389/fphar.2020.569651
Mi, X., Xu, R., Hong, S., Xu, T., Zhang, W., and Liu, M. (2020). M2 Macrophage-Derived Exosomal lncRNA AFAP1-AS1 and MicroRNA-26a Affect Cell Migration and Metastasis in Esophageal Cancer. Mol. Ther. - Nucleic Acids 22, 779–790. doi:10.1016/j.omtn.2020.09.035
Miao, H., Lu, J., Guo, Y., Qiu, H., Zhang, Y., Yao, X., et al. (2021). LncRNA TP73-AS1 Enhances the Malignant Properties of Pancreatic Ductal Adenocarcinoma by Increasing MMP14 Expression through miRNA -200a Sponging. J. Cell. Mol. Med. 25 (7), 3654–3664. doi:10.1111/jcmm.16425
Mu, M., Niu, W., Zhang, X., Hu, S., and Niu, C. (2020). LncRNA BCYRN1 Inhibits Glioma Tumorigenesis by Competitively Binding with miR-619-5p to Regulate CUEDC2 Expression and the PTEN/AKT/p21 Pathway. Oncogene 39, 6879–6892. doi:10.1038/s41388-020-01466-x
Nagy, Á., Munkácsy, G., and Győrffy, B. (2021). Pancancer Survival Analysis of Cancer Hallmark Genes. Sci. Rep. 11 (1), 6047. doi:10.1038/s41598-021-84787-5
Oliveros, J. C. (2007-2015). Venny. An Interactive Tool for Comparing Lists with Venn's Diagrams. https://bioinfogp.cnb.csic.es/tools/venny/index.html.
Paulson, A. S., Tran Cao, H. S., Tempero, M. A., and Lowy, A. M. (2013). Therapeutic Advances in Pancreatic Cancer. Gastroenterology 144, 1316–1326. doi:10.1053/j.gastro.2013.01.078
Peng, W., Zhang, C., Peng, J., Huang, Y., Peng, C., Tan, Y., et al. (2020). lnc-FAM84B-4 Acts as an Oncogenic lncRNA by Interacting with Protein hnRNPK to Restrain MAPK Phosphatases-DUSP1 Expression. Cancer Lett. 494, 94–106. doi:10.1016/j.canlet.2020.08.036
Philip, P. A., Mooney, M., Jaffe, D., Eckhardt, G., Moore, M., Meropol, N., et al. (2009). Consensus Report of the National Cancer Institute Clinical Trials Planning Meeting on Pancreas Cancer Treatment. Jco 27 (33), 5660–5669. doi:10.1200/JCO.2009.21.9022
Ping, Y., Zhou, Y., Hu, J., Pang, L., Xu, C., and Xiao, Y. (2020). Dissecting the Functional Mechanisms of Somatic Copy-Number Alterations Based on Dysregulated ceRNA Networks across Cancers. Mol. Ther. - Nucleic Acids 21, 464–479. doi:10.1016/j.omtn.2020.06.012
Pipas, J. M., Zaki, B. I., McGowan, M. M., Tsapakos, M. J., Ripple, G. H., Suriawinata, A. A., et al. (2012). Neoadjuvant Cetuximab, Twice-Weekly Gemcitabine, and Intensity-Modulated Radiotherapy (IMRT) in Patients with Pancreatic Adenocarcinoma. Ann. Oncol. 23, 2820–2827. doi:10.1093/annonc/mds109
Qiu, Y. L., Zheng, H., Devos, A., Selby, H., and Gevaert, O. (2020). A Meta-Learning Approach for Genomic Survival Analysis. Nat. Commun. 11, 6350. doi:10.1038/s41467-020-20167-3
Raman, P., Maddipati, R., Lim, K. H., and Tozeren, A. (2018). Pancreatic Cancer Survival Analysis Defines a Signature that Predicts Outcome. PLoS One 13 (8), e0201751. doi:10.1371/journal.pone.0201751
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007
Salmena, L., Poliseno, L., Tay, Y., Kats, L., and Pandolfi, P. P. (2011). A ceRNA Hypothesis: the Rosetta Stone of a Hidden RNA Language? Cell. 146, 353–358. doi:10.1016/j.cell.2011.07.014
Shang, M., Zhang, L., Chen, X., and Zheng, S. (2019). Identification of Hub Genes and Regulators Associated with Pancreatic Ductal Adenocarcinoma Based on Integrated Gene Expression Profile Analysis. Discov. Med. 28 (153), 159–172.
Shao, Y., Guo, X., Zhao, L., Shen, Y., Niu, C., Wei, W., et al. (2020). A Functional Variant of the miR-15 Family Is Associated with a Decreased Risk of Esophageal Squamous Cell Carcinoma. DNA Cell. Biol. 39 (9), 1583–1594. doi:10.1089/dna.2020.5606
Shen, R., Wang, X., Wang, S., Zhu, D., and Li, M. (2021). Long Noncoding RNA CERS6-AS1 Accelerates the Proliferation and Migration of Pancreatic Cancer Cells by Sequestering MicroRNA-15a-5p and MicroRNA-6838-5p and Modulating HMGA1. Pancreas 50 (4), 617–624. doi:10.1097/MPA.0000000000001806
Siegel, R. L., Miller, K. D., and Jemal, A. (20182018)., 68. CA, 7–30. doi:10.3322/caac.21442Cancer Statistics, 2018CA A Cancer J. Clin.1
Siegel, R. L., Miller, K. D., and Jemal, A. (20202020). Cancer Statistics, 2020. CA A Cancer J. Clin. 70, 7–30. doi:10.3322/caac.21590
Tang, Z., Li, C., Kang, B., Gao, G., Li, C., and Zhang, Z. (2017). GEPIA: a Web Server for Cancer and Normal Gene Expression Profiling and Interactive Analyses. Nucleic Acids Res. 45 (W1), W98–W102. doi:10.1093/nar/gkx247
Teng, F., Zhang, J.-X., Chang, Q.-M., Wu, X.-B., Tang, W.-G., Wang, J.-F., et al. (2020). lncRNA MYLK-AS1 Facilitates Tumor Progression and Angiogenesis by Targeting miR-424-5p/E2F7 axis and Activating VEGFR-2 Signaling Pathway in Hepatocellular Carcinoma. J. Exp. Clin. Cancer Res. 39, 235. doi:10.1186/s13046-020-01739-z
Thul, P. J., Åkesson, L., Wiking, M., Mahdessian, D., Geladaki, A., Ait Blal, H., et al. (2017). A Subcellular Map of the Human Proteome. Science 356 (6340). doi:10.1126/science.aal3321
Tripathi, M., Nandana, S., Yamashita, H., Ganesan, R., Kirchhofer, D., and Quaranta, V. (2008). Laminin-332 Is a Substrate for Hepsin, a Protease Associated with Prostate Cancer Progression. J. Biol. Chem. 283, 30576–30584. doi:10.1074/jbc.M802312200
Tseng, L.-L., Cheng, H.-H., Yeh, T.-S., Huang, S.-C., Syu, Y.-Y., Chuu, C.-P., et al. (2020). Targeting the Histone Demethylase PHF8-Mediated PKCα-Src-PTEN axis in HER2-Negative Gastric Cancer. Proc. Natl. Acad. Sci. U.S.A. 117, 24859–24866. doi:10.1073/pnas.1919766117
Turner, C. A., Calvo, N., Frost, D. O., Akil, H., and Watson, S. J. (2008). The Fibroblast Growth Factor System Is Downregulated Following Social Defeat. Neurosci. Lett. 430, 147–150. doi:10.1016/j.neulet.2007.10.041
Wang, L., and Liu, J. (2017). Research Progress of Competing Endogenous RNA. Sheng Wu Yi Xue Gong Cheng Xue Za Zhi 34, 967–971. doi:10.7507/1001-5515.201606080
Wang, X., Li, X., Wei, X., Jiang, H., Lan, C., Yang, S., et al. (2020). PD-L1 Is a Direct Target of Cancer-FOXP3 in Pancreatic Ductal Adenocarcinoma (PDAC), and Combined Immunotherapy with Antibodies against PD-L1 and CCL5 Is Effective in the Treatment of PDAC. Sig Transduct. Target Ther. 5, 38. doi:10.1038/s41392-020-0144-8
Wu, J., Liu, J., Wei, X., Yu, Q., Niu, X., Tang, S., et al. (2019). A Feature-Based Analysis Identifies COL1A2 as a Regulator in Pancreatic Cancer. J. Enzyme Inhibition Med. Chem. 34 (1), 420–428. doi:10.1080/14756366.2018.1484734
Wu, M., Li, X., Zhang, T., Liu, Z., and Zhao, Y. (2019). Identification of a Nine-Gene Signature and Establishment of a Prognostic Nomogram Predicting Overall Survival of Pancreatic Cancer. Front. Oncol. 9, 996. doi:10.3389/fonc.2019.00996
Yan, W., Liu, X., Wang, Y., Han, S., Wang, F., Liu, X., et al. (2020). Identifying Drug Targets in Pancreatic Ductal Adenocarcinoma through Machine Learning, Analyzing Biomolecular Networks, and Structural Modeling. Front. Pharmacol. 11, 534. doi:10.3389/fphar.2020.00534
Yang, C., Liu, Z., Zeng, X., Wu, Q., Liao, X., Wang, X., et al. (2019). Evaluation of the Diagnostic Ability of Laminin Gene Family for Pancreatic Ductal Adenocarcinoma. Aging 11 (11), 3679–3703. doi:10.18632/aging.102007
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS A J. Integr. Biol. 16, 284–287. doi:10.1089/omi.2011.0118
Zhang, J., Sun, H., Liu, S., Huang, W., Gu, J., Zhao, Z., et al. (2020). Alteration of Tumor-Associated Macrophage Subtypes Mediated by KRT6A in Pancreatic Ductal Adenocarcinoma. aging 12 (22), 23217–23232. doi:10.18632/aging.104091
Zhang, Q., Song, G., Yao, L., Liu, Y., Liu, M., Li, S., et al. (2018). miR-3928v Is Induced by HBx via NF-Κb/egr1 and Contributes to Hepatocellular Carcinoma Malignancy by Down-Regulating VDAC3. J. Exp. Clin. Cancer Res. 37, 14. doi:10.1186/s13046-018-0681-y
Zhang, Y., Huang, Y.-X., Wang, D.-L., Yang, B., Yan, H.-Y., Lin, L.-H., et al. (2020). LncRNA DSCAM-AS1 Interacts with YBX1 to Promote Cancer Progression by Forming a Positive Feedback Loop that Activates FOXA1 Transcription Network. Theranostics 10, 10823–10837. doi:10.7150/thno.47830
Zhou, R., Sun, H., Zheng, S., Zhang, J., Zeng, D., Wu, J., et al. (2020). A Stroma‐related lncRNA Panel for Predicting Recurrence and Adjuvant Chemotherapy Benefit in Patients with Early‐stage Colon Cancer. J. Cell. Mol. Med. 24, 3229–3241. doi:10.1111/jcmm.14999
Keywords: pancreatic adenocarcinoma, ceRNA network, immune infiltration, therapeutic targets, prognosis
Citation: Wang Y, Wang Z, Li K, Xiang W, Chen B, Jin L and Hao K (2022) lncRNAs Functioned as ceRNA to Sponge miR-15a-5p Affects the Prognosis of Pancreatic Adenocarcinoma and Correlates With Tumor Immune Infiltration. Front. Genet. 13:874667. doi: 10.3389/fgene.2022.874667
Received: 17 February 2022; Accepted: 26 May 2022;
Published: 11 July 2022.
Edited by:
Manal S. Fawzy, Suez Canal University, EgyptReviewed by:
Rongying Ou, First Affiliated Hospital of Wenzhou Medical University, ChinaLi Xueling, Inner Mongolia University, China
Copyright © 2022 Wang, Wang, Li, Xiang, Chen, Jin and Hao. 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: LiQin Jin, bGlxaW5qaW5AMTI2LmNvbQ==; Ke Hao, aGFva2VAaG1jLmVkdS5jbg==