We organized the 1689 MEs from the HMDB database (Supplementary Table S1) (Wishart et al., 2018). We used 180 cases RNA-seq of pancreatic ductal adenocarcinoma (PAAD and 4 cases of normal tissue), with clinical information for WGCNA analysis in TCGA (http://www.cancer.gov/tcga). Meanwhile, we also organized 36 normal samples and 36 tumor samples (GSE15471) of microarray from the GEO database.

The microarray data GSE15471 were downloaded from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2005). GSE15471 dataset contained 36 PC tissues and 36 paired adjacent normal tissues. For further analysis, the data was profiled on the GPL570 platform (Affymetrix, Human Genome U133 plus 2.0 Array). Differentially expressed genes (DEGs) were identified using the R software (version 4.0.5) package limma. The cutoff criteria were p_value <0.05, and |log2(fold change)(FC)| > 0.586, while RNAs with low expression values was filtered out.

The R software (version 4.0.5) package DEGseq (Wang et al., 2010) were used to analyze the RNA-seq from the TCGA-PAAD datasets and to identify DEGs. The cutoff criteria were p_value <0.05 and |log2(fold change)(FC)|>0.586, and filtered out RNAs with low expression values (the number of samples with normalized gene expression <1 was more than half of the total number of samples).

RNA-seq and microarray datasets used above was analyzed by WGCNA to construct a co-expression network (Langfelder and Horvath, 2008). The expression profiles of MEs were obtained from the TCGA dataset and GSE15471 datasets, respectively. Then the samples were clustered through the systematic cluster tree to determine any outliers.

The soft-thresholding function was used calculate the power parameter. The dynamic tree cut method was used to identify the modules of co-expressed gene. Then a dendrogram of genes was obtained using a hierarchical clustering approach, which was based on dissimilarity of the unsigned topological overlap matrix (TOM). Lastly, genes with similar expression profiles were grouped and network modules was produced.

We chose two modules from the TCGA and GSE15471 by screening edges and the criterion of weight value was >0.02 for the brown module (TCGA) and the black module (GSE15471) respectively. To visualize the co-expression network and identify the nodes and hub genes, the results were input into Cytoscape to (v3.8.0; downloaded from the website https://cytoscape.org/) (Shannon et al., 2003).

We used the KEGG Orthology-Based Annotation System (KOBAS) (Xie et al., 2011) for GO and KEGG pathway enrichment analysis. A p_value <0.05 was set as the cutoff criterion.

CIBERSORT (https://cibersort.stanford.edu/index.php) (Newman et al., 2015) was used to analyze the mRNA expression data for immune microenvironment components determination.

We used the HMDB database and identified a total of 1689 metabolism-related genes (Supplementary Table S1). Then we extracted data for 180 samples of PAAD and 4 samples of normal tissue from the TCGA dataset and used the GSE15471 dataset comprising 36 PAAD samples and 36 normal samples. We found that 92.5% (1563) of the metabolism-related genes were expressed in both datasets (Figure 1A). The expression levels of metabolism-related genes were significantly higher than those of the other expressed genes in these datasets (p_value <0.01) (Figures 1B,C). Next, we identified the significantly DEGs between the PAAD and normal tissue data, which showed that 279 and 421 metabolism-related genes were significantly expressed in the TCGA and GSE15471 datasets, respectively (|log2(FC)| > 0.586 and p_value <0.05) (Figure 1D). KEGG enrichment analysis of both datasets showed that metabolic pathways are significantly disordered in tumor tissues of PC, including glutathione metabolism pathway (Figure 1E).

We performed WGCNA analysis on metabolic genes in GSE15471. We found that tumor tissue could be well differentiated from normal tissue using the created metabolic gene clusters (Figure 2A). We detected 7 modules in the network (Figure 2B). Genes in the seven modules showed a high correlation with each other (Figure 2C). The black module was also found to have the strongest correlation with pancreatic cancer patients, with a coefficient of correlation of 0.56 (p = 1 × 10⁻⁴) (Figure 2D). It shows that in the black module, gene expression was higher in most tumor samples than in normal tissue (Figure 2E). Therefore, this module was selected as the most clinically important screening module for further analysis.

We also subjected metabolic genes in the TCGA to WGCNA analysis. Since there were only four cases of normal samples of pancreatic cancer in the RNA-seq data of TCGA, we used other clinical indicators, such as sex, age, neoplasm histologic grade (NHG), pathologic N, and pathologic stage in the analysis (Figure 3A). Subsequently, we detected nine modules in the network (Figure 3B). Genes in the nine modules showed a high correlation with each other (Figure 3C). The results indicated that the brown module had the most significant correlation with pathologic_stage (correlation coefficient = −0.21, p = 0.008, Figure 3D). We examined the modules’ average gene significance related to pathologic stage, and the brown module was also found to have the most negetive association with the pathologic stage of PC (Figure 3E). Therefore, the brown module was selected as the most clinically important screening module for further analysis.

Then the black module of the GSE15471 data was used to construct a PPI network (Figure 4A), in which GPX2 and GSTP1 (encoding glutathione S-Transferase Pi 1) were identified as the hub genes. The black module was analyzed by KEGG, which showed that the glutathione pathway was significantly enriched (Figure 4B). Then we also used the brown module of TCGA to create a PPI network in which the hub genes was MAN1C1 (mannosidase alpha class 1C member 1) gene (Figure 4C). The brown module was then analyzed by KEGG, which showed that the glutathione pathway was also significantly enriched (Figure 4D). We then compared the expression levels of the metabolic enzyme genes and found that the gene expression levels of members of the glutathione pathway were significantly higher than those of other metabolic genes (Figure 4E). These results indicated that the glutathione pathway are significantly dysregulated in pancreatic cancer.

Since glutathione pathway has a great impact on the clinical evaluation of PC, we examined the expression level of genes involved in glutathione pathway. In the GSE15471 data, 14 genes in the glutathione pathway were differentially expressed (|log2(FC)| >0.586, p_value <0.05) (Figure 5A). In the TCGA, 12 glutathione pathway genes were found to be significantly differentially expressed (Figure 5B). We found that glutathione pathway genes, including GPX2, GSPT1, RRM2, are differentially expressed in both data (Figure 5C). We also verified the expression of GPX2, GSTP1 and RRM2 in glutathione pathway in our own RNA sequencing data and the results was consistent with that of public data we used (Figure 5D). Among them, GPX2, GSTP1 and RRM2 (Ribonucleotide Reductase Regulatory Subunit M2) showed different expression levels in different stages of PAAD (Figure 5E) (**p_value < 0.01). In stage II, the expression levels of these three genes were relatively high, suggesting they could be clinical indicators and providing a preliminary basis for later tumor heterogeneity. Then we performed the survival analysis of the glutathione pathway genes in patients with PC. The results indicated that patients with higher expression of glutathione pathway genes, including GSTP1 and RRM2, have shorter survial term (Figure 5F). Collectively, these data indicated that glutathione pathway genes may play critical roles in PC.

Since glutathione pathway genes, including GPX2, showed different expression levels in different stages of PAAD, we tried to predict immune heterogeneity based on gene expression profiles in PC at different stages. To comprehensively determine the cellular components of the tumor microenvironment across different stages of pancreatic cancer, we used the CIBERSORT in silico cytometry method (Newman et al., 2015) to evaluate 22 different immune cell types in 36 normal tissues and 36 tumor tissues in GSE15471, and 153 tumor tissues at different stages in TCGA. In the results for the GSE15471 data, PC appeared to promote an M2 macrophage gene signature compared with normal pancreatic tissue (Figure 6A). In the results of the TCGA analysis, M2 macrophages were more abundant in patients with stage IV disease compared with the patients at other stages (Figure 6B). We evaluated the correlation between the expression of GPX2 and CD163, the marker of M2, and it showed they are significantly correlated (Figure 6C). These results suggested that GPX2 may play important roles in cancer immune heterogeneity.