Integrated transcriptome meta-analysis of pancreatic ductal adenocarcinoma and matched adjacent pancreatic tissues

Selection of microarray datasets

NCBI Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) and ArrayExpress Archive of Functional Genomics Data were systematically searched for eligible datasets using the keyword “pancreatic ductal adenocarcinoma”. The inclusion criteria were: (i) gene expression microarray data, (ii) human-derived pancreatic ductal adenocarcinoma tissues and matched adjacent non-tumor tissue samples. When tumor-adjacent tissue pairs were not specified clearly in the overall design of the study or the sample description, the authors of the relevant paper were consulted for the confirmation of sample pairs. Non-confirmed or unspecified sample pairs were excluded from the meta-analysis.

Generation of gene expression matrix files and evaluation of data quality

All data processing and integration procedures were performed using ExAtlas, which is an online software tool for meta-analysis and visualization of gene expression data (Sharov, Schlessinger & Ko, 2015). Briefly, the datasets that have been selected for the meta-analysis were uploaded to ExAtlas. Unpaired samples were removed from the sample files, and then the gene expression matrix file was generated from each dataset separately. All the extracted data had been normalized using the RMA algorithm. Individual sample quality was evaluated by checking the correlation of log10-transformed expression level with other data for a set of pre-selected housekeeping genes and the level of the global standard deviation. Samples, where the correlation of expression of housekeeping genes in the range from 0.5 to 0.95, and the level of standard deviation from the global mean for each set of genes grouped by the average expression is less than 0.3 were considered of good quality. Samples of low quality were removed from the datasets.

Standard meta-analysis

In the Pairwise comparison section of ExAtlas (Sharov, Schlessinger & Ko, 2015), one of the tumor gene expression profiles was added as a sample for examination, and its adjacent non-tumor tissue pair was added for baseline control. Then, the meta-analysis section was used to add more gene expression profile pairs. The random-effects method (DerSimonian & Laird, 1986), which takes into account the variance of heterogeneity between studies, was used to perform the meta-analysis. False discovery rate (FDR) is less than 0.05 and the change of gene expression is ≥2-fold were considered significant. The analysis was performed for each gene symbol, and the effects were presented as combined fold changes and combined log-ratios (log10).

External validation of the identified DEGs

The external validation was done using GEPIA Database (http://gepia.cancer-pku.cn/index.html) (Tang et al., 2017) by comparing transcriptomic data from the TCGA PAAD (pancreatic adenocarcinoma), and the TCGA normal and GTEx data. For genes, whose probes are not found in the validation dataset, Logsdon pancreas dataset (available at the Oncomine Database; Rhodes et al., 2004) including data from pancreatic ductal adenocarcinoma and healthy pancreatic tissues was used. The external validation in a study including PDAC and matched adjacent tissues could not be performed because all studies that passed the inclusion criteria were included in the meta-analysis.

Transcription factor binding site enrichment analysis

TFactS Database (http://www.tfacts.org) (Essaghir & Demoulin, 2012; Essaghir et al., 2010) was used as a tool to predict which transcription factors (TF) are regulated, inhibited, or activated based on the list of DEGs. Briefly, the up-and down-regulated genes were uploaded to the TFact database, and the DEGs were then compared with the sign-sensitive catalog of validated target genes of TFs. Transcription factors whose target genes show a significant overlap with the DEGs were reported. The value of <0.05 for all four indexes (P-value, q-value, E-value, and FDR) in the TFacts Database was considered statistically significant.

Finding prognostic genes for PDAC

Kaplan Meier Plotter (http://kmplot.com/analysis/), an open-access database that provides log-rank P-value and Hazard ratio with 95% confidence intervals for Kaplan–Meier analysis of the correlation between mRNA expression level and patient overall survival (Nagy et al., 2018), was employed to predict the prognostic importance of the DEGs detected in this study. The TCGA Pan-Cancer pancreatic ductal adenocarcinoma cohort was used for the analysis. Auto-select best cutoff value was used to split the patients in survival analysis. P ≤ 0.05 and FDR ≤ 0.05 were used as the cutoff for significance. Then, four independent external validation datasets (GSE21501, GSE50827, GSE57495, and GSE71729) including clinical and transcriptomic data from pancreatic ductal adenocarcinoma patients were analyzed using the PROGgenev2 Prognostic Database (Goswami & Nakshatri, 2014) to evaluate the prognostic relevance of the identified candidate prognostic genes by Km-Plotter. The patient cohorts were divided into two equal groups based on median expression for each gene. P ≤ 0.05 was accepted as statistically significant.

The correlation between mRNA expression of the identified prognostic genes and pathological stages of the disease was evaluated in TCGA PAAD data using the GEPIA Database. P < 0.05 was accepted as statistically significant.

Protein-protein interaction analysis

Search Tool for the Retrieval of Interacting Genes (STRING, version 11.0, https://string-db.org/) is an online database designated to evaluate physical and functional associations of proteins (Szklarczyk et al., 2014). STRING app in Cytoscape Software (version 3.6.1) was used to detect the interactions among the DEGs with a confidence score cut-off >0.9. Then, the protein-protein interaction (PPI) network was constructed using Cytoscape Software (Shannon et al., 2003). The Molecular Complex Detection plug-in (MCODE) was used to identify clusters in the PPI network with the parameters false degree cutoff 2 and K-Core 2 (Bader & Hogue, 2003). Finally, the Gene Ontology and KEGG pathway enrichment analyzes of the DEGs in the clusters with MCODE score >5 were performed with the STRING Enrichment app in Cytoscape Software by retrieving functional enrichment for selected clusters only.

Gene ontology and pathway enrichment analysis of the DEGs

Gene Ontology Enrichment Analysis was performed using the Functional Enrichment Analysis Tool (FunRich) (Pathan et al., 2015) to identify the biological processes, molecular functions, and cellular components that are shared by the up-or down-regulated genes separately. The Gene Ontology Database (The Gene Ontology Consortium, 2017; Ashburner et al., 2000; The Gene Ontology Consortium, 2018) was selected as a background for the analysis, and a P-value of <0.05 was used as the cut-off for significance.

The list of differentially up-and down-regulated genes in PDAC was explored for functionally enriched pathways using a WEB-based GEne SeT AnaLysis Toolkit (WebGestalt 2017) (Wang et al., 2017) with a cut-off criterion of FDR < 0.05. Enrichment categories were selected as Kyoto Encyclopedia for Genes and Genomes (KEGG) (Kanehisa & Goto, 2000), Reactome Pathway Knowledgebase (Fabregat et al., 2018), PANTHER (Thomas et al., 2003), and Wikipathways (Kelder et al., 2012). The category size was calculated based on the number of overlapping genes between the annotated genes in the category and the ranked gene list for the GSEA method. Categories with sizes smaller than 15 and greater than 1,000 were removed during the analysis.

Microarray datasets

The workflow and summary of the results of the presented meta-analysis is shown in Fig. 1. A systematic search of the studies was carried out up to January 2019. The search with the keyword Pancreatic Ductal Adenocarcinoma’ resulted in 18 microarray gene expression datasets which had been submitted only to the ArrayExpress Archive of Functional Genomics Data. Among these, none of the datasets was found to contain adjacent normal tissue as a control group and therefore these datasets were not included in the meta-analysis. Gene Expression Omnibus Datasets were searched with the term “pancreatic ductal adenocarcinoma” and the search results were filtered by selecting organism as homo sapiens, study type as expression profiling by array and attribute name as tissue. As a result, 54 studies were found in the GEO Database. These studies were carefully evaluated and eight studies (GSE62452, GSE15471, GSE62165, GSE56560, GSE60646, GSE55643, GSE18670, and GSE58561) which contain data from 144 PDAC tissues and matched adjacent non-tumor tissue samples were found to meet the inclusion criteria. After exclusion of non-confirmed sample pairs and samples with low quality, 105 PDAC and matched adjacent non-tumor tissue pairs from four datasets (GSE62452, GSE15471, GSE62165, and GSE56560) were decided to be eligible for the meta-analysis (Table 1).

Differentially expressed genes in PDAC

The meta-analysis identified 344 (342 protein-coding genes) differentially over-expressed and 168 (157 protein-coding genes) repressed genes in PDAC compared to the matched adjacent non-tumor samples (Table S1). A list of the top 20 differentially expressed genes (FDR < 0.05) with at least a 2-fold differential expression between groups is shown in Table 2. Along with the genes reported to be associated with pancreatic cancer in previous studies, the results of the meta-analysis revealed candidate over-expressed genes for PDAC, many of which have previously received little or no attention, such as KYNU (kynureninase), ITGBL1 (integrin beta-like 1), and ADGRF1 (adhesion G protein-coupled receptor).

The differentially expressed genes identified in PDAC was validated in TCGA combined GTEx data containing data from PDAC tissues and pancreatic tissues from healthy individuals. This confirmation was also aimed at determining whether the DEGs identified could also distinguish PDAC tissues from healthy pancreatic tissues. Gene expression analysis of TCGA combined GTEx data validated the up-regulation of 336 and the down-regulation of 116 protein-coding genes in PDAC (P ≤ 0.001, Fold Change ≥2). Additionally, the elevated expression of AK4 was validated in the Logsdon Pancreas Dataset in Oncomine Database due to the absence of the probes for this gene in TCGA combined GTEx data (P = 1.16E−5). The identified DEGs and the results of the validation analysis are shown in Table S1.

Taken together, the expression of a significant portion of the identified DEGs was consistent with the gene expression profile from TCGA, implying that the identified genes may have the potential to discriminate PDAC from both adjacent and healthy pancreatic tissues.

Activated and inhibited transcription factors in PDAC

To predict which transcription factors (TF) are activated or inhibited to regulate the transcription of the identified DEGs, TFs whose target genes show a significant overlap with the DEGs were identified using TFactS Database. Table S2 shows the results of the TF enrichment analysis. A total of five transcription factors were found to reach statistical significance (P-value, q-value, E-value, and FDR < 0.05). Among these, Transcription Factor 7 (TCF7), Catenin Beta 1 (CTNNB1), Smad Family Member 3 (SMAD3) and Jun Proto-Oncogene (JUN) were found to be significantly activated in PDAC while Smad Family Member 7 (SMAD7) was detected as significantly inhibited.

Prognostic values of the identified DEGs in PDAC

The prognostic values of the identified DEGs genes were evaluated by analyzing the TCGA Pan-Cancer pancreatic ductal adenocarcinoma data and survival outcomes of patients included in the Km-Plotter. Sixty-one up-regulated and one down-regulated genes (P ≤ 0.05 and FDR ≤ 0.05) were found to be significantly associated with the overall survival rate of patients with PDAC. TCGA combined GTEx data validated the up- and down regulation of the identified candidate prognostic genes in PDAC. Subsequently, four external datasets (GSE21501, GSE50827, GSE57495, and GSE71729) were analyzed to validate the prognostic significance of the identified 62 candidate genes for PDAC. The candidate prognostic gene list identified by Km-Plotter and the results of the validation analysis are shown in Table S3.

Genes whose prognostic potentials were confirmed in one of the four validation sets were included in the “prognostic gene list” and their predicted prognostic values were ranked based on (i) validation status (validated in 1-2-3-4 dataset(s) =validation status 1-2-3-4), (ii) the average of p-values, and (iii) the average of hazard ratios. None of the genes could reach validation status 4. However, two genes, LAMC2 and ITGB6, were found to be significantly correlated with worse overall survival in three of the four validation datasets. SERPINB5, KRT7 TGM2, IGF2BP3, INPP4B, IL1RN, DCBLD2, and MPZL2 were validated in two datasets. Three genes, TMPRSS4, SEMA3C, and MGLL, were associated with increased overall survival in one of the validation datasets, therefore excluded from the prognostic gene list. Gene expression data of CDK1, PKM, AK4, AREG were not available in some validation datasets. While AK4 and PKM, which were found to be correlated with overall survival in one of the four external datasets, were included in the prognostic gene list, CDK1 and MGLL were excluded from the list due to insufficient gene expression data.

Consequently, a total of 28 up-regulated genes were individually validated in at least one of the four external validation datasets. The ranked prognostic gene list for PDAC and the results of the Kaplan–Meier survival analysis in five datasets are shown in Table 3. Kaplan–Meier survival plots for the identified prognostic genes are shown in Fig. 2.

Further analysis of TCGA pancreatic ductal adenocarcinoma data in GEPIA showed that, twenty-one of the identified prognostic genes (ITGB6, LAMC2, KRT7, SERPINB5, IGF2BP3, IL1RN, MPZL2, SFTA2, MET, LAMA3, ARNTL2, SLC2A1, LAMB3, COL17A1, EPSTI1, IL1RAP, AK4, ANXA2, S100A16, KRT19, and GPRC5A) were also correlated with the pathological stages of the disease, underlying their prognostic value for PDAC (Fig. 3).

Protein-protein interaction (PPI) network

The PPI analysis was performed to evaluate the physical and functional associations of proteins encoded by the identified differentially expressed genes in PDAC. The PPI network was constructed by Cytoscape based on the PPI correlations from the STRING database. PPIs among the DEGs with a confidence score cut-off >0.9 were selected to construct the PPI network (Fig. 4). Then the degree of each node in the network was calculated by using the NetworkAnalyzer tool of Cytoscape to identify hub proteins in the PPI network. The degree of a node is the number of edges connected to the node, and it has been stated that nodes with higher degrees -which correspond to hub proteins in the PPI network-, play an essential role in the organization of the PPI network, therefore might be more crucial and relevant than non-hub genes (Vallabhajosyula et al., 2009). In this study, nodes with degrees >15 are considered to indicate “hub proteins” and are presented in Table 4. Next, the PPI network was analyzed by MCODE to identify clusters in the network. Each clustered protein group was then analyzed to predict its biological function in PDAC.

Among the 27 hub proteins, ITGA2 and COL17A1 were found to be associated with unfavorable prognosis for PDAC based on the TCGA Pan-Cancer pancreatic ductal adenocarcinoma dataset in Km-Plotter (P ≤ 0.05 and FDR ≤ 0.05). However, this finding could not be validated in the four external datasets, except the significant correlation between the up-regulation of COL17A1 and decreased overall survival found in GSE57495. Moreover, Fibronectin 1 (FN1), tissue inhibitor of metalloproteinases 1 (TIMP1), and fibrillin 1 (MSLN) constituted the top three proteins with degrees exceeding 20. All these three proteins were members of CLUSTER1, which includes proteins whose functions mostly associated with extracellular matrix signaling pathways and structural organization, underlying the importance of extracellular dynamics in the pathogenesis of PDAC. The most significantly up-regulated gene in this study, POSTN, was found to be located close to the collagen sub-cluster included in CLUSTER1, emphasizing their close interactions. CLUSTER1 also included up-regulated genes with unknown importance in pancreatic cancer, such as MATN3, SERPINA1, and IGFBP5.

CLUSTER3 covered known prognostic biomarkers such as TOP2A (Tsiambas et al., 2007), CDK1 (Piao et al., 2019), and MKI67 (Striefler et al., 2016). This cluster was surrounded by genes (CENPK, ANLN, ECT2, and FAM83D) related to the cell cycle, but it appears that they were not included in the cluster due to the stringent confidence score cut-off criteria of the analysis. Amongst, the function of FAM83D in PDAC is not yet known. However, there are previous reports that describe this protein as a probable proto-oncogene that regulates cell proliferation, growth, migration, and epithelial to mesenchymal transition (Wang et al., 2015; Wang et al., 2013), suggesting that it may have similar functions in PDAC.

Furthermore, CLUSTER2 and 4, were mapped to the chemokine signaling pathway and type I interferon signaling pathway, respectively, and located close to each other in the PPI network. Although the importance of chemokine signaling in PDAC has been previously reported (Geismann et al., 2019), the exact role of cellular innate antiviral response in the pathogenesis of PDAC remains elusive. In this study, OAS1, OAS2, and RSAD2, which play critical roles in cellular innate antiviral responses induced by type I and type II interferon, were found to be up-regulated in PDAC. Additionally, Kaplan–Meier survival analysis of the TCGA data using Km-Plotter indicated an association of high OAS1 expression with the low survival rate of PDAC patients, however, could not be validated in other four external validation datasets. The roles of these proteins in PDAC are still unknown, and further research is needed for the evaluation of their potential as diagnostic or prognostic biomarkers for PDAC.

While CLUSTER6 shows the known interleukin-laminin-EGF crosstalk in the focal adhesion pathway in PDAC, CLUSTER9 indicated a novel down-regulated gene in PDAC, P2RX1, a ligand-gated ion channel with relatively high calcium permeability, that may be linked to apoptosis by increasing the intracellular concentration of calcium in the presence of ATP (Uniprot Database, by similarity). CLUSTER8 included cell surface antigens such as CD109, CD66e, and CD90. However, as shown in the PPI network, the interaction between these proteins and the proteins that have roles in the regulation of glycolysis such as PKM, PFKP, and ENO2 was intriguing and may indicate stemness associated interactions in PDAC. The other clusters were significantly mapped to those specific KEGG pathways: “Mucin-Type O-glycan Biosynthesis” (CLUSTER5) and “Estrogen Signaling Pathway” (CLUSTER7), indicating once more the high expression of keratins and mucin-type O-glycans in pancreatic tumors. The proteins assigned to the pathways are shown in Table S4.

Gene ontology analysis of the DEGs

Gene ontology analysis of the DEGs was performed to identify enriched biological processes, molecular functions and cellular locations associated with differential gene expression in PDAC. As shown in Fig. 5A, the analysis identified nine significantly overrepresented individual categories of GO Molecular Function for the up-regulated genes including “extracellular matrix structural constituents” (P < 0.001), “collagen binding” (P < 0.001) and “integrin-binding” (P < 0.001). “Extracellular matrix organization” and “cell adhesion” were the most significantly enriched biological processes (P < 0.001, Fig. 5B), while most of the proteins encoded by the up-regulated genes were found to be located in the “collagen-containing extracellular matrix” (P < 0.001, Fig. 5C).

For the down-regulated genes, “serine-type endopeptidase activity” was the only GO molecular function category that reached statistical significance (P = 0.008, Fig. 5D). Three GO Biological Processes were found to be significantly enriched for the down-regulated genes: “proteolysis” (P = 0.001), “cellular zinc ion homeostasis” (P = 0.002), and “cellular response to copper ion” (P = 0.042), (Fig. 5E). Moreover, most of the over-represented GO cellular locations were associated with extracellular space for the down-regulated genes (Fig. 5F).

Pathway enrichment analysis of the DEGs

Significantly dysregulated pathways in PDAC were identified by using four different pathway databases in WebGestalt (Table S5). The obtained results from the pathway enrichment analysis of the identified DEGs underlined the importance of the crosstalk between tumor cells and extracellular matrix, including integrin and PI3K-Akt-mTOR signaling pathways, in PDAC pathogenesis. Pathways related to metabolism and pancreatic secretion were identified as negative related categories by KEGG and Reactome Pathway Knowledgebase.