Oncomine online database (www.oncomine.org) (Rhodes et al., 2004; Rhodes et al., 2007) and TCGA cohort (https://portal.gdc.cancer.gov) (Tomczak et al., 2015) were used to analyze the expression differences of SLC2A1 in varied cancers. Student’s t-test was used to compare the expression difference of SLC2A1 in tumor and control samples in Oncomine database, and data with fold change >2 and P < 1E-8 were selected for display. In the TCGA cohort, the CRC cohort used included the colon adenocarcinoma (COAD) and rectum adenocarcinoma (READ) datasets. We also analyzed CRC cohort from the TCGA and GEO (www.ncbi.nlm.nih.gov/geo; GSE9348, GSE110223, and GSE110224) (Barrett et al., 2012) databases to investigate differences in SLC2A1 expression level between tumor and control samples. By analyzing the clinical datasets of TCGA CRC cohort, the relation between SLC2A1 expression and clinicopathological characteristics of patients with CRC was investigated, and the diagnostic value of SLC2A1 for CRC was evaluated by ROC curve. Eventually, we certified the difference in SLC2A1 expression between CRC and normal samples by used qRT-PCR and IHC staining. Refer to our previous work (Liu et al., 2021a) for specific steps, including specific methods of qRT-PCR and IHC staining, as detailed in the Supplementary Material.

The TCGA CRC cohort was analyzed utilizing the R statistical computing language (version 3.6.3) to investigate co-expressed genes related with SLC2A1 expression. The Pearson’s correlation coefficient was used for statistical analysis, and correlations were considered as significant when |cor| > 0.3 and p < 0.05. Ggplot2 software package of R language (https://ggplot2.tidyverse.org) was utilized to draw volcano and heat map for demonstration. Top 200 genes that were positively correlated with SLC2A1 expression were selected as hub genes. The Genome Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg) pathway investigation of hub genes was carried by using the clusterProfiler (Yu et al., 2012) software package of R language, and the data were analyzed by using the ggplot2 software package.

To study the potential function of SLC2A1 in CRC, we divided the samples in the CRC cohort into high and low groups based on the SLC2A1 median expression in CRC and carried Gene Set Enrichment Analysis (GSEA) (www.gsea-msigdb.org/gsea/index.jsp) (Subramanian et al., 2005) to research whether genes in the two groups were enriched with meaningful biological functions. The annotated gene set c2.cp.v7.2.symbols.gmt (Curated) was selected as the reference gene set. The FDR (q-value) < 0.25 and p < 0.05 were considered statistically significant.

To investigate the possible regulatory mechanism of SLC2A1 in the regulation of CRC immune cells, we utilized the TIMER database (https://cistrome.shinyapps.io/timer) to evaluate the relationship between SLC2A1 expression and immunofiltration cells in TCGA CRC samples. Immunofiltration cells include B cell, neutrophil, CD4+ T cell, CD8+ T cell, macrophage, and dendritic cell (DC). We utilized the somatic copy number alteration module of the TIMER tool to correlate genetic copy number variation (CNV) of SLC2A1 with the relative proportion of immune infiltrating cells. The GSVA (Hänzelmann et al., 2013) software package of R language was utilized to analyze the relative abundance of 24 immune cells in CRC samples with high and low SLC2A1 expression, and the specific algorithm was ssGSEA. Furthermore, we analyzed the relation between SLC2A1 and immunofiltration cell marker genes in CRC samples utilizing Tumor IMmune Estimation Resource (TIMER), Gene Expression Profiling Interactive Analysis (GEPIA), and TCGA databases. Immunofiltration marker genes refer to the previous study, in which immune markers of different immune cells are described (Liu et al., 2021b).

The association between the expression level of SLC2A1 and the m6A relative genes, including YTHDF1, WTAP, RBM15, FTO, ALKBH5, ZC3H13, YTHDF3, HNRNPC, YTHDC2, METTL14, METTL3, IGF2BP3, IGF2BP2, RBMX, RBM15B, IGF2BP1, YTHDC1, HNRNPA2B1, VIRMA, and YTHDF2 (Li et al., 2019), in GSE110224 and TCGA CRC cohort was analyzed utilizing the R statistical computing language. In the TCGA cohort, the CRC cohort used included the COAD and READ datasets. The ratio of m6A relative genes in CRC samples with high and low SLC2A1 expression was analyzed by R statistical computing language. The Kaplan–Meier curve indicated the relation between the expression level of relative genes and the prognosis of patients with CRC. Ggplot2 software package was utilized for visual analysis of the cohort.

The R statistical computing language was utilized to analyze the association between the expression level of SLC2A1 and ferroptosis relative genes in the GSE110224 and TCGA CRC cohort, including CDKN1A, HSPA5, EMC2, SLC7A11, NFE2L2, MT1G, HSPB1, GPX4, FANCD2, CISD1, FDFT1, SLC1A5, SAT1, TFRC, RPL8, NCOA4, LPCAT3, GLS2, DPP4, CS, CARS, ATP5MC3, ALOX15, ACSL4, and AIFM2 (Doll et al., 2019; Liu et al., 2020). In the TCGA cohort, the CRC cohort used included the COAD and READ datasets. R statistical computing language was utilized to analyze the ratio of ferroptosis relative genes in CRC samples with high and low SLC2A1 expression. The Kaplan–Meier curve indicated the relation between the expression level of relative genes and the prognosis of patients with CRC. Ggplot2 software package was used for visualizing analysis of the cohort.

The mirDIP (www.ophid.utoronto.ca/mirDIP/) (Tokar et al., 2018), micorT CDS (www.microrna.gr/microT-CDS) (Paraskevopoulou et al., 2013), and miRNet (https://www.mirnet.ca/miRNet/home.xhtml) (Chang et al., 2020) tools were used to predict the miRNAs of target SLC2A1, and the correlation between SLC2A1 and these miRNAs was analyzed in the TCGA CRC cohort to screen negatively correlated miRNAs as target miRNAs. miRNet and starBase (www.starbase.sysu.edu.cn/starbase2/index.php) (Li et al., 2014) tools were used to predict the long non-coding RNA (lncRNAs) of target microRNA (miRNAs), and the correlation between the target miRNAs and these lncRNAs was analyzed in the TCGA CRC dataset, and the negatively correlated lncRNAs were selected as target lncRNAs. To further narrow the prediction range, the expression level of target lncRNA in CRC and its correlation with SLC2A1 were further analyzed according to the ceRNA theory. The igraph software package of R language (http://igraph.org) is used to analyze the data visually.

All statistical analysis was conducted through Xiantao platform (www.xiantao.love). Xiantao platform is a database integrating TCGA tumor microarray data, which is mainly used for gene expression analysis, correlation analysis, enrichment analysis, interactive network analysis, clinical significance analysis, and related plotting.

The differences in SLC2A1 mRNA expression level between CRC and control group were analyzed by using Oncomine database and TCGA cohort. Oncomine database analysis indicated that SLC2A1 expression was higher in bladder cancer (Sanchez-Carbayo et al., 2006), breast cancer (Zhao et al., 2004; Curtis et al., 2012), CRC (Ki et al., 2007; Skrzypczak et al., 2010), esophageal cancer (Su et al., 2011), kidney cancer (Jones et al., 2005; Beroukhim et al., 2009), leukemia (Andersson et al., 2007), lung cancer (Su et al., 2007; Landi et al., 2008; Hou et al., 2010; Okayama et al., 2012; Selamat et al., 2012), ovarian cancer (Yoshihara et al., 2009), and pancreatic cancer (Pei et al., 2009) than in normal tissues. Some studies also found that the expression of SLC2A1 in breast cancer (Finak et al., 2008), esophageal cancer (Kim et al., 2010), and leukemia (Haferlach et al., 2010) was lower than that in normal samples (Figure 2A). Table 1 summarizes the details of SLC2A1 expression levels in pan-cancers.

We further analyzed SLC2A1 mRNA expression level in different human cancers utilizing TCGA cohort. Figure 2B indicates the differences between SLC2A1 in varied tumor samples and control samples. SLC2A1 expression level in breast invasive carcinoma, cervical squamous cell carcinoma and endocervical adenocarcinoma, cholangiocarcinoma, COAD, esophageal carcinoma, glioblastoma multiforme, head and neck squamous cell carcinoma, kidney renal clear cell carcinoma, kidney renal papillary cell carcinoma, liver hepatocellular carcinoma, lung adenocarcinoma, lung squamous cell carcinoma, READ, stomach adenocarcinoma, thyroid carcinoma, and uterine corpus endometrial carcinoma was remarkably higher than that in control tissues and significantly decreased in kidney chromophobe and prostate adenocarcinoma.

By analyzing CRC datasets from TCGA and GEO, we further determined the difference in SLC2A1 expression between CRC and normal samples. The results of TCGA and GEO cohort analysis demonstrated that the expression level of SLC2A1 remarkably increased in CRC compared to the control group (Figures 2C,D). To further confirm the precision of the analysis results, we performed the qRT-PCR and IHC staining experiments for further verification. As shown in Figure 2F, qRT-PCR results indicated that SLC2A1 mRNA expression remarkably increased in human CRC cell lines SW480 (2.197 ± 0.127 vs. 0.954 ± 0.058) and HCT116 (2.360 ± 0.061 vs. 0.954 ± 0.058) compared with normal human colorectal mucosa cells. The IHC staining results indicated that SLC2A1 was primarily expression in the CRC cell membrane. SLC2A1 IHC scores in tumor tissues were remarkably higher than those in paracancerous tissues (2.375 ± 0.606 vs. 0.896 ± 0.592, Figures 2G–I). These results point out that overexpression of SLC2A1 may conduce to the development of CRC. To estimate the diagnostic ability of SLC2A1 in CRC, we carried the ROC curve analysis. The ROC curve indicated that SLC2A1 had high accuracy in predicting the outcomes of normal and tumor (Figure 2E), and the area under the ROC curve was 0.901 (95% CI: 0.869–0.934).

To confirm the significance of SLC2A1 in clinical environment, we analyzed clinical data from TCGA CRC cohort. As shown in Figure 3, the results showed that the SLC2A1 expression in Stage II group was lower than that in Stage III and Stage IV groups. SLC2A1 expression in N0 group was lower than that in N2 group. SLC2A1 expression in M0 group was lower than that in M1 group. During progression-free survival events, SLC2A1 expression levels in patients who died were significantly higher than those in the surviving group.

We analyzed co-expressed genes related to SLC2A1 expression in the TCGA CRC dataset using the R statistical computing language. Only the protein-coding genes were kept. The main results presented in Figure 4A genes were positively related with SLC2A1 expression, whereas 3,911 genes were negatively related with SLC2A1 expression (p < 0.05). Under the conditions of |cor| > 0.3 and p < 0.05, a total of 515 genes were obtained, including 513 positively related genes and two negatively related genes. When the threshold values were cor > 0.5 and p < 0.05, the three genes had the strong association: EPHA2 (cor = 0.538, p = 6.36E-50), KRT80 (cor = 0.538, p = 7.33E-50), and KRT19 (cor = 0.528, p = 9.93E-48), respectively. The results, as shown in Figures 4B,C, indicated that the top 50 genes are positively and negatively associated with SLC2A1 expression, respectively. The detailed description of co-expressed genes is indicated in Supplementary Table S1.

The GO term and KEGG pathway investigation of the first 200 co-expressed genes mainly associated with SLC2A1 expression were carried through using the clusterProfiler software package of R language. SLC2A1 co-expressed genes involved 242 biological processes, 80 cell components, 31 molecular functions, and 18 KEGG under the condition of p.adj <0.1 and q-value <0.2. The bubble graph shows the top five messages of biological processes, cell components, molecular functions, and KEGG, respectively. The GO term annotation indicated that these genes were primarily participated in epidermis development, cell–cell junction, and cadherin binding (Figures 4D–F). The KEGG pathway analysis indicated that these genes were chiefly participated in Shigellosis and Wnt signaling pathway (Figure 4G). Supplementary Table S2 summarizes the GO term and KEGG pathway details of SLC2A1 co-expression enrichment analysis.

To investigate the possible mechanism of SLC2A1 in CRC, the GSEA analysis was carried on the differential genes. A total of 194 gene sets were found, among which the top three enrichment pathways with the strongest correlation were REACTOME G ALPHA I SIGNALING EVENTS (FDR = 0.031, p = 0.001), REACTOME M PHASE (FDR = 0.031, p = 0.001), and REACTOME NEURONAL SYSTEM (FDR = 0.031, p = 0.001). At the same time, we also found that the different genes were involved in WP LNCRNA INVOLVEMENT IN CANONICAL WNT SIGNALING AND COLORECTAL CANCER (FDR = 0.200, p = 0.013), PID HIF1 TFPATHWAY (FDR = 0.232, p = 0.017), and REACTOME GLYCOLYSIS (FDR = 0.244, p = 0.019), respectively (Figure 5). Detailed enrichment analysis information is indicated in Supplementary Table S3.

The association between SLC2A1 expression level and CRC samples immunofiltration cells was analyzed by TIMER database. The results indicated that the expression level of SLC2A1 was positively associated with the infiltrating level of CD4+ T cell (r = 0.101, p = 4.31E-2), neutrophil (r = 0.181, p = 2.69E-4), and DC (r = 0.124, p = 1.32E-2) and negatively associated with the infiltrating level of B cell (r = –0.151, p = 2.41E-3) (Figure 6A). In addition, SLC2A1 CNV was found to be significantly correlated with the infiltration levels of B cell, neutrophil, CD8+ T cell, and DC in COAD cohort (Figure 6B).

The infiltration abundance of immune cells between the high and low SLC2A1 expression group was analyzed by ssGSEA algorithm, and the results showed that there were differences in the expressions of various immune cells between the two groups (Figure 6C), including DC (p = 0.012), eosinophils (p = 0.001), immature DC (iDC; p = 0.032), neutrophils (p = 0.001), natural killer (NK) CD56 bright cells (p < 0.001), NK cells (p = 0.002), plasmacytoid DC (pDC; p = 0.005), T helper cells (p = 0.003), and Th2 cells (p = 0.035), respectively.

To further evaluate the association between SLC2A1 and CRC samples immunofiltration cells, we analyzed the TIMER, GEPIA database, and TCGA CRC cohort to investigate the relationship between SLC2A1 and immunofiltration cells marker genes in various immunofiltration cells. As shown in Table 2, all investigations showed that SLC2A1 expression was correlated with neutrophil immune marker genes, including FCGR3A (r = 0.12, p = 4.6E-2), CD55 (r = 0.26, p = 1.1E-5) and ITGAM (r = 0.13, p = 2.6E-2). The scatter plots showed the correlation between SLC2A1 expression level and neutrophil immune marker genes, respectively (Figure 6D).

The modification of m6A has a remarkable effect in the make progress of CRC. We analyzed GSE110224 and TCGA CRC cohort to explore the association between expression level of SLC2A1 and 20 m6A relative genes in CRC. As shown in Figure 7A, analysis indicated that, in the GSE110224 and TCGA CRC cohort, the SLC2A1 expression was remarkably positively associated with IGF2BP3 and YTHDF1 (p < 0.05). Furthermore, in the TCGA CRC cohort, SLC2A1 expression level was remarkably positively associated with ALKBH5, FTO, IGF2BP1, IGF2BP2, METTL3, and RBM15B (p < 0.05) but negatively associated with METTL14, RBM15, and YTHDC2 (p < 0.05).

As shown in Figures 7B,C, the scatter plot indicated the relationship between SLC2A1 and m6A relative genes. Furthermore, TCGA CRC samples were divided into high expression group and low expression group according to the SLC2A1 expression level. As shown in Figure 7D, we attempt to analyze the differential expression level of m6A relative gene between the high SLC2A1 expression group and low SLC2A1 expression group to confirm whether there is a difference in m6A modification between the two groups in CRC. The analyses indicated that, compared with the low expression group, the expression of ALKBH5, IGF2BP1, IGF2BP3, RBM15B, VIRMA, YTHDF1, and ZC3H13 increased in the SLC2A1 high expression group, whereas the expression of METTL4 decreased (p < 0.05). As shown in Figure 7E, the analysis indicated the key genes (IGF2BP3 and YTHDF1) that have expression correlation and differential expression relationship with SLC2A1. The Kaplan–Meier curves indicated that high IGF2BP3 expression was remarkably related with poor prognosis of CRC (p < 0.05, Figure 7F), whereas YTHDF1 was not related with poor prognosis of patients with CRC. These analyses point out that SLC2A1 may be strongly associated to the m6A modification of CRC, particularly through the regulating of IGF2BP3, and eventually impact the development and prognosis of CRC.

We analyzed GSE110224 and TCGA CRC cohort to study the association between SLC2A1 and 25 ferroptosis relative genes in CRC in terms of expression level. The correlation results showed that in GSE110224 and TCGA, the expression of SLC2A1 was significantly correlated with FDFT1, GPX4, RPL8, and SLC1A5 (p < 0.05, Figure 8A). In addition, in the TCGA CRC cohort, SLC2A1 expression was positively correlated with AIFM2, ATP5MC3, CARS1, CDKN1A, CS, HSPA5, HSPB1, LPCAT3, MT1G, and TFRC (p < 0.05). In GSE110224 cohort, SLC2A1 expression was greatly negatively correlated with MT1G (p < 0.05).

The scatter graph demonstrated the correlation between SLC2A1 and ferroptosis of relevant genes (Figure 8B). In addition, according to the expression level of SLC2A1, TCGA CRC samples were set to high expression group and low expression group according to the SLC2A1 expression level. To analyze the differential expression levels of ferroptosis relative genes between the two groups in CRC (Figure 8C). The analyses indicated that, compared with the low expression group, the expressions of AIFM2, ALOX15, CDKN1A, DPP4, GPX4, HSPA5, LPCAT3, RPL8, SLC1A5, and TFRC increased in SLC2A1 high expression group, whereas the expression of CISD1 decreased (p < 0.05). As shown in Figure 8D, the analysis indicated the key genes (GPX4, RPL8, and SLC1A5) that have expression correlation and differential expression relationship with SLC2A1. The Kaplan–Meier curves indicated that high GPX4 expression was significantly correlated with poor prognosis in CRC (p < 0.05), whereas RPL8 and SLC1A5 expression were not (Figure 8E). These results point out that SLC2A1 may be significantly correlated to ferroptosis of CRC, especially through the regulation of GPX4, and eventually affect the development of CRC.

Studies have shown that the regulatory mechanism of ceRNA plays a key role in CRC, so we attempted to predict and screen the ceRNA regulatory network of SLC2A1 in CRC. We predicted 47, 110, and 89 miRNAs targeting SLC2A1 using mirDIP, micorT CDS, and miRNet tools, respectively. Venn diagrams show the results predicted by the mirDIP, micorT CDS, and miRNet tools. As shown in Figures 5, 9A, miRNAs were predicted by the three tools, including hsa-miR-148a-3p, hsa-miR-148b-3p, hsa-miR-328-3p, hsa-miR-132-3p, and hsa-miR-330-5p. Furthermore, we analyzed the association between these miRNAs and SLC2A1 expression levels and found that only hsa-miR-148a-3p (r = −0.130, p = 0.002) was negatively correlated with SLC2A1 expression levels in CRC (Figure 9C).

We further predicted that the lncRNAs might bind to hsa-miR-148a-3p using miRNet and starBase tools, and the predicted results were the same as previous studies (Liu et al., 2021a). Venn diagrams show the results predicted by the miRNet and starBase tools. A total of 14 lncRNAs were predicted by both tools, including KIF9-AS1, HCG18, HOTAIRM1, HOXA11-AS, LINC00174, NUTM2A-AS1, H19, KCNQ1OT1, HOTAIR, DHRS4-AS1, OIP5-AS1, LINC00667, A1BG-AS1, and XIST (Figure 9B). Further analysis of the correlation between these lncRNAs and hsa-miR-148a-3p expression levels revealed that AKIF9-AS1 (r = −0.240, p < 0.001), DHRS4-AS1 (r = −0.120, p = 0.002), H19 (r = −0.250, p < 0.001), HCG18 (r = −0.200, p < 0.001), HOTAIR (r = −0.200, p < 0.001), KIF9-AS1 (r = −0.140, p < 0.001), NUTM2A-AS1 (r = −0.280, p < 0.001), and OIP5-AS1 (r = −0.200, p < 0.001) were negatively correlated with hsa-miR-148a-3p expression levels in CRC (Figure 9C). The network diagram shows the relationship of the predicted ceRNA network (Figure 9D).

To further narrow the prediction range, we analyzed the expression of target lncRNA and its correlation with SLC2A1 expression in the CRC cohort. The results showed that only H19 was upregulated in CRC (p < 0.001), and there was a significant positive correlation with SLC2A1 expression (r = 0.280, p < 0.001, Figure 9E). Therefore, we predict that the H19-hsa-miR-148a-3p-SLC2A1 ceRNA network may play a significant role in the development of CRC. The network diagram shows the relationship of the final ceRNA network (Figure 9F).