Skip to main content

Novel insight into pancreatic adenocarcinoma pathogenesis using liquid association analysis



Pancreatic ductal adenocarcinoma (PDAC) is a lethal malignancy associated with a poor prognosis. High-throughput disease-related-gene expression data provide valuable information on gene interaction, which consequently lead to deeper insight about pathogenesis. The co-expression analysis is a common approach that is used to investigate gene interaction. However, such an approach solely is inadequate to reveal the complexity of the gene interaction. The three-way interaction model is known as a novel approach applied to decode the complex relationship between genes.


In the current study, the liquid association method was used to capture the statistically significant triplets involved in the PDAC pathogenesis. Subsequently, gene set enrichment and gene regulatory network analyses were performed to trace the biological relevance of the statistically significant triplets.


The results of the current study suggest that “response to estradiol” and “Regulation of T-cell proliferation” are two critical biological processes that may be associated with the PDAC pathogenesis. Additionally, we introduced six switch genes, namely Lamc2, Klk1, Nqo1, Aox1, Tspan1, and Cxcl12, which might be involved in PDAC triggering.


In the current study, for the first time, the critical genes and pathways involved in the PDAC pathogenesis were investigated using the three-way interaction approach. As a result, two critical biological processes, as well as six potential biomarkers, were suggested that might be involved in the PDAC triggering. Surprisingly, strong evidence for the biological relevance of our results can be found in the literature.

Peer Review reports


It is well established that adenocarcinoma (PDAC) is responsible for more than 90% of all pancreatic cancer cases [1]. The median survival duration of PDAC is less than 6 months, and only 6% of patients survive by passing 5 years from their diagnosis. Unfortunately, due to late diagnosis and resistance to the currently available chemotherapy drugs, the proposed treatments have not improved significantly [2, 3]. Therefore, obtaining a better understanding on molecular mechanisms of the disease’s progression and discovery of novel therapeutic targets is an extremity needed to improve the treatments of PDAC. Although the pathogenesis of PDAC is extensively studied, it remains unknown yet.

Gene interactions provide critical information on disease pathogenesis. In this regard, high-throughput gene expression data is widely adopted for elucidating gene interactions [4, 5]. Accordingly, such data have been extensively investigated to explain the pathogenesis of PDAC [6]. Giulietti et al. have suggested several central genes involved in PDAC using the weighted gene co-expression network analysis on the PDAC microarray-based gene expression data. Additionally, they proposed several biological processes, including lipid metabolism and transmembrane transport, which are critical in the PDAC pathogenesis [7]. On the other hand, Zhou et al. have introduced several hub genes associated with both the progression and prognosis of PDAC by constructing a co-expression network. Besides, they have suggested that the cell cycle plays a vital role in PDAC [8].

To the best of our knowledge, a large amount studies to detect essential genes and pathways involved in PDAC were performed using the two-way interaction approach [7,8,9]. Such an approach can capture co-expressed genes, may possibly encode some functionally-related proteins. However, sometimes the two-way interaction approach is inefficient to trace functionally-related genes. A reason for this fact is the molecular state of a cell is susceptible to intra- and inter-cellular alteration. Therefore, the co-expression relationship of two functionally-related genes may change based on the cellular condition. In other words, the co-expression relationships may have a dynamic nature depending on the cellular state. Hence, the two-way interaction approach is too simplistic to describe the complex relationships among genes [10,11,12].

The three-way interaction approach is a more efficient strategy to trace functionally-related genes. Such an approach is inferred from the cross-shaped co-expression pattern, which in the expression levels of two genes is directly correlated under a certain condition, while these are inversely correlated under another condition. More precisely, the three-way interaction approach captures the dynamic nature of the co-expression pattern using introducing the third gene, called a switch gene. In other words, depending on the expression level (or genotype) of the switch gene, the expression levels of a gene pair are either directly or inversely correlated [13,14,15]. Therefore, determining the three-way interaction would be considered a pivotal step to decode the complex biological systems.

Moreover, the three-way interaction approach can be a competent compromise between the reality of the complexity of biological interactions and the discoverability of interactions. Because the two-way interaction approach is too simplistic for explaining complex molecular interactions; on the other hand, any four-way or higher-level models are less utilizable because, in the transcriptomics level, these models require a huge computational demand (as the number of combinations increases exponentially with the number of interaction levels) [13, 16].

A well-known representative of the three-way interaction is the relationship among thyroid hormone expression, growth hormone, and thyroid hormone receptor (TR-RXR). TR-RXR plays an essential role as a suppressor or as an activator depending on the absence or the presence (amount) of thyroid hormone, respectively. Therefore, it can be stated that the expression levels of both growth hormone and TR-RXR are directly correlated in the presence of thyroid hormone, but they are inversely correlated in the absence or with the low amount of thyroid hormone, as a switch gene [17].

To the best of our knowledge, gene expression relationships in PDAC have not been studied using the three-way interaction model so far. Therefore, the present study aimed to investigate three-way interaction in the PDAC-related microarray gene expression data, in order to find the potential biomarkers and biological processes involved in their pathogenesis (Fig. 1).

Fig. 1
figure 1

Flowchart for identification biologically relevant three-way interaction

Materials and methods

Choosing the dataset and reducing the dimension

In order to gather PDAC associated gene expression profiling datasets, the publicly available transcriptome data repositories include ArrayExpress [18], NCBI gene Expression Omnibus (GEO) [19] and The Cancer Genome Atlas (TCGA) [20], were comprehensively screened. Among numerous transcriptome datasets such as GDS4336 [21] and GSE15471 [22], carried out with PDAC and adjacent non-tumor tissue samples, we preferred to employ GDS4336 considering its large sample size and sampling characteristic.

The employed databases included pancreatic ductal adenocarcinoma tumour (45 samples) and adjacent non-tumour tissue gene expression (45 samples). There were 35,556 probe sets on Affymatrix Human Gene 1.0 ST Array [23]. Thereafter, we used the genefilter package [24] in Bioconductor to remove the duplicate probes. Finally, the 18,232 genes remained for further statistical analyses.

Selecting the candidate switching genes

According concept the three-way interaction, a {X1, X2} gene pair have a contradictory co-expression behavior depending on the expression level of corresponding switch gene. Therefore, it is expected a disease related-switch gene exhibit the following traits (1) differentially expressed in a disease samples compare to healthy ones and consequently, (2) its gene expression distribution be bimodal.

Thereafter, fold change and standard errors were estimated by fitting a linear model (using the lmFit function in Limma package [25] for each gene in the groups. Genes with empirical Bayes t-test and p values of 0.05 were firstly selected and then corrected by calculating the false discovery rate p < 0.05. Furthermore, the bimodality of each one of the gene’s expression was calculated by diptest function in diptest package [26] in the selected genes by considering p < 0.05.

In this way, the two above-mentioned constraints were applied to define the candidate switching genes.

Determining the liquid association triplets

The three-way interactions between the candidate switching genes as well as all the possible combinations were calculated using the fastMLA function in the fastLiquidAssociation package [27]. Accordingly, this package works using a modified liquid association algorithm for the determination of changes in the co-expression relationships of the other two genes X1 and X2, in terms of the expression level of the third gene (X3). Notably, before running the liquid association algorithm, performing two pre-processing stages for each variable on the data are essential, which are as follows: the first is transforming into normal distribution [28], and the second one is standardizing the mean to 0 and standard deviation to 1 [29]. In the present study, the first and second steps were performed by an in-house script and CTT package [30], respectively. Next, the Bonferroni correction approach was used to estimate the false discovery rate (FDR) [31]. As well, the statistically significant triplets were considered the three-way interactions with FDR < 0.05.

Constructing the gene regulatory network

One of the fundamental properties of the gene regulatory network (GRN) is providing a systematic proficiency regarding the complicated molecular mechanisms, which can be used to check gene expression under diverse biological conditions [32, 33]. The geWorkbench_2.6.0 software was then applied to reconstruct GRN among all the studied genes in the dataset as hub markers. Moreover, the Algorithm for the Reconstruction of Accurate Cellular Networks (ARACNE) was used to construct the GRN [32]. Accordingly, this algorithm is able to make mutual information (MI) as a measure to discover the direct regulatory interactions between each transcriptional regulator and its potential target(s).

Gene set enrichment analysis

Gene set enrichment analysis (GSEA), which is a popular statistical method, can establish the validity of the shared association among a set of genes using a predefined Gene Ontology [34]. Afterward, GSEA was done in this study using the Biological Process of Gene ontology database [35]. These analyses were applied using the ClueGO plugin [35] within the Cytoscape v.3.3.0 environment [36]. Additionally, the two-sided hypergeometric test along with the Bonferroni step-down correction approach and a kappa score level of 0.4 were regarded as the statistically significant analyses.

In silico validation in other related and unrelated gene expression datasets

The SurvExpress [37] web-based biomarker validation tool was used to assess the suggested switch genes as potential biomarkers in various cancer types. Validation analyses were carried out through related as well as unrelated gene expression datasets obtained from TCGA [20], ICGC [38] and NCBI-GEO [19].

The suggested switch genes; namely, Lamc2, Klk1, Nqo1, Aox1, Tspan1, and Cxcl12, were surveyed through pancreatic ductal adenocarcinoma (ICGC, n = 189), cervical squamous cell carcinoma (TCGA, n = 191), lymphoma (NCBI-GEO, GSE10846, n = 420), glioblastoma (TCGA, n = 148), hepatocellular carcinoma (NCBI-GEO, GSE10143, n = 162) and stomach adenocarcinoma (TCGA, n = 352).

For each utilized dataset, samples according to their prognostic index were classified into low and high-risk groups. Subsequently, in each dataset, the prognostic performance of suggested switch genes was determined using Kaplan–Meier plots, log-rank test p values, hazard ratios (HR), and their confidence intervals (CI). Additionally, the correlation of the survival analysis with gene expression levels was presented by Heatmaps, which in samples were sorted by their prognostic index, and genes were clustered based on Euclidean distance. Moreover, differences in gene expression levels between high and low-risk groups were determined using the t-test; furthermore, box plots were created to illustrate such differences.

Ethical approval statements

Since the current investigation that was performed using a public GEO dataset (accession number: GDS4336) can be considered a computational survey, obtaining both ethics approval and consent statements was not required in this research.


Data pre-processing and candidate switch gene selection

After pre-processing and normalization of the obtained data using the robust multi-array analysis method, the final dataset consisted of 90 samples, including 45 pancreatic ductal adenocarcinoma tumour samples and 45 adjacent non-tumour samples. Accordingly, each Array included 28,896 probes. After removing duplicate probes, this number dropped to 18,232 probes.

After restricting the dataset to include both the differentially expressed genes and bimodality expressed genes, the 53 genes were remained. Such genes considered as the candidate switch genes to be analysed using fastLA. The list of the candidate switch genes is reported in the Additional file 1: Table S1.

Tracing statistically significant cross-shaped triplets

The liquid association analysis was performed for each combination consisting of a candidate switch gene (X3) and each possible pair of genes in the dataset {X1, X2}. Subsequently, we selected the top 300,000 triplets with the highest significance levels in terms of their p-values, which were defined as the outputs of this analysis. Figure 2 shows the changes in FDR versus − log (p-value) for the first 300,000 triplets. In addition, a set of significant cross-shaped triplets comprising of 220 triplet combinations were chosen for conducting more analyses, in terms of FDR < 0.05. The list of the statistically significant triplets is reported in the Additional file 1: Table S2.

Fig. 2
figure 2

FDR versus − log (p-value). The changes in FDR (Bonferroni-corrected p-value) versus − log (p-value) for the first 10,000 results of fastLA [21]. As shown FDR = 0.05 corresponds to − log (p-value) = 6.78

Identifying the biological process involved in PDAC

GSEA was applied to detect the biologically-relevant triplets as well as finding the BP involved in PDAC. Correspondingly, GSEA was performed with a p value < 0.05 and an FDR < 0.1 for all the genes involved in the statistically significant triplets. Since the terms in lower levels of gene ontology are more general, so the enriched terms in levels lower than five were not reported. Figure 3 shows all the significantly enriched terms, including “Response to estradiol”, “Embryonic limb morphogenesis”, “Positive regulation of cellular amide metabolic process”, “Establishment of the protein localization to plasma membrane”, “Type I interferon signalling pathway”, “Positive regulation of the developmental growth”, and “Regulation of T cell proliferation”.

Fig. 3
figure 3

Gene set enrichment analysis (GSEA). Enriched terms based on biological process for all genes involved in statistically significant triplets. The biological relevance of two statistically significant triplets was confirmed by GSEA

According to the definition of the three-way interaction, gene expressions of both X1 and X2 were found to be directly correlated under a specific condition. So, it is expected that both X1 and X2 be in a common pathway or biological process. Therefore, in the current study, we detected the statistically significant triplets in the enriched terms to identify the biologically relevant triplets among them. In this regard, we identified two biologically relevant triplets in terms of GSEA, including Klk, and {Cst3, Gh1}, which are involved in the "Response to estradiol" process, as well as Lamc2 and {Ccr7, Dhps}, which are involved in the "Regulation of T cell proliferation" process.

Additionally, Fig. 4 presents the scatter plot of these triplets. This plot also reveals a remarkable change in the correlation between X1 and X2 resulting from the changes in X3 expression level.

Fig. 4
figure 4

Scatter plot of two biologically relevant triplets. Based on the fastLA algorithm, the samples are divided into three-bin according to the expression of the X3 gene. In each case, a considerable change in the correlation of X1and X2 occurs as a result of the change in X3

Detecting the cross-shaped triplets in GRN

In order to further analyse the functional relevance of the statistically significant triplets, a GRN was constructed using ARACNE. The detailed information of such network is reported in the Additional file 1: Table S3. According to the concept of the three-way interaction, the expression level of the switch gene (X3) affects the correlation direction of the {X1, X2} gene pair. In this regard, it is expected to observe a regulatory relationship between these genes.

By detecting the statistically significant triplets in GRN, we found that seven triplets may be biologically relevant to each other. The position of such triplets is presented as a sub-network of GRN in Fig. 5. The information on the liquid association analysis of such triplets is reported in Table 1.

Fig. 5
figure 5

Regulatory relationships within triplets. The regulatory relationships of significant triplets obtained from liquid association analysis were traced in the GRN. The Mutual information (MI) value identifies the magnitude of each relationship

Table 1 The liquid association analysis information of the seven triplets that the regulatory relationships of them were traced in gene regulatory network

Altogether, the biological relevance of nine triplets was approved by either GSEA or GRN methods.

The specificity of the switch genes to PDAC

In order to in-silico validation of suggested switch genes as PDAC's potential biomarker, several analyses were carried over various human cancers.

First of all, the prognostic performance of suggested switch genes was assessed in a related gene expression dataset (i.e., PDAC). As indicated in Fig. 6a, the suggested switch genes show significant prognostic performance in such a PDAC-related dataset (HR = 1.93, p = 7.8 × 10−5). Then, the prognostic performance of suggested switch genes was assessed in five gene expression datasets of unrelated but highly prevalent cancers. The survival analyses pointed out insignificant performance in such cancers, cervical squamous cell carcinoma (p = 0.10), lymphoma (p = 0.22), glioblastoma (p = 0.06), hepatocellular carcinoma (p = 0.87) and stomach adenocarcinoma (p = 0.10). (Additional file 1: Fig. S4). As a consequence, our suggested switch genes could be considered as prognostic for PDAC.

Fig. 6
figure 6

The prognostic power of suggested switch genes through related and unrelated datasets. (A) Pancreatic ductal adenocarcinoma as a related dataset, (B) cervical squamous cell carcinoma and lymphoma as two exemplary unrelated datasets

The prognostic power of suggested switch genes in cervical squamous cell carcinoma and lymphoma as two exemplary datasets shown in Fig. 6b. The survival analysis results of other PDAC-unrelated datasets were presented in Additional file 1: Fig. S4.


The improvements are seen in our descriptive comprehension of pancreatic ductal adenocarcinoma; however, some considerable cancer biology features like the entire involved genes, still remain poorly understood. It is obvious that obtaining a better understanding on PDAC biology may lead to provide more effective treatments in this regard [39]. Herein, we used a novel bio-computational approach to identify new candidate genes involved in this aggressive cancer process. To the best of our knowledge, the present study introduced a novel bio-computational approach for the first time, named as the three-way interaction model, to investigate new candidate genes as well as the possible mechanisms of the gene switching in PDAC.

Gene set enrichment analysis

The enriched biological processes (BP) are “response to estradiol”, “embryonic limb morphogenesis pathway”, “positive regulation of cellular amide metabolic process”, “establishment of protein localization to plasma membrane”, “type I interferon signalling pathway”, “positive regulation of the developmental growth”, and “regulation of T-cell proliferation”. As well, all the seven common enriched terms are biologically relevant to PDAC. For more information, see below.

In this regard, the first enriched biological process is the "response to estradiol". Pancreatic cancer cell proliferation is highly estrogen-sensitive in vitro, and estrogen receptor alpha and beta are frequently expressed in cancerous cells [40]. ERbeta/ERalpha ratio may even affect the estrogen-mediated growth stimulation and then reduce cytotoxicity at physiological concentrations that may possess some clinical implications for pancreatic cancer therapy [40]. A recent study demonstrated the critical roles of estrogen receptors and the efficacies of anti-estrogen therapy in both pancreatic cancer cell’s proliferation as well as cancer progression’s inhibition [41]. In some previous studies, it was suggested that T-box factors, especially Tbx2, as transcription factors, play important roles in controlling cell cycle progression, embryonic development, and cancer genesis. Tbx2 is implicated in various developmental processes like the morphogenesis of various organs and tissues such as limbs, kidneys, lungs, heart. Moreover, it was shown that it is overexpressed in various types of cancers such as pancreatic, liver, bladder, and breast cancers and can overwhelm senescence, which is a cellular process serving as a cancer development inhibitor [42].

Another enriched biological procedure is "positive regulation of cellular amide metabolic process". Nitrogen-based nucleotides synthesis and nonessential amino acids (NEAAs) are known as the essential metabolic stages in tumour cell growth. Accordingly, by donating its amide group and converting it into glutamic acid, it can provide an essential nitrogen donor in 3 autonomous enzymatic reactions in charge of purine synthesis. In addition, it is responsible for pyrimidine synthesis in two genes [43, 44]. Hence, the amide metabolism pathway and the change of glutamine from an NEAA in normal cells into an essential amino acid within cancerous cells play the roles of the main substrates feeding the tumorigenic cells’ anabolic growth procedures.

Interestingly, a previous study [44] indicated that in PDAC cells, both amide and glutamine metabolism pathways support tumorigenic growth through the noncanonical metabolic pathway, which is oriented by the general biochemical position adopted by several types of cancer cells [45].

"Establishment of the protein localization to plasma membrane" is another enriched biological process. Such a biological process is considered an underlying mechanism for increasing the aggressive performance of pancreatic cancer. Accordingly, it was identified that some localized proteins, like the protein placenta-specific 8 (PLAC8, Onzin), are related to the pancreatic tumour progression. Furthermore, these can be expressed invasively and ectopically within the human PDAC and advanced preneoplastic lesions. Nevertheless, the precise molecular function of the localized proteins is still unclear, and some pieces of evidence highly proposed its role, which greatly is based on both physiologic and cellular contexts. It was demonstrated that such proteins are localized to the plasma membrane in pancreatic cancer cells, within which they interact with particular membranous structures in a spatially and temporally stable mode [46].

As mentioned earlier, "Type I interferon signalling pathway" is another enriched biological process. It was shown that Type I interferons (e.g., IFNα/β) possess various antitumor activities. Some clinical studies evaluated the effects of adjuvant IFNα treatment on the obtained equivocal findings regarding pancreatic cancer [47]. Importantly, unlike permissive PDAC cell lines, the resistant PDA cell lines can mainly respond and secrete type I interferon (IFN). Correspondingly, this indicates the contribution of the intact type I IFN responses to their resistance phenotype [48]. However, it is thought that cancer cells are mostly defective in type 1 IFN responses and also in production [49, 50]. As well, since IFN responses are normally pro-apoptotic, anti-proliferative, and anti-angiogenic [51], they are undesirable circumstances for tumour formation. Nevertheless, some cancer cells, such as some mesotheliomas [52], melanomas [53], lymphomas [54], bladder cancers [55], renal cancers [56], and possibly other types of cancers [47, 57], create or respond to type I IFN [58]. Therefore, despite the fact that human pancreatic cancer cell lines variably respond to IFNα and β, it seems that both type I interferon signalling pathway and the expression level of the type I IFN receptor play an important role in pancreatic cancer progression.

The "positive regulation of developmental growth" is also another enriched biological process. Scientifically, the idea of cancer cells sharing some features of their embryonic predecessors is historical [59]. Tumour cells, or at least a subset of them, similar to embryonic development, can preserve indefinite growth and cellular plasticity. It was indicated that the developmental paths could direct the start of PDA precursors from their cellular ancestors; however, embryonic signalling paths such as TGFβ, Hedgehog, and Wnt-β-catenin alone are insufficient for the PDA initiation [60]. Recent data displayed the role of Sox9, which is expressed throughout the development, in the PDA initiation [61, 62]. Therefore, the positive regulation of the developmental growth pathway is associated with pancreatic cancer progression.

The last enriched biological process is "regulation of T-cell proliferation". It was shown that immune disorders are one of the most significant hallmarks of cancer, and pancreatic cancer is characterized by abnormal immune cell infiltration. Accordingly, this is related to the immunosuppressive circumstances facilitating the tumour cells’ escape from immune cells [63, 64]. In pancreatic cancer, tumour-infiltrating regulatory T-cells (Tregs) are considered independent prognostic factors, which increase the secretion of different immunosuppressive cytokines, prevent cytotoxic lymphocyte function, and play a role in the immune escape [65,66,67].

Biologically relevant triplets

In this study, among the statistically significant triplets, two triplets whose X1 and X2 genes are involved in the same biological procedure are known as the biologically-relevant triplets.

The relationship between the involved genes in triplet Klk1 and {Gh1, Cst3}

One of the biologically relevant triplets contains Kallikrein-1 (Klk1), as the switch gene, and a gene pair {growth hormone 1 (Gh1), cystatin C (Cst3)} (Fig. 4). According to Fig. 3, genes in this triplet are involved in the cancer’s progression via "response to estradiol" biological process. It was also demonstrated that estrogens, like estradiol, are a major factor in the cancer’s progression, like breast cancer [68]. Additionally, these are tightly linked to the growth hormone (Gh/Gh1) [68].

Gh1 is a pituitary peptide hormone with some stimulating pleiotropic effects. The secretion of Gh1 occurs through two systems, including autocrine and paracrine. In human beings, Gh1 signal transduction not only is primarily mediated via binding to the GH receptor (GHR), but it could also be done by the prolactin receptor (PRLR) [69]. In several types of cancers like triple-negative breast cancers (TNBC), a correlation exists between a higher expression of the membrane-bound G protein‑coupled estrogen receptor (GPER) and a worse consequence caused by that. On the other hand, a strong connection potentially exists between GPER expression and growth hormone receptor (GHR). Previously, Girgert et al. in their study have reported that the inhibition of GH receptor reduces the expression of G protein‑coupled estrogen receptor (GPER) and inhibits growth stimulation of breast cancer indeed via the inhibition of estrogen (estradiol) signal transduction [68].

The Cst3 gene belongs to the cystatin superfamily that encompasses some proteins containing multiple cystatin-like sequences, which are considered endogenous inhibitors of lysosomal cysteine proteinases. Moreover, it was found that this inhibitor contributes to various pathological and normal processes and seems to be implicated in the malignant progression of different human tumours. The increased Cst3 levels have also been reported in patients with malignant diseases [70, 71]; however, the exact role of Cst3 in malignant diseases is still debatable [70]. It has been previously reported that hormones, including Gh1 [72, 73] and steroidal agents (estradiol), can increase Cst3 levels in the circulation [70]. Therefore, the increased expression levels of both Gh1 and Cst3 (as a positive co-expression) and cross talk between the factors, especially via responding to the estradiol pathway, may play a key role in the progression of cancer [73]. Therefore, as displayed in the above-mentioned pieces of evidence, a direct co-expression relationship exists between the Gh1 and Cst3 genes. So, those pieces of evidence confirmed our results.

The human kallikreins, as a cluster of 15 serine protease genes, are located in the chromosomal band 19q13.4. It is considered a non-randomly reset area in numerous solid tumours like pancreatic cancer [74]. Such family genes are thought to play a role in the tumorigenesis process due to their high expression profile in hormone-dependent cancers (HDCs), in their regulation by the steroid hormones like estradiol, and in their processing abilities that have been found to be associated with tumour’s progression [75]. In comparison to pancreatic juices obtained from PDAC patients with noncancerous controls, a previous study [76] recognized some differential proteins in cancer samples like kallikrein 1 (Klk1). Accordingly, this indeed indicates the importance of this factor in the PDAC progression [77]. Previously, Jones et al. in their study proved that the release of kinins is catalyzed by rat Klk1 that is able to stimulate both prolactin and Gh1 secretion [78]. The Klk1 gene was also observed to be associated with prolactin-secreting cells within human Gh1-secreting adenomas [79]. Thus, as we found in the current study, in hormone-dependent cancers, Klk1, as the switching gene, can indirectly regulate both Gh1 and Cst3 and also provide a target for promotor-specific therapeutics in PDAC.

The relationship between the involved genes in triplet Lamc2 and {Dhps, CCr7}

Another biologically relevant triplet contains the laminin subunit gamma 2 (Lamc2) gene, as the switch gene, as well as a gene pair [deoxyhypusine synthase (Dhps), chemokine receptor 7 (Ccr7)] (Fig. 4). All the genes involved in such triplet were shown to participate in the cancer’s progression via the regulation of the "T-cell proliferation" pathway (Fig. 3). Recently, the relationship between T-cell conversion and some mutations contributing to an immunosuppressive tumour’s microenvironment within the immune escape of pancreatic cancer, has been proved [67, 80].

It was previously indicated that the binding of the chemokines to the chemokine receptor affects the signal transduction of the T lymphocytes activation pathway. This signal pathway plays a crucial role in the cancer’s progression because it will lead to an affinity increment of a surface adhesion molecule, which is known as integrins, and enhance T lymphocyte extravasation via capillary walls [81]. The C–C chemokine receptor type 7 is a protein, encoded by the Ccr7 gene in human beings. For this receptor, two ligands were recognized, including the (C–C motif) ligand 21 (CCL21) and chemokines (C–C motif) ligand 19 (CCL19/ELC). In this regard, concerning Ccr7, the relationships between this chemokine receptor expression in lymph node and tumour cells metastasis were represented in various types of cancers such as colorectal [82], gastric [83], esophageal [84], hepatocellular [85], thyroid cancer [86], malignant melanoma [87], breast cancer [88], cervical cancer [89], and pancreatic cancer [90]. Moreover, it is well-known that by ligating both the chemokines CCL21 and CCL19 to Ccr7, the homeostatic homing of T-cells is regulated into secondary lymphoid organs (SLOs) [91]. Recent studies performed on cancer revealed the essential roles of Ccl21-Ccr7 in neoangiogenesis, proliferation, activation, retention, and accumulation of T-cells at the persistent inflammation sites [92, 93].

Deoxyhypusine synthase (Dhps) through hypusine formation leads to the activation of eukaryotic initiation factor 5A (eIF-5A). The eIF-5A expression is significantly upregulated during the processes of dendritic cells (DCs) maturation and DC-mediated T lymphocyte activation. Hence, a new possibility may be presented by pharmacological interferences with hypusine formation for the modulation of DC function, which causes both T lymphocyte proliferation and activation [94]. Previously, Cheon Kim et al. in their study have reported the differentially expressed Dhps gene, which was observed to be significantly associated with canonical tumorigenesis and tumour’s progression in sporadic colorectal cancers [95].

Therefore, both Ccr7 and Dhps could affect the T lymphocyte proliferation pathway that plays a critical role in the cancer’s progression. Previous studies also indicated that both Ccr7 and Dhps play an important role in the recruitment of T-cells like diabetogenic Th1 cells, into inflamed islets, and thus in the pathogenesis of type 1 diabetes [96, 97]. So, the above-mentioned studies indeed confirmed the results of this study.

The protein encoded by the laminin subunit gamma 2 (Lamc2) gene is a crucial protein in the basal lamina affecting both the transformed and normal cell’s differentiation, adhesion, migration, survival, and phenotype. The basement membranes play a role as a mechanical barrier to tumour growth. However, these molecules, like laminins, are also considered the important autocrine factors created by various cancers, in order to promote tumorigenesis. The accumulating evidence in this regard revealed that the signalling network mediated by Lamc2 plays key roles in the progression, invasion, and migration of multiple kinds of cancers such as pancreas, stomach, tongue, bladder, and colorectal cancers [98]. Moreover, it is suggested as a potential therapeutic anticancer target to inhibit tumorigenesis [98]. Recently, Kosanam et al. have conducted a comprehensive PDAC tissue proteomic study that demonstrated both the prognostic and diagnostic potentials of Lamc2, and its behaviour in prospective biomarker panels for the improved pancreatic cancer diagnosis [99].

In another study, Jiang et al. have recently reported that there is a significant correlation between Lamc2 and Ccr7 of neutrophils in head and neck squamous cell carcinoma after adjusting for tumour purity [100]. The result of this study was in line with ours.

The gene regulatory relationship between the genes involved in the statistically significant triplets

The gene regulatory relationship between the genes involved in the statistically significant triplets was traced at this stage. As shown in Fig. 5, the regulatory relationships between Fhl5 → Peli2, Aox1 → Cnksr1, Gucy1a3 → Ephx1, Svep1 → Peli2, and Slc9a9 → Man1a2 were found with the maximum of seven intermediates. Previous studies have ascertained some of our findings. See below.

Chan et al. in their study demonstrated that Peli2 and Fhl5 genes, as the transcriptional targets of the myocyte enhancer factor 2 (MEF2), are simultaneously up-regulated in human neural progenitor/stem cells (hNPCs) [101]. Previously, it was shown that the MEF2 family of transcription factors is highly expressed in the brain, constituting a key cause of neurodegeneration [102]. In another research, Lee et al. in gene expression dataset GSE29174, have found some breast cancer-related genes after filtering the significantly and differentially expressed genes by fold change. Importantly, those genes, such as Peli2 and Fhl5, were found to be simultaneously down-regulated in breast cancers [103]. Recently, Namani et al. in their study have discovered the genes related to the Keap1 mutations, as a result, prognostic genes were observed to be greatly related to the upregulation of the nuclear factor erythroid-2-related factor 2 (NFE2L2/NRF2( path, among the Keap1-mutated Lung adenocarcinoma patients. Surprisingly, using the Position Weight Matrix (PWM) scan and ChIPseek web instruments, it was demonstrated that not only Peli2 and Fhl5, but also Nqo1 (the corresponding switch gene), as NRF2 binding sites, are involved in this pathway [104].

In this study, two other regulatory relationships were found, which are Aox1 → Cnksr1 and Svep1 → Peli2. In a previous study, Namani et al. have found that not only AOX1 serves as an NRF2 binding site, but Cnksr1, Peli2, and Svep1 also are highly correlated with the NRF2 pathway upregulation observed among the Keap1-mutated Lung adenocarcinoma patients [104].

It was demonstrated that Aox1 serves as a transcriptional target of MEF2 in human neural progenitor/stem cells (hNPCs) [101]. Furthermore, it was shown that a relationship exists between Aox1 and breast cancer, and a previous study reported that this factor is down-regulated in breast cancer cells [103]. As well, the Aox1 genewas reported as an NRF2 binding site in the Keap1-mutated Lung adenocarcinoma patients [104].

Another regulatory relationship was found between Gucy1a3 and Ephx1genes. Accordingly, both Ephx1 and Gucy1a3 genes play roles in some biological procedures such as cell cycle/apoptosis androgen signalling, transcriptional regulation, and signal transduction.

The last regulatory relationship was observed between Slc9a9 and Man1a2 genes. Namani et al. have also found that not only Cxcl12 (as a switching gene) serves as an NRF2 binding site, but also both Slc9a9 and Man1a2 are highly correlated with the NRF2 pathway upregulation among the KEAP1-mutated Lung adenocarcinoma patients [104].

As indicated in Fig. 5, such gene regulatory relationships are involved in the seven statistically significant triplets. Additionally, such triplets are controlled by Nqo1, Aox1, Tspan1, and Cxcl12 genes as the switch genes. Interestingly, the central roles of Nqo1, Aox1, and Cxcl12 are demonstrated in pancreatic cancer. See below.

Crnogorac-Jurcevic et al. [105] have conducted a large-scale immunoblotting analysis by including more than 900 primary antibodies on pancreatic cancer tissue, chronic pancreatitis, and normal pancreas. Numerous proteins like Aqx1, have been shown with different expressions between chronic pancreatitis and pancreatic cancer, indicating their potentially key roles in pancreatic carcinogenesis and as a biomarker for early recognition of pancreatic cancer [106]. The elevated Nqo1 expression in the pancreatic tumour is also represented in numerous early kinds of cancers like pancreatic intraepithelial neoplasia (PanINs) [107]. As mentioned earlier, Nqo1 is overexpressed in pancreatic tumour versus the associated normal tissue; however, catalase expression reduces in comparison to the levels related to normal pancreas tissues. Previously, Shaalan Beg et al. have conducted a clinical trial and as a result, explored a novel Nqo1 bioactivatable drug, named as ARQ761, to target pancreatic cancer treatment and to improve the chemotherapeutic impacts by metabolic modulation in pancreatic cancer [108]. Recently, Garg et al. in their study demonstrated that NFkB activity in pancreatic stellate cells could promote tumour growth by increasing the expression of Cxcl12, thereby preventing cytotoxic T-cells from infiltrating the tumour and killing cancer cells. In addition, they suggested some approaches for blocking Cxcl12 in pancreatic tumour cells, in order to increase antitumor immunity [109].

Although the current study provided new insight into the prognosis of PDAC, more effort should be performed to validate them clinically and extend such findings. We investigated the prognosis of PDAC through suggested switch genes captured in the context of the three-way interaction model. Results were in silico validated through various cancer-related gene expression datasets, and it can be concluded that the suggested switch genes might be significant for PDAC prognosis. Our findings are based on public gene expression datasets (switch genes investigation and validation), so a crucial extension of current work would be to learn whether clinical trials can reiterate the patterns.

Conclusion and future work

Nowadays, several high-throughput disease-related "omics" datasets are freely available. These datasets comprise valuable information on the disease-related pathways as well as their corresponding gene interactions. In this work, for the first time, we used the three-way interaction model to trace the switch genes that include PC. This approach has some advantages in comparison with the pairwise co-expression analyses; more precisely, the three-way interaction model can deal with the co-expression relationships’ dynamic nature. Therefore, the three-way interaction model can potentially shed light on some cellular alterations causes. In the current research, we suggested several switch genes as diagnostic biomarkers or potential therapeutic anticancer targets for anti-pancreatic cancer therapy. Surprisingly, in previous studies, it has been proved that all of our defined switch genes, including Lamc2 [99], Klk1 [77], Nqo1 [110], Aox1 [105], Tspan1 [111], and Cxcl12 [112], are closely associated with the pancreatic cancer’s progression. In this regard, some biological processes such as "Response to estradiol", "Embryonic limb morphogenesis pathway", "Positive regulation of cellular amide metabolic process", "The pathway of establishment of the protein localization to plasma membrane", "Type I interferon signalling pathway", "Positive regulation of the developmental growth", and "Regulation of T-cell proliferation" were also indicated to be associated with PDAC. Besides, we have introduced two biologically relevant triplets, namely Klk1 and {Gh1, Cst3} and Lamc2 and {Dhps, Ccr7}, which may be specifically involved in the onset of pancreatic cancer. Our investigation was performed in the Drug Bank database and showed that some drugs, including Peroxisomal acyl-coenzyme A oxidase 1 and Methocarbamol, are designed for Aox1 and Tspan1, as therapeutic anti pancreatic cancer targets, respectively. The next step may be studying these kinds of drugs on pancreatic cancer cell lines, in order to investigate their exact effect on PDAC suppression.

Availability of data and materials

The authors confirm that the data supporting the findings of this study are available within the article and its supplementary materials.


  1. Sarantis P, et al. Pancreatic ductal adenocarcinoma: treatment hurdles, tumor microenvironment and immunotherapy. World J Gastrointest Oncol. 2020;12(2):173.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Yeh JJ. Prognostic signature for pancreatic cancer: are we close? Future Oncol. 2009;5:313–21.

    Article  CAS  PubMed  Google Scholar 

  3. Collisson EA, et al. Subtypes of pancreatic ductal adenocarcinoma and their differing responses to therapy. Nat Med. 2011;17(4):500–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Shabani S, Khayer N, Motalebzade J. Characterization of pathways involved in colorectal cancer using real-time RT-PCR gene expression data. Gastroenterol Hepatol From Bed Bench. 2021;14(2):123.

    Google Scholar 

  5. Khayer N, et al. Rps27a might act as a controller of microglia activation in triggering neurodegenerative diseases. PLoS ONE. 2020;15(9):e0239219.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Iacobuzio-Donahue CA, et al. Exploration of global gene expression patterns in pancreatic adenocarcinoma using cDNA microarrays. Am J Pathol. 2003;162(4):1151–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Giulietti M, et al. Weighted gene co-expression network analysis reveals key genes involved in pancreatic ductal adenocarcinoma development. Cell Oncol. 2016;39(4):379–88.

    Article  CAS  Google Scholar 

  8. Zhou Z, et al. Ten hub genes associated with progression and prognosis of pancreatic carcinoma identified by co-expression analysis. Int J Biol Sci. 2018;14(2):124.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Skoda J, et al. Co-expression of cancer stem cell markers corresponds to a pro-tumorigenic expression profile in pancreatic adenocarcinoma. PLoS ONE. 2016;11(7):e0159255.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Stuart JM, et al. A gene-coexpression network for global discovery of conserved genetic modules. Science. 2003;302(5643):249–55.

    Article  CAS  PubMed  Google Scholar 

  11. Lee HK, et al. Coexpression analysis of human genes across many microarray data sets. Genome Res. 2004;14(6):1085–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Majd HA, et al. Two-way gene interaction from microarray data based on correlation methods. Iran Red Crescent Med J. 2016;18(6):e24373.

    Google Scholar 

  13. Khayer N, et al. Three-way interaction model with switching mechanism as an effective strategy for tracing functionally-related genes. Expert Rev Proteom. 2019;16(2):161–9.

    Article  CAS  Google Scholar 

  14. Khayer N, et al. Three-way interaction model to trace the mechanisms involved in Alzheimer’s disease transgenic mice. PLoS ONE. 2017;12(9):e0184697.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Khayer N, et al. Nkx3-1 and Fech genes might be switch genes involved in pituitary non-functioning adenoma invasiveness. Sci Rep. 2021;11(1):20943.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Zhang J, Ji Y, Zhang L. Extracting three-way gene interactions from microarray data. Bioinformatics. 2007;23(21):2903–9.

    Article  CAS  PubMed  Google Scholar 

  17. Lazar MA. Thyroid hormone action: a binding contract. J Clin Investig. 2003;112(4):497–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Brazma A, et al. ArrayExpress—a public repository for microarray gene expression data at the EBI. Nucleic Acids Res. 2003;31(1):68–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Barrett T, et al. NCBI GEO: mining millions of expression profiles—database and tools. Nucleic Acids Res. 2005;33(suppl_1):D562–6.

    CAS  PubMed  Google Scholar 

  20. Tomczak K, Czerwińska P, Wiznerowicz M. The cancer genome atlas (TCGA): an immeasurable source of knowledge. Contemp Oncol. 2015;19(1A):A68.

    Google Scholar 

  21. Zhang G, et al. Integration of metabolomics and transcriptomics revealed a fatty acid network exerting growth inhibitory effects in human pancreatic cancer. Clin Cancer Res. 2013;19(18):4983–93.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Badea L, et al. Combined gene expression analysis of whole-tissue and microdissected pancreatic ductal adenocarcinoma identifies genes specifically overexpressed in tumor epithelia. Hepatogastroenterology. 2008;55(88):2016–27.

    CAS  PubMed  Google Scholar 

  23. Zhang S, et al. Mast cell tryptase induces microglia activation via protease-activated receptor 2 signaling. Cell Physiol Biochem. 2012;29(5–6):931–40.

    Article  PubMed  Google Scholar 

  24. Gentleman R et al. Genefilter: methods for filtering genes from high-throughput experiments. R package version. 2015;1(1).

  25. Ritchie ME, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47–e47.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Maechler M. Package ‘diptest’. R Package Version 0.75–5. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2013.

    Google Scholar 

  27. Gunderson T. The fastLiquidAssociation Package. 2016.

  28. Li K-C. Genome-wide coexpression dynamics: theory and application. Proc Natl Acad Sci. 2002;99(26):16875–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Ho YY, et al. Modeling liquid association. Biometrics. 2011;67(1):133–41.

    Article  PubMed  Google Scholar 

  30. Willse JT, Willse MJT. Package ‘CTT’. 2018.

  31. Weisstein EW. Bonferroni correction. 2004.

  32. Schlitt T, Brazma A. Current approaches to gene regulatory network modelling. BMC Bioinform. 2007;8(6):1–22.

    Google Scholar 

  33. Zarnegarnia Y, et al. Application of fuzzy clustering in analysis of included proteins in esophagus, stomach and colon cancers based on similarity of Gene Ontology annotation. Koomesh. 2010;12(1):14–21.

    Google Scholar 

  34. Subramanian A, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. 2005;102(43):15545–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Consortium GO. The gene ontology (GO) database and informatics resource. Nucleic Acids Res. 2004;32(suppl_1):D258–61.

    Article  Google Scholar 

  36. Shannon P, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Aguirre-Gamboa R, et al. SurvExpress: an online biomarker validation tool and database for cancer gene expression data using survival analysis. PLoS ONE. 2013;8(9):e74250.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Zhang J, et al. International cancer genome consortium data portal—a one-stop shop for cancer genomics data. Database. 2011;2011:bar026.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Hezel AF, et al. Genetics and biology of pancreatic ductal adenocarcinoma. Genes Dev. 2006;20(10):1218–49.

    Article  CAS  PubMed  Google Scholar 

  40. Konduri S, Schwarz RE. Estrogen receptor β/α ratio predicts response of pancreatic cancer cells to estrogens and phytoestrogens. J Surg Res. 2007;140(1):55–66.

    Article  CAS  PubMed  Google Scholar 

  41. Xue J, et al. Important roles of estrogen receptor alpha in tumor progression and anti-estrogen therapy of pancreatic ductal adenocarcinoma. Life Sci. 2020;260:118302.

    Article  CAS  PubMed  Google Scholar 

  42. Abrahams A, Parker MI, Prince S. The T-box transcription factor Tbx2: its role in development and possible implication in cancer. IUBMB Life. 2010;62(2):92–102.

    CAS  PubMed  Google Scholar 

  43. Erickson JW, Cerione RA. Glutaminase: A hot spot for regulation of cancer cell metabolism? Oncotarget. 2010;1(8):734.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Son J, et al. Glutamine supports pancreatic cancer growth through a KRAS-regulated metabolic pathway. Nature. 2013;496(7443):101–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Blum R, Kloog Y. Metabolism addiction in pancreatic cancer. Cell Death Dis. 2014;5(2):e1065–e1065.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Kaistha BP, et al. PLAC8 localizes to the inner plasma membrane of pancreatic cancer cells and regulates cell growth and disease progression through critical cell-cycle regulatory pathways. Cancer Res. 2016;76(1):96–107.

    Article  CAS  PubMed  Google Scholar 

  47. Moerdyk-Schauwecker M, et al. Resistance of pancreatic cancer cells to oncolytic vesicular stomatitis virus: role of type I interferon signaling. Virology. 2013;436(1):221–34.

    Article  CAS  PubMed  Google Scholar 

  48. Murphy AM, et al. Vesicular stomatitis virus as an oncolytic agent against pancreatic ductal adenocarcinoma. J Virol. 2012;86(6):3073–87.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Barber GN. Vesicular stomatitis virus as an oncolytic vector. Viral Immunol. 2004;17(4):516–27.

    Article  CAS  PubMed  Google Scholar 

  50. Lichty BD, et al. Vesicular stomatitis virus: re-inventing the bullet. Trends Mol Med. 2004;10(5):210–6.

    Article  CAS  PubMed  Google Scholar 

  51. Wang BX, Rahbar R, Fish EN. Interferon: current status and future prospects in cancer therapy. J Interferon Cytokine Res. 2011;31(7):545–52.

    Article  PubMed  Google Scholar 

  52. Saloura V, et al. Evaluation of an attenuated vesicular stomatitis virus vector expressing interferon-β for use in malignant pleural mesothelioma: heterogeneity in interferon responsiveness defines potential efficacy. Hum Gene Ther. 2010;21(1):51–64.

    Article  CAS  PubMed  Google Scholar 

  53. Linge C, et al. Interferon system defects in human malignant melanoma. Cancer Res. 1995;55(18):4099–104.

    CAS  PubMed  Google Scholar 

  54. Sun WH, et al. Interferon-α resistance in a cutaneous T-cell lymphoma cell line is associated with lack of STAT1 expression. Blood J Am Soc Hematol. 1998;91(2):570–6.

    CAS  Google Scholar 

  55. Matin SF, et al. Impaired α-interferon signaling in transitional cell carcinoma: lack of p48 expression in 5637 cells. Cancer Res. 2001;61(5):2261–6.

    CAS  PubMed  Google Scholar 

  56. Pfeffer LM, et al. Human renal cancers resistant to IFN’s antiproliferative action exhibit sensitivity to IFN’s gene-inducing and antiviral actions. J Urol. 1996;156(5):1867–71.

    Article  CAS  PubMed  Google Scholar 

  57. Stojdl DF, et al. VSV strains with defects in their ability to shutdown innate immunity are potent systemic anti-cancer agents. Cancer Cell. 2003;4(4):263–75.

    Article  CAS  PubMed  Google Scholar 

  58. Naik S, Russell SJ. Engineering oncolytic viruses to exploit tumor specific defects in innate immune signaling pathways. Expert Opin Biol Ther. 2009;9(9):1163–76.

    Article  CAS  PubMed  Google Scholar 

  59. Virchow R. Krankheitswesen und Krankheitsursachen. Arch Pathol Anat Physiol Klin Med. 1880;79(2):185–228.

    Article  Google Scholar 

  60. Morris JP, Wang SC, Hebrok M. KRAS, Hedgehog, Wnt and the twisted developmental biology of pancreatic ductal adenocarcinoma. Nat Rev Cancer. 2010;10(10):683–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Kopp JL, et al. Identification of Sox9-dependent acinar-to-ductal reprogramming as the principal mechanism for initiation of pancreatic ductal adenocarcinoma. Cancer Cell. 2012;22(6):737–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Reichert M, et al. Developmental pathways direct pancreatic cancer initiation from its cellular origin. Stem Cells Int. 2016;2016:1–8.

    Article  Google Scholar 

  63. Neesse A, et al. Stromal biology and therapy in pancreatic cancer. Gut. 2011;60(6):861–8.

    Article  PubMed  Google Scholar 

  64. Bauer C, et al. Prevailing over T cell exhaustion: new developments in the immunotherapy of pancreatic cancer. Cancer Lett. 2016;381(1):259–68.

    Article  CAS  PubMed  Google Scholar 

  65. Cheng H, et al. The combination of systemic inflammation-based marker NLR and circulating regulatory T cells predicts the prognosis of resectable pancreatic cancer patients. Pancreatology. 2016;16(6):1080–4.

    Article  CAS  PubMed  Google Scholar 

  66. Hanke T, et al. High intratumoral FOXP3+ T regulatory cell (Tregs) density is an independent good prognosticator in nodal negative colorectal cancer. Int J Clin Exp Pathol. 2015;8(7):8227.

    PubMed  PubMed Central  Google Scholar 

  67. Cheng H, et al. KrasG12D mutation contributes to regulatory T cell conversion through activation of the MEK/ERK pathway in pancreatic cancer. Cancer Lett. 2019;446:103–11.

    Article  CAS  PubMed  Google Scholar 

  68. Girgert R, Emons G, Gründker C. Inhibition of growth hormone receptor by Somavert reduces expression of GPER and prevents growth stimulation of triple-negative breast cancer by 17β-estradiol. Oncol Lett. 2018;15(6):9559–66.

    PubMed  PubMed Central  Google Scholar 

  69. Perry JK, et al. Growth hormone and cancer: an update on progress. Curr Opin Endocrinol Diabetes Obes. 2013;20(4):307–13.

    Article  CAS  PubMed  Google Scholar 

  70. Leto G, Crescimanno M, Flandina C. On the role of cystatin C in cancer progression. Life Sci. 2018;202:152–60.

    Article  CAS  PubMed  Google Scholar 

  71. Staun-Ram E, Miller A. Cathepsins (S and B) and their inhibitor Cystatin C in immune cells: modulation by interferon-β and role played in cell migration. J Neuroimmunol. 2011;232(1–2):200–6.

    Article  CAS  PubMed  Google Scholar 

  72. Sze L, et al. Impact of growth hormone on cystatin C. Nephron Extra. 2013;3(1):118–24.

    Article  PubMed  PubMed Central  Google Scholar 

  73. Zhu X-R, et al. Corticosteroids significantly increase cystatin C levels in the plasma by promoting cystatin C production in rats. Renal Fail. 2019;41(1):698–703.

    Article  CAS  Google Scholar 

  74. Yousef GM, et al. In-silico analysis of kallikrein gene expression in pancreatic and colon cancers. Anticancer Res. 2004;24(1):43–52.

    CAS  PubMed  Google Scholar 

  75. Myers SA. Kallikrein gene regulation in hormone-dependent cancer cell lines. Brisbane City: Queensland University of Technology; 2003.

    Google Scholar 

  76. Chen R, et al. Quantitative proteomic profiling of pancreatic cancer juice. Proteomics. 2006;6(13):3871–9.

    Article  CAS  PubMed  Google Scholar 

  77. Pan S, Brentnall TA, Chen R. Proteomics analysis of bodily fluids in pancreatic cancer. Proteomics. 2015;15(15):2705–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Jones T, Figueroa C, Bhoola K. Bioregulatory role of the kallikrein–kinin system in the normal pituitary gland and its tumours. Eur J Endocrinol. 1992;127(6):481–4.

    Article  CAS  Google Scholar 

  79. Komatsu N, et al. Proteolytic processing of human growth hormone by multiple tissue kallikreins and regulation by the serine protease inhibitor Kazal-Type5 (SPINK5) protein. Clin Chim Acta. 2007;377(1–2):228–36.

    Article  CAS  PubMed  Google Scholar 

  80. Yaguchi T, et al. The mechanisms of cancer immunoescape and development of overcoming strategies. Int J Hematol. 2011;93(3):294–300.

    Article  PubMed  Google Scholar 

  81. Slaney CY, Kershaw MH, Darcy PK. Trafficking of T cells into tumors. Cancer Res. 2014;74(24):7168–74.

    Article  CAS  PubMed  Google Scholar 

  82. Günther K, et al. Prediction of lymph node metastasis in colorectal carcinoma by expressionof chemokine receptor CCR7. Int J Cancer. 2005;116(5):726–33.

    Article  PubMed  Google Scholar 

  83. Mashino K, et al. Expression of chemokine receptor CCR7 is associated with lymph node metastasis of gastric carcinoma. Cancer Res. 2002;62(10):2937–41.

    CAS  PubMed  Google Scholar 

  84. Ding Y, et al. Association of CC chemokine receptor 7 with lymph node metastasis of esophageal squamous cell carcinoma. Clin Cancer Res. 2003;9(9):3406–12.

    CAS  PubMed  Google Scholar 

  85. Schimanski CC, et al. Chemokine receptor CCR7 enhances intrahepatic and lymphatic dissemination of human hepatocellular cancer. Oncol Rep. 2006;16(1):109–13.

    CAS  PubMed  Google Scholar 

  86. Sancho M, et al. Expression and function of the chemokine receptor CCR7 in thyroid carcinomas. J Endocrinol. 2006;191(1):229–38.

    Article  CAS  PubMed  Google Scholar 

  87. Takeuchi H, et al. CCL21 chemokine regulates chemokine receptor CCR7 bearing malignant melanoma cells. Clin Cancer Res. 2004;10(7):2351–8.

    Article  CAS  PubMed  Google Scholar 

  88. Cabioglu N, et al. CCR7 and CXCR4 as novel biomarkers predicting axillary lymph node metastasis in T1 breast cancer. Clin Cancer Res. 2005;11(16):5686–93.

    Article  CAS  PubMed  Google Scholar 

  89. Kodama J, et al. Association of CXCR4 and CCR7 chemokine receptor expression and lymph node metastasis in human cervical cancer. Ann Oncol. 2007;18(1):70–6.

    Article  CAS  PubMed  Google Scholar 

  90. Nakata B, et al. Chemokine receptor CCR7 expression correlates with lymph node metastasis in pancreatic cancer. Oncology. 2008;74(1–2):69–75.

    Article  CAS  PubMed  Google Scholar 

  91. Moschovakis GL, Förster R. Multifaceted activities of CCR7 regulate T-cell homeostasis in health and disease. Eur J Immunol. 2012;42(8):1949–55.

    Article  CAS  PubMed  Google Scholar 

  92. Pickens SR, et al. Characterization of CCL19 and CCL21 in rheumatoid arthritis. Arthritis Rheum. 2011;63(4):914–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  93. Pickens SR, et al. Role of the CCL21 and CCR7 pathways in rheumatoid arthritis angiogenesis. Arthritis Rheum. 2012;64(8):2471–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  94. Kruse M, et al. Inhibition of CD83 cell surface expression during dendritic cell maturation by interference with nuclear export of CD83 mRNA. J Exp Med. 2000;191(9):1581–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  95. Kim JC, et al. Gene expression profiling: canonical molecular changes and clinicopathological features in sporadic colorectal cancers. World J Gastroenterol WJG. 2008;14(43):6662.

    Article  CAS  PubMed  Google Scholar 

  96. Shan Z, et al. CCR7 directs the recruitment of T cells into inflamed pancreatic islets of nonobese diabetic (NOD) mice. Immunol Res. 2014;58(2–3):351–7.

    Article  CAS  PubMed  Google Scholar 

  97. Cui Y, et al. Differential expression network analysis for diabetes mellitus type 2 based on expressed level of islet cells. Ann d’endocrinol. 2016;77:22–9.

    Article  Google Scholar 

  98. Garg M, Braunstein G, Koeffler HP. LAMC2 as a therapeutic target for cancers. Expert Opin Ther Targets. 2014;18(9):979–82.

    Article  CAS  PubMed  Google Scholar 

  99. Kosanam H, et al. Laminin, gamma 2 (LAMC2): a promising new putative pancreatic cancer biomarker identified by proteomic analysis of pancreatic adenocarcinoma tissues. Mol Cell Proteom. 2013;12(10):2820–32.

    Article  CAS  Google Scholar 

  100. Jiang P, et al. Identification of therapeutic and prognostic biomarkers of lamin C (LAMC) family members in head and neck squamous cell carcinoma. Med Sci Monit Int Med J Exp Clin Res. 2020;26:e925735–41.

    CAS  Google Scholar 

  101. Chan SF, et al. Transcriptional profiling of MEF2-regulated genes in human neural progenitor cells derived from embryonic stem cells. Genom Data. 2015;3:24–7.

    Article  PubMed  Google Scholar 

  102. Chen H, et al. Neuroprotective and neurogenic effects of novel tetramethylpyrazine derivative T-006 in Parkinson’s disease models through activating the MEF2-PGC1α and BDNF/CREB pathways. Aging (Albany NY). 2020;12(14):14897.

    Article  CAS  Google Scholar 

  103. Lee C-H, et al. MicroRNA-regulated protein-protein interaction networks and their functions in breast cancer. Int J Mol Sci. 2013;14(6):11560–606.

    Article  PubMed  PubMed Central  Google Scholar 

  104. Namani A, et al. Systematic identification of multi omics-based biomarkers in KEAP1 mutated TCGA lung adenocarcinoma. J Cancer. 2019;10(27):6813.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  105. Crnogorac-Jurcevic T, et al. Proteomic analysis of chronic pancreatitis and pancreatic adenocarcinoma. Gastroenterology. 2005;129(5):1454–63.

    Article  CAS  PubMed  Google Scholar 

  106. Jenkinson C, et al. Biomarkers for early diagnosis of pancreatic cancer. Expert Rev Gastroenterol Hepatol. 2015;9(3):305–15.

    Article  CAS  PubMed  Google Scholar 

  107. DeNicola GM, et al. Oncogene-induced Nrf2 transcription promotes ROS detoxification and tumorigenesis. Nature. 2011;475(7354):106–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  108. Beg MS, et al. Using a novel NQO1 bioactivatable drug, beta-lapachone (ARQ761), to enhance chemotherapeutic effects by metabolic modulation in pancreatic cancer. J Surg Oncol. 2017;116(1):83–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  109. Garg B, et al. Nfκb in pancreatic stellate cells reduces infiltration of tumors by cytotoxic T cells and killing of cancer cells, via up-regulation of CXCL12. Gastroenterology. 2018;155(3):880–91.

    Article  CAS  PubMed  Google Scholar 

  110. Lewis AM, et al. Targeting NAD (P) H: quinone oxidoreductase (NQO1) in pancreatic cancer. Mol Carcinog Publ Cooper Univ Texas MD Anderson Cancer Center. 2005;43(4):215–24.

    CAS  Google Scholar 

  111. Zhang X, et al. TSPAN1 upregulates MMP2 to promote pancreatic cancer cell migration and invasion via PLCγ. Oncol Rep. 2019;41(4):2117–25.

    CAS  PubMed  PubMed Central  Google Scholar 

  112. Sleightholm RL, et al. Emerging roles of the CXCL12/CXCR4 axis in pancreatic cancer progression and therapy. Pharmacol Ther. 2017;179:158–70.

    Article  CAS  PubMed  Google Scholar 

Download references


We would like to thank Institute for Research in Fundamental Sciences (IPM) for providing us with research facilities.


This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

Author information

Authors and Affiliations



ZSE performed the experiments. NK, as the corresponding author, designed and managed the study. AT and AA analyzed and interpreted the data. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Nasibeh Khayer.

Ethics declarations

Ethics approval and consent to participate

Since the current investigation that was performed using a public GEO dataset (accession number: GDS4336), can be considered a computational survey, obtaining both ethics approval and consent statements was not mandatory in this research.

Consent for publication

Since the current investigation that was performed using a public GEO dataset (accession number: GDS4336), can be considered a computational survey, obtaining consent for publication was not mandatory in this research.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

S1 Table. A list of candidate switch genes. S2 Table. A list of statistically significant triplets. S3 Table. Detailed information of Gene Regulatory Network. S4 Figure.The prognostic power of suggested switch genes through glioblastoma, hepatocellular carcinoma and stomach adenocarcinoma as PDAC-unrelated datasets.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Shokati Eshkiki, Z., Khayer, N., Talebi, A. et al. Novel insight into pancreatic adenocarcinoma pathogenesis using liquid association analysis. BMC Med Genomics 15, 30 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Pancreatic ductal adenocarcinoma
  • Liquid association analysis
  • Three-way gene interaction
  • Gene set enrichment analysis
  • Therapeutic targets