The RNA sequence transcriptome data of COAD patients and relevant clinical information were extracted from The Cancer Genome Atlas (TCGA; https://cancergenome.nih.gov/) database (24). Patients without clinical information were excluded. A total of 429 COAD samples and 37 normal samples were included in the study. The clinical characteristics of the patients are shown in Table 1 . The corresponding gene set of 210 glycosyltransferases was downloaded from GlycoGene DataBase (GGDB; https://acgg.asia/ggdb2/) (25). The GEO dataset (https://www.ncbi.nlm.nih.gov/geo/) GSE39582, which contains 579 samples, was used to verify the prognostic value of the risk model (26).

We used the limma software package in R version (3.6) to identify the differentially expressed GT genes (p value < 0.05 and | log2 (fold change) | > 1.5) between COAD and normal colorectal tissues (27). Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were conducted to analyze the enrichment pathways associated with differentially expressed genes (DEGs). Protein–protein interaction (PPI) network of DEGs was analyzed by the Search Tool for Interaction Genes (STRING) database (28). Single nucleotide variant (SNV) mutation frequency of DEGs was downloaded from GSCALite (http://bioinfo.life.hust.edu.cn/web/GSCALite/) (29). Copy number variation (CNV) mutation frequencies of GT DEGs in COAD patients were extracted from the TCGA database and analyzed by R software.

We used “limma” R package to identify GT-related lncRNAs (27). Pearson correlation analysis was conducted to detect the association between the expression levels of GT genes and lncRNAs. A total of 1,140 GT-related lncRNAs were identified based on the correlation coefficient >0.3 and p < 0.001.

To identify GT-related lncRNAs with prognostic value, we conducted the univariate and multivariate COX analysis. We first screened GT-related lncRNAs with potential survival impact on COAD patients based on the standard of p < 0.05 and subsequently incorporated those lncRNAs into LASSO Cox regression analysis and multivariate COX regression analysis (30). The following formula was used to calculate the risk score of each patient: risk score = (coefficient lncRNA₁ × expression of lncRNA₁) + (coefficient lncRNA₂ × expression of lncRNA₂) +…+ (coefficient lncRNA_n × expression of lncRNA_n). The COAD patients were divided into high- and low-risk groups based on the median value of risk score. R package “timeROC” was used to evaluate the accuracy of risk model (31). Additionally, we used the chi-square test to analyze the relationship between the model and clinical features in order to evaluate the prognostic role of this model.

Univariate Cox and multivariate Cox regression analyses were used to evaluate whether the risk score was an independent prognostic factor. Receiver operating characteristic (ROC) curve was made to evaluate the risk score and different clinical characteristics in predicting outcomes (32). Kaplan–Meier (KM) plot was used to analyze the outcome of high- and low-risk groups in subgroups according to clinical characteristics. A p-value <0.05 was considered statistically significant.

With “rms” R package, all significantly independent prognostic factors, including age, stage, and risk score were used to build up a nomogram for the 1-, 3-, and 5-year OS (33). The calibration plot was applied to evaluate whether the prediction outcome of the nomogram showed good consistency with practical application.

With curated gene set (kegg.v7.4.symbols.gmt) and GSEA software (https://www.gsea-msigdb.org/gsea/login.jsp), GSEA analysis was conducted in the COAD cohort to explore the significantly enriched pathways between the high- and low-risk groups based on the criterion: p<0.05 and false discovery rate (FDR) <0.25 (34, 35).

Immune Cell Abundance Identifier (ImmuCellAI) is an analytical tool that provides the quantitative infiltration of immune cells and predict the response of immune checkpoint blockade (ICB) therapy by using gene expression matrix data (36). We estimated the relative abundance of immune cell subtypes for each sample based on gene expression data. The abundance of immune cell subtypes and response to ICB therapy were compared between the high- and low-risk groups based on p < 0.05. The correlation between risk scores and immune infiltration was calculated by Pearson correlation analysis. Besides, the expression of immune-specific genes was compared between high- and low-risk groups.

The COAD cell lines (HCT116, DLD1, and HT-29) and normal human colon mucosal epithelial cell line NCM460 were acquired from the American Type Culture Collection (ATCC, USA). All cells were cultured in Roswell Park Memorial Institute (RPMI)-1640 medium (Gibco, CA, USA) with 10% fetal bovine serum (FBS, Gibco, CA, USA), 100 U/ml penicillin, and 100 mg/ml streptomycin. To evaluate the expression levels of GT-related lncRNAs, we used RNA Trizol reagent (Invitrogen, USA) to extract total cellular RNA, which were then reverse transcribed to cDNA using PrimeScript™ RT reagent Kit (Takara, Japan). Real-time fluorescent quantitative PCR was performed by using SYBR Green Master Mix (Yeasen, China). The β-actin mRNA was chosen as an endogenous control. The expression level of GT-related lncRNAs was analyzed using 2^−ΔΔCT. Each PCR reaction was performed in triplicate. The primer sequences are listed in Supplementary Table S1 .

The data analysis of this study was performed by R software (version 3.6; https://www.r-project.org/) and GraphPad Prism 8 software. One-way analysis of variance (ANOVA) was used to compare the expression level of GT in 429 COAD tissues and 37 normal tissues. The two-tailed t-test was performed to analyze continuous variables between the two subgroups, and the Chi-square test was conducted to evaluate categorical data. The DEGs co-expression network and GT-lncRNA-mRNA co-expression network were conducted by the Cytoscape software (version 3.6.0; https://cytoscape.org/). The Kaplan–Meier method was used to compare the OS of each group based on the median value, and log-rank tests were applied to evaluate the significance of differences. The RT-qPCR results were analyzed using one-way ANOVA, a p < 0.05 on both sides was considered statistically significant.

To identify the differentially expressed GTs and explore the biological functions of GTs in COAD, we used the TCGA database to analyze 429 COAD and 37 normal tissues. We uncovered 46 differentially expressed GTs (18 upregulated and 28 downregulated, Figures 1A, B ; Supplementary Table S2 ). The GO analyses revealed that these DEGs were most enriched in glycoprotein biosynthetic and metabolic processes ( Figures 1C, D ). KEGG analysis showed that 46 DEGs mainly involved in the mucin type O-glycan biosynthesis and glycosphingolipid biosynthesis including lacto and neolacto series, globo and isoglobo series, and ganglio series. We also clarified the relationship between these differentially expressed GTs via STRING database ( Figure 1E ). The incidence of single nucleotide variant (SNV) and copy number variations (CNVs) of differentially expressed GTs were analyzed in COAD. As shown in Figure 1F , SNV mutations were found in all 46 DEGs, among which WBSCR17 had the highest mutation frequency (31%), followed by UGGT2 (29%) and HS6ST3 (14%). We also found that CNV alterations were common in DEGs. Nearly half of the DEGs had copy number deletion, while the CNV amplification frequencies of PIGU, DPM1, PIGZ, ALG3, ST6GALNAC1, UGGT2, LFNG, HS6ST3, B4GALNT2, UST, HS3ST4, and GALNT6 were widespread ( Figure 1G ).

From the TCGA database, we identified 1,140 GT-related lncRNAs. We further combined these lncRNAs with COAD survival data in univariate COX regression analysis to identify 33 GT-related lncRNAs with significant prognostic value. As a result, 32 GT-related lncRNAs were found to be associated with increased cancer risk (p < 0.05 and HR > 1), while only 1 lncRNA, AC124067.4, was a protective factor for COAD (p < 0.05 and HR < 1) ( Table 2 ; Supplementary Figure S1A ). In correlation analysis, almost all prognostic lncRNAs have a weak to moderate correlation with other lncRNAs, whereas LINC00174 and AC026471.4, and LINC00174 and AC107375.1 had the strongest correlation (r =0.58, and 0.57, respectively) ( Supplementary Figure S1B ).

To better clarify the prognostic potential of GT-related lncRNAs, we conducted the LASSO and multivariate COX analysis on 33 GT-related lncRNAs ( Figures 2A, B ). Ultimately, seven lncRNAs, namely, LINC02381, AC002310.1, ZEB1AS1, AC020558.2, AC105219.1, MIR210HG, and AC009237.14, were identified to construct the prediction model ( Figure 2C ). The prognostic risk model was established based on the following formulation: risk score = (0.2552 × expression value of LINC02381) + (0.5370 × expression value of AC002310.1) + (0.4937 × expression value of ZEB1AS1) + (0.3940 × expression value of AC020558.2) + 0.1416 × expression value of AC105219.1) + (0.1187 × expression value of MIR210HG) + (0.1524 × expression value of AC009237.14). The risk score of patients was calculated, and COAD patients were divided into high- and low-risk groups according to the median risk score. The OS in high-risk patients was significantly lower than that in low-risk patients (p < 0.001, Figure 2D ). The distribution of the risk score and survival status are shown in Figure 2E . The risk score had a good prognostic predictive ability with the area under the ROC curve (AUC) values at 1, 3, and 5 years of 0.726, 0.748, and 0.844, respectively ( Figure 2F ). Furthermore, we investigated the prognostic implication of the expression of lncRNAs. The OS of patients in the high-expression group of seven lncRNAs was significantly lower than that in the low-expression group (all p < 0.05, Figures 3A–G ). We conducted a lncRNA–mRNA co-expression network to visualize the relationship among seven prognostic lncRNAs and GTs ( Figure 3H ). According to the risk score, a Sanky diagram was produced to exhibit the association among GTs, GT-related lncRNA, and outcome types of lncRNAs, indicating a single lncRNA corresponding to one or more mRNA, and these seven GT-related lncRNAs were all risk factors for a poor prognosis ( Figure 3I ). To verify the prognostic value of the risk model, we acquired the GSE39582 dataset from GEO database and then extracted the expression levels of seven GT-related lncRNAs and the survival data of patients. The results showed that patients with low-risk score had a better prognosis compared to those with high-risk score ( Supplementary Figure S2 ).

In the next step, we explored the relationship between the risk model and clinical variables. Significant differences were found between high- and low-risk groups in the pathological stage, T stage, M stage, and N stage, whereas the two groups showed no significant differences in age and gender ( Figures 4A–D ; Supplementary Figure S3 ). The high-risk group underwent significant upstaging compared with those in the low-risk group (52.2% vs. 2.1%, p=0.004). With the respect to the distribution of the American Joint Committee on Cancer (AJCC) stage, the high-risk group also comprised much poorer AJCC stage (T3/4, M1, and N1/2) than the low-risk group (85.1% vs. 74.7%, p = 0.078; 14.9% vs. 9.9%, p = 0.021; 51.9% vs. 28.4%, p < 0.001, respectively). Next, we performed an in-depth analysis of the correlation between seven lncRNAs and clinical variables ( Figures 4E–H ). In the pathological stage, the expression level of four lncRNAs (ZEB-AS1, AC105219.1, MIR210HG, and AC009237.14) was gradually elevated as the tumor stage progressed. For T stage, only one lncRNA was identified with significant differences across subgroups. Regarding M and N stage, three and five lncRNAs, respectively, were found with different expressions among different subgroups. Interestingly, we noticed that the expression level of ZEB1-AS1 was significantly increased with the progression of all stages (T, M, N, and S), suggesting that this lncRNA might be a core prognostic factor in COAD.

We used univariate and multivariate COX regression analyses to evaluate whether the risk model had an independent prognosis value. The results of univariate and multivariate analysis revealed that age (HR=1.622 and 2.519, respectively; 95%CI, 0.930–2.829 and 1.382–4.592, respectively; p = 0.088 and 0.003, respectively), stage (HR=3.163 and 2.383, respectively; 95%CI, 1.939–5.158 and 1.176–4.831, respectively; p < 0.001 and 0.016, respectively), and risk score (HR=3.064 and 2.271, respectively; 95%CI, 1.843–5.092 and 1.341–3.847, respectively; p < 0.001 and 0.002, respectively) can be independent prognostic factors ( Figures 5A, B ). Time-dependent ROC curves were conducted to evaluate the prognostic and predictive ability of clinical features and risk score ( Figures 5C–E ). Results showed that the AUC value of risk scores at 1, 3, and 5 years was 0.694, 0.741, and 0.834, respectively. More importantly, the risk score had higher precision compared with other predictors. As we all know, age, gender, and clinical stage are identified as independent prognostic characteristics of COAD patients. Then, patients were stratified according to age (≤60 and >60), gender (female and male), stage (I/II and III/IV), T stage (T1/2 and T3/4), M stage (M0 and M1), and N stage (N0 and N1/2) ( Figures 5F–H ; Supplementary Figure S4 ). Among the subgroups, patients in the high-risk group exhibited much poorer survivability than those in the low-risk group.

To evaluate the clinical characters and risk model for colorectal cancer prognosis, we integrated age, stage, and risk score to build a nomogram ( Figure 6A ). Additionally, a time-dependent ROC analysis was performed to evaluate whether risk score and clinical traits were capable for survival prediction ( Figure 6B ). The AUC value of risk score in time of 5-year follow-up was 0.726, 0.703, 0.729, 0.757, and 0.825, respectively. The AUC value of stage in time of 5-year follow-up was 0.755, 0.720, 0.758, 0.704, and 0.656, respectively. For age, the average AUC was below 0.6 during a follow-up of 5 years. These results showed that risk score exhibited much more powerful capacity of survival prediction compared with other clinical traits. Besides, calibration plots in 1, 3, and 5 years are shown in Figures 6C–E , and it revealed that the observed vs. predicted rates of 1-, 3-, and 5-year overall survival were in perfect concordance.

To further explore the abnormally activated signal pathways between the high- and low-risk groups, we performed GSEA enrichment analysis. The results showed that the low-risk group was involved in several important pathways, including citrate cycle TCA cycle, fatty acid metabolism, glutathione metabolism, glycolysis gluconeogenesis, oxidative phosphorylation, P53 signaling pathway, proteasome, pyrimidine metabolism, ribosome, and valine leucine and isoleucine degradation ( Figure 7A ). Given the important role of tumor immune microenvironment (TIM), we looked into the expression levels of immune-related genes between the high- and low-risk score groups. The result showed that CD200R1, ARG1, CD160, ADORA2A, IL4, VEGFA, TNFRSF4, TNFRSF14, TNFRSF25, TBX2, VSIR, TGFB1, and NOS3 are increased, while CXCL8, NOS2, EZH2, and HHLA2 are decreased in the high-risk group compared with the low-risk group ( Figure 7B ). In addition, Pearson correlation analysis was conducted to analyzed the correlation between expression levels of immune-related genes and seven GT-related lncRNAs ( Supplementary Table S3 ). The results turn out that AC020558.2, AC105219.1, and MIR210HG were positively correlated with TNFRSF25 and TNFRSF14, which are immune suppression genes, while LINC02381, AC002310.1, and ZEB1-AS1 were negatively associated with immune response gene NOS2 and HHLA2. ImmuCellAI is an online tool to estimate the response of immune checkpoint blockade (ICB) therapy and the infiltration of immune cells based on the gene expression dataset (36). We used ImmuCellAI to predict the response of ICB therapy between the high- and low-risk groups. The result revealed that patients in the low-risk group tend to have a better response to ICB therapy than those in the high-risk group, even though the difference was not significant (p=0.0723, Figure 7C ). The abundance of immune cells with significant differences between the high- and low-risk groups is shown in Figure 7D . The levels of B cell and CD4+ T cell in the high-risk group were higher than those in the low-risk group, while the levels of T-helper 1 cells (Th1), neutrophil, and natural killer T (NKT) in the high-risk group were lower than those in the low-risk group (p<0.05). To determine the correlation between the risk score and tumor-infiltrating immune cells, we found that the risk score was positively correlated with B cell and CD4+ T cell and negatively correlated with Th1 and neutrophil (p<0.05) ( Figure 7E ). Taken together, these data suggest that our risk model can distinguish different characteristics of tumor immune cells in COAD.

To further verify our results, we used RT-qPCR assay to analyze the expression level of seven GT-related lncRNAs in COAD cell lines ( Figures 8A–G ). Three human COAD cell lines, namely, HCT116, DLD1, and HT-29, were used, and NCM460, a normal human colon mucosal epithelial cell line, was used as control. Almost all seven prognostic GT-related lncRNAs were upregulated in COAD cell lines, indicating that these GT-related lncRNAs play key roles during progression of COAD. The results were generally consistent with our work.