Glycerophosphodiester phosphodiesterase 1 (GDE1) acts as a potential tumor suppressor and is a novel therapeutic target for non-mucin-producing colon adenocarcinoma

Data collection

We downloaded 289 COAD tumor samples (transcripts per million, TPM) from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/repository) and 41 peritumoral normal tissues of COAD (TPM) from the UCSC Xena website (https://xena.ucsc.edu/). The data were analyzed with a unifying pipeline, which used a CutAdapt trimming adapter, STAR for alignment, and RSEM and Kallisto as quantifiers. The COAD clinical data were downloaded from the TCGA database using the ‘cgdsr’ package (Null, Team & Null, 2009; Jacobsen, 2015).

Differentially expressed gene (DEG) screening

We downloaded metabolism-associated genes from RECON3, the most comprehensive human metabolic network model which contains 13543 reaction, 4140 metabolites, and 3288 metabolic genes (Brunk et al., 2018), and identified 3,267 genes that corresponded to the TPM data. We removed lowly expressed genes (TPM <1) from the COAD and normal samples, and thus retained 2,747 genes. We used principal component analysis (PCA). To obtain significantly changed genes in COAD, we used the ‘aov’ function to analyze the DEGs between the COAD and normal samples. We obtained 2,114 significant DEGs in COAD with a p < 0.01 cutoff.

Co-expression network construction by WGCNA

WGCNA can be used to analyze global gene expression profiling and identify highly co-expressed genes. We installed WGCNA with Bioconductor (http://bioconductor.org/biocLite.R) for co-expression analysis of metabolism-related genes in COAD. The soft threshold method was used for correlation analysis of the expression profiles. We used average linkage hierarchical clustering to group transcripts based on topological overlap dissimilarity in network connection strengths. We set the minimum gene number of each module to 30 and cutHeight parameter to 0.25 to merge similar modules.

Identification of significant histological type modules

To incorporate external information into the co-expression network, WGCNA uses gene significance measures GS (minus log of a p-value). A gene significance measure indicates (the absolute value of) the correlation between individual gene in a given module and the trait. The average GS across the module genes was used to characterize the correlation between modules and clinical traits. For each module, we also define a quantitative measure of module membership MM as the correlation of the module eigengene MEs (defined as the first principal component of a given module) and the gene expression profile. The relationships between MEs and clinical traits were analyzed, and the relevant modules were identified.

Functional and pathway enrichment analysis

We used the Database for Annotation Visualization and Integrated Discovery (DAVID) (http://david.abcc.ncifcrf.gov/) for gene functional annotation to determine biological functions. The genes in the significant histological type modules were uploaded in DAVID for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. We set a cutoff of adjusted p < 0.05 for significant biological processes and significant pathways.

Survival analysis

The prognosis profiles of the genes in significant histological type modules were analyzed using the R survival package from Bioconductor (https://www.bioconductor.org). Samples were classified as two subtypes based on gene expression. Samples with larger than median gene expression levels were defined as the “high expression” subtype; all other samples were defined as the “low expression” subtype. The log-rank p-value and survival visualization were performed using the R “survminer” package.

Metabolism-related DEGs in COAD

We obtained gene expression TPM data from 289 COAD samples and 41 adjacent normal tissue samples, which included 60,498 genes. The QC analyzing result showed that the median and the quartile log 2 TPM value are similar between adjacent normal samples and cancer samples, indicating they are comparable (Fig. S1). We used PCA to identify the possible outlier samples, and found that all the samples are located in the confidential ellipse (Fig. 1A). We thus remained all of the samples to further analysis. To determine if metabolism-related genes influenced COAD, we first obtained 3,288 metabolism genes in RECON3, a most comprehensive human metabolic network model (Brunk et al., 2018). We removed lowly expressed genes (TPM <1) from the COAD and normal samples, and thus acquired 2,722 genes. We then determined DEGs in COAD compared with the normal tissue samples using ANOVA in R. We identified 2,114 DEGs with a cutoff of p < 0.01 (Fig. 1B), which included 1224 highly expressed genes and 890 lowly expressed genes (Table S1). The profiles of the 2,114 DEGs are shown with heatmap in Fig. 1C.

Metabolism gene co-expression network constructed by WGCNA

To explore the hub metabolism-related genes in COAD, we analyzed the co-expression network of the 2,114 DEGs using WGCNA. The power value is the most critical parameter that can affect the independence and average connectivity degree of co-expression modules. We screened network topology using different soft thresholding powers for later analysis (Figs. 2A, 2B). The scale-free R2 equated to 0.90 at the power value of six. Therefore, the power value of six was used to construct the co-expression network with WGCNA, which was based on the hierarchical clustering of the calculated dissimilarities. The distance between clusters was defined as the average distance between all inter-cluster pairs. With a cutHeight of 0.25 (Fig. 2C) and minimum gene number of 30, we merged the modules, with 12 then obtained (Fig. 2D). We randomly selected 500 genes in the hierarchical clustering results, with the rows and columns representing each gene in each module. The deep red color indicated high topological overlap, suggesting that the co-expressed genes had similarity at the network topology level (Fig. 2E, Fig. S2).

Identifying genes in red module associated with histological typing

To investigate the modules associated with clinical traits, we analyzed the module trait relationship. We first convert clinical trait data into variables which expressed in quantitative form, next we quantify gene relationship to clinical traits and important modules by using the gene significance and module membership (MM) measures. The results showed that the red module was significantly but negatively correlated with COAD clinical histological typing (r =  − 0.33, p = 4e−08) (Fig. 3A). We also tested the relevance between each module and COAD clinical traits. The red module demonstrated much higher module significance than any other module, suggesting its stronger connection with COAD histological typing (Fig. 3B), viz., non-mucin-producing colon adenocarcinoma (N-COAD) and mucin-producing colon adenocarcinoma (M-COAD). Red module membership was also significantly correlated with gene significance (Fig. 3C). Taken together, our results indicated that the red module was significantly correlated to COAD histological typing.

The red module contained 199 genes (Table S2). To clarify the function of these genes, GO and KEGG pathway analyses were conducted using DAVID (v6.8). Genes in the red module were significantly enriched in ATP-binding (Benjamini–Hochberg method adjusted p-value = 4.83e−13) molecular function (Fig. 3D) and in N-Glycan biosynthesis pathways (Benjamini–Hochberg method adjusted p-value = 6.23e−04) (Fig. 3E).

Significant association of GDE1 gene with histological typing and prognosis

To identify which genes in the red module were associated with different histological typing and prognosis of COAD, we performed survival analysis (Table 1). We found GDE1 to be the most significant prognosis-related gene. This gene demonstrated lower expression in COAD at the RNA level (Fig. 4A) and its expression level was significantly associated with survival time (p = 0.001; Table 1). Interestingly, we also observed that GDE1 expression in N-COAD was significantly lower than that in normal tissues (p < 1e−6) but higher than that in M-COAD (p = 0.003) (Fig. 4A). Furthermore, GDE1 was the most significant gene correlated with N-COAD subtype prognosis (Table 1).

Although different histological subtypes of COAD showed similar prognosis profiles (Fig. 4B), we found that the correlation of GDE1 with prognosis was distinct in the N-COAD and M-COAD subtypes (Figs. 4C and 4D). Specifically, GDE1 was significantly associated with N-COAD prognosis (p = 0.0017), and N-COAD with high GDE1 expression demonstrated significantly longer survival time (Fig. 4C). In M-COAD, however, GDE1 was not correlated with prognosis (p = 0.92) (Fig. 4D). Taken together, our results showed that GDE1 was specially associated with N-COAD subtype prognosis, suggesting that it could be a potential novel cancer target for N-COAD.