Next Article in Journal
TERT Expression Induces Resistance to BRAF and MEK Inhibitors in BRAF-Mutated Melanoma In Vitro
Next Article in Special Issue
Phase II Trial Evaluating Olaparib Maintenance in Patients with Metastatic Castration-Resistant Prostate Cancer Responsive or Stabilized on Docetaxel Treatment: SOGUG-IMANOL Study
Previous Article in Journal
Unraveling the Role of Epithelial–Mesenchymal Transition in Adenoid Cystic Carcinoma of the Salivary Glands: A Comprehensive Review
Previous Article in Special Issue
c-MYC-Induced AP4 Attenuates DREAM-Mediated Repression by p53
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Cancers
Volume 15
Issue 11
10.3390/cancers15112887
cancers-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Systematic Review

A Comprehensive Transcriptional Signature in Pancreatic Ductal Adenocarcinoma Reveals New Insights into the Immune and Desmoplastic Microenvironments

by
Irene Pérez-Díez
1,2,3,†,
Zoraida Andreu
3,†,
Marta R. Hidalgo
1,3,
Carla Perpiñá-Clérigues
1,3,4,
Lucía Fantín
1,
Antonio Fernandez-Serra
3,5,
María de la Iglesia-Vaya
2,3,
José A. Lopez-Guerrero
3,5,6,* and
Francisco García-García
1,3,*
1
Bioinformatics and Biostatistics Unit, Principe Felipe Research Center (CIPF), 46012 Valencia, Spain
2
Biomedical Imaging Unit FISABIO-CIPF, Fundación para el Fomento de la Investigación Sanitaria y Biomédica de la Comunidad Valenciana, 46012 Valencia, Spain
3
IVO-CIPF Joint Research Unit of Cancer, Príncipe Felipe Research Center (CIPF), 46012 Valencia, Spain
4
Department of Physiology, School of Medicine and Dentistry, University of Valencia, 46010 Valencia, Spain
5
Laboratory of Molecular Biology, Fundación Instituto Valenciano de Oncología, 46009 Valencia, Spain
6
Department of Pathology, Medical School, Catholic University of Valencia, 46001 Valencia, Spain
*
Authors to whom correspondence should be addressed.
†
These authors contributed equally to this work.
Cancers 2023, 15(11), 2887; https://doi.org/10.3390/cancers15112887
Submission received: 7 April 2023 / Revised: 11 May 2023 / Accepted: 23 May 2023 / Published: 24 May 2023
(This article belongs to the Special Issue Oncogenes in Cancer 2.0)
Browse Figures
Review Reports Versions Notes

Abstract

:

Simple Summary

Pancreatic ductal adenocarcinoma (PDAC) is a highly lethal disease with a few curative options. Desmoplastic stroma and immune system evasion in PDAC represent challenges to the success of therapeutic strategies are used to suitably treat other tumor types. Characterizing the PDAC microenvironment (including the immune environment) remains critical to developing safe and efficient therapies. Here, we present a comprehensive meta-analysis identifying 1153 significantly dysregulated genes, which mainly impact extracellular matrix remodeling and the immune system. We identify two signatures of twenty-eight immune-related genes and eleven stroma-related genes influencing PDAC patients’ survival. Additionally, five immune genes are associated with PDAC prognosis for the first time.

Abstract

Pancreatic ductal adenocarcinoma (PDAC) prognoses and treatment responses remain devastatingly poor due partly to the highly heterogeneous, aggressive, and immunosuppressive nature of this tumor type. The intricate relationship between the stroma, inflammation, and immunity remains vaguely understood in the PDAC microenvironment. Here, we performed a meta-analysis of stroma-, and immune-related gene expression in the PDAC microenvironment to improve disease prognosis and therapeutic development. We selected 21 PDAC studies from the Gene Expression Omnibus and ArrayExpress databases, including 922 samples (320 controls and 602 cases). Differential gene enrichment analysis identified 1153 significant dysregulated genes in PDAC patients that contribute to a desmoplastic stroma and an immunosuppressive environment (the hallmarks of PDAC tumors). The results highlighted two gene signatures related to the immune and stromal environments that cluster PDAC patients into high- and low-risk groups, impacting patients’ stratification and therapeutic decision making. Moreover, HCP5, SLFN13, IRF9, IFIT2, and IFI35 immune genes are related to the prognosis of PDAC patients for the first time.

1. Introduction

Pancreatic ductal adenocarcinoma (PDAC) is the most common type of pancreatic cancer, representing over 80% of all diagnosed pancreatic neoplasms. This highly lethal cancer has a poor prognosis, with a median survival rate of fewer than six months, although its five-year survival rate increased to 12% in recent years [1]. While it is currently the third leading cause of cancer-related deaths worldwide [1], the yearly increase in its incidence may make PDAC the second leading cause of cancer-related deaths by 2030 [1]. The absence of reliable biomarkers for effective screening and early diagnosis at the pre-symptomatic stages when treatments function most effectively represents a primary reason why most PDAC cases remain incurable. Currently, most patients present locally advanced (30–35%) or metastatic (50–55%) PDAC at diagnosis [2].
In advanced-stage PDAC patients, curative surgery remains impossible, and systemic therapeutic options (including immunotherapy) remain limited and ineffective [3]. Among the solid tumors, PDAC represents an immunologically “cold” tumor characterized by sparse T cell infiltration [4,5]; in contrast, immunologically “hot” tumors (such as melanoma) suffer from a high neoantigen load and immune cell infiltration [6]. PDAC tumors possess distinctive features such as an extracellular matrix (ECM) composition and a fibrotic stroma, which make it highly desmoplastic and significantly influence immune responses [7]. PDAC cells strongly interact with the surrounding microenvironment, which includes components such as immune cells, cytokines, metabolites, fibroblasts, and hyaluronan. These interactions create a highly fibrotic and active organized stroma (desmoplastic stroma) and an immunosuppressive environment that makes PDAC invasive and highly resistant to immunotherapy [5,8]; therefore, the characterization of the stroma and tumor immune microenvironment in PDAC patients represents a critical step in developing more effective therapeutic strategies. In the last few years, several investigations have focused on studying gene expression in PDAC to better understand the molecular composition of this devastating cancer and identify different molecular subtypes of pancreatic cancer that improve the stratification of patients for clinical strategies [9,10]. Bailey and colleagues defined four molecular subtypes of pancreatic cancer: squamous, pancreatic progenitor, immunogenic and aberrantly differentiated endocrine exocrine (ADEX) [10], while Moffitt’s group identified two stromal subtypes that were defined as “normal” and “activated” [9]. Nevertheless, in clinical practice, it is difficult to perform this broad molecular test on each patient. Therefore, and despite these new insights in pancreatic cancer, the diagnostic and prognostic outcomes of PDAC patients are extremely poor compared to those of other types of cancers. Additionally, new studies need to be conducted to understand the extreme complexity of PDAC and find simpler genetic signatures that can be incorporated into clinical practice and improve the clinical setting for PDAC patients and families.
We aimed to understand the stroma and tumor immune microenvironments of PDAC patients by retrieving and analyzing transcriptomic data from 21 different studies (representing a population of 922 samples; 320 controls and 602 cases) from the Gene Expression Omnibus (GEO)-NCBI and ArrayExpress data repositories. Through meta-analysis, we identified a series of gene signatures with survival prognostic value that may play a significant role in therapeutic decision making for PDAC patients, including five genes not previously related to PDAC survival. We also provide a friendly user web tool with detailed and interactive visualization of our comprehensive meta-analysis results.

2. Materials and Methods

For all bioinformatics and statistical analyses, we employed R software v. 4.1.3 [11] (Supplementary Table S1 details the R packages and versions).

2.1. Study Search and Selection

Publicly available datasets were collected from GEO-NCBI [12] and ArrayExpress databases [13]. Data available in the Cancer Genome Atlas (TCGA) [14] were excluded from the original search with the purpose of using this dataset as an external cohort for survival analysis. Following the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) guidelines [15], a systematic search of published studies was conducted in 2021 (period: 2002–2021). The protocol has not been registered. Three researchers in the study conducted the literature search (C.P.C., L.F., and I.P.D.), and the consistency of the review and selection procedures used was evaluated and confirmed. A broad search was performed using the MeSH (Medical Subject Headings) thesaurus keyword “pancreatic cancer”, after which stringent filters were applied. The final inclusion criteria were:
  • Normal and PDAC samples available.
  • RNA extracted directly from human pancreas biopsies.
  • Patients had not undergone treatment before biopsy.
  • Sample size > 4 for PDAC and control groups.
Finally, normalized gene expression from twenty-seven microarray studies (GSE86436, GSE71989, GSE62452, GSE62165, GSE60979, GSE56560, GSE55643, GSE46234, GSE43795, GSE43288, GSE41368, GSE32676, GSE28735, GSE27890, GSE22780, GSE19650, GSE18670, GSE16515, GSE15471, GSE1542, GSE11838, GSE102238, GSE101448, E-MTAB-3365, E-MTAB-1791, E-MEXP-950, and E-EMBL-6) and the count matrices of two RNA-sequencing (RNA-seq) (GSE119794 and GSE136569) datasets were retrieved for further analysis.

2.2. Individual Preprocessing and Analysis

Datasets were individually analyzed in two steps: preprocessing and differential expression analysis.
The nomenclature of clinical variables included in each study was standardized for data preprocessing, and then, exploratory analysis was performed. Prior to exploratory analysis, RNA-seq raw count matrices were normalized using the trimmed mean of m values from the edgeR package [16,17]. The normalization method performed by the original authors for each microarray dataset was assessed, and the matrices were log2 transformed when necessary. Exploratory analysis included expression boxplots, unsupervised clustering, and principal component analysis (PCA) to detect patterns of expression between samples and genes and the presence of batch effects in each study.
Differential gene expression analyses were performed in R using limma (v. 3.48.3) [18], and a paired sample design was implemented in those datasets where applicable. Differentially expressed genes were identified using p values with Benjamini-Hochberg correction [19] for a false discovery rate (FDR) at a significance level of 0.05.

2.3. Gene Expression Meta-Analysis

Gene expression analysis results were integrated into a meta-analysis using the DerSimonian & Laird random effects model [20], considering individual study heterogeneity. This model considers the variability in individual studies by increasing the weights of studies with less variability when meta-analysis results are computed.
A total of 24,365 genes were evaluated. p values, FDR-corrected p values, the logarithm of Fold Change (log2FC), and 95% confidence intervals of log2FC were calculated for each evaluated gene, and both funnel and forest plots were computed for each gene. These representations were assessed for possible biased results, where log2FC represents the effect size of a function, and the standard error of the log2FC serves as a study precision measure [21]. Genes were considered significant when FDR < 0.05, absolute log2FC > 0.6, and were measured in at least eleven studies. Sensitivity analysis (leave-one-out cross-validation [22]) was conducted for each significant gene to verify alterations in the results, owing to the inclusion of any study.
Statistically significant results from the gene expression meta-analysis were functionally enriched by over-representation analysis (ORA) using clusterProfiler [23,24] and ReactomePA [25]. Gene Ontology (GO) terms [26,27] and Reactome pathway [28] enrichment were performed following this approach. Only those functions and pathways with more than ten differentially expressed genes found in the gene set were considered. Functional enrichment was explored and visualized with the rrvgo package [29].

2.4. Web Tool

To make the data and results of our research widely accessible, a web tool was developed using the shiny package in R. The tool was developed in a user-friendly manner, allowing users to navigate and interact with the data. Users can then select different variables and parameters to visualize the data in numerous ways. The tool also includes interactive plots and tables to display the analysis results. The web tool is hosted on a secure server and is regularly maintained to ensure stability and performance. The source code for the tool is also publicly available and can be accessed through our GitHub repository: https://github.com/ipediez/ShinyReport (accessed on 1 May 2023).

2.5. Survival Analysis

RNA-seq expression data and metadata from patients in the Pancreatic adenocarcinoma (PAAD) TCGA cohort were downloaded from cBioPortal [30]. Z-scores of RNA-seq expression were used for survival analysis. For each analyzed gene, samples were divided into two groups based on their expression levels. Samples with expression Z-scores below the lower quartile were classified as having low levels of expression, whereas samples exceeding the upper quartile were classified as having high levels of expression. Forty-five samples with high levels of expression and forty-five samples with low levels of expression were included for survival analysis. Gene-wise Kaplan–Meier survival analysis compared the low-level and high-level expression groups. This method estimates the probability of survival over time based on the expression levels of the gene of interest. The log-rank test was used to compare the survival curves between distinct groups of samples.
For risk-score-based survival, genes were tagged as highly expressed for a given sample when the expression levels were above the upper quartile. Then, samples were clustered into “high-risk” and “low-risk” groups based on the number of highly expressed genes. The cutoff was set as the median of the highly expressed genes in each sample. Furthermore, a proportional hazard model using Cox regression was implemented to study the impact of clinicopathological variables on survival and evaluate the contribution of the risk score in a multivariate model.

3. Results

We performed a systematic review and differential gene expression analysis of PDAC transcriptomic studies from the GEO-NCBI [12] and ArrayExpress [13] databases to explore the stroma and immune environments in PDAC patients. We then integrated the results of each differential gene expression analysis into a meta-analysis. The biological context of the meta-analysis results was explored via functional enrichment using an ORA of GO terms and pathways (Figure 1). Finally, we conducted survival analysis to explore the impact of specific candidate genes on patients’ outcomes.

3.1. Systematic Review

The systematic review identified 143 non-duplicated studies. Then, we excluded studies with samples from patients under cancer treatment and studies where the sample size was less than four in the PDAC or the control group, resulting in a subset of twenty-nine studies (Figure 2). We discarded eight studies after exploratory analysis, giving a final set of twenty-one homogeneous and comparable studies for further analysis. The selected studies included 922 samples (320 controls and 602 cases). Although most studies did not include relevant sample metadata, we assessed the clinical characteristics when they were available. Supplementary Tables S2 and S3 contain further information regarding the selected studies and clinicopathological characteristics of the study population.

3.2. Integration of Differential Expression Profiles

Exploratory analysis found abnormal normalization or a lack of annotation in eight studies, which we excluded from further analysis (listed in Supplementary Table S4). Then, we performed the independent differential gene expression analysis of each study and meta-analysis for 24,365 genes evaluated in the different datasets, including every gene found in at least two studies. We considered results with an FDR < 0.05, an absolute log2FC > 0.6, and those evaluated in at least eleven studies to be significant; overall, 1153 genes met these criteria (Figure 3; further details are given in Supplementary Table S4).
We noted the presence of genes encoding ECM components (e.g., collagens, fibronectin, laminin, and stratifin), proteoglycans (e.g., versican), cell adhesion molecules, integrins, matrix metallopeptidases, and additional peptidases and enzymes that impact mechano-contractility, epithelial tension, and the stiffness of the tumoral stroma, which can promote tumor progression and resistance to therapy (Figure 4). Table 1 displays the twenty genes with the highest and lowest log2FC values from the meta-analysis; these genes mainly play roles in ECM remodeling, desmoplasia, metabolism, and the immune system. Supplementary Table S4 reports a complete list of significantly affected genes.
We performed ORA using GO biological process terms to identify the possible implications of 1153 significantly differentially expressed genes in the PDAC samples. We considered only those biological processes with at least ten associated genes and an adjusted p value under 0.05. We found 546 over-represented biological processes among the over-expressed genes and 40 biological processes over-represented among the under-expressed genes (Supplementary Table S5). ORA revealed the enrichment of terms related to the tumor microenvironment (Figure 5), with GO terms related to the immune system, cell adhesion, and ECM remodeling/degradation. Of note, additional over-represented functions were related to metastasis (vascularization, cell migration, collagen, mesenchymal transition, cell proliferation, and peptidyl modifications) [7,31].

3.3. Interactive Tool for Results Visualization

The web tool contains comprehensive information regarding the data and results of the meta-analysis of gene expression. The application includes tables and plots for the differential expression results of twenty-one datasets included in the study and meta-analysis results. Statistical indicators, such as the log odds ratio, confidence intervals, and adjusted p values, are provided to estimate each study’s global expression and specific contribution. The web tool is available online: https://bioinfo.cipf.es/MetaPDAC/ (accessed on 1 January 2020).

3.4. Immune System: A Functional Overview in PDAC

To focus our analysis on the tumor immune microenvironment, we extracted a consensus list of genes related to the immune system and inflammation from NCBI and GO databases (mainly framed in the categories of HLA, interleukin, CD, interferon, chemokine, and S100 genes, Supplementary Table S6). Considering an FDR threshold of 0.05 and an absolute fold change greater than 0.6, we discovered the significant differential expression of 322 immune genes in our meta-analysis results. To explore the functional involvement of these results, we performed ORA on this group of genes using GO biological process terms and Reactome pathways. We considered significant functional terms with at least ten associated genes and an adjusted p value < 0.05. We discovered the over-representation of thirty-three GO terms and twenty-seven pathways among the over-expressed immune-related genes and none when under-expressed genes were analyzed. The enriched terms suggest the increased activity of neutrophil-related immune response, the negative regulation of cell killing, interferon signaling, and an antigen presentation via major histocompatibility complex II.

3.5. Immune and Stromal Survival Signatures Impact PDAC Prognosis

We explored the 322 differentially expressed immune-related genes and identified a set of 70 genes of particular interest in our experimental research (Table 2). We performed survival analysis using the TCGA PAAD cohort for each of these genes and found statistically significant differences in twenty-eight genes (IFI27, IL1R2, IL1RN, IL1RAP, IL18, IL22RA1, HCP5, SLFN13, CD58, CD109, IFI44L, IFI16, IFITM1, IFIT1, IFIT3, IRF9, IFIT2, IFI35, CXCL10, CXCL5, CXCL9, S100P, S100A6, S100A2, S100A16, S100A11, S100A14, and S100A10), which shared a pattern: a higher expression in patients associated with a lower rate of survival. As far as we are aware, this is the first time that HCP5, SLFN13, IRF9, IFIT2, and IFI35 have been related to prognosis value in PDAC patients (Supplementary Figure S1).
We analyzed genes that displayed statistical significance as a “signature,” dividing the samples into high-risk and low-risk groups based on the number of highly expressed genes (above the upper quartile). We set the median (six highly expressed genes) as the cutoff value to divide the samples into groups. Interestingly, patients in the high-risk group possessed shorter survival times than those in the low-risk group did (p value < 0.0001, Figure 6A). Furthermore, we studied the effect of this signature in a multivariate Cox model including age, alcoholic history, the presence of chronic pancreatitis, diabetes diagnostic, tumor grade, and the AJCC classification of a metastatic tumor and a residual tumor as covariates. The proposed signature was the only variable with p value < 0.05 and showed a hazard ratio of 2.36 (Supplementary Figure S2). We then analyzed the co-occurrence of highly expressed genes in the samples, finding two main co-occurrence groups that related to high-risk patients: i) the interferon gene family (IFN genes) and ii) the S100 and IL genes (S100A14, S100A16, S100A6, S100A11, IL1R2, IL1RN, and S100P) (Figure 6B).
To explore how a desmoplastic environment can affect patients’ survival, we employed an homologous approach using genes related to ECM remodeling (Table 1). We discovered eleven genes whose survival analysis showed statistically significant differences (CEACAM5, CEACAM6, FN1, GJB2, GPRC5A, LAMB3, LAMC2, SFN, SLC6A14, TSPAN1, and VCAN). Again, we divided the samples into high-risk and low-risk groups using the median of the number of highly expressed genes as the cutoff value (median = 3). Patients with high levels of expression in three or more genes from the signature presented lower survival times than those with fewer highly expressed genes did (p value = 0.00012, Figure 7A). Of note, we distinguished a cluster of co-occurrence of patients with high levels of GJB2, FN1, and VCAN at the same time (Figure 7B).
Finally, we performed comparative analysis between the immune and stromal survival signatures identified in our work and other signatures generated in previous works for patient stratification [9,10]. These results provided insight into the level of intersection between this group of signatures (Supplementary Tables S7–S9).

4. Discussion

Using comprehensive meta-analysis, we explored the immune environment and desmoplastic stroma of PDAC tumors to contribute to a deeper understanding of tumorigenesis and the design of effective therapeutic strategies, such as immunotherapies. ECM components from the desmoplastic stroma tightly interact with the immune environment and contribute to immune evasion by modulating immune cell infiltration, thus influencing cell proliferation, tumor progression, and overall survival [32,33]. The meta-analysis and ORA results characterized differences in the gene-expression landscape of PDAC tumors and identified more than 1000 dysregulated genes, most of them with immune system- and desmoplasia-related roles. We discovered thirty-nine genes (twenty-eight immune-related genes and eleven stroma-related genes) that impact PDAC patients’ survival.
Among the top forty dysregulated genes (Table 1), we observed the upregulation of collagens (COL11A1 and COL10A1), which influence immune infiltration and chemoresistance and confer a poor prognosis [34,35,36]. PDAC patients also presented with upregulated periostin expression, which has been linked to a shorter overall survival [37], and cystatin SN, which contributes to pancreatic cancer cell proliferation and may represent a potential biomarker for the early detection of pancreatic cancer [38]. Stratifin and matrix metallopeptidase 1 also appeared to be upregulated in PDAC patients; stratifin stimulates matrix metallopeptidase 1 expression in fibroblasts, contributing to remodel ECM [39]. The increased expression of fibronectin in the PDAC stroma has also been reported. The observed upregulation of cathepsin E and sulfatase 1 expression in the PDAC microenvironment might also benefit the development of therapeutic strategies with polymer drug conjugates since they may contribute to drug release [40,41,42].
The analysis of the top forty dysregulated genes also provided evidence for the downregulation of genes coding for proteolytic enzymes released by the pancreas (e.g., chymotrypsin, chymotrypsinogen, lipases, and phospholipases). Pancreatic cancer cells express around 20% of chymotrypsin C normal cells expression, with this enzyme participating in cancer cell apoptosis and migration [43]. A recent report suggested that a combination of trypsinogen and chymotrypsinogen displayed an anti-tumorigenic potential [44].
Focusing on the immune environment, PDAC tumors develop a wide range of mechanisms to evade the immune system (e.g., a low level of expression of HLA antigens, immunosuppressive signals that inhibit natural killer and T cell functions, and the presence of immunosuppressive cells). This creates an immunotolerant environment in which the immune system of PDAC patients does not robustly recognize and target cancer cells [45]. We explored the expression of seventy genes of particular interest, including those from the HLA, interleukin, CD, interferon, chemokine, and S100 categories. The survival analysis of these genes in the TCGA PAAD cohort identified a twenty-eight immune-related gene signature with a prognostic value that was used to cluster PDAC patients into high-risk and low-risk groups.
The proposed signature possessed significance in univariate and multivariate Cox models with clinicopathological variables, significantly adding statistical power to the survival analysis. This signature could aid in the stratification of patients (Figure 8) who could benefit from immunotherapeutic strategies, given that it could contribute to distinguishing “cold” PDAC tumors (characterized by the low presence of T cells (CD8+) and natural killer cells, high presence of immunosuppressive cell populations, and poor prognoses and responses to immunotherapy) from “hot tumors” (with an opposite profile) [46,47]. We uncovered two high gene-expression co-occurrence patterns, one composed of IFN genes and the other of S100/IL genes. The IFN signaling pathways participate in PDAC development, while the over-expression of S100 genes blocks the infiltration and cytotoxic activity of CD8+ T cells, and the low level of expression of IL1RN and IL1R2 has been associated with increase survival in PDAC patients [48,49,50].
To the best of our knowledge, this is the first report of data suggesting a link between the HCP5, SLFN13, IRF9, IFIT2, and IFI35 immune genes and PDAC prognosis, presenting the discriminatory power of clustering PDAC patients. The remaining genes of the immune gene signature have been individually associated with PDAC or other cancers, with data suggesting that their overexpression could impact patients’ diagnosis, prognosis, and response to treatment [51,52,53,54,55,56]; however, we report that a joint gene expression signature of these genes impacts PDAC patients’ survival.
Focusing on the PDAC stroma, the altered genes include several types of collagens, fibronectins, and proteolytic enzymes, such as metalloproteases and peptidases (Table 1 and Supplementary Table S4), which significantly contribute to ECM composition and stromal remodeling and support desmoplasia and immunosuppression [57]. The survival analysis of significantly dysregulated stromal gene expression from the meta-analysis of the TCGA PAAD cohort revealed a gene signature with prognostic capacity that clustered PDAC patients into high-risk and low-risk groups. We observed a co-occurrence pattern in high-risk patients, indicating a subgroup of PDAC patients with a high level of expression of GJB2, FN1, and VCAN genes. These results indicate stromal heterogeneity in PDAC [58] and the need to characterize it to stratify patients (Figure 8).
With respect to other dysregulated genes, the upregulation of CEACAM5 and CEACAM6 represents an early event in pancreatic carcinogenesis, with these genes being candidates for immunotherapies [59,60,61]. Furthermore, laminins LAMBC2 and LAMB3 support cancer progression and resistance to gemcitabine—one of the main chemotherapeutics used in PDAC patients [62,63]. In general, the association of the stroma signature with a poor prognosis is consistent with the one described in previous studies for each gene: CEACAM5 [64], CEACAM6 [65], FN1 [66], GJB2 [67], GPRC5A [68], LAMB3 [69,70], LAMC2 [69,70], SFN [71], SLC6A14 [72], TSPAN1 [73], and VCAN [66].
With respect to other similar approaches, we are aware of two additional studies in which expression datasets were integrated to explore the nature of the PDAC in depth: one by Gooneskere and colleagues, who integrated six PDAC and three other pancreatic carcinomas datasets [74], and one by Irigoyen and colleagues, who integrated two peripheral blood datasets [75]. Both approaches integrate different datasets at the gene level to increase the number of samples and perform unique DGE analysis. In contrast, our approach analyzed each dataset independently, and then integrated the results, evaluating their robustness. From the experimental design point of view, both studies differ greatly from ours, since Grooneskere et al.’s one is not specifically focused on PDAC, and Irigoyen et al.’s one does not analyze pancreatic tissue. From a methodological point of view, our study contributes to a more profound and robust analysis of the PDAC expression landscape by integrating data after DGE analysis had been performed, thus avoiding the necessity to control heterogeneity among studies and retaining the full potential of biological differences.
Other molecular studies based on whole transcriptome and genomic analyses of pancreatic tumors have found specific gene signatures that identify different molecular subtypes [9,10]. However, the aim of this study was not to identify molecular subtypes, such as in the cited works. The immune signature or the stromal signature presented in this work establishes patient survival groups (high-risk group and low-risk group), which could help practitioners to decide if the patient could benefit from immunotherapy, for example, or not. Intersection analysis indicated that there is hardly any overlap between the gene signatures found in our study and the signatures described by Bailey et al. [10] or Moffitt et al. [9], as shown in the supplemental analysis (Supplementary Tables S7–S9). Therefore, the proposed gene signatures show subtype-independent survival value and display a reasonable number of genes for them to be translated to clinics. Nevertheless, more and deeper studies are needed for this purpose. Additionally, the works by Moffitt et al. and Bailey et al. are enormously rich and provide comprehensive molecular stratification to facilitate personalized treatment and the identification of therapeutic targets. Unfortunately, extensive molecular analyzes are difficult to translate to clinical practice for individual patients.
A potential limitation of our study has been the relative heterogeneity among the sample sizes and sequencing platforms used. The meta-analysis methodology, which integrates data groups and provides results with higher statistical power and precision [76,77], addresses this issue by independently comparing each study and combining the results. A lack of clinical and/or molecular information in most studies, such as survival time, stage condition, or molecular pattern, represents an additional limitation. We employed TCGA data for survival analysis, but additional analyses should integrate other covariates of interest in the study.
Finally, we provided an interactive web tool that allows users to explore our results, facilitating the accessibility, transparency, and reusability of our research. Overall, the web tool provides a detailed and interactive visualization of the meta-analysis results, allowing users to further explore and understand the gene expression patterns identified in the studies. Other functionalities include the capability to customize and filter the data to further investigate specific aspects of the analysis in more detail. In this manner, we aim to align our research with the FAIR principles to share our data in a way that can be of further use to the scientific community who studies this aggressive and lethal tumor.

5. Conclusions

Therapeutic strategies to overcome the immune microenvironment and the desmoplastic stroma barriers remain limited and generally unsuccessful. This study performs a comprehensive transcriptional meta-analysis of the molecular PDAC environment. The results highlight the relevance of the interaction between the immune system and stroma, revealing an impact on patients’ survival. The identified gene signatures provide new insights into the potential therapeutic targets for this deadly disease that can help to stratify its heterogeneity. Future studies are needed to explore the benefits of targeting the immune and stromal microenvironments as a treatment strategy for PDAC.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/cancers15112887/s1. Table S1: Software version. Table S2: Dataset inclusion. Table S3: Clinical characteristics. Table S4: Gene meta-analysis results. Table S5: ORA results. Table S6: NCBI and GO Immune system genes. Table S7: Gene signature comparisons. Table S8: Moffit-Basal-like. Table S9. Moffit-Classical.

Author Contributions

Conceptualization, F.G.-G.; data curation, I.P.-D., L.F., C.P.-C. and Z.A.; formal analysis, I.P.-D.; funding acquisition, J.A.L.-G. and F.G.-G.; investigation, I.P.-D., Z.A., M.R.H., J.A.L.-G., M.d.l.I.-V. and F.G.-G.; methodology, M.R.H., I.P.-D., Z.A. and F.G.-G.; project administration, F.G.-G.; software, I.P.-D.; supervision, J.A.L.-G. and F.G.-G.; validation, I.P.-D. and Z.A.; visualization, I.P.-D.; writing—original draft, I.P.-D., Z.A., M.R.H., J.A.L.-G., A.F.-S., M.d.l.I.-V. and F.G.-G.; writing—review and editing, I.P.-D., Z.A., M.R.H., J.A.L.-G., A.F.-S., M.d.l.I.-V. and F.G.-G. I.P.-D. and Z.A. contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by and partially funded by the Institute of Health Carlos III (project IMPaCT-Data, exp. IMP/00019), co-funded by the European Union, European Regional Development Fund (ERDF, “A way to make Europe”), and PID2021-124430OA-I00 funded by MCIN/AEI/10.13039/501100011033/FEDER, UE (“A way to make Europe”). This work has been carried out under the framework of the ULISES project: H2020-FETOPEN-2018-2019-2020-01 Contract nº: 899708.

Data Availability Statement

Publicly available datasets were analyzed in this study. Data are openly available in GEO and ArrayExpress, with the following accession numbers: GSE71989, GSE62452, GSE62165, GSE60979, GSE56560, GSE55643, GSE43795, GSE41368, GSE32676, GSE28735, GSE22780, GSE18670, GSE16515, GSE15471, GSE1542, GSE11838, GSE101448, GSE119794, GSE136569, E-MEXP-950, and E-EMBL-6.

Acknowledgments

The authors thank the Principe Felipe Research Center (CIPF) for providing access to the cluster, co-funded by European Regional Development Funds (ERDF) in Valencian Community 2014–2020. The authors also thank Stuart P. Atkinson for reviewing the manuscript. The results published here are in whole or part based upon data generated by the TCGA Research Network: https://www.cancer.gov/tcga (accessed on 1 January 2020).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Siegel, R.L.; Miller, K.D.; Wagle, N.S.; Jemal, A. Cancer statistics, 2023. CA Cancer J. Clin. 2023, 73, 17–48. [Google Scholar] [CrossRef]
  2. Park, W.; Chawla, A.; O’Reilly, E.M. Pancreatic Cancer: A Review. JAMA 2021, 326, 851–862. [Google Scholar] [CrossRef] [PubMed]
  3. Henriksen, A.; Dyhl-Polk, A.; Chen, I.; Nielsen, D. Checkpoint Inhibitors in Pancreatic Cancer. Cancer Treat. Rev. 2019, 78, 17–30. [Google Scholar] [CrossRef] [PubMed]
  4. Nezhad Shamohammadi, F.; Yazdanifar, M.; Oraei, M.; Kazemi, M.H.; Roohi, A.; Mahya Shariat Razavi, S.; Rezaei, F.; Parvizpour, F.; Karamlou, Y.; Namdari, H. Controversial Role of Γδ T Cells in Pancreatic Cancer. Int. Immunopharmacol. 2022, 108, 108895. [Google Scholar] [CrossRef] [PubMed]
  5. Ullman, N.A.; Burchard, P.R.; Dunne, R.F.; Linehan, D.C. Immunologic Strategies in Pancreatic Cancer: Making Cold Tumors Hot. J. Clin. Oncol. 2022, 40, 2789–2805. [Google Scholar] [CrossRef]
  6. Dong, C.; Dang, D.; Zhao, X.; Wang, Y.; Wang, Z.; Zhang, C. Integrative Characterization of the Role of IL27 In Melanoma Using Bioinformatics Analysis. Front. Immunol. 2021, 12, 713001. [Google Scholar] [CrossRef]
  7. Ostios-Garcia, L.; Villamayor, J.; Garcia-Lorenzo, E.; Vinal, D.; Feliu, J. Understanding the Immune Response and the Current Landscape of Immunotherapy in Pancreatic Cancer. World J. Gastroenterol. 2021, 27, 6775–6793. [Google Scholar] [CrossRef]
  8. Di Federico, A.; Mosca, M.; Pagani, R.; Carloni, R.; Frega, G.; De Giglio, A.; Rizzo, A.; Ricci, D.; Tavolari, S.; Di Marco, M.; et al. Immunotherapy in Pancreatic Cancer: Why Do We Keep Failing? A Focus on Tumor Immune Microenvironment, Predictive Biomarkers and Treatment Outcomes. Cancers 2022, 14, 2429. [Google Scholar] [CrossRef]
  9. Moffitt, R.; Marayati, R.; Flate, E.; Volmar, K.E.; Herrera Loeza, S.G.; Hoadley, K.A.; Rashid, N.U.; Williams, L.A.; Eaton, S.C.; Chung, A.H.; et al. Virtual microdissection identifies distinct tumor- and stroma-specific subtypes of pancreatic ductal adenocarcinoma. Nat Genet. 2015, 47, 1168–1178. [Google Scholar] [CrossRef]
  10. Bailey, P.; Chang, D.K.; Nones, K.; Johns, A.L.; Patch, A.; Gingras, M.; Miller, D.K.; Christ, A.N.; Bruxner, T.J.C.; Quinn, M.C.; et al. Genomic analyses identify molecular subtypes of pancreatic cancer. Nature 2016, 531, 47–52. [Google Scholar] [CrossRef]
  11. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2021. [Google Scholar]
  12. Edgar, R.; Domrachev, M.; Lash, A.E. Gene Expression Omnibus: NCBI Gene Expression and Hybridization Array Data Repository. Nucleic Acids Res. 2002, 30, 207–210. [Google Scholar] [CrossRef] [PubMed]
  13. Athar, A.; Füllgrabe, A.; George, N.; Iqbal, H.; Huerta, L.; Ali, A.; Snow, C.; Fonseca, N.A.; Petryszak, R.; Papatheodorou, I.; et al. ArrayExpress Update—From Bulk to Single-Cell Expression Data. Nucleic Acids Res. 2019, 47, D711–D715. [Google Scholar] [CrossRef]
  14. Cancer Genome Atlas Research Network; Weinstein, J.N.; Collisson, E.A.; Mills, G.B.; Shaw, K.R.M.; Ozenberger, B.A.; Ellrott, K.; Shmulevich, I.; Sander, C.; Stuart, J.M. The Cancer Genome Atlas Pan-Cancer Analysis Project. Nat. Genet. 2013, 45, 1113–1120. [Google Scholar] [CrossRef] [PubMed]
  15. Moher, D.; Liberati, A.; Tetzlaff, J.; Altman, D.G. PRISMA Group Preferred Reporting Items for Systematic Reviews and Meta-Analyses: The PRISMA Statement. PLoS Med. 2009, 6, e1000097. [Google Scholar] [CrossRef] [PubMed]
  16. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics 2010, 26, 139–140. [Google Scholar] [CrossRef]
  17. Robinson, M.D.; Oshlack, A. A Scaling Normalization Method for Differential Expression Analysis of RNA-Seq Data. Genome Biol. 2010, 11, R25. [Google Scholar] [CrossRef]
  18. Ritchie, M.E.; Phipson, B.; Wu, D.; Hu, Y.; Law, C.W.; Shi, W.; Smyth, G.K. Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 2015, 43, e47. [Google Scholar] [CrossRef]
  19. Benjamini, Y.; Hochberg, Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J. R. Stat. Soc. Ser. B Methodol. 1995, 57, 289–300. [Google Scholar] [CrossRef]
  20. DerSimonian, R.; Laird, N. Meta-Analysis in Clinical Trials. Control. Clin. Trials 1986, 7, 177–188. [Google Scholar] [CrossRef]
  21. Sterne, J.A.; Egger, M. Funnel Plots for Detecting Bias in Meta-Analysis: Guidelines on Choice of Axis. J. Clin. Epidemiol. 2001, 54, 1046–1055. [Google Scholar] [CrossRef]
  22. Viechtbauer, W. Conducting Meta-Analyses in R with the Metafor Package. J. Stat. Softw. 2010, 36, 1–48. [Google Scholar] [CrossRef]
  23. Yu, G.; Wang, L.-G.; Han, Y.; He, Q.-Y. ClusterProfiler: An R Package for Comparing Biological Themes among Gene Clusters. Omics J. Integr. Biol. 2012, 16, 284–287. [Google Scholar] [CrossRef] [PubMed]
  24. Wu, T.; Hu, E.; Xu, S.; Chen, M.; Guo, P.; Dai, Z.; Feng, T.; Zhou, L.; Tang, W.; Zhan, L.; et al. ClusterProfiler 4.0: A Universal Enrichment Tool for Interpreting Omics Data. Innovation 2021, 2, 100141. [Google Scholar] [CrossRef] [PubMed]
  25. Yu, G.; He, Q.-Y. ReactomePA: An R/Bioconductor Package for Reactome Pathway Analysis and Visualization. Mol. Biosyst. 2016, 12, 477–479. [Google Scholar] [CrossRef]
  26. Ashburner, M.; Ball, C.A.; Blake, J.A.; Botstein, D.; Butler, H.; Cherry, J.M.; Davis, A.P.; Dolinski, K.; Dwight, S.S.; Eppig, J.T.; et al. Gene Ontology: Tool for the Unification of Biology. The Gene Ontology Consortium. Nat. Genet. 2000, 25, 25–29. [Google Scholar] [CrossRef] [PubMed]
  27. Gene Ontology Consortium. The Gene Ontology Resource: Enriching a GOld Mine. Nucleic Acids Res. 2021, 49, D325–D334. [Google Scholar] [CrossRef] [PubMed]
  28. Gillespie, M.; Jassal, B.; Stephan, R.; Milacic, M.; Rothfels, K.; Senff-Ribeiro, A.; Griss, J.; Sevilla, C.; Matthews, L.; Gong, C.; et al. The Reactome Pathway Knowledgebase 2022. Nucleic Acids Res. 2022, 50, D687–D692. [Google Scholar] [CrossRef]
  29. Sayols, S. rrvgo: A Bioconductor Package to Reduce and Visualize Gene Ontology Terms. microPubl. Biol. 2023. [Google Scholar] [CrossRef]
  30. Cerami, E.; Gao, J.; Dogrusoz, U.; Gross, B.E.; Sumer, S.O.; Aksoy, B.A.; Jacobsen, A.; Byrne, C.J.; Heuer, M.L.; Larsson, E.; et al. The CBio Cancer Genomics Portal: An Open Platform for Exploring Multidimensional Cancer Genomics Data. Cancer Discov. 2012, 2, 401–404. [Google Scholar] [CrossRef]
  31. Hessmann, E.; Buchholz, S.M.; Demir, I.E.; Singh, S.K.; Gress, T.M.; Ellenrieder, V.; Neesse, A. Microenvironmental Determinants of Pancreatic Cancer. Physiol. Rev. 2020, 100, 1707–1751. [Google Scholar] [CrossRef]
  32. Whatcott, C.J.; Diep, C.H.; Jiang, P.; Watanabe, A.; LoBello, J.; Sima, C.; Hostetter, G.; Shepard, H.M.; Von Hoff, D.D.; Han, H. Desmoplasia in Primary Tumors and Metastatic Lesions of Pancreatic Cancer. Clin. Cancer Res. 2015, 21, 3561–3568. [Google Scholar] [CrossRef] [PubMed]
  33. Zhou, X.; Liu, Y.; Hu, M.; Wang, M.; Liu, X.; Huang, L. Relaxin Gene Delivery Modulates Macrophages to Resolve Cancer Fibrosis and Synergizes with Immune Checkpoint Blockade Therapy. Sci. Adv. 2021, 7, eabb6596. [Google Scholar] [CrossRef] [PubMed]
  34. Liu, Q.; Zhao, H.; Guo, Y.; Zhang, K.; Shang, F.; Liu, T. Bioinformatics-Based Analysis: Noncoding RNA-Mediated COL10A1 Is Associated with Poor Prognosis and Immune Cell Infiltration in Pancreatic Cancer. J. Healthc. Eng. 2022, 2022, 7904982. [Google Scholar] [CrossRef]
  35. Wang, H.; Ren, R.; Yang, Z.; Cai, J.; Du, S.; Shen, X. The COL11A1/Akt/CREB Signaling Axis Enables Mitochondrial-Mediated Apoptotic Evasion to Promote Chemoresistance in Pancreatic Cancer Cells through Modulating BAX/BCL-2 Function. J. Cancer 2021, 12, 1406–1420. [Google Scholar] [CrossRef] [PubMed]
  36. Zheng, X.; Liu, X.; Zheng, H.; Wang, H.; Hong, D. Integrated Bioinformatics Analysis Identified COL11A1 as an Immune Infiltrates Correlated Prognosticator in Pancreatic Adenocarcinoma. Int. Immunopharmacol. 2021, 90, 106982. [Google Scholar] [CrossRef] [PubMed]
  37. Neuzillet, C.; Nicolle, R.; Raffenne, J.; Tijeras-Raballand, A.; Brunel, A.; Astorgues-Xerri, L.; Vacher, S.; Arbateraz, F.; Fanjul, M.; Hilmi, M.; et al. Periostin- and Podoplanin-Positive Cancer-Associated Fibroblast Subtypes Cooperate to Shape the Inflamed Tumor Microenvironment in Aggressive Pancreatic Adenocarcinoma. J. Pathol. 2022, 258, 408–425. [Google Scholar] [CrossRef]
  38. Jiang, J.; Liu, H.-L.; Liu, Z.-H.; Tan, S.-W.; Wu, B. Identification of Cystatin SN as a Novel Biomarker for Pancreatic Cancer. Tumour Biol. 2015, 36, 3903–3910. [Google Scholar] [CrossRef]
  39. Chavez-Muñoz, C.; Morse, J.; Kilani, R.; Ghahary, A. Primary Human Keratinocytes Externalize Stratifin Protein via Exosomes. J. Cell. Biochem. 2008, 104, 2165–2173. [Google Scholar] [CrossRef]
  40. Mohamed, M.M.; Sloane, B.F. Cysteine Cathepsins: Multifunctional Enzymes in Cancer. Nat. Rev. Cancer 2006, 6, 764–775. [Google Scholar] [CrossRef]
  41. Berquin, I.M.; Sloane, B.F. Cathepsin B Expression in Human Tumors. Adv. Exp. Med. Biol. 1996, 389, 281–294. [Google Scholar] [CrossRef]
  42. Atkinson, S.P.; Andreu, Z.; Vicent, M.J. Polymer Therapeutics: Biomarkers and New Approaches for Personalized Cancer Treatment. J. Pers. Med. 2018, 8, 6. [Google Scholar] [CrossRef] [PubMed]
  43. Wang, H.; Sha, W.; Liu, Z.; Chi, C.-W. Effect of Chymotrypsin C and Related Proteins on Pancreatic Cancer Cell Migration. Acta Biochim. Biophys. Sin. 2011, 43, 362–371. [Google Scholar] [CrossRef] [PubMed]
  44. González-Titos, A.; Hernández-Camarero, P.; Barungi, S.; Marchal, J.A.; Kenyon, J.; Perán, M. Trypsinogen and Chymotrypsinogen: Potent Anti-Tumor Agents. Expert Opin. Biol. Ther. 2021, 21, 1609–1621. [Google Scholar] [CrossRef] [PubMed]
  45. Makkouk, A.; Weiner, G.J. Cancer Immunotherapy and Breaking Immune Tolerance: New Approaches to an Old Challenge. Cancer Res. 2015, 75, 5–10. [Google Scholar] [CrossRef]
  46. Liu, Y.-T.; Sun, Z.-J. Turning Cold Tumors into Hot Tumors by Improving T-Cell Infiltration. Theranostics 2021, 11, 5365–5386. [Google Scholar] [CrossRef]
  47. Rubin, S.J.S.; Sojwal, R.S.; Gubatan, J.; Rogalla, S. The Tumor Immune Microenvironment in Pancreatic Ductal Adenocarcinoma: Neither Hot nor Cold. Cancers 2022, 14, 4236. [Google Scholar] [CrossRef]
  48. Zhuang, H.; Chen, X.; Dong, F.; Zhang, Z.; Zhou, Z.; Ma, Z.; Huang, S.; Chen, B.; Zhang, C.; Hou, B. Prognostic Values and Immune Suppression of the S100A Family in Pancreatic Cancer. J. Cell. Mol. Med. 2021, 25, 3006–3018. [Google Scholar] [CrossRef]
  49. Fujisawa, M.; Kanda, T.; Shibata, T.; Sasaki, R.; Masuzaki, R.; Matsumoto, N.; Nirei, K.; Imazu, H.; Kuroda, K.; Sugitani, M.; et al. Involvement of the Interferon Signaling Pathways in Pancreatic Cancer Cells. Anticancer Res. 2020, 40, 4445–4455. [Google Scholar] [CrossRef]
  50. Herremans, K.M.; Szymkiewicz, D.D.; Riner, A.N.; Bohan, R.P.; Tushoski, G.W.; Davidson, A.M.; Lou, X.; Leong, M.C.; Dean, B.D.; Gerber, M.; et al. The Interleukin-1 Axis and the Tumor Immune Microenvironment in Pancreatic Ductal Adenocarcinoma. Neoplasia 2022, 28, 100789. [Google Scholar] [CrossRef]
  51. Yuan, B.; Guan, Q.; Yan, T.; Zhang, X.; Xu, W.; Li, J. LncRNA HCP5 Regulates Pancreatic Cancer Progression by MiR-140-5p/CDK8 Axis. Cancer Biother. Radiopharm. 2020, 35, 711–719. [Google Scholar] [CrossRef]
  52. Liu, Y.; Wang, J.; Dong, L.; Xia, L.; Zhu, H.; Li, Z.; Yu, X. Long Noncoding RNA HCP5 Regulates Pancreatic Cancer Gemcitabine (GEM) Resistance By Sponging Hsa-MiR-214-3p To Target HDGF. OncoTargets Ther. 2019, 12, 8207–8216. [Google Scholar] [CrossRef] [PubMed]
  53. Xu, J.; Chen, S.; Liang, J.; Hao, T.; Wang, H.; Liu, G.; Jin, X.; Li, H.; Zhang, J.; Zhang, C.; et al. Schlafen Family Is a Prognostic Biomarker and Corresponds with Immune Infiltration in Gastric Cancer. Front. Immunol. 2022, 13, 922138. [Google Scholar] [CrossRef]
  54. Rodolosse, A.; Chalaux, E.; Adell, T.; Hagège, H.; Skoudy, A.; Real, F.X. PTF1alpha/P48 Transcription Factor Couples Proliferation and Differentiation in the Exocrine Pancreas [Corrected]. Gastroenterology 2004, 127, 937–949. [Google Scholar] [CrossRef] [PubMed]
  55. Hu, Y.; Wang, B.; Yi, K.; Lei, Q.; Wang, G.; Xu, X. IFI35 Is Involved in the Regulation of the Radiosensitivity of Colorectal Cancer Cells. Cancer Cell Int. 2021, 21, 290. [Google Scholar] [CrossRef] [PubMed]
  56. Shen, H.; Zhan, M.; Zhang, Y.; Huang, S.; Xu, S.; Huang, X.; He, M.; Yao, Y.; Man, M.; Wang, J. PLZF Inhibits Proliferation and Metastasis of Gallbladder Cancer by Regulating IFIT2. Cell Death Dis. 2018, 9, 71. [Google Scholar] [CrossRef]
  57. Ho, W.J.; Jaffee, E.M.; Zheng, L. The Tumour Microenvironment in Pancreatic Cancer—Clinical Challenges and Opportunities. Nat. Rev. Clin. Oncol. 2020, 17, 527–540. [Google Scholar] [CrossRef]
  58. Hosein, A.N.; Brekken, R.A.; Maitra, A. Pancreatic Cancer Stroma: An Update on Therapeutic Targeting Strategies. Nat. Rev. Gastroenterol. Hepatol. 2020, 17, 487–505. [Google Scholar] [CrossRef]
  59. Zińczuk, J.; Zaręba, K.; Romaniuk, W.; Kamińska, D.; Nizioł, M.; Baszun, M.; Kędra, B.; Guzińska-Ustymowicz, K.; Pryczynicz, A. Expression of Chosen Carcinoembryonic-Related Cell Adhesion Molecules in Pancreatic Intraepithelial Neoplasia (PanIN) Associated with Chronic Pancreatitis and Pancreatic Ductal Adenocarcinoma (PDAC). Int. J. Med. Sci. 2019, 16, 583–592. [Google Scholar] [CrossRef]
  60. Rizeq, B.; Zakaria, Z.; Ouhtit, A. Towards Understanding the Mechanisms of Actions of Carcinoembryonic Antigen-Related Cell Adhesion Molecule 6 in Cancer Progression. Cancer Sci. 2018, 109, 33–42. [Google Scholar] [CrossRef]
  61. Han, Z.-W.; Lyv, Z.-W.; Cui, B.; Wang, Y.-Y.; Cheng, J.-T.; Zhang, Y.; Cai, W.-Q.; Zhou, Y.; Ma, Z.-W.; Wang, X.-W.; et al. The Old CEACAMs Find Their New Role in Tumor Immunotherapy. Investig. New Drugs 2020, 38, 1888–1898. [Google Scholar] [CrossRef]
  62. Okada, Y.; Takahashi, N.; Takayama, T.; Goel, A. LAMC2 Promotes Cancer Progression and Gemcitabine Resistance through Modulation of EMT and ATP-Binding Cassette Transporters in Pancreatic Ductal Adenocarcinoma. Carcinogenesis 2021, 42, 546–556. [Google Scholar] [CrossRef]
  63. Zhang, H.; Pan, Y.-Z.; Cheung, M.; Cao, M.; Yu, C.; Chen, L.; Zhan, L.; He, Z.-W.; Sun, C.-Y. LAMB3 Mediates Apoptotic, Proliferative, Invasive, and Metastatic Behaviors in Pancreatic Cancer by Regulating the PI3K/Akt Signaling Pathway. Cell Death Dis. 2019, 10, 230. [Google Scholar] [CrossRef] [PubMed]
  64. Lu, Y.; Li, D.; Liu, G.; Xiao, E.; Mu, S.; Pan, Y.; Qin, F.; Zhai, Y.; Duan, S.; Li, D.; et al. Identification of Critical Pathways and Potential Key Genes in Poorly Differentiated Pancreatic Adenocarcinoma. OncoTargets Ther. 2021, 14, 711–723. [Google Scholar] [CrossRef] [PubMed]
  65. Johnson, B.; Mahadevan, D. Emerging Role and Targeting of Carcinoembryonic Antigen-Related Cell Adhesion Molecule 6 (CEACAM6) in Human Malignancies. Clin. Cancer Drugs 2015, 2, 100–111. [Google Scholar] [CrossRef] [PubMed]
  66. Lei, X.; Chen, G.; Li, J.; Wen, W.; Gong, J.; Fu, J. Comprehensive Analysis of Abnormal Expression, Prognostic Value and Oncogenic Role of the Hub Gene FN1 in Pancreatic Ductal Adenocarcinoma via Bioinformatic Analysis and in Vitro Experiments. PeerJ 2021, 9, e12141. [Google Scholar] [CrossRef] [PubMed]
  67. Fukuhisa, H.; Seki, N.; Idichi, T.; Kurahara, H.; Yamada, Y.; Toda, H.; Kita, Y.; Kawasaki, Y.; Tanoue, K.; Mataki, Y.; et al. Gene Regulation by Antitumor MiR-130b-5p in Pancreatic Ductal Adenocarcinoma: The Clinical Significance of Oncogenic EPS8. J. Hum. Genet. 2019, 64, 521–534. [Google Scholar] [CrossRef]
  68. Jahny, E.; Yang, H.; Liu, B.; Jahnke, B.; Lademann, F.; Knösel, T.; Rümmele, P.; Grützmann, R.; Aust, D.E.; Pilarsky, C.; et al. The G Protein-Coupled Receptor RAI3 Is an Independent Prognostic Factor for Pancreatic Cancer Survival and Regulates Proliferation via STAT3 Phosphorylation. PLoS ONE 2017, 12, e0170390. [Google Scholar] [CrossRef]
  69. Islam, S.; Kitagawa, T.; Baron, B.; Abiko, Y.; Chiba, I.; Kuramitsu, Y. ITGA2, LAMB3, and LAMC2 May Be the Potential Therapeutic Targets in Pancreatic Ductal Adenocarcinoma: An Integrated Bioinformatics Analysis. Sci. Rep. 2021, 11, 10563. [Google Scholar] [CrossRef]
  70. Lin, H.; Yang, P.; Li, B.; Chang, Y.; Chen, Y.; Li, Y.; Liu, K.; Liang, X.; Chen, T.; Dai, Y.; et al. S100A10 Promotes Pancreatic Ductal Adenocarcinoma Cells Proliferation, Migration and Adhesion through JNK/LAMB3-LAMC2 Axis. Cancers 2022, 15, 202. [Google Scholar] [CrossRef]
  71. Robin, F.; Angenard, G.; Cano, L.; Courtin-Tanguy, L.; Gaignard, E.; Khene, Z.-E.; Bergeat, D.; Clément, B.; Boudjema, K.; Coulouarn, C.; et al. Molecular Profiling of Stroma Highlights Stratifin as a Novel Biomarker of Poor Prognosis in Pancreatic Ductal Adenocarcinoma. Br. J. Cancer 2020, 123, 72–80. [Google Scholar] [CrossRef]
  72. Schniers, B.K.; Wachtel, M.S.; Sharma, M.; Korac, K.; Rajasekaran, D.; Yang, S.; Sniegowski, T.; Ganapathy, V.; Bhutia, Y.D. Deletion of Slc6a14 Reduces Cancer Growth and Metastatic Spread and Improves Survival in KPC Mouse Model of Spontaneous Pancreatic Cancer. Biochem. J. 2022, 479, 719–730. [Google Scholar] [CrossRef] [PubMed]
  73. Zhou, C.; Liang, Y.; Zhou, L.; Yan, Y.; Liu, N.; Zhang, R.; Huang, Y.; Wang, M.; Tang, Y.; Ali, D.W.; et al. TSPAN1 Promotes Autophagy Flux and Mediates Cooperation between WNT-CTNNB1 Signaling and Autophagy via the MIR454-FAM83A-TSPAN1 Axis in Pancreatic Cancer. Autophagy 2021, 17, 3175–3195. [Google Scholar] [CrossRef] [PubMed]
  74. Goonesekere, N.C.W.; Wang, X.; Ludwig, L.; Guda, C. A Meta Analysis of Pancreatic Microarray Datasets Yields New Targets as Cancer Genes and Biomarkers. PLoS ONE 2014, 9, e93046. [Google Scholar] [CrossRef] [PubMed]
  75. Irigoyen, A.; Jimenez-Luna, C.; Benavides, M.; Caba, O.; Gallego, J.; Ortuño, F.M.; Guillen-Ponce, C.; Rojas, I.; Aranda, E.; Torres, C.; et al. Integrative Multi-Platform Meta-Analysis of Gene Expression Profiles in Pancreatic Ductal Adenocarcinoma Patients for Identifying Novel Diagnostic Biomarkers. PLoS ONE 2018, 13, e0194844. [Google Scholar] [CrossRef]
  76. Normand, S.L. Meta-Analysis: Formulating, Evaluating, Combining, and Reporting. Stat. Med. 1999, 18, 321–359. [Google Scholar] [CrossRef]
  77. Higgins, J.P.T.; Thomas, J.; Chandler, J.; Cumpston, M.; Li, T.; Page, M.J.; Welch, V.A. Cochrane Handbook for Systematic Reviews of Interventions; John Wiley & Sons: Hoboken, NJ, USA, 2019; ISBN 978-1-119-53661-1. [Google Scholar]
Cancers 15 02887 g001 550
Figure 1. Workflow and analysis design. Relevant studies from the GEO-NCBI and ArrayExpress databases were retrieved, and data exploration and preprocessing were then performed. After DGE analysis, the results from different studies were integrated into a gene meta-analysis. Functional profiling methodologies were applied to explore the biological implications of the results.
Figure 1. Workflow and analysis design. Relevant studies from the GEO-NCBI and ArrayExpress databases were retrieved, and data exploration and preprocessing were then performed. After DGE analysis, the results from different studies were integrated into a gene meta-analysis. Functional profiling methodologies were applied to explore the biological implications of the results.
Cancers 15 02887 g001
Cancers 15 02887 g002 550
Figure 2. Flow of information through the distinct phases of the systematic review, following PRISMA Statement guidelines.
Figure 2. Flow of information through the distinct phases of the systematic review, following PRISMA Statement guidelines.
Cancers 15 02887 g002
Cancers 15 02887 g003 550
Figure 3. Volcano plot summarizing the gene expression meta-analysis. Significantly over-expressed genes are shown in red, and significantly under-expressed genes are shown in blue (FDR < 0.05; absolute log2FC > 0.6). Genes that do not show significant differential expression are represented in black. Only genes found in at least eleven studies are shown.
Figure 3. Volcano plot summarizing the gene expression meta-analysis. Significantly over-expressed genes are shown in red, and significantly under-expressed genes are shown in blue (FDR < 0.05; absolute log2FC > 0.6). Genes that do not show significant differential expression are represented in black. Only genes found in at least eleven studies are shown.
Cancers 15 02887 g003
Cancers 15 02887 g004 550
Figure 4. Overview of PDAC microenvironment. Meta-analysis results indicated an overexpression of several ECM components, e.g., stratifin, fibronectin 1, different laminin subtypes (gamma2 and beta3), collagens, and proteoglycans that characterize the dense and desmoplastic stroma of PDAC tumors. Additionally, the results highlight the presence of immune components such as IFN27, which contribute to an increase in the number of M2 macrophages and a decrease in the number of CD8+ T cells. Therefore, the desmoplastic stroma and the immune system favor immune tolerance and poor prognosis in PDAC. The red upward-pointing arrows denote genes exhibiting significant overexpression in the conducted meta-analysis. IFN27: interferon alpha inducible protein; MMP1: matrix metallopeptidase 1; NK cells: natural killer cells; T cells: T effector lymphocytes; Tregs: T regulatory lymphocytes T.
Figure 4. Overview of PDAC microenvironment. Meta-analysis results indicated an overexpression of several ECM components, e.g., stratifin, fibronectin 1, different laminin subtypes (gamma2 and beta3), collagens, and proteoglycans that characterize the dense and desmoplastic stroma of PDAC tumors. Additionally, the results highlight the presence of immune components such as IFN27, which contribute to an increase in the number of M2 macrophages and a decrease in the number of CD8+ T cells. Therefore, the desmoplastic stroma and the immune system favor immune tolerance and poor prognosis in PDAC. The red upward-pointing arrows denote genes exhibiting significant overexpression in the conducted meta-analysis. IFN27: interferon alpha inducible protein; MMP1: matrix metallopeptidase 1; NK cells: natural killer cells; T cells: T effector lymphocytes; Tregs: T regulatory lymphocytes T.
Cancers 15 02887 g004
Cancers 15 02887 g005 550
Figure 5. Scatter plot of ORA results. The scatterplot reports the GO biological process representative terms after redundancy reduction in a two-dimensional space derived from the semantic similarities between GO terms. The dot size represents the number of biological processes related to a GO term. The parent terms of the main clusters are labeled.
Figure 5. Scatter plot of ORA results. The scatterplot reports the GO biological process representative terms after redundancy reduction in a two-dimensional space derived from the semantic similarities between GO terms. The dot size represents the number of biological processes related to a GO term. The parent terms of the main clusters are labeled.
Cancers 15 02887 g005
Cancers 15 02887 g006 550
Figure 6. Survival analysis of immune system genes. A twenty-eight gene signature clustered patients into high-risk or low-risk groups based on the number of highly expressed signature genes in their transcriptomic profile. Patients with at least six highly expressed genes were classified as having a high risk, whereas those with five or fewer were classified as having a low risk. (A) Kaplan–Meier curve. Patients from the high-risk group (red) had shorter survival times than patients from the low-risk group did (blue). Below, the number of still alive patients and percentage in each group at 0, 25, 50, 75, and 100 months, and the censored events. (B) Heatmap demonstrating the patterns of high expression between genes and samples. Gene expression was coded as 1 for a sample above the upper quartile.
Figure 6. Survival analysis of immune system genes. A twenty-eight gene signature clustered patients into high-risk or low-risk groups based on the number of highly expressed signature genes in their transcriptomic profile. Patients with at least six highly expressed genes were classified as having a high risk, whereas those with five or fewer were classified as having a low risk. (A) Kaplan–Meier curve. Patients from the high-risk group (red) had shorter survival times than patients from the low-risk group did (blue). Below, the number of still alive patients and percentage in each group at 0, 25, 50, 75, and 100 months, and the censored events. (B) Heatmap demonstrating the patterns of high expression between genes and samples. Gene expression was coded as 1 for a sample above the upper quartile.
Cancers 15 02887 g006
Cancers 15 02887 g007 550
Figure 7. Survival analysis of ECM remodeling genes. An eleven-gene signature clustered patients into high-risk or low-risk groups based on the number of highly expressed signature genes in their transcriptomic profile. Patients with at least three highly expressed genes were classified as having a high risk, whereas those with five or fewer were classified as having a low risk. (A) Kaplan–Meier curve. Patients from the high-risk group (red) had shorter survival times than patients from the low-risk group did (blue). Below, the number of still alive patients and percentage in each group at 0, 25, 50, 75, and 100 months, and the censored events. (B) Heatmap demonstrating the patterns of high expression between genes and samples. Gene expression was coded as 1 for a sample above the upper quartile.
Figure 7. Survival analysis of ECM remodeling genes. An eleven-gene signature clustered patients into high-risk or low-risk groups based on the number of highly expressed signature genes in their transcriptomic profile. Patients with at least three highly expressed genes were classified as having a high risk, whereas those with five or fewer were classified as having a low risk. (A) Kaplan–Meier curve. Patients from the high-risk group (red) had shorter survival times than patients from the low-risk group did (blue). Below, the number of still alive patients and percentage in each group at 0, 25, 50, 75, and 100 months, and the censored events. (B) Heatmap demonstrating the patterns of high expression between genes and samples. Gene expression was coded as 1 for a sample above the upper quartile.
Cancers 15 02887 g007
Cancers 15 02887 g008 550
Figure 8. Patient stratification based on PDAC molecular features. Meta-analysis from transcriptomic studies allows a better understanding of the PDAC environment. In this study, the found gene signatures might contribute to the stratification of PDAC patients. In a first step, the immune or the stroma gene signatures can divide patients into high- and low-risk populations. After, with a focus on the immune signature co-occurrence, patients could be divided into those with more S100/IL genes and those with a more IFN expressed genes. The knowledge about these molecular features of PDAC tumors may guide the design of more effective therapeutic strategies.
Figure 8. Patient stratification based on PDAC molecular features. Meta-analysis from transcriptomic studies allows a better understanding of the PDAC environment. In this study, the found gene signatures might contribute to the stratification of PDAC patients. In a first step, the immune or the stroma gene signatures can divide patients into high- and low-risk populations. After, with a focus on the immune signature co-occurrence, patients could be divided into those with more S100/IL genes and those with a more IFN expressed genes. The knowledge about these molecular features of PDAC tumors may guide the design of more effective therapeutic strategies.
Cancers 15 02887 g008
Table
Table 1. Top twenty genes up- and down-regulated in PDAC patients.
Table 1. Top twenty genes up- and down-regulated in PDAC patients.
Gene SymbolGene NameExpression LevelFunction
CEACAM6CEA cell adhesion molecule 6UPEMR
SLC6A14Solute carrier family 6 member 14UPEMR
S100PS100 calcium-binding protein PUPEMR
CTSECathepsin EUPEMR
SULF1Sulfatase 1UPEMR
POSTNPeriostinUPEMR
GJB2Gap junction protein beta 2UPEMR
GPRC5AG protein-coupled receptor class C group 5 member AUPEMR
SFNStratifinUPEMR
FN1Fibronectin 1UPEMR
LAMC2Laminin subunit gamma 2UPEMR
CEACAM5CEA cell adhesion molecule 5UPEMR
MMP1Matrix metallopeptidase 1UPEMR
COL11A1Collagen type XI alpha 1 chainUPEMR
TSPAN1Tetraspanin 1UPEMR
IFI27Interferon alpha inducible Protein 27UPIS
CST1Cystatin SNUPEMT
LAMB3Laminin subunit beta 3UPEMR
COL10A1Collagen type X alpha 1 chainUPEMR
VCANVersicanUPEMR
CTRB2Chymotrypsinogen B2DOWNEMR
PLA2G1BPhospholipase A2 group IBDOWNMetabolism
CTRCChymotrypsin CDOWNEMR
GNMTGlycine N-methyltransferaseDOWNMetabolism
AQP8Aquaporin 8DOWNH2O2 transport
SYCNSyncolinDOWNExocytosis
CPA2Carboxypeptidase A2DOWNMetabolism
CELA2AChymotrypsin-like elastase 2ADOWNEMR
GP2Glycoprotein 2DOWNMetabolism
KLK1Kallikrein 1DOWNSerine protease
ALBAlbuminDOWNOncotic pressure
CTRB1Chymotrypsinogen B1DOWNEMR
ERP27Endoplasmic reticulum protein 27DOWNLipid and protein synthesis
TMED6Transmembrane p24 trafficking protein 6DOWNInsulin secretion
PNLIPRP1Pancreatic lipase-related protein 1DOWNMetabolism
CUZD1CUB and zona pellucida-like domain 1DOWNEMR and IS
CELA2BChymotrypsin-like elastase 2BDOWNEMR
PNLIPRP2Pancreatic lipase-related protein 2DOWNMetabolism
CTRLChymotrypsin-likeDOWNEMR
SERPINI2Serpin family I member 2DOWNProtease inhibitor
EMR = ECM remodeling; IS = immune system; EMT = epithelial–mesenchymal transition.
Table
Table 2. Subset of immune-related genes.
Table 2. Subset of immune-related genes.
Functional GroupGenes
HLAHLA-F, HLA-DRB5, HLA-B, HLA-A, HCP5, HLA-DRA, HLA-DPA1, HLA-DQB1, HLA-DQA1, HLA-DMB, HLA-DRB1, HLA-G, HLA-DPB1, SLFN12, SLFN13, and SLFN11
InterleukinIL1R2, IL1RN, IL1RAP, IL7R, IL2RG, IRAK3, IL18, LIF, and IL22RA1
CDCD58, CD109, CD52, CD53, CD74, CD14, CCDC80, CCDC141, CCDC69, DCDC2, and PDCD4
InterferonIFI27, IFI44L, IFI6, STING1, IFI16, IFITM1, ISG20, IFIT1, IFIT3, IFITM2, IRF9, IFIT2, IFNGR2, IFITM3, and IFI35
ChemokineCCL20, CCL18, CXCL10, CXCL5, CXCL8, CXCR4, CKLF, CXCL9, CXCL3, CXCL14, and CXCL12
S100S100P, S100A6, S100A2, S100A16, S100A11, S100A4, S100A14, and S100A10
Genes in bold possess statistically significant differences according to survival analysis.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Pérez-Díez, I.; Andreu, Z.; Hidalgo, M.R.; Perpiñá-Clérigues, C.; Fantín, L.; Fernandez-Serra, A.; de la Iglesia-Vaya, M.; Lopez-Guerrero, J.A.; García-García, F. A Comprehensive Transcriptional Signature in Pancreatic Ductal Adenocarcinoma Reveals New Insights into the Immune and Desmoplastic Microenvironments. Cancers 2023, 15, 2887. https://doi.org/10.3390/cancers15112887

AMA Style

Pérez-Díez I, Andreu Z, Hidalgo MR, Perpiñá-Clérigues C, Fantín L, Fernandez-Serra A, de la Iglesia-Vaya M, Lopez-Guerrero JA, García-García F. A Comprehensive Transcriptional Signature in Pancreatic Ductal Adenocarcinoma Reveals New Insights into the Immune and Desmoplastic Microenvironments. Cancers. 2023; 15(11):2887. https://doi.org/10.3390/cancers15112887

Chicago/Turabian Style

Pérez-Díez, Irene, Zoraida Andreu, Marta R. Hidalgo, Carla Perpiñá-Clérigues, Lucía Fantín, Antonio Fernandez-Serra, María de la Iglesia-Vaya, José A. Lopez-Guerrero, and Francisco García-García. 2023. "A Comprehensive Transcriptional Signature in Pancreatic Ductal Adenocarcinoma Reveals New Insights into the Immune and Desmoplastic Microenvironments" Cancers 15, no. 11: 2887. https://doi.org/10.3390/cancers15112887

APA Style

Pérez-Díez, I., Andreu, Z., Hidalgo, M. R., Perpiñá-Clérigues, C., Fantín, L., Fernandez-Serra, A., de la Iglesia-Vaya, M., Lopez-Guerrero, J. A., & García-García, F. (2023). A Comprehensive Transcriptional Signature in Pancreatic Ductal Adenocarcinoma Reveals New Insights into the Immune and Desmoplastic Microenvironments. Cancers, 15(11), 2887. https://doi.org/10.3390/cancers15112887

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop