Xiao-Xi Li,Hong Li,Li-Quan Jin,Yun-Bo Tan
Exploration and Validation of Pancreatic Cancer Hub Genes Based on Weighted Gene Co-Expression Network Analysis and Immune Infiltration Score Analysis Li et al
DPYD,pancreatic cancer,tumor immunity,WGCNA
Objective To find pancreatic cancer (PC)-related hub genes based on weighted gene co-expression network analysis (WGCNA) construction and immune infiltration score analysis and validate them immunohistochemically by clinical cases, to generate new concepts or therapeutic targets for the early diagnosis and treatment of PC. Material and Methods In this study, WGCNA and immune infiltration score were utilized to identify the relevant core modules of PC and the hub genes within these core modules. Results Using WGCNA analysis, data from PC and normal pancreas integrated with TCGA and GTEX were analyzed and brown modules were chosen from the six modules. Five hub genes, including DPYD, FXYD6, MAP6, FAM110B, and ANK2, were discovered to have differential survival significance via validation tests utilizing survival analysis curves and the GEPIA database. The DPYD gene was the only gene associated with PC survival side effects. Validation of the Human Protein Atlas (HPA) database and immunohistochemical testing of clinical samples showed positive results for DPYD expression in PC. Conclusion In this study, we identified DPYD, FXYD6, MAP6, FAM110B, and ANK2, as immune-related candidate markers for PC. Only the DPYD gene had a negative impact on the survival of PC patients. Through validation of the HPA database and immunohistochemical testing of clinical cases, we believe that the DPYD gene brings novel ideas and therapeutic targets in the diagnosis and treatment of PC.
Exploration and Validation of Pancreatic Cancer Hub Genes Based on Weighted Gene Co-Expression Network Analysis and Immune Infiltration Score Analysis Li et al

To find pancreatic cancer (PC)-related hub genes based on weighted gene co-expression network analysis (WGCNA) construction and immune infiltration score analysis and validate them immunohistochemically by clinical cases, to generate new concepts or therapeutic targets for the early diagnosis and treatment of PC.

In this study, WGCNA and immune infiltration score were utilized to identify the relevant core modules of PC and the hub genes within these core modules.

Using WGCNA analysis, data from PC and normal pancreas integrated with TCGA and GTEX were analyzed and brown modules were chosen from the six modules. Five hub genes, including DPYD, FXYD6, MAP6, FAM110B, and ANK2, were discovered to have differential survival significance via validation tests utilizing survival analysis curves and the GEPIA database. The DPYD gene was the only gene associated with PC survival side effects. Validation of the Human Protein Atlas (HPA) database and immunohistochemical testing of clinical samples showed positive results for DPYD expression in PC.

In this study, we identified DPYD, FXYD6, MAP6, FAM110B, and ANK2, as immune-related candidate markers for PC. Only the DPYD gene had a negative impact on the survival of PC patients. Through validation of the HPA database and immunohistochemical testing of clinical cases, we believe that the DPYD gene brings novel ideas and therapeutic targets in the diagnosis and treatment of PC.

Pancreatic cancer (PC) is a malignant tumor with a very high mortality rate for which there is no effective treatment other than surgery.1,2 PC is prone to metastasis and invasion in the early stage, and there is no effective means for early detection. Most patients often have cancer metastasis at the time of diagnosis and lose the opportunity for surgical resection.3 For more than 40 years, with the advancement of technology and increasing research on pancreatic cancer, some treatments such as targeted therapy, neoadjuvant therapy, and translational therapy have gradually emerged in people’s vision, but their therapeutic effects on pancreatic cancer are often limited, and the survival rate of pancreatic cancer patients remains very low.4 The unique microenvironment of PC promotes cancer cell proliferation and immune escape. Two main features characterize the immune microenvironment of PC: the dense fibrous tissue proliferation and extensive immunosuppression. Both of these features can promote immune surveillance escape of PC cells and lead to their proliferation by directly suppressing anti-tumor immunity.5 Over the last decade, substantial progression in molecular subtyping for PDAC has facilitated the understanding of molecular pathogenesis and provided clues for advanced therapy designing.6–11 However, to date, there is still lacking an immune feature- based molecular subtyping for better understand the molecular mechanism of pancreatic cancer invasion and progression in order to facilitate early diagnosis and effective treatment. PC treatment should emphasize on early diagnosis, prevention of metastasis, and prevention of immune evasion. Moreover, the current understanding of PC is incomplete; consequently, it is essential to comprehend the molecular mechanisms of PC and to find novel prognostic biomarkers and potential therapeutic targets associated with immunity in PC. Weighted gene co-expression network analysis (WGCNA) is a promising emerging bioinformatics analysis method.12 It is primarily a systems biology approach to describe gene association patterns between different samples, to identify highly synergistic sets of genes, and ultimately screen essential hub genes for a specific trait as alternate biomarker genes or therapeutic targets based on the inline nature of the gene set and its association with clinical performance. In this paper, we screened for immune signature-related hub genes in PC using the WGCNA method to select target genes from among them. This approach provides potential biomarkers and new therapeutic targets for PC treatment.

The patient RNA sequencing and clinical information data were downloaded from the PC cohort in the TCGA database (https://www.cancer.gov/aboutnci/organization/ccg/research/structural-genomics/tcga) on July 7, 2022. The expression data was 182. The clinical data was 196. At the same time, the RNA sequencing data and clinical information of normal pancreatic tissues were downloaded from the GTEX database (https://www.g texportal.org/home/), with 167 gene expression data and 204 clinical data. Combining the two data using R language, the PC tissues that could not be distinguished between the two were then screened out from normal tissues using principal component analysis (PCA). We used PCA, a multifactorial statistical method, to reduce the dimensionality of the data. The DPYD gene immunohistochemical test results for the clinical case were obtained from the First Affiliated Hospital of Dali University, Yunnan, China. The specific workflow is as follows (Figure 1).

Limma13 (linear models for microarray) data is a differential expression screening method based on generalized linear models. In this study, we utilized the R package “Limma (version 3.40.6)” to differentially analyze the differentially expressed genes (DEG) of the combined gene expression profiles14 to identify the differentially up-regulated and down-regulated genes between normal pancreatic tissues and PC tissues. We analyzed the functional enrichment of differentially expressed genes (DEGs) in R using Gene Ontology (GO) for up-regulated and down-regulated genes.

A total of 22,515 differential genes were evaluated and a gene co-expression network was constructed using the “WGCNA” R language package for these differential genes. Co-expression modules are groups of genes that share a high degree of topological overlap similarity and genes within the same module typically co-express to a greater extent. We constructed a hierarchical clustering tree by clustering genes according to a dissimilarity matrix representing gene linkage, and used a dynamic tree-cutting algorithm to identify gene co-expression modules. The goodSamplesGenes method of the R package WGCNA15,16 was used to remove outlier genes and samples. A scale-free co-expression network was subsequently constructed using WGCNA to obtain six co-expression modules, with the gray module representing the set of genes that could not be assigned to a module.

Immune scores were assigned to CD3+ T cells, CD8+ T cells, cytotoxic lymphocytes, NK cells, B lymphocytes, monocytes, bone marrow dendritic cells, neutrophils and endothelial cells, and fibroblasts in each of the six co-expression modules. Using the R language package immune-oncology biological research (IOBR),17 this score was used to quantify the degree of immune infiltration in PC and its response to therapy. This score quantifies the degree of immune infiltration in PC and represents the immune properties of the tumor.18 The immune score was used as the immune signature of PC to correlate with the modules screened. The hub genes with the highest association and variability with the immune signature and modules were identified.

Gene Module Membership (MM) indicates the relationship between genes and modules. GeneTrait Significance (GS) refers to the relationship between the clinical characteristics or immune information of the gene and the gene itself. The criteria for hub genes in the immune-related gene module network are “cor.geneTraitSignificance” > 0.6 and “cor.gene ModuleMembership” > 0.8. Screening and filtering under conditions satisfying this criterion yields hub genes with high relevance for the desired immune profile.15

Using the R package SURVIVAL, Kaplan-Meier survival curves and Log rank tests for differences in overall survival between high and low-expression groups of hub genes were utilized to identify hub genes having prognostic importance for survival. Log < 0.01 and P < 0.05 were judged statistically significant.

The GEPIA online database (http://gepia.cancer-pku.cn/) is a website that incorporates data from the TCGA and GTEX databases, and has interactive and customizable features.19 The screened hub genes were then uploaded to the GEPIA website so that their expression and survival curves could be analyzed. The verification of the dependability of the hub genes was investigated in this study.

We verified the hub genes that were selected at the protein level. The Human Protein Atlas (HPA) database employs transcriptomics and proteomics technologies to examine protein expression levels at the RNA and protein levels in various human tissues and organs.20 It has nearly 26,000 antibodies, and all results have immunohistochemical (IHC) staining in various tissues. In normal pancreatic and PC tissue, we found IHC markers for the expressed proteins of hub genes in HPA and downloaded them for comparison and localization.

In this study, IHC staining was applied to the cancer tissues of three PC patients from the First Affiliated Hospital of Dali University.21,22 The cancer tissues were formalin-fixed and embedded in paraffin wax and then three different PC tissue sections of the patients were taken. The tissue sections were sequentially treated with alcohol in the order of xylene, 100%X2, 95%, 90%, 80%, and 70%. Tissue antigen was repaired by boiling with sodium citrate buffer and adding endogenous peroxidase inhibitors. It was then blocked in 5% goat serum for 1 hour to prevent non-specific binding. Subsequently, anti-DPYD antibodies (dilution ratio: 1:300; Proteintech:27662-1-AP) were dropped onto the tissue sections separately and incubated overnight at 4 °C. Secondary antibodies were then dropped and incubated. Then, the sections were blocked and photographed after adding diaminobenzidine for color development and then hematoxylin re-staining. Finally, the staining results were evaluated by experienced researchers.

There were 182 samples from the TCGA database with 4 normal or para-cancerous tissue samples and 178 (182–4) PC samples. Then, 167 normal pancreatic tissue samples from the TCGA and GTEX database and 4 normal tissue samples were enrolled, for a total of 171 (167+4) samples (Figure 2A). Five normal pancreatic tissue samples were excluded after filtering by PCA, which can better distinguish tumor samples from normal ones. The first component, represented in the red color was significantly distinguished from the second component, represented in the blue color (Figure 2B), and further analysis and screening was conducted on the remaining 344 (178+171-5) samples.

Using the “Limma (version 3.40.6)” function of the R program, samples of 173 normal pancreatic tissues and 171 PC tissues were detected, and DEGs were defined as those with |log2 (fold change) | > 1 and P < 0.01. In all, 22,515 differential genes were identified, with 333 up-regulated and 22,182 down-regulated genes (Figure 2C). It is assumed that the abnormally high number of down-regulated genes is due to the integration of the normal pancreas gene expression profile into the GTEX data. At the same time, DEGs were also analyzed by GO enrichment to understand their potential biological functions. As can be seen from Figure 2D and E, the up-regulated DEGs in cancer are mainly involved in small molecule ribonucleic acid, aminoacyl tRNA biosynthesis, interactions between viral proteins and cytokines and cytokine receptors, chemokine signaling pathway, IL-17 signaling pathway and P53 signaling pathway (Figure 2D). Most down-regulated DEGs were involved in the macromolecular synthesis, cytosolic nitrogenation biosynthesis, metal ion binding, and biosynthesis of base-containing chemicals (Figure 2E). The biofunctional annotations of these genes with significant alterations, such as the chemokine signaling route and the IL-17 signaling pathway, are comparable with some of the recognized dysfunctions of PC. On the other hand, this demonstrates the excellent reliability of our findings.

By running the MCP-counter algorithm, we brought 178 PC samples with gene expression information into the R package, thus quantifying the transcriptome data of the samples into an absolute abundance of 8 immune cells and 2 stromal cell populations. This algorithm has been used in several gene expression studies to quantify the abundance of immune cells in different samples.23,24 We used the scores of individual immune cells of PC calculated by MCPcounter as the immune signature of PC.

WGCNA was used to construct an expression matrix of 22,515 DEGs and a co-expression network from the immune infiltrating cell score data calculated by MCPcount. We clustered the samples of DEGs using Pearson correlation coefficients, and after removing outliers, drew a clustering tree of the samples. First, we calculated the soft threshold power and raised the co-expression similarity to calculate the neighborhood. We used the pickSoftThreshold function in WGCNA to analyze its network topology. The figure shows that the average connectivity is best when the plotted scalar independence reaches 0.86 and the soft threshold power is set to 8. Therefore, the soft threshold β is set to 8 (Figure 3A and B). After establishing a soft threshold value of 8, the created gene network conforms to the scale-free distribution, allowing the hub modules and genes to be searched. We used the WGCNA R program to create the gene network further and identify the modules. We set the minimum module size for its clustering analysis to 30 and the sensitivity to 3, and after clustering DEGs with similar expression patterns into the same modules, modules set to show a difference in cut height < 0.25 were merged (Figure 3C). We built the connection matrix and the topological overlap matrix as an extension. Six modules were discovered based on the average hierarchical clustering and dynamic tree cropping, whose distribution spans from black, blue, brown, gray, red, and yellow modules (Figure 3D and E). Among them, the gray modules were genes that could not be integrated into other modules. The modules screened in this study were correlated with the immune features represented by the immune infiltration score of PC calculated by the MCPcounter. This analysis showed (Figure 4A) that the brown modules showed a high correlation with the immune features (cor=0.89, p=2.2x10-105) and the immune features represented by T cell and Monocytic-lineage. The highest correlation was found with the immune features represented by T cell and Monocytic-lineage. The findings clearly show that genes within the brown module contribute to pancreatic carcinogenesis and development and are closely associated with immunity.

The analysis shows that MM and GS of the two modules representing T cell and Monocytic-lineage in the brown module are strongly positively correlated with each other, with correlation coefficients of r=0.94 and r=0.91, respectively (Figure 4B and C). Under the condition that “cor.geneModuleMembership” > 0.8 threshold and satisfying “cor.geneTraitSignificance” > 0.6 thresholds, we screened 330 hub genes under the more strongly correlated T cell module.

Based on the expression data and clinical information of PC samples in TCGA, we examined the potential association between the expression of 330 genes screened in the brown module and patient survival. After univariate Cox proportional risk regression analysis of the 330 genes by the survival package in R, 26 genes with prognostic relevance for survival were screened (log < 0.01, P < 0.05) (Figure 5A). For the 26 screened genes, validation using the GEPIA website resulted in survival differences for only 5 genes, including DPYD, FXYD6, MAP6, FAM110B, and ANK2. We defined these 5 genes as the final hub genes. As shown by the survival analysis curves (Figure 5B), only the DPYD gene has a negative impact on the survival of PC patients. That is, DPYD may be able to stimulate the proliferation and differentiation of PC cells, hence shortening patient survival time.

To get additional insight into the protein expression level of the hub gene DPYD, we compared the IHC profiles of PC and normal pancreatic tissues from the HPA database. The comparison revealed that DPYD was substantially expressed in the cytoplasm of PC cells, but it was nearly absent in normal pancreatic tissue cells (Figure 6). Therefore, we propose that the DPYD gene has the potential to become a novel PC-specific marker.

We took pathological sections from surgically confirmed PC patients from the First Affiliated Hospital of Dali University in Yunnan, China. The heterogeneous expression of the DPYD gene in PC cells was detected clinically by IHC analysis in all cases (Figure 7), thus further validating the clinical application of the DPYD gene as an immune-related hub gene. We did not use normal pancreatic tissues for comparison due to the harsh conditions of obtaining normal pancreatic tissues in the clinic, which is a shortcoming of this study.

The mortality rate of PC has been high and the cure rate and prognosis are unsatisfactory. The main reason is that PC is prone to metastasis and immune escape due to its unique immune microenvironment. In addition, early PC metastasis is associated with Epithelial-Mesenchymal Transition (EMT). The immune infiltrating cells within the EMT play a key role.25,26 There is an inextricable link between immunity and PC. The molecular mechanisms involved in the immune microenvironment of PC are still unclear, so it is crucial to explore the genes associated with immunity in patients with PC.

According to GO enrichment analysis, the 333 up-regulated DEGs in our study were connected to small molecule RNAs, amyl tRNA production, viral protein-cytokine interactions, chemokine signaling pathways, IL-17 signaling pathways, and P53 signaling pathways in cancer. Pancreatic cancer is also linked to a chemokine signaling pathway with substantial heterogeneity in enrichment analysis. Chemokines are members of a family of tiny molecules of cytokines whose primary function is to draw immune cells to the site of inflammation. Their downstream signaling pathways, JAK/STAT, Ras, ERK, and Akt, are also linked to PC, encouraging cancer proliferation, differentiation, and metastasis,27–30 and primarily activate these pathways.31 According to a recent study, the inflammatory chemokine CXCL3 stimulates the spread of PC via a brand-new cancer escape mechanism hijacked by myofibroblasts.32 Additionally, research has revealed that the IL-17 and p53 pathways are involved in the initiation, progression, and prognosis of PC and that the action of IL-17 promotes PC metastasis and invasion as well as predicts PC prognosis and gemcitabine therapy effectiveness.33,34 The high P53 mutation level in the pancreas also demonstrates that our projections generally match the facts.35

By WGCNA construction, we identified five crucial genes—DPYD, FXYD6, MAP6, FAM110B, and ANK2—linked to immunity and survival. In prior investigations and research, FXYD6, MAP6, FAM110B, and ANK2 were not significantly related to PC. FXYD6 is significantly altered in hepatocellular carcinoma. Its inhibition reduces the growth of mouse hepatocellular carcinoma cells, suggesting that it may serve as a novel therapeutic target for hepatocellular carcinoma.36 MAP6 is an intraluminal protein that stimulates the coiling of neuronal microtubules.37 Recent research indicates that MAP6 may influence the malignant course of osteosarcoma via the bax/bcl and Wnt/-catenin signaling pathways.38 FAM110B is a member of the gene family FAM110. This family protein is localized to the centrosome and accumulates at the center of microtubule organization in interphase and at the mitotic spindle pole; FAM110B expression inhibits cell cycle progression in the G1 phase.39 FAM110B has also been demonstrated to decrease the proliferation and differentiation of lung small-cell carcinoma.40 Ank2 is mainly expressed in genes associated with human breast cancer and autism risk.41

Five genes, including DPYD, FXYD6, MAP6, FAM110B, and ANK2, were discovered to have survival differences after the 26 genes examined in this study were validated in the GEPIA website. Survival analysis curves revealed that only the DPYD gene negatively impacts the survival of PC patients, ie, DPYD can stimulate the proliferation and differentiation of PC cells, reducing patient survival time. The protein that DPYD expresses is dihydropyrimidine dehydrogenase. This protein is a catabolic enzyme for pyrimidine. It is a rate-limiting factor and an initiator in the uracil and thymidine catabolic pathway.42,43 DPYD is currently of interest primarily in cancer patients receiving 5-fu uracil chemotherapy due to the increased risk of chemotherapy toxicity caused by DPYD deficiency resulting from the deficiency of this gene.44–46 It has also been shown that elevated DPYD expression in hepatocellular carcinoma cells promotes the formation and progression of hepatocellular carcinoma and increases the likelihood of cancer cell infiltration and spread.47 In a study using mucuna proteins, a chemopreventive drug for PC, reduction of DPYD led to an elevation of STAT3 expression, along with PC cell proliferation, and inhibited the capacity of cancer cells to penetrate and metastasize.48 This is consistent with our finding that DPYD can stimulate the proliferation and differentiation of PC cells. In addition, there appears to be a relationship between DPYD and the invasive capacity of cancer cells, which requires additional investigation.

Comparative analysis was performed using the IHC profiles of PC and normal pancreatic tissues from the HPA database. The comparison revealed that DPYD was highly expressed in the cytoplasm of PC cells but nearly nonexistent in normal pancreatic tissue cells; heterogeneity of the DPYD gene in PC cells was also observed in clinical PC cases. Consequently, it is evident that the DPYD gene has the potential to become a potential novel marker for PC and may become a new therapeutic target for PC research.

We used a combination of bioinformatics analysis, experimental verification, and dataset cross-validation to study PC-related DEGs and obtained one novel hub gene (DPYD) that may be associated with the prognosis of PC. Based on these results, the present study provides new ideas for studying the molecular mechanisms of PC. Previous studies on WGCNA usually only focused on clinical data of genes to discover clinically relevant hub genes. However, in this study we used WGCNA and immune-related data to identify hub genes with immune relevance. The results of this study may provide new ideas for studying the molecular mechanism of PC and provide valuable clues for its treatment.