Three datasets were involved in this study. Clinical information and survival data of COAD patients from The Cancer Genome Atlas (TCGA-COAD) were downloaded from the UCSC XENA database (https://xenabrowser.net/), as well as the gene expression matrix (HTSeq-FPKM), somatic mutation, copy number variation (CNV, gene-level), and methylation (Illumina Human Methylation 450) sequencing results. Clinical information and survival data of colon cancer patients in the GSE39582 dataset were retrieved from the Gene Expression Omnibus (GEO) dataset (http://www.ncbi.nlm.nih.gov/geo). Microarray RAW “CEL” data files of colon cancer patients in the GSE39582 dataset were also downloaded, and normalization was conducted through the Robust Multiarray Average (RMA) method in R package “Affy.” A single-cell transcriptomic dataset of COAD patients in GSE146771 (Smart-seq2) was also obtained from the GEO database. The KEGG_GLYCEROLIPID_METABOLISM gene set and 50 cancer hallmark gene sets were downloaded from the Molecular Signature Database v. 7.4 (MSigDB). A total of 49 glycerolipid metabolism-related genes were obtained from the KEGG_GLYCEROLIPID_METABOLISM pathway.

The following preprocessing steps were applied to both TCGA-COAD and GSE39582 datasets. (1) Patients who had no clinical information were excluded from the datasets. (2) Patients whose survival time was 0 or survival status was unknown were excluded from the analysis. (3) The type of gene ID in the gene expression matrix of both datasets was synchronized to Gene Symbol. When multiple gene IDs were synchronized to one Gene Symbol, the median value was selected to represent the expression level.

The 49 glycerolipid metabolism-related genes were mapped to the gene expression matrix and gene methylation matrix. The gene expression level and methylation level of these genes were compared between tumor and normal tissues by the Mann–Whitney test, and P-value <0.05 was considered significant. The CNV amplification and deletion frequency, as well as the somatic single-nucleotide variation, of the glycerolipid metabolism-related genes were calculated from relative sequence results of tumor samples.

Single-sample gene set enrichment analysis (ssGSEA) is an unsupervised single-sample gene set enrichment analysis method (11). ssGSEA calculates the enrichment score of each sample by the cumulative distribution function difference between gene expression ranks inside the settled gene set and outside. The enrichment score was further normalized by the range of total gene expression matrix for all genes and samples (11). The ssGSEA analysis method included in R package “GSVA” was used to calculate the KEGG_GLYCEROLIPID_METABOLISM gene set enrichment score (GLMS) based on the gene expression matrix of TCGA-COAD and GSE39582 (12). Samples were categorized into high- and low-GLMS groups based on the median of GLMS for both TCGA-COAD and GSE39582 datasets.

The consensus molecular subtypes (CMSs) of TCGA-COAD datasets were estimated via the R package “CMScaller” (13). The clinical characteristics of TCGA-COAD datasets such as age, sex, BMI (classified into high and low groups based on median), stage, microsatellite instability, venous invasion, and tumor location were included for assessing the association between the GLMS and clinical characteristics, together with the CMS subtypes.

The ssGSEA analysis method was also applied to calculate the 50-cancer hallmark gene set enrichment score based on the gene expression matrix of TCGA-COAD dataset. Spearman correlation analysis was used to quantify the correlation between the GLMS and 50 cancer hallmark gene sets based on their ssGSEA scores of TCGA-COAD samples. The correlation between every single glycerolipid metabolism-related gene in the KEGG gene set and 50 cancer hallmark gene sets was also calculated via R package “WGCNA” corAndPvalue function based on gene expression and the hallmark ssGSEA score (14).

The proportions of 64 cell types in the tumor microenvironment of TCGA-COAD tumor samples were estimated by R package “xCell” (15). The expressions of 14 immune checkpoints were compared between high- and low-GLMS groups including CD274, CD276, CD40, CTLA4, HAVCR2, ICOS, IDO1, LAG3, PDCD1, TIGIT, TNFRSF18, TNFRSF4, TNFRSF9, and VTCN1.

The pRRophetic algorithm was originally developed for drug response prediction based on gene expression microarray data and the half-maximal inhibitory concentration (IC50) of 138-drug response data in almost 700 cell lines from the Cancer Genome Project (CGP) database (16, 17). The drug sensitivity of those drugs in TCGA-COAD dataset, as indicated by IC50 prediction, was calculated based on the R package “pRRophetic” and gene expression matrix of these samples. The difference of estimated drug sensitivity between high- and low-GLMS groups was evaluated for every drug via the Wilcoxon test, and P-value <0.05 was considered significant. Subsequently, the correlation between every single gene in the KEGG_GLYCORLIPID_METABOLISM gene set and drug sensitivity was calculated via R package “WGCNA” corAndPvalue function based on expression and IC50 (14). Furthermore, potential drug response-related genes were selected according to P-value and correlation coefficient (P-value< 0.05, |cor| >0.5).

GLMS of single-cell transcriptomic dataset GSE146771 was calculated by the ssGSEA method as well. The cell annotation information corresponding to this dataset was downloaded from the Tumor Immune Single-cell Hub (TISCH) database (http://tisch.comp-genomics.org). The GLMS of every single cell was shown by uniform manifold approximation and projection (UMAP) plots at two resolutions: “Celltype malignancy” and “Celltype major lineage.” The Kruskal–Wallis test was applied to evaluate whether GLMS varied among different cell types, and P-value <0.05 was considered significant. Finally, the cell subpopulations on which glycerolipid metabolism affected were estimated by the log2FoldChange (log2FC) value of each cell types via “one vs. rest” calculation.

Data preprocessing and R package analyses were completed in R (version 4.1). Statistical analysis and figures were completed in GraphPad Prism (version 9). Continuous variables were displayed as median (range), and categorical variables were displayed as number (percentage). The Wilcoxon test was applied to calculate the significance of the difference in continuous variables between two groups, and the chi-square test or Fisher’s exact test was used to calculate the significance of the difference in categorical variables between groups. Kaplan–Meier analysis was performed to estimate the overall survival (OS) rates, and the difference of survival curves between high- and low-GLMS groups was compared by the log-rank test through R package “survival.” Univariable and multivariable Cox proportional hazards regression analyses were conducted via R package “survival” coxph function to determine whether the GLMS could be an independent prognostic factor for colon cancer. Statistical significance was considered as P-value <0.05.

Based on the 430 cancer samples and 39 normal samples in TCGA-COAD dataset, the expression levels of 43 glycerolipid metabolism-related genes were available out of 49 genes from the KEGG_GLYCEROLIPID_METABOLISM gene set, while the methylation information was provided for 42 genes. The baseline characteristics of the samples in TCGA-COAD dataset are described in Supplementary Table . Differences in expression and methylation levels were compared between tumor and normal tissues ( Figure 1A ), which identified nine genes with a negative correlation between expression and methylation levels, including CEL, DGAT2, LCLAT1, and LIPG (hypomethylated and upregulated in tumor tissues), and AGPAT3, AKR1B1, ALDH2, DGKD, and MBOAT1 (hypermethylated and downregulated in tumor tissues), indicating an inhibitory function of methylation in the mRNA expression of those genes. However, the methylation of six genes was positively correlated with expression levels, including DGKH and DGKZ (hypermethylated and upregulated in tumors), and AGPAT1, DGKB, DGKQ, and LIPC (hypomethylated and downregulated in tumors), suggesting the potential regulatory mechanisms of mRNA expression in those genes beyond methylation. The copy number variations and gene mutations were also analyzed based on cancer samples. Copy number amplification appeared in AGPAT1, AKR1B1, DGKB, DGAT2, and DGKH, while AGPAT3, ALDH2, DGKD, MBOAT1, and LIPC exhibited significant copy number deletions ( Figure 1 ). A total of 120 (27.9%) tumor samples harbored mutations in glycerolipid metabolism-related genes, with two genes identified as high frequently mutated, including DGKB (15%) and DGKD (14%) ( Figure 1 ).

Each sample in TCGA-COAD dataset was assigned with a glycerolipid metabolism-related gene enrichment score (GLMS) according to the single-sample gene set enrichment analysis (ssGSEA) algorithm and divided into high- and low-GLMS groups by the cutoff of the median value. We explored whether the GLMS could predict the prognosis of colon cancer in TCGA-COAD dataset. Kaplan–Meier analysis showed that the patients with higher GLMS exhibited improved overall survival (OS) compared with those with lower GLMS (log-rank P = 0.03; HR, 0.63; 95% CI, 0.42–0.94, Figure 2 ). Univariable analysis was also conducted in other clinical factors, including age, sex, stage, tumor location, MSI status, and CMS, among which only clinical stage and GLMS exhibited statistical significance (P < 0.05, Figure 2 ). These clinical factors along with GLMS were included in the multivariable Cox regression, and the association between GLMS and OS remained significant (P = 0.04; HR, 0.64; 95% CI, 0.42–0.98, Figure 2 ), suggesting that the GLMS was a potential independent prognostic factor of colon cancer.

To further verify the robustness of the GLMS in predicting the prognosis of colon cancer, another dataset (GSE39582) was used as a validation cohort. A total of 562 tumor samples were included in the GSE39582 dataset, with the baselined characteristics described in Supplementary Table . With the same algorithm of GLMS, patients in the high-GLMS group also had a better OS (log-rank P < 0.001; HR, 0.57, 95% CI, 0.43–0.76, Figure 2 ). Similar results were also observed in the univariable ( Figure 2 ) and multivariable Cox regression ( Figure 2 ). Univariable analysis has identified tumor stage and GLMS that were associated with OS ( Figure 2 ). The GLMS was further confirmed to be an independent prognostic factor of OS in the multivariable analysis (P < 0.001; HR, 0.60; 95% CI, 0.45–0.80, Figure 2 ). Taken together, these results suggested that the GLMS was an independent prognostic factor for colon cancer.

The potential association between GLMS and clinical features was further analyzed including age, sex, body mass index (BMI), tumor stage by American Joint Committee on Cancer (AJCC) TNM Classification, microsatellite status, vascular invasion status, tumor location, and consensus molecular subtypes (CMS) ( Figure 3 ). Results showed that left colonic carcinomas exhibited significantly higher GLMS compared with right ones (P = 0.008, Figure 3 ), suggesting the capacity of the GLMS in differentiating the tumor location. Besides, CMS subtypes among tumors have revealed significant differences in GLMS distribution (P < 0.05, Figure 3 ), while patients in different groups by age, sex, BMI, tumor stage, microsatellite instability status, and venous invasion status were observed to have a similar distribution of the GLMS ( Figure 3C ).

To explore the potential pathways glycerolipid metabolism was involved in, the correlation between the GLMS and cancer hallmarks was analyzed. A significant positive association was observed between GLMS and pathways including bile acid metabolism, xenobiotic metabolism, and peroxisome, while pathways such as interferon gamma response, allograft rejection, apoptosis, and inflammatory response were found to be negatively correlated with the GLMS ( Figure 4 ). The 42 genes involved in the GLMS were analyzed as well. Samples with an upregulated expression of AGPAT4, AKR1B1, and DGKI showed positive enrichment in pathways including cholesterol homeostasis, bile acid metabolism, inflammatory response, and NOTCH signaling ( Figure 4 ), while expressions of AGPAT4, AKR1B1, DGKI, and MGLL were negatively associated with pathways such as protein secretion and IL2-STAT5 signaling ( Figure 4 ), suggesting the core genes and pathways involved in glycerolipid metabolism.

The immune microenvironment was revealed by utilizing R package “xCell” to estimate the infiltration of the proportion of different immune cells, which was further compared between high- and low-GLMS groups. As presented in Figure 4 , a total of 33 immune cells or stromal cells exhibited a significant infiltration difference between samples with high- and low-GLMS. Most of the immune-promoting cells, including B cells, macrophages, activated dendritic cells (aDC), classical dendritic cells (cDC), and immature dendritic cells (iDC), accounted for higher proportions in the low-GLMS group, while cells negatively regulating immune response such as NKT cells were significantly infiltrated in high-GLMS group samples.

Furthermore, the expression of 14 genes involved in immune checkpoints was compared between the two groups. As a result, a total of 12 immune checkpoint molecules were upregulated in low-GLMS group samples ( Figure 4 ), suggesting a potential benefit of immune therapies in patients with low GLMS.

Considering that lower GLMS was associated with poor prognosis, we conducted further analysis exploring the relationship between GLMS and drug sensitivity. Based on the Genomics of Drug Sensitivity in Cancer (GDSC) database, the half-maximal inhibitory concentration (IC50) information was acquired to predict the treatment response. Among the 96 drugs which exhibited a significant response difference between groups with low and high GLMS, those that were currently used in the treatment of colon cancer were paid extra attention. Low-GLMS group samples were observed to be more sensitive to chemotherapeutics including cisplatin ( Figure 5 ), gemcitabine ( Figure 5 ), and camptothecin ( Figure 5 ), which had analogues of oxaliplatin, capecitabine, and irinotecan, respectively. Besides, ponatinib ( Figure 5 ) as an anti-angiogenesis inhibitor targeting vascular endothelial growth factor receptors (VEGFRs) and epidermal growth factor receptor (EGFR) tyrosine kinase inhibitors (TKIs) including gefitinib ( Figure 5 ) and afatinib ( Figure 5 ) also exhibited significant increased sensitivity associated with low-GLMS group samples, suggesting potential benefits utilizing these therapies in patients with lower GLMS who had poorer survival.

In addition, further analysis also explored the association between IC50 of the drugs and 42-glycerolipid metabolism-related gene expression. According to P-value <0.05 and the absolute value of correlation >0.5, a total of 14 genes were identified including AKR1B1, DGKI, AGPAT2, GPAM, LCLAT, DGKZ, MGLL, DGKQ, AGPAT4, DGAT1, AKR1A1, DGKD, DGKH, and AGK ( Figure 5 ), suggesting that these genes might play a crucial role in drug response.

The single-cell sequencing dataset GSE146771 was introduced in this current study to investigate the origination of glycerolipid metabolic disorders. Annotations of cell samples were obtained from the TISCH database to identify the subtype of the cells, which were presented by a rough classification of “malignant-immune-stromal cells” ( Figure 6 ) and detailed subsets ( Figure 6 ). Each cell sample was assigned with a GLMS according to the same ssGSEA algorithm ( Figure 6 ). Results showed that the scores were generally higher in malignant cells and stromal cells, as well as in CD4Tconv, CD8Tex, and plasma cells when it came to the sophisticated categories. Further statistical analysis revealed a significant difference of the score distribution among various cell types (P < 2.2e-16) ( Figure 6 ). Enrichment of the scores was analyzed based on the log2FC value of the scores in each cell subset, which showed significant enrichment in malignant cells and Tprolif cells ( Figure 6 ), suggesting potential targets and originations of the glycerolipid metabolic disorders.