Wei Yang1†
Shuaibing Lu1†
Liangqun Peng1†Zhandong Zhang1Yonglei Zhang1Dandan Guo2
Fei Ma1
Yawei Hua1
Xiaobing Chen3,4*- 1Department of General Surgery, The Affiliated Cancer Hospital of Zhengzhou University, Henan Cancer Hospital, Zhengzhou, China
- 2Department of Radiology, The Third Affiliated Hospital of Zhengzhou University, Zhengzhou, China
- 3Department of Medical Oncology, The Affiliated Cancer Hospital of Zhengzhou University, Henan Cancer Hospital, Zhengzhou, China
- 4Zhengzhou Key Laboratory for Precision Therapy of Gastrointestinal Cancer, Zhengzhou, China
Background: Colon cancer (CC) is the second most common gastrointestinal malignancy. About one in five patients have already developed distant metastases at the time of initial diagnosis, and up to half of patients develop distant metastases from initial local disease, which leads to a poor prognosis for CC patients. Necroptosis plays a key role in promoting tumor growth in different tumors. The purpose of this study was to construct a prognostic model composed of necroptosis-related genes (NRGs) in CC.
Methods: The Cancer Genome Atlas was used to obtain information on clinical features and gene expression. Gene expression differential analysis, weighted gene co-expression network analysis, univariate Cox regression analysis and the least absolute shrinkage and selection operator regression algorithm were utilized to identify prognostic NRGs. Thereafter, a risk scoring model was established based on the NRGs. Biological processes and pathways were identified by gene ontology and gene set enrichment analysis (GSEA). Further, protein-protein interaction and ceRNA networks were constructed based on mRNA-miRNA-lncRNA. Finally, the effect of necroptosis related risk score on different degrees of immune cell infiltration was evaluated.
Results: CALB1, CHST13, and SLC4A4 were identified as NRGs of prognostic significance and were used to establish a risk scoring model. The time-dependent receiver operating characteristic curve analysis revealed that the model could well predict the 1-, 3-, and 5-year overall survival (OS). Further, GSEA suggested that the NRGs may participate in biological processes, such as the WNT pathway and JAK-Stat pathway. Eight key hub genes were identified, and a ceRNA regulatory network, which comprised 1 lncRNA, 5 miRNAs and 3 mRNAs, was constructed. Immune infiltration analysis revealed that the low-risk group had significantly higher immune-related scores than the high-risk group. A nomogram of the model was constructed based on the risk score, necroptosis, and the clinicopathological features (age and TNM stage). The calibration curves implied that the model was effective at predicting the 1-, 3-, and 5-year OS of CC.
Conclusion: Our NRG-based prognostic model can assist in the evaluation of CC prognosis and the identification of therapeutic targets for CC.
1. Introduction
Colon cancer (CC) is a deadly tumor that affects individuals worldwide. The incidence of CC is increasing, especially in cities and regions with rapid economic development in the United States (1). With the popularization of cancer screening and advancements in treatment-related medical technology, patient outcomes have improved significantly. But as the onset of CC is insidious, there’s still a lot patients diagnosed in the advanced stage, where the condition is severe and difficult to treat, and palliative care is the only available treatment option (2, 3). Only few biomarkers are available for the diagnosis and therapy of CC. Circulating tumor DNA (ctDNA) is a subset of circulating free DNA (cfDNA) from tumor cells. In many studies, ctDNA has been found to be of great value in the early diagnosis, efficacy evaluation, drug resistance monitoring and prognosis prediction of tumors. Among them, targeted drug therapy guided by ctDNA is the most important clinical application at present (4, 5). At present, ctDNA as biological markers have been found to be associated with the prognosis of colon cancer, but they have not been widely applied in clinic (6–8). Therefore, exploring potential biomarkers of CC remains the focus of CC-related research.
Necroptosis is a lytic manner of programmed cell death that prevents the self-destruction of activated cells that are blocked by apoptosis. In some degenerative or inflammatory diseases, necrotizing apoptosis plays a role in destroying infected cells or damaged cells (9). Unlike apoptosis, the activation of necroptosis does not depend on caspase kinase activation. Under caspase inhibition, the binding of death receptor and ligand can trigger necrotizing apoptosis (10). Necrotizing apoptosis plays a dual role in tumorigenesis and development, which can not only enhance cellular immunity (11) and play an anti-tumor role, but also stimulate the tumor to form an immunosuppressive microenvironment and promote tumor progression (12). According to previous studies (13) and owing to the activity of intracellular RIP-1 and MLKL, the combination of 5-FU and ZVAD (caspase inhibitor) can promote necroptosis of colorectal cancer (CRC) cells, highlighting the important value of necroptosis in the study of tumor drug resistance. However, a prognostic scoring system for CC based on the tags of genes associated with necroptosis has not been established.
Recently, high-throughput sequencing and the gene chip technology have been widely used in the field of life science (14, 15). Bioinformatics is an important tool for analyzing large volumes of existing biological data. By analyzing the potentially important core genes or prognostic factors within the data (16, 17), potential tumor markers or therapeutic targets can be explored (18). Several previous studies focused on single genes as diagnostic and prognostic indicators (19, 20). However, these biomarkers, especially individual gene expression levels that may be influenced by multiple factors, are insufficient to accurately and independently predict patient outcomes. As a result, these markers cannot be used as reliable and independent prognostic indicators. Therefore, in this study, statistical models composed of multiple prognostic necroptosis-related markers were employed to improve the predictive power of CC.
In recent years, the ceRNA hypothesis has attracted attention and has become one of the hot spots in the study of RNA interaction. The regulatory mechanisms among mRNA, miRNA, lncRNA, or circRNA are extremely complex and have important biological significance. LncRNA or circRNA can compete with mRNA to bind to miRNA, thereby forming a complex lncRNA-miRNA-mRNA network or circRNA-miRNA-mRNA network. However, an imbalance in the ceRNA regulatory network can lead to the initiation and progression of tumors (21). To date, the function of most mRNAs as ceRNAs in the progression and prognosis of CC has not been thoroughly defined.
In this study, the prognostic risk model of necroptosis-related genes (NRGs) was established, and the diagnostic and predictive significance of the model was evaluated. Thereafter, a ceRNA network was constructed based on mRNA-miRNA-lncRNA, and the effects of necroptosis-related risk score on different degrees of immune cell infiltration were evaluated. Overall, the findings of this study provide a theoretical foundation for further assessments of the diagnosis, treatment, and molecular mechanism of CC.
2. Materials and methods
2.1. Data download
The Cancer Genome Atlas (TCGA) Genomic Data Commons (GDC) website (https://portal.gdc.cancer.gov/) was used to obtain the expression spectrum data for colon cancer (colon adenocarcinoma, COAD) patients (n = 514), such as the count, and FPKM and TPM values; and patient clinical data (n=430), such as gender, age, TMN stage, and survival prognosis. “Masked somatic mutation” was selected as the somatic mutation data (n=420) and downloaded. The somatic mutations were visualized using maftools package (22) in R to obtain tumor mutation burden (TMB) for per patient. Additionally, the MSI data in the TCGA-COAD patients’ dataset was obtain the tumor mutation burden (TMB) per patient. Additionally, the MSI data in TCGA-COAD patient dataset were obtained from the cBioPortal database (https://www.cbioportal.org). The baseline information of TCGA-COAD patients is provided in Table 1.
Table 1 Baseline data table of patients in TCGA-COAD dataset.
Gene expression data and the clinical characteristics of patients of GSE17536 (23) and GSE39582 (24) were downloaded from the GEO database. The data samples were obtained from Homo Sapiens. The chip platforms were grounded in the GPL570 [HG-U133_plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. After deleting patients lacking clinical information, 177 COAD tissue samples were included in GSE17536, all of which were used in the analysis. GSE39582 included 585 COAD tissue samples. As survival information was not available for 5 samples, 580 were included in the analysis. R’s limma package (25) was used to standardize the two data sets separately. Table S1 shows the information from GEO.
2.2. Calculation of the necroptosis score based on gene expression matrix
For all samples in the combined dataset and based on 36 necrotizing apoptosis-related genes from previously published literature (26), the necroptosis score (NPs) of every sample in the TCGA-COAD dataset was determined using the R package, GSVA (27), and the ssGSEA method according to the gene expression matrix of the respective sample.
2.3. Screening of differentially expressed genes (DEGs)
Using the above method, the NPs of each sample was obtained, and the optimal cut-off value was selected in combination with the patient’s survival data. Based on the NPs score, the NPs group was separated into high and low NPs. To identify genes associated with NPs, the DEGs between high NPs and low NPs in TCGA-COAD samples were identified using the R package, limma (25). The screening threshold of the DEGs was set to |log2 fold change (FC)|>1 and adjusted P < 0.05. The results of the difference analysis are presented as a heatmap and a volcano plot.
2.4. Weighted gene co-expression network analysis (WGCNA)
WGCNA was implemented using the R package, WGCNA (28). First, the weighted values of the calculated correlation coefficients between any two genes were used to generate connections between genes in the network to assemble a scale-free network. A hierarchical clustering tree was then established according to the correlation coefficients. The branches of the cluster tree highlighted various genetic modules, and various colors signified different modules. The module saliency was then calculated. All mRNAs of the sample were input into WGCNA to measure the associations between the two NPs groups and different modules. All genes were recorded in their respective modules. The genes in the respective modules were considered as modular characteristic genes (MEs). The correlation between the NPs values and genes was determined based on the significance of genes. Module membership was determined according to the relevance between module genes and DEG expression profile. The modules of interest were selected using the MS score, and all genes in these modules were recognized to have a high correlation with NPs.
2.5. Subtype analysis of patients with CC based on necroptotic genes
Based on the necroptosis characteristic genes and TCGA-COAD expression data, the k-means method in “ConsensusClusterPlus” R package (29) was used to perform unsupervised cluster analysis to identify the necroptosis subtypes. The concordant clustering algorithm was used to discern the cluster number. The analysis was repeated 1000 times to ensure stability of the category. Principal component analysis (PCA) was conducted for patients with grouped subtypes to determine the differences between samples. Survival analysis was implemented after grouping to determine the influence of various subtypes on prognosis.
2.6. Establishment and testing of the prognostic risk models based on NPs characteristic genes
The results of differential expression analysis combined with WGCNA analysis were used to acquire the NPs-related characteristic genes. The significantly differentially expressed NPs-related characteristic genes were involved in the model, and the prognostic genes were screened using univariate Cox regression analysis, with a cut-off P value of 0.1. Subsequently, the selected genes were regularized and dimensionally reduced using the least absolute shrinkage and selection operator regression (LASSO) algorithm to further identify prognostic-related genes. Thereafter, the weighted normalized gene expression value of the penalty coefficient acquired by multivariate Cox analysis (STEP method) was used to establish a risk score formula. Using the median risk score, patients were divided into high-risk and low-risk groups.
The above dataset based on TCGA-COAD served as a training set, and the internal test was conducted using the bootstrap method with 1000 re-sampling. Thereafter, the coefficient based on model variables was used to calculate the risk score for each sample in test sets, GSE17536 and GSE39582, using the predict function in the “survival” R package (https://CRAN.R-project.org/package=survival>). Finally, a time-dependent receiver operating characteristic (ROC) curve was plotted. The area under the curve (AUC) was used to reflect the performance of the model.
2.7. Analysis of DEGs in the NPs-related metabolic model
To acquire genes relevant to the NPs model, DEGs between the high-risk group and low-risk group of TCGA-COAD patients were analyzed using the R limma package. The screening threshold of the significantly different DEGs was defined as |LogFC| > 1 and adj. P value < 0.05. The DEGs were visualized using volcano maps and heat maps.
2.8. Functional enrichment analysis
GO (30) analysis is an approach adopted for massive functional enrichment research, including biological processes (BP), molecular functions (MF), and cellular components (CC). Kyoto Encyclopedia of Genes and Genomes (KEGG) (31) is an extensively used database that stores data regarding genomes, biological pathways, diseases, and drugs. The ClusterProfiler package of R (32) was used for GO and KEGG analyses of significant DEGs. The critical value of FDR less than 0.05 indicated significant difference.
To explore the discrepancies in biological processes among different subgroups, GSEA was performed according to the gene expression profiling dataset of COAD patients. The gene set “c2.cp.v7.2.Symbols. gmt” obtained from the MSigDB (33) database was used for GSEA. An FDR < 0.25 indicated statistical significance.
2.9. Identification and correlation analysis of the tumor immune infiltrating cells
To quantitatively analyze the relative tumor infiltration degree of various immunocytes in COAD, the ssGSEA algorithm was employed to differentiate between highly sensitive and specific phenotypes of various human immune cells in the tumor microenvironment (TME). The algorithm revealed 28 gene sets for labeling various tumor-infiltrating immunocyte types based on a study by Bindea et al. (34). The gene sets comprised various human immunocyte subtypes, such as macrophages, mast cells, etc. Enrichment scores obtained using ssGSEA in R’s GSVA package indicated the degree of infiltration of various immune cell types in every sample. Meanwhile, R’s ESTIMATE package (35) was employed to evaluate the immunological activity of the tumor. ESTIMATE quantitatively analyzes the immune activity of a tumor sample according to its gene expression profile to obtain an immune score per tumor sample. Herein, the discrepancies in immune infiltration features between the two groups of patients with COAD were compared.
2.10. Analysis of copy number variation (CNV)
To compare the copy number differences between the two groups of TCGA-COAD patients, the TCGAbiolinks package (36) of R was used to obtain the Masked Copy Number Segment information. The downloaded CNV fragments were subjected to GISTIC 2.0 analysis by GenePattern (https://cloud.genepattern.org); default parameters were used for the analysis.
2.11. Establishment of the prognostic model according to the NPs risk score
The predictive power of the NPs risk score in combination with clinicopathological characteristics on OS based on univariate and multivariate Cox analysis was used to demonstrate that the NPs risk score in combination with clinicopathological features can be used to estimate patient prognosis. The risk scoring model was then combined with clinicopathological features to establish a nomogram, and the accuracy of the model was reflected by the AUC values under the time-ROC curve. The performance of the rosette was assessed using a calibration curve that compared the predicted values of the rosette with the observed actual values. Testing of the model was carried out using the bootstrap method, and internal re-sampling was performed 1000 times.
2.12. Establishment of the PPI network and screening of hub-genes
The STRING (37) online tool was applied to establish the PPI network. Genes with scores > 0.7, which indicates high credibility, were selected from the STRING database to construct the network model visualized using Cytoscape (version3.7.2) (38). The Maximal Clique Centrality (MCC) of each node was calculated using the cytoHubba plug-in (39) in Cytoscape software.
2.13. Establishment of the ceRNA network according to mRNA-miRNA-lncRNA
Information on the miRNA-mRNA interactions was collected from the miRTarBase database (40). The core mRNAs acquired from the PPI analysis were used to predict the miRNAs that might be regulated. The relevant lncRNAs were further predicted based on evidence from the luciferase reporter gene assay. The results of ceRNA analysis were visualized using Cytoscape software. The P values of all hypothetical tests were two-sided, and a p value of less than 0.05 was considered to indicate statistical significance.
2.14. Statistical analysis
Data analysis was performed using R software. version 4.0.2. Independent Student’s t test and Mann–Whitney U test (namely Wilcoxon rank-sum test) were used to estimate the differences between two groups of normally distributed and two groups of non-normally distributed continuous variables, respectively. The χ2 test or Fisher exact test was carried out to determine the difference between the two groups of categorical variables. Survival analysis was carried out using R’s survival package. The Kaplan–Meier survival curve was applied to display survival differences, and log-rank test was performed to compare the differences in survival. Univariate and multivariate Cox analyses were based on the R survival package. LASSO analysis was carried out using the glmnet R package (41).
3. Results
3.1. Expression and mutation of necroptotic genes in CC patients
The whole research design was illustrated in Figure 1. First, 36 necroptotic genes were extracted from the RNA-seq data of TCGA-COAD and their expression differences were compared between the normal group and tumor group. A total of 30 necroptotic genes were found to be differentially expressed, and only 6 genes were not differentially expressed. Such results suggest that necroptosis may play crucial role in COAD (Figures 2A, B). According to the somatic mutation data of TCGA-COAD samples, mutation information was obtained for the 36 necroptotic genes using the maftools package. The necroptotic genes were not found to mutate significantly in COAD patients, except TP53, in which the mutation frequency of TNF, CASP6, and TNFSF10 was less than 1% (Figure 2C). In addition, the prognostic status of the 36 genes was analyzed. Only TRAF2, RIPK3, and IPMK genes were found to have significant prognostic differences, and may thus serve as potential prognostic markers (Figure S1).
Figure 1 The flow chart.
Figure 2 Differential expression and mutation information of necroptotic genes in TCGA-COAD datasets. (A, B) Necroptotic genes were compared between the normal and tumor groups in TCGA-COAD dataset using the Wilcoxon test. *P < 0.05, **P < 0.01, ***P < 0.001. (C) Mutation information of the necroptotic genes in TCGA-COAD. "ns" represents "no significance".
3.2. Calculation of the necroptosis score and screening of characteristic genes
Based on the above differential necroptosis genes, the NPs of each COAD patient was obtained using the ssGSEA algorithm to represent the necroptosis level of the patient. The optimal cut-off value was determined using prognostic analysis based on the necroptosis score. Thereafter, TCGA-COAD samples were divided into groups with high NPs and low NPs. Survival analysis revealed that patients with low necrosis apoptosis scores had worse prognosis than patients with high necroptosis score (log-rank P < 0.024; Figure 3A). Subsequently, 766 DEGs were acquired through rigorous analysis, including 380 significantly upregulated and 386 downregulated genes, respectively (Figures 3B, C). A gene co-expression network was also established to identify biologically significant gene modules through WGCNA and further identify genes closely related to COAD necroptosis. In this study, 4 modules (except grey module) were obtained for subsequent analysis (Figures 3D, E). As shown in Figure 3F, we integrated the difference analysis results with the MEturquoise, MEblue, and MEbrown modules to obtain a total of 209 necrotizing apoptotic characteristic genes.
Figure 3 Screening of genes associated with necroptosis. (A) Survival analysis results revealed marked difference in survival status between groups with high and low necrotizing apoptosis scores (log-rank P =0.024); (B, C) Volcano maps and heat maps revealing DEG expression among COAD samples in the groups with high and low necrotizing apoptosis. (D) Quality control result selected by WGCNA softpower as 4. (E) Set of genes associated with the necroptosis phenotype analyzed and screened using WGCNA. The heat map demonstrated the correlation and significant difference between different gene modules and necroptosis score, where the P values are shown in parentheses. (F) Intersection of DEGs and genes in the significant module of WGCNA.
3.3. Identification of necroptotic subtypes
Based on the above genes with necroptosis characteristics, consistent clustering was employed to cluster LIHC samples. Here, K=2 was selected and two subgroups, subgroup 1 and subgroup 2, were obtained (Figure S2A). Dimension reduction analysis was conducted via PCA and the PCA results of the two groups were plotted. Based on the results, the degree of differentiation between the two groups was not obvious, which may be due to the insignificant clustering gene characteristics, resulting in insignificant grouping differences (Figure S2B). The prognostic characteristics of subgroup 1 and subgroup 2 was subsequently analyzed using the KM curve; however, no distinct difference in prognosis was found between the two groups, suggesting that clustering subtypes based on all necroptotic characteristic genes could not distinguish differences in the prognosis of patients (Figure S2C). However, when the expression distributions of characteristic genes of type 1 and type 2 were compared, necroptotic characteristic genes were found to be significantly differentiated in the two subtypes, suggesting that the classification is of guiding significance for assessing the mechanism of necroptosis, but not suitable for identifying clinical prognostic markers (Figure S2D).
3.4. Construction and evaluation of risk models related to necroptosis
Based on the necroptosis characteristic genes, a necroptosis-related risk score system was constructed to quantitatively evaluate the prognostic information of each COAD patient by risk score. First, univariate Cox regression analysis revealed that 23 genes met the screening criteria (P < 0.05). Dimension reduction was analyzed using LASSO (Figure 4A). When 5 variables were present, the most stable model was obtained. Multivariate Cox regression analysis revealed that Calbindin 1 (CALB1), Carbohydrate sulfotransferases (CHST13), and solute carrier family 4 member 4 (SLC4A4) were independent prognostic factors (Figure 4B). Multivariate Cox analysis was also carried out to obtain the model coefficients of important characteristic genes. Thereafter, the gene expression was multiplied and summed with its coefficients to construct a risk score. The final risk score (necroptosis risk score related to prognosis) was calculated for each sample. In terms of the risk score and gene expression values of patients, a heat map of the risk factors was plotted to show the distribution of the risk score (Figure 4C). The time-dependent ROC curve analysis revealed AUC values of 0.684, 0.657, and 0.710 for the 1-, 3-, and 5-year OS, respectively, which indicated that risk score was an ideal predictor of OS in COAD patients (Figure 4D). For the external dataset test, the GSE39582 and GSE17536 datasets were employed. After data normalization, the model was tested. For GSE39582, the ROC curve revealed AUC values of 0.682, 0.623, and 0.708 for the 1-, 3-, and 5-year OS, respectively. For GSE17536, the AUC values were 0.665, 0.712, and 0.758 for the 1-, 3-, and 5-year OS, respectively, indicating that the model had been well tested in external datasets (Figure S3).
Figure 4 Prognostic model and model test based on necroptotic characteristic genes. (A) LASSO regression analysis; the number of variables corresponding to the optimal λ value is 5. (B) Three genes identified as independent prognostic factors through multivariate Cox stepwise regression analysis. *P < 0.05. (C) Risk score distribution and survival status of COAD patients; (D) TimeROC curve of TCGA-COAD (training set). Internal validation was performed using the Bootstrap method, with 1000 iterations.
3.5. Analysis of DEGs and functional enrichment in patients with high- and low-risk necroptosis score
To determine the role of the necroptosis-related risk model on the evolution of COAD samples, TCGA-COAD patients were divided into high-risk group and low-risk group based on the expression median value of the COAD patient risk model score in TCGA dataset. Subsequently, the DEGs in the two groups of patients were identified. Overall, 317 genes were significantly differentially expressed in COAD patients, among which 185 and 132 genes were significantly upregulated and downregulated, respectively (Figures 5A, B).
Figure 5 DEG analysis and functional enrichment analysis based on the necroptosis-related risk model. (A, B) Volcano map and heat map revealing DEG expression between the high- and low- risk groups in TCGA-COAD dataset. (C) GO analysis revealed that the differential genes were correlated with GO:0044421 extracellular region part, GO:0005576 extracellular region, GO:0042588 zymogen granule, GO:0071752 secretory dimeric IgA immunoglobulin complex, and other functions. (D) KEGG results revealed that these DEGs participated in nitrogen metabolism, bile secretion, rheumatoid arthritis, and other pathways. (E) GSEA results suggested that the KEGG results were similar to those of differential gene enrichment, and the main enrichment pathways were the WNT pathway, JAK-STAT pathway, immune-related pathway, etc.
Functional enrichment analysis was performed using 317 genes identified as significantly DEGs. The GO analysis results revealed that the significant DEGs were related to GO:0044421 extracellular region part, GO:0005576 extracellular region, GO:0042588 zymogen granule, GO:0071752 secretory dimeric IgA immunoglobulin complex, and other functions (Figure 5C). KEGG functional analysis suggested that significant DEGs mainly had an impact on nitrogen metabolism, bile secretion, and rheumatoid arthritis, and other pathways (Figure 4D). Many pathways were found to be related to immunity, such as the chemokine signaling pathway, WNT signaling pathway, etc. Detailed GO and KEGG results are presented in Tables S2 and S3.
Based on the results of expression analysis, we continued GSEA and summarized the related results of the pathway database based on C2. KEGG pathway results and GSEA results revealed distinct differences in the activity of the JAK-STAT signaling pathway, WNT signaling pathway, fructose and mannose metabolism, primary immunodeficiency, and nitrogen metabolism (Figure 5). The detailed results of GSEA and the metabolism-related pathways are provided in Table S4. These findings coincide with those obtained from the KEGG database, and suggest the activation and inhibition characteristics of the high- and low-risk groups.
3.6. Protein interaction and regulatory network analysis
In terms of the DEGs in the high- and low- risk groups, we aimed to identify the hub gene that played a key role, and its potential molecular interaction mechanism. First, the STRING database was used to analyze the protein interaction mechanism. As shown in Figure 6A, after screening with a confidence of 0.700, the number of PPI nodes (protein) was 233. Further, 81 edges were identified, with an average connection degree of 0.695 for each node. The enrichment statistic P value of the whole PPI network was less than 1.0e-16.
Figure 6 PPI and regulatory network analysis. (A) PPI regulation network, detailed display of the network node information, connection line information, and the composition of the different sub-network information. (B) Hub gene regulation network based on cytoHubba calculation. (C) CeRNA regulatory network predicted using the miRTarBase database. Blue represents miRNA; green represents LncRNA; and Brown represents mRNA.
The interacting proteins were further identified as hub genes using the cytoHubba plugin in Cytoscape. After the calculations, MUC5AC, MUC5B, WNT16, WIF1, and other interacting proteins that had the top 8 scores were found (Figure 6B). Subsequently, miRNA molecules and lncRNAs that potentially regulate these hub genes were analyzed using miRTarBase database. Finally, a ceRNA regulatory network was established using Cytoscape (Figure 6C).
3.7. Differential expression analysis of immunocyte infiltration
The influence of necroptosis-associated risk scores in patients with TCGA-COAD on their holistic immune characteristics and varying degrees of immunocyte infiltration was analyzed. Patients in the high-risk group were found to have significantly lower immune-related scores (P < 0.001); however, no distinct difference was found in the matrix score (Figures 7A, B). ssGSEA was used to appraise the changes and effects of immunological characteristics of COAD tissues during pathogenesis. Through ssGSEA, the relative enrichment scores of 28 different subtypes of immunocytes in the high- and low-risk groups of COAD patients was obtained. Heat maps were generated to illustrate their expression in different patients (Figure 7C). Based on the results, the expression abundance of immunocytes in the low-risk group was lower than that in the high-risk group. The correlation analysis results revealed that most immune cell infiltration levels were positively correlated (Figure 7D). Differential analysis also revealed distinctions in the infiltration levels of various immune cells between COAD samples in the high and low necrotizing apoptosis-related risk groups. Only effector memory CD8+ T cells, immature dendritic (iDC) cells, and other cells were not found to significantly differ between the groups (Figure 7E). A distinct difference was found between the HLA family expression levels and the various immunological targets of the high- and low-risk groups. Further, the immunoactive genes were almost all elevated in the low-risk group (Figures 7F, G).
Figure 7 Correlation between necrosis risk scores and the infiltrates of different immunocytes (A, B) Immune score and stromal score between the low- and high- risk groups; (C) Heat map showing the invasion degrees of 28 different immune cells in TCGA and GEO database; (D) Association heat map showing the association between various levels of immunocyte infiltration. (E) Differential analysis of 28 different immunocyte infiltration levels between the two groups; (F) Analysis of differences in the expression of multiple members of the HLA family between the high- and low-risk related subgroups; (G) Differential expression analysis of multiple immunotherapy-related targets between the high- and low-risk related groups. *P < 0.05, **P < 0.01, ***P < 0.001. "ns" represents "no significance".
3.8. Effects of the necroptosis-related risk score on genomic changes in COAD samples
The influence of the necroptosis-related risk score on changes in genetic variation levels, including single nucleotide polymorphism and CNV, in COAD patients was evaluated. The analysis of single-nucleotide mutations in common tumor-driven genes revealed that the high mutation levels were similar or close between patients with high and low scores in the necroptosis-related model (Figure 8A). Based on assessments of the frequency of CNV changes, CNV was found to be widely present in high- and low-risk samples. However, no distinct discrepancy in CNV was found between the two groups (Figures 8B, C). When TMB and MSI were compared between the two patient groups, no distinct discrepancy in TMB and MSI was found between the high and low risk groups. This result suggests that changes at the genomic level were not significant in the two groups (Figures 8D, E).
Figure 8 Impact of necroptosis-related risk grouping on genetic variation and immunotherapy in COAD samples. (A) Mutation map of common tumorigenic driver genes in patients in the high- and low-risk groups. Mutation information per gene per sample is presented as a waterfall plot, and different colors represent different types of mutation. The subsection above the legend shows the sudden change load; (B, C) Changes in the copy number levels of different genes in the high-risk and low-risk groups, where genes with significant copy number increase in red and genes with significant copy number deletions in blue; (D, E) Comparison of the difference in MSI level and TMB level between patients in the high- and low-risk groups, respectively; (F) Discrepancy between the high-risk and low-risk groups based on the tide score calculated from the tide database. ***P < 0.001. "ns" represents "no significance".
Based on the significant role of immunotherapy in tumors, the TIDE algorithm was employed to calculate the sensitivity of patients in the high- and low-risk groups to immunotherapy. The TIDE score in the low-risk group was higher than that in the high-risk group, suggesting that the immunotherapy response of the high-risk group might be better than that of the low-risk group (Figure 8F).
3.9. Establishment of a prognostic model according to the necroptosis-related risk score
To further probe the clinical value of necroptosis-related risk score, the clinical characteristics related to the high-risk and low-risk groups, such as the discrepancy in age and TNM stage, were analyzed. Notably, no distinct discrepancy was found in the age of patients in the high-risk group (Figure 9A). For gender, the proportion of women in the high-risk group increased (Figure 9B). In terms of stage, a significantly higher proportion of advanced patients was identified in the high-risk group (Figure 9C). Subsequently, based on the risk scores associated with necroptosis and clinicopathologic features (age and TNM stage), we established a prognostic model for COAD patients (Figure 9D) and analyzed the model via 1000 resampling using the bootstrap method. Based on timeROC, the AUC values were 0.798, 0.772, and 0.741 for 1-, 3-, and 5- years, respectively (Figure 9E). Calibration curves were generated to present the consistency of the model. A good consistency was found between the model’s estimated 1-, 3-, and 5-year OS and the actual observed OS of patients (Figure 9F).
Figure 9 Performance of the necroptosis risk scores in the prediction of prognosis for patients with COAD. (A-C) Superimposed histogram showing the proportion of age, sex, and stage in patients in the high- and low-risk groups. The effect of age was similar in both groups, with an increased proportion of women in the high-risk group and significantly more advanced patients in the high-risk group. (D) Nomogram of the model. (E) Time-dependent ROC curve of the clinical prediction model based on risk score. (F) For the calibration curve of the nomogram, the bootstrap method was adopted, and resampling was performed 1000 times. The abscissa is the survival predicted by the nomogram, and the ordinate is the actual observed survival. The calibration plot revealed that the bias-corrected line for 1-, 3-, and 5- years OS was close to the ideal line, indicating good consistency between the predicted value and the actual value.
4. Discussion
In recent years, the incidence of CC among young people has gradually increased (42). Necroptosis has different functions in diverse tumors, including promoting tumor progression in lung cancer, pancreatic cancer, and glioblastoma (43–45), or inhibiting tumor growth in gastric cancer (GC), head and neck squamous cell carcinoma, melanoma and CRC (46–49). Necroptosis also has a two-way effect of promoting cancer and suppressing cancer in breast cancer (50, 51). As a result, we cannot appraise the prognosis of CC according to the expression of individual necrosis regulators alone. Targeting NRGs is regarded as one of the effective methods for reducing tumor chemotherapy resistance, opening up a new approach for cancer treatment (52). A prior study revealed the construction of a prognostic model of lncRNA associated with GC necroptosis to differentiate hot and cold tumors of gastric carcinoma, to ultimately predict prognosis and the effectiveness of immunotherapy (53). Nevertheless, the theory of necroptosis in CC remains indistinct. In this study, prognostic risk models based on NPs characteristic genes was constructed to predict prognosis and immunotherapy, and systematically analyze the correlation between immune cell infiltration, immune checkpoints, and CC.
A total of 30 necroptotic genes were found to be differentially expressed. Thereafter, an analysis of DEGs revealed 766 DEGs in the high NPs group and low NPs group. A gene co-expression network was also established to identify biologically significant gene modules through WGCNA. Finally, a total of 209 necrotizing apoptotic characteristic genes were identified. The results of univariate Cox regression analysis, LASSO, and multivariate Cox regression analysis revealed that CALB1, CHST13, and SLC4A4 are independent prognostic marks. The final risk score was then calculated for each sample.
CALB1 is a vitamin D-dependent calcium-binding protein with six EF hands on the long arm of chromosome 8 at position 21.3 (54). CALB1, a component of Calbindin, has been confirmed to restrain tumor cell apoptosis. A prior study suggested that CALB1 may exert carcinogenic effects in ovarian cancer by inhibiting the p53 pathway (55). CALB1 is overexpressed in nonsmall cell lung cancer (NSCLC) tissues, and has a significant connection with lymph node metastasis and prediction of worse survival (56). In osteosarcoma, the downregulation of CALB1 gene expression resulted in reduced cell proliferation and cell clonal formation (57). In this study, CALB1 was verified to be an independent risk factor for prognosis. Further, its expression was found to increase, indicating poor prognosis of patients. Previous studies did not directly explain the relationship between CALB1 and CC. However, this study provides ideas for future diagnosis and treatment using CALB1 as an oncogene.
Chondroitin sulfate (CS) is a glycosaminoglycans (GAGs) that participates in multiple biological processes and exerts crucial function in the interaction among stromal tumor cells (58). CHST13 gene is located on chromosome 3q21.3. A prior study suggested that CHST13 may serve as a negative regulator of HCC cell invasion and chemotherapy sensitivity by modulating Mitogen-Activated Protein Kinase (MAPK) activity (59). The mRNA expression of CHST13 was found to be significantly higher in in ovarian cancer specimens than in non-malignant tumor specimens (60). The results of this study indicate that CALB1 is an independent prognostic marker that plays the role of an oncogenic gene in the occurrence and development of CC. Thus, CALB1 could serve as an original biomarker for the diagnosis and prognosis evaluation of CC. However, no prior study has revealed that CHST13 can serve as a high-risk independent prognostic factor for OS. Accordingly, the present study is an important supplement to this field.
Homo sapiens solute vector family Member 4 (SLC4A4) is a member of the solute vector family which encodes an electrogenic Na+/HCO3− cotransporter (61). A previous study showed that SLC4A4 is increasingly expressed in prostate cancer tissues and cell lines. Further, the SLC4A4 expression level in cancer tissues was significantly associated with the degree of disease progression. SLC4A4 promotes prostate cancer progression through the Akt-mediated signaling pathway (62). Mir-222-3p expression was increased in PTC, while that of SLC4A4 was low. SLC4A4 could reverse the promoting function of Mir-222-3p on the proliferation, invasion, and migration of PTC cells (63). SLC4A4 had a lower expression in CRC than normal tissue, indicating that SLC4A4 was associated with poor prognosis (64). This study revealed that SLC4A4 may be an individual prognostic factor for CC patients and may exert a protective function in the tumorigenesis and progression of CC, which aligned well with the proposals from existing studies.
TCGA-COAD samples were divided into high- and low-risk groups based on the median expression value of the COAD patient risk model score in TCGA dataset. DEG analysis and functional enrichment analysis were then performed. The main enrichment pathways included the WNT pathway, JAK-STAT pathway, primary immunodeficiency pathway, chemokine pathway, fructose and mannose metabolism, nitrogen metabolism, etc.
The abnormal WNT signaling pathway is highly relevant to tumorigenesis and progression of multiple tumors, including CC (65–67). A previous study confirmed that activation of Wnt/β-catenin signaling contributes to the aberrant expression of several oncogenes that regulate the dedifferentiation phenotype and EMT in CC cells (68). Another study demonstrated that RBBP4 activates the Wnt/β-catenin pathway to accelerate the progression of CC (69). Based on this analysis, the WNT signaling pathway was identified to be significantly differentially enriched in the high-risk group phenotype.
The JAK/STAT pathway plays an increasingly vital role in regulating immune function, cell proliferation, differentiation, and death (70–72). Fibroblast growth factor receptor was reported to mediate PD-L1 expression in CC by activating the JAK2/STAT3 signaling pathway (73). Notably, activation of the JAKs structure promotes phosphorylation of the STAT family (74). The STAT3 signaling pathway is in close contact with the construction of a tumorigenic inflammatory microenvironment (75). The proliferation and viability of macrophages were reported to be enhanced by STAT3 activation, the immune tolerance of CC cells, and inhibition of extracellular matrix remodeling, thereby playing tumor-promoting roles (76). Based on this analysis, the JAK/STAT signaling pathway was identified to be significantly differentially enriched in the low-risk group phenotype. Consistent with the results of this study, necroptosis may promote tumor progression by inhibiting the JAK/STAT pathway.
A systematic review of all cases of clinically diagnosed primary immunodeficiency and early-onset gastrointestinal (GI) cancer in three publicly available databases (MEDLINE, SCOPUS, and EMBASE) was previously conducted. Based on the results, primary immunodeficiency may be linked with potential risk factors for GI tumor. Adenocarcinomas of the stomach and colon were identified as the most common GI tumor (77). A previous literature revealed the involvement of chemokine (CC theme) ligand 7 (CCL7) in the progression of CRC (78). Another literature revealed the ectopic expression of the novel chemokine, CXCL17, in primary CC. The expression of CXCL17 might inform the prognosis of CC patients as CXCL17 enhances angiogenesis and attracts immune cells (79).
Based on the differential genes in the high and low risk groups, we opted to identify the key hub genes and their underlying molecular interaction mechanisms. These hub genes were MUC5AC, MUC5B, WNT16, WNT11, WIF1, SFRP5, B3GNT6, and GALNT12.
Mucin is a type of high molecular weight glycoprotein that is mainly involved in protecting epithelial cells of different organs from physical, chemical, and pathogenic damage (80). Mucin has abnormal expression in many malignant tumors, which is correlated with the proliferation, migration, invasion, adhesion and metastasis of tumor cells (81, 82). Changes in mucin expression have been reported to have a high correlation with the occurrence of CRC (81). Normally, expression of the secreted mucin, MUC5AC, is restricted to the stomach, lung, ear, conjunctiva, nasopharynx, and gallbladder. Several studies revealed that secreted MUC5AC is overexpressed in pancreatic cancer, lung cancer, and breast cancer (83–85). In fact, the secreted mucin, MUC5AC, was not identified in normal colonic mucosa, but was present in benign and malignant colon (80, 86). Prior literature confirmed that MUC5AC across the membrane protein, CD44, mediated the initiation and progression of CC, and provided resistance to chemotherapy in CRC through the β-catenin/p53/p21 signaling pathways (87). Secreted MUC5B mucin is generally not expressed in normal adult gastrointestinal mucosa, but has been proven to be differentially overexpressed in some subtypes of GC and CRC (88–90).
Numerous studies proved that over-activation of the Wnt signaling pathway is the main culprit in the onset of most human malignant tumors (91, 92). The Wnt signaling pathway plays a crucial role in multiple biological processes, such as embryogenesis and tissue homeostasis, exerting significant functions in the tumorigenesis and progression of CRC (93). A previous study found one or more mutations downstream in the Wnt signaling pathway, especially adenomatous polyposis coli (APC), in more than 90% of patients with CRC (94). WNT16 is one of the most impressive members of the WNT pathway (95, 96). The Wnt signaling pathway consists of canonical signals and noncanonical signals. Transmembrane proteins and their receptors mediate canonical Wnt signaling. Atypical Wnt signaling involves two pathways: the Wnt/Ca2+ pathway and the Wnt/c-jun N-terminal kinase (JNK)(planar cell polarity) pathway (97, 98). Wnt11 exerts its role via the noncanonical WNT pathway (97). Studies have confirmed that Wnt11 has a vital effect in the regulation of CRC cell proliferation, migration, and invasion (99, 100).
Wnt inhibitory factor 1 (WIF1) can interact with the Wnt protein to inhibit the canonical and non-canonical Wnt pathways to exert tumor inhibitory effect. WIF1 silenced by methylation has been found to participate in CC progression (101).
Secreted frizzled-related protein 5 (SFRP5) is a new type of adipocytokine, belonging to the SFRP family. Plasma SFRP5 levels were found to be distinctly decreased in obese patients and patients with diabetes, coronary artery disease, and other related diseases (102, 103). SFRP5 is underexpressed in moderate tumor tissues including lung cancer, ovarian cancer, GC, and breast cancer tissues, and is associated with poor prognosis (104–107). GALNT12 has been revealed to be a strong candidate for CRC susceptibility (108).
The B3GNT protein family is differentially expressed in multiple cancers, such as GI cancer, pancreatic carcinoma, and prostate cancer (109–111). The expression of B3GNT was found to be significantly decreased in GC and CRC (112). Although the 8 hub genes are relevant to tumorigenesis and progression, relevant studies on CRC are insufficient.
Based on increasing evidence, the ceRNA regulatory network plays a key role in the progression of various common cancers (113, 114). Shang et al. (115) found that the tumor-derived exosome, circPACRGL, acts as a sponge molecule of miR-142-3p/miR-506-3p, promoting the propagation, diversion, invasion, and adhesion of CC cells and N1 to N2 neutrophil differentiation. Wu et al. (116) showed that the LNC473-MIR574/miR15B-APAF1 IRES signaling axis could manipulate the propagation and apoptosis of CRC cells to influence the initiation and progression of CRC. In this study, a ceRNA regulatory network was constructed with 1 lncRNA, 5 miRNAs, and 3 mRNAs, revealing the potential regulatory mechanism of lncRNA-miRNA-mRNA in CC, and indicating the direction for further exploration of the pathogenesis of CC.
The immune microenvironment of CC and immunotherapy for CC patients should be explored. Immune cells in the TME perform vital functions in tumor progression (117). Based on prior studies, immune checkpoint inhibitors (ICIs) have great potential in immunotherapy of CC (118).
Immunocheckpoint inhibitor therapy of CC is in the “MSI era” because microsatellite instability (MSI) or mismatch repair gene status (MMR) is the best predictor of efficacy. Based on MSI status, CC patients can be divided into two groups according to the efficacy of immunotherapy: “advantaged population” – MSI-H/dMMR type (MSI-H type for short); “Invalid population” – MSS/pMMR type cancer (MSS type cancer for short) (119). However, only about 5% of metastatic colorectal cancer (mCRC) is MSI-H, and about 95% is MSS type (120). How to turn “cold tumor” into “hot tumor” effective for immunotherapy has been a hot research direction. EYNOTE 016 Phase II clinical trial results showed that the objective response rate (ORR) of MSI-H mCRC patients was 40%, while the ORR of MSS mCRC patients was 0 (121). Immunotherapy has enriched the treatment modalities of multiple malignancies, including cytotoxic T lymphocyte-associated protein 4 (CTLA-4) and inhibitor of programmed death-1 (P D-1)/programmed cell death ligand 1(PD-L1). In patients with MSS mCRC, single-agent immunotherapy has failed, and multiple clinical trials of immunotherapy in combination with other therapies are being actively explored, including combination immunotherapy and immunotherapy combined with targeted therapy, radiotherapy, oncolytic virus, bisspecific antibodies, etc. Pd-1/PD-L1 inhibitors in combination with other immunotherapies may play a synergistic role in enhancing the antitumor effect. In a Phase II trial, durvalumab, a PD-L 1 inhibitor, combined with CTLA-4 inhibitor tremelimumab in refractory mCRC patients (92% pMMR/MSS) showed significant benefits in overall survival (OS) (122). In addition, in 2019, 24 patients with pMMR/MSS colon cancer who had failed standard treatment were included in the REGONIVO study. The ORR reached 33.3% after treatment with regorafenib combined with navulizumab, which significantly improved progression free survival (PFS) and OS (123).
In this study, the low-risk group was found to have higher levels of infiltration of multiple immunocytes, several HLA family members, and multiple immunotherapy targets. In addition, based on the significant role of immunotherapy in tumors, the TIDE algorithm was used to assess the sensitivity of both groups to immunotherapy. The TIDE score was lower in the high-risk group than the low-risk group, suggesting that the immunotherapy response of the high-risk group might be better than that of the low-risk group.
To further enhance the prediction accuracy of the model, a nomogram model based on the risk score prognostic model and clinical indicators (including age and pathological stage), which markedly improved the precision of the model, was established. The time-dependent ROC curve suggested that the risk score had favorable predictive performance for the OS of COAD patients. Calibration curves revealed a good consistency between the model’s estimated 1-, 3-, and 5-year OS and the actual observed OS of patients.
This study had some limitations. First, clinical information and basic experimental verification are lacking. Further, the reliability of the results is dependent on the accuracy of TCGA dataset. In the future, the results of this study should first be verified through clinical trials and basic experiments. Prospective studies are also needed as retrospective studies may be subject to bias. Finally, clinical follow-up data are lacking to prove the accuracy of our prognostic model.
In this study, a prognostic risk model based on NRGs was established, and the diagnostic and predictive significance of the risk model was evaluated. The results of this study will help to reveal the pathogenesis of CRC, enabling the development of new diagnostic ideas, and facilitate the search for new therapeutic targets and prognostic molecular markers.
Data availability statement
The datasets for this study can be found in the Cancer Genome Atlas (TCGA) Genomic Data Commons (GDC) website (https://portal.gdc.cancer.gov/). The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.
Author contributions
WY designed the ideas for this paper. WY, DG, and FM contributed to the writing of the manuscript. SL, LP, ZZ, and YZ contributed to data collation and data analysis. WY, SL, YH, and XC analyzed and interpreted the data. All authors contributed to the article and approved the submitted version.
Funding
This study was supported by Henan Province Medical Science and Technology Tackling Program Joint Co-construction Project (No. LHGJ20200188).
Acknowledgments
We appreciate all participants and the staff of Fusheng Ge, SL, LP, ZZ, YH, and XC, and Affiliated Tumor Hospital of Zhengzhou University.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.1085038/full#supplementary-material
References
1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA: Cancer J Clin (2022) 72(1):7–33. doi: 10.3322/caac.21708
2. Biller LH, Schrag D. Diagnosis and treatment of metastatic colorectal cancer: A review. Jama (2021) 325(7):669–85. doi: 10.1001/jama.2021.0106
3. Dekker E, Tanis PJ, Vleugels JLA, Kasi PM, Wallace MB. Colorectal cancer. Lancet (London England) (2019) 394(10207):1467–80. doi: 10.1016/S0140-6736(19)32319-0
4. Cohen JD, Li L, Wang Y, Thoburn C, Afsari B, Danilova L, et al. Detection and localization of surgically resectable cancers with a multi-analyte blood test. Sci (New York NY) (2018) 359(6378):926–30. doi: 10.1126/science.aar3247
5. Jiang T, Li X, Wang J, Su C, Han W, Zhao C, et al. Mutational landscape of cfDNA identifies distinct molecular features associated with therapeutic response to first-line platinum-based doublet chemotherapy in patients with advanced NSCLC. Theranostics (2017) 7(19):4753–62. doi: 10.7150/thno.21687
6. Osumi H, Shinozaki E, Yamaguchi K, Zembutsu H. Clinical utility of circulating tumor DNA for colorectal cancer. Cancer Science (2019) 110(4):1148–55. doi: 10.1111/cas.13972
7. Osumi H, Shinozaki E, Takeda Y, Wakatsuki T, Ichimura T, Saiura A, et al. Clinical relevance of circulating tumor DNA assessed through deep sequencing in patients with metastatic colorectal cancer. Cancer Med (2019) 8(1):408–17. doi: 10.1002/cam4.1913
8. Holm M, Andersson E, Osterlund E, Ovissi A, Soveri LM, Anttonen AK, et al. Detection of KRAS mutations in liquid biopsies from metastatic colorectal cancer patients using droplet digital PCR, idylla, and next generation sequencing. PLoS One (2020) 15(11):e0239819. doi: 10.1371/journal.pone.0239819
9. Snyder AG, Oberst A. The antisocial network: Cross talk between cell death programs in host defense. Annu Rev Immunol (2021) 39:77–101. doi: 10.1146/annurev-immunol-112019-072301
10. Newton K, Wickliffe KE, Dugger DL, Maltzman A, Roose-Girma M, Dohse M, et al. Cleavage of RIPK1 by caspase-8 is crucial for limiting apoptosis and necroptosis. Nature (2019) 574(7778):428–31. doi: 10.1038/s41586-019-1548-x
11. Yatim N, Jusforgues-Saklani H, Orozco S, Schulz O, Barreira da Silva R, Reis e Sousa C, et al. RIPK1 and NF-κB signaling in dying cells determines cross-priming of CD8+ T cells. Sci (New York NY) (2015) 350(6258):328–34. doi: 10.1126/science.aad0395
12. Wang W, Marinis JM, Beal AM, Savadkar S, Wu Y, Khan M, et al. RIP1 kinase drives macrophage-mediated adaptive immune tolerance in pancreatic cancer. Cancer Cell (2020) 38(4):585–90. doi: 10.1016/j.ccell.2020.09.020
13. Oliver Metzig M, Fuchs D, Tagscherer KE, Gröne HJ, Schirmacher P, Roth W. Inhibition of caspases primes colon cancer cells for 5-fluorouracil-induced TNF-α-dependent necroptosis driven by RIP1 kinase and NF-κB. Oncogene (2016) 35(26):3399–409. doi: 10.1038/onc.2015.398
14. Angelescu R, Dobrescu R. MIDGET:Detecting differential gene expression on microarray data. Comput Methods Programs Biomed (2021) 211:106418. doi: 10.1016/j.cmpb.2021.106418
15. Wang Y, Zhao Y, Bollas A, Wang Y, Au KF. Nanopore sequencing technology, bioinformatics and applications. Nat Biotechnol (2021) 39(11):1348–65. doi: 10.1038/s41587-021-01108-x
16. Czech L, Huerta-Cepas J, Stamatakis A. A critical review on the use of support values in tree viewers and bioinformatics toolkits. Mol Biol Evol (2017) 34(6):1535–42. doi: 10.1093/molbev/msx055
17. Ma J, Chen X, Lin M, Wang Z, Wu Y, Li J. Bioinformatics analysis combined with experiments predicts CENPK as a potential prognostic factor for lung adenocarcinoma. Cancer Cell Int (2021) 21(1):65. doi: 10.1186/s12935-021-01760-y
18. Wang YC, Tian ZB, Tang XQ. Bioinformatics screening of biomarkers related to liver cancer. BMC Bioinf (2021) 22(Suppl 3):521. doi: 10.1186/s12859-021-04411-1
19. Long X, Deng Z, Li G, Wang Z. Identification of critical genes to predict recurrence and death in colon cancer: Integrating gene expression and bioinformatics analysis. Cancer Cell Int (2018) 18:139. doi: 10.1186/s12935-018-0640-x
20. Jia W, He YF, Qian XJ, Chen J. TPMT mRNA expression: A novel prognostic biomarker for patients with colon cancer by bioinformatics analysis. Int J Gen Med (2022) 15:151–60. doi: 10.2147/IJGM.S338575
21. van Eijck MM, Schoonman GG, van der Naalt J, de Vries J, Roks G. Diffuse axonal injury after traumatic brain injury is a prognostic factor for functional outcome: A systematic review and meta-analysis. Brain Injury (2018) 32(4):395–402. doi: 10.1080/02699052.2018.1429018
22. Mayakonda A, Koeffler HP. Maftools: Efficient analysis, visualization and summarization of MAF files from large-scale cohort based cancer studies. (2016). doi: 10.1101/052662
23. Smith JJ, Deane NG, Wu F, Merchant NB, Zhang B, Jiang A, et al. Experimentally derived metastasis gene expression profile predicts recurrence and death in patients with colon cancer. Gastroenterology (2010) 138(3):958–68. doi: 10.1053/j.gastro.2009.11.005
24. Marisa L, de Reyniès A, Duval A, Selves J, Gaub MP, Vescovo L, et al. Gene expression classification of colon cancer into molecular subtypes: characterization, validation, and prognostic value. PLoS Med (2013) 10(5):e1001453. doi: 10.1371/journal.pmed.1001453
25. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47–e. doi: 10.1093/nar/gkv007
26. Qi L, Xu R, Ren X, Zhang W, Yang Z, Tu C, et al. Comprehensive profiling reveals prognostic and immunogenic characteristics of necroptosis in soft tissue sarcomas. Front Immunol (2022) 13:877815. doi: 10.3389/fimmu.2022.877815
27. Hänzelmann S, Castelo R, Guinney JJBb. GSVA. Gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics (2013) 14(1):1–15. doi: 10.1186/1471-2105-14-7
28. Langfelder P, Horvath S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559
29. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinf (Oxford England) (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170
30. Ashburner M, Ball CA, Blake JA, Botstein D, Cherry JMJNG. Gene ontology: Tool for the unification of biology. Gene Ontol Consortium (2000) 25(1):25–9. doi: 10.1038/75556
31. Ogata H, Goto S, Sato K, Fujibuchi W, Kanehisa MJNAR. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res (1999) 27(1):29–34. doi: 10.1093/nar/27.1.29
32. Yu G, Wang L-G, Han Y, He Q-Y. JOajoib. clusterProfiler: an r package for comparing biological themes among gene clusters. OMICS (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118
33. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Systems (2015) 1(6):417–25. doi: 10.1016/j.cels.2015.12.004
34. Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity (2013) 39(4):782–95. doi: 10.1016/j.immuni.2013.10.003
35. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4(1):1–11. doi: 10.1038/ncomms3612
36. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: An R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res (2016) 44(8):e71. doi: 10.1093/nar/gkv1507
37. Mering CV, Huynen M, Jaeggi D, Schmidt S, Bork P, Snel B. STRING: a database of predicted functional associations between proteins. Nucleic Acids Res (2003) 31(1):258–61. doi: 10.1093/nar/gkg034
38. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res (2003) 13(11):2498–504. doi: 10.1101/gr.1239303
39. Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: Identifying hub objects and sub-networks from complex interactome. BMC Syst Biol (2014) 8 Suppl 4(Suppl 4):S11. doi: 10.1186/1752-0509-8-S4-S11
40. Huang HY, Lin YC, Li J, Huang KY, Shrestha S, Hong HC, et al. miRTarBase 2020: Updates to the experimentally validated microRNA-target interaction database. Nucleic Acids Res (2020) 48(D1):D148–d54. doi: 10.1093/nar/gkz896
41. Simon N, Friedman J, Hastie T, Tibshirani RJ. Regularization paths for cox’s proportional hazards model via coordinate descent. J Stat Softw (2011) 39(5):1. doi: 10.18637/jss.v039.i05
42. Islami F, Ward EM, Sung H, Cronin KA, Tangka FKL, Sherman RL, et al. Annual report to the nation on the status of cancer, part 1: National cancer statistics. J Natl Cancer Institute (2021) 113(12):1648–69. doi: 10.1093/jnci/djab131
43. Wang Q, Chen W, Xu X, Li B, He W, Padilla MT, et al. RIP1 potentiates BPDE-induced transformation in human bronchial epithelial cells through catalase-mediated suppression of excessive reactive oxygen species. Carcinogenesis (2013) 34(9):2119–28. doi: 10.1093/carcin/bgt143
44. Seifert L, Werba G, Tiwari S, Giao Ly NN, Alothman S, Alqunaibit D, et al. The necrosome promotes pancreatic oncogenesis via CXCL1 and mincle-induced immune suppression. Nature (2016) 532(7598):245–9. doi: 10.1038/nature17403
45. Park S, Hatanpaa KJ, Xie Y, Mickey BE, Madden CJ, Raisanen JM, et al. The receptor interacting protein 1 inhibits p53 induction through NF-kappaB activation and confers a worse prognosis in glioblastoma. Cancer Res (2009) 69(7):2809–16. doi: 10.1158/0008-5472.CAN-08-4079
46. Ke H, Augustine CK, Gandham VD, Jin JY, Tyler DS, Akiyama SK, et al. CYLD inhibits melanoma growth and progression through suppression of the JNK/AP-1 and β1-integrin signaling pathways. J Invest Dermatol (2013) 133(1):221–9. doi: 10.1038/jid.2012.253
47. Feng X, Song Q, Yu A, Tang H, Peng Z, Wang X. Receptor-interacting protein kinase 3 is a predictor of survival and plays a tumor suppressive role in colorectal cancer. Neoplasma (2015) 62(4):592–601. doi: 10.4149/neo_2015_071
48. McCormick KD, Ghosh A, Trivedi S, Wang L, Coyne CB, Ferris RL, et al. Innate immune signaling through differential RIPK1 expression promote tumor progression in head and neck squamous cell carcinoma. Carcinogenesis (2016) 37(5):522–9. doi: 10.1093/carcin/bgw032
49. Ertao Z, Jianhui C, Kang W, Zhijun Y, Hui W, Chuangqi C, et al. Prognostic value of mixed lineage kinase domain-like protein expression in the survival of patients with gastric caner. Tumour Biol J Int Soc Oncodevelopmental Biol Med (2016) 37(10):13679–85. doi: 10.1007/s13277-016-5229-1
50. Koo GB, Morgan MJ, Lee DG, Kim WJ, Yoon JH, Koo JS, et al. Methylation-dependent loss of RIP3 expression in cancer represses programmed necrosis in response to chemotherapeutics. Cell Res (2015) 25(6):707–25. doi: 10.1038/cr.2015.56
51. Liu X, Zhou M, Mei L, Ruan J, Hu Q, Peng J, et al. Key roles of necroptotic factors in promoting tumor growth. Oncotarget (2016) 7(16):22219–33. doi: 10.18632/oncotarget.7924
52. Wu Y, Dong G, Sheng C. Targeting necroptosis in anticancer therapy: mechanisms and modulators. Acta Pharm Sin B (2020) 10(9):1601–18. doi: 10.1016/j.apsb.2020.01.007
53. Zhao Z, Liu H, Zhou X, Fang D, Ou X, Ye J, et al. Necroptosis-related lncRNAs: Predicting prognosis and the distinction between the cold and hot tumors in gastric cancer. J Oncol (2021) 2021:6718443. doi: 10.1155/2021/6718443
54. Kang CM, Lee JH, Jeon M, Song JS, Kim SO. The effect of MMP-13, MMP-12, and AMBN on gingival enlargement and root deformation in a new type of gingival fibromatosis. J Clin Pediatr Dentistry (2018) 42(1):50–4. doi: 10.17796/1053-4628-42.1.9
55. Cao LQ, Wang YN, Liang M, Pan MZ. CALB1 enhances the interaction between p53 and MDM2, and inhibits the senescence of ovarian cancer cells. Mol Med Rep (2019) 19(6):5097–104. doi: 10.3892/mmr.2019.10212
56. Jin C, Lin T, Shan L. Downregulation of calbindin 1 by miR-454-3p suppresses cell proliferation in nonsmall cell lung cancer in vitro. Cancer Biother Radiopharm (2019) 34(2):119–27. doi: 10.1089/cbr.2018.2598
57. Huang Z, Fan G, Wang D. Downregulation of calbindin 1, a calcium-binding protein, reduces the proliferation of osteosarcoma cells. Oncol Lett (2017) 13(5):3727–33. doi: 10.3892/ol.2017.5931
58. Sugahara K, Kitagawa H. Recent advances in the study of the biosynthesis and functions of sulfated glycosaminoglycans. Curr Opin Struct Biol (2000) 10(5):518–27. doi: 10.1016/S0959-440X(00)00125-1
59. Zhou H, Li Y, Song X, Zhao Y, Cheng L, Zhao L, et al. CHST11/13 regulate the metastasis and chemosensitivity of human hepatocellular carcinoma cells Via mitogen-activated protein kinase pathway. Digest Dis Sci (2016) 61(7):1972–85. doi: 10.1007/s10620-016-4114-5
60. Oliveira-Ferrer L, Heßling A, Trillsch F, Mahner S, Milde-Langosch K. Prognostic impact of chondroitin-4-sulfotransferase CHST11 in ovarian cancer. Tumour Biol J Int Soc Oncodevelopmental Biol Med (2015) 36(11):9023–30. doi: 10.1007/s13277-015-3652-3
61. Huynh KW, Jiang J, Abuladze N, Tsirulnikov K, Kao L, Shao X, et al. CryoEM structure of the human SLC4A4 sodium-coupled acid-base transporter NBCe1. Nat Commun (2018) 9(1):900. doi: 10.1038/s41467-018-03271-3
62. Liu Z, Wang Q, Zhai G, Ke S, Yu X, Guo J. SLC4A4 promotes prostate cancer progression in vivo and in vitro via AKT-mediated signalling pathway. Cancer Cell Int (2022) 22(1):127. doi: 10.1186/s12935-022-02546-6
63. Zhang C, Chang Q, Hu Y, Chang W, Guo X, Fu L, et al. MiR-222-3p promotes the proliferation, migration and invasion of papillary thyroid carcinoma cells through targeting SLC4A4. Histol Histopathol (2021) 36(11):1199–207. doi: 10.14670/HH-18-387
64. Chen X, Chen J, Feng Y, Guan W. Prognostic value of SLC4A4 and its correlation with immune infiltration in colon adenocarcinoma. Med Sci Monit Int Med J Exp Clin Res (2020) 26:e925016. doi: 10.12659/MSM.925016
65. Carreira-Barbosa F, Nunes SC. Wnt signaling: Paths for cancer progression. Adv Exp Med Biol (2020) 1219:189–202. doi: 10.1007/978-3-030-34025-4_10
66. Bienz M. The subcellular destinations of APC proteins. Nat Rev Mol Cell Biol (2002) 3(5):328–38. doi: 10.1038/nrm806
67. Dietinger V, García de Durango CR, Wiechmann S, Boos SL, Michl M, Neumann J, et al. Wnt-driven LARGE2 mediates laminin-adhesive O-glycosylation in human colonic epithelial cells and colorectal cancer. Cell Communication Signaling CCS (2020) 18(1):102. doi: 10.1186/s12964-020-00561-6
68. Cheng R, Sun B, Liu Z, Zhao X, Qi L, Li Y, et al. Wnt5a suppresses colon cancer by inhibiting cell proliferation and epithelial-mesenchymal transition. J Cell Physiol (2014) 229(12):1908–17. doi: 10.1002/jcp.24566
69. Li YD, Lv Z, Zhu WF. RBBP4 promotes colon cancer malignant progression via regulating wnt/β-catenin pathway. World J Gastroenterol (2020) 26(35):5328–42. doi: 10.3748/wjg.v26.i35.5328
70. Salas A, Hernandez-Rocha C, Duijvestein M, Faubion W, McGovern D, Vermeire S, et al. JAK-STAT pathway targeting for the treatment of inflammatory bowel disease. Nat Rev Gastroenterol Hepatol (2020) 17(6):323–37. doi: 10.1038/s41575-020-0273-0
71. Villarino AV, Kanno Y, O'Shea JJ. Mechanisms and consequences of jak-STAT signaling in the immune system. Nat Immunol (2017) 18(4):374–84. doi: 10.1038/ni.3691
72. Yu H, Lee H, Herrmann A, Buettner R, Jove R. Revisiting STAT3 signalling in cancer: new and unexpected biological functions. Nat Rev Cancer (2014) 14(11):736–46. doi: 10.1038/nrc3818
73. Li P, Huang T, Zou Q, Liu D, Wang Y, Tan X, et al. FGFR2 promotes expression of PD-L1 in colorectal cancer via the JAK/STAT3 signaling pathway. J Immunol (Baltimore Md 1950) (2019) 202(10):3065–75. doi: 10.4049/jimmunol.1801199
74. Zhao H, Wu L, Yan G, Chen Y, Zhou M, Wu Y, et al. Inflammation and tumor progression: Signaling pathways and targeted intervention. Signal Transduct Target Ther (2021) 6(1):263. doi: 10.1038/s41392-021-00658-5
75. Wu J, Liu W, Williams JP, Ratner N. EGFR-Stat3 signalling in nerve glial cells modifies neurofibroma initiation. Oncogene (2017) 36(12):1669–77. doi: 10.1038/onc.2016.386
76. Yu H, Pardoll D, Jove R. STATs in cancer inflammation and immunity: A leading role for STAT3. Nat Rev Cancer (2009) 9(11):798–809. doi: 10.1038/nrc2734
77. Zheng B, Artin MG, Chung H, Chen B, Sun S, May BL, et al. Immunogenetics of gastrointestinal cancers: A systematic review and retrospective survey of inborn errors of immunity in humans. J Gastroenterol Hepatol (2022) 37(6):973–82. doi: 10.1111/jgh.15848
78. Kurzejamska E, Sacharczuk M, Landázuri N, Kovtonyuk O, Lazarczyk M, Ananthaseshan S, et al. Effect of chemokine (C-c motif) ligand 7 (CCL7) and its receptor (CCR2) expression on colorectal cancer behaviors. Int J Mol Sci (2019) 20(3):686. doi: 10.3390/ijms20030686
79. Ohlsson L, Hammarström ML, Lindmark G, Hammarström S, Sitohy B. Ectopic expression of the chemokine CXCL17 in colon cancer cells. Br J Cancer (2016) 114(6):697–703. doi: 10.1038/bjc.2016.4
80. Niv Y, Rokkas T. Mucin expression in colorectal cancer (CRC): Systematic review and meta-analysis. J Clin Gastroenterol (2019) 53(6):434–40. doi: 10.1097/MCG.0000000000001050
81. Krishn SR, Kaur S, Smith LM, Johansson SL, Jain M, Patel A, et al. Mucins and associated glycan signatures in colon adenoma-carcinoma sequence: Prospective pathological implication(s) for early diagnosis of colon cancer. Cancer Lett (2016) 374(2):304–14. doi: 10.1016/j.canlet.2016.02.016
82. Reynolds IS, Fichtner M, McNamara DA, Kay EW, Prehn JHM, Burke JP. Mucin glycoproteins block apoptosis; promote invasion, proliferation, and migration; and cause chemoresistance through diverse pathways in epithelial cancers. Cancer Metastasis Rev (2019) 38(1-2):237–57. doi: 10.1007/s10555-019-09781-w
83. Kaur S, Smith LM, Patel A, Menning M, Watley DC, Malik SS, et al. A combination of MUC5AC and CA19-9 improves the diagnosis of pancreatic cancer: A multicenter study. Am J Gastroenterol (2017) 112(1):172–83. doi: 10.1038/ajg.2016.482
84. Lakshmanan I, Rachagani S, Hauke R, Krishn SR, Paknikar S, Seshacharyulu P, et al. MUC5AC interactions with integrin β4 enhances the migration of lung cancer cells through FAK signaling. Oncogene (2016) 35(31):4112–21. doi: 10.1038/onc.2015.478
85. Pereira MB, Dias AJ, Reis CA, Schmitt FC. Immunohistochemical study of the expression of MUC5AC and MUC6 in breast carcinomas and adjacent breast tissues. J Clin Pathol (2001) 54(3):210–3. doi: 10.1136/jcp.54.3.210
86. Krishn SR, Ganguly K, Kaur S, Batra SK. Ramifications of secreted mucin MUC5AC in malignant journey: A holistic view. Carcinogenesis (2018) 39(5):633–51. doi: 10.1093/carcin/bgy019
87. Pothuraju R, Rachagani S, Krishn SR, Chaudhary S, Nimmakayala RK, Siddiqui JA, et al. Molecular implications of MUC5AC-CD44 axis in colorectal cancer progression and chemoresistance. Mol Cancer (2020) 19(1):37. doi: 10.1186/s12943-020-01156-y
88. Perrais M, Rousseaux C, Ducourouble MP, Courcol R, Vincent P, Jonckheere N, et al. Helicobacter pylori urease and flagellin alter mucin gene expression in human gastric cancer cells. Gastric Cancer (2014) 17(2):235–46. doi: 10.1007/s10120-013-0267-5
89. Pinto-de-Sousa J, Reis CA, David L, Pimenta A, Cardoso-de-Oliveira M. MUC5B expression in gastric carcinoma: Relationship with clinico-pathological parameters and with expression of mucins MUC1, MUC2, MUC5AC and MUC6. Virchows Archiv an Int J Pathol (2004) 444(3):224–30. doi: 10.1007/s00428-003-0968-y
90. Mino-Kenudson M, Tomita S, Lauwers GY. Mucin expression in reactive gastropathy: An immunohistochemical analysis. Arch Pathol Lab Med (2007) 131(1):86–90. doi: 10.5858/2007-131-86-MEIRGA
91. van Neerven SM, Vermeulen L. The interplay between intrinsic and extrinsic Wnt signaling in controlling intestinal transformation. Differentiation; Res Biol Diversity (2019) 108:17–23. doi: 10.1016/j.diff.2019.02.002
92. Nie X, Liu Y, Chen WD, Wang YD. Interplay of miRNAs and canonical wnt signaling pathway in hepatocellular carcinoma. Front Pharmacol (2018) 9:657. doi: 10.3389/fphar.2018.00657
93. Pandurangan AK, Divya T, Kumar K, Dineshbabu V, Velavan B, Sudhandiran G. Colorectal carcinogenesis: Insights into the cell death and signal transduction pathways: A review. World J Gastrointest Oncol (2018) 10(9):244–59. 10.4251/wjgo.v10.i9.244
94. Cancer Genome Atlas Network. Comprehensive molecular characterization of human colon and rectal cancer. Nature (2012) 487(7407):330–7. doi: 10.1038/nature11252
95. Ara H, Takagishi M, Enomoto A, Asai M, Ushida K, Asai N, et al. Role for daple in non-canonical wnt signaling during gastric cancer invasion and metastasis. Cancer Science (2016) 107(2):133–9. doi: 10.1111/cas.12848
96. Mao J, Fan S, Ma W, Fan P, Wang B, Zhang J, et al. Roles of wnt/β-catenin signaling in the gastric cancer stem cells proliferation and salinomycin treatment. Cell Death Dis (2014) 5(1):e1039. doi: 10.1038/cddis.2013.515
97. Cohen ED, Tian Y, Morrisey EE. Wnt signaling: an essential regulator of cardiovascular differentiation, morphogenesis and progenitor self-renewal. Dev (Cambridge England) (2008) 135(5):789–98. doi: 10.1242/dev.016865
98. Veeman MT, Axelrod JD, Moon RT. A second canon. functions and mechanisms of beta-catenin-independent wnt signaling. Dev Cell (2003) 5(3):367–77. doi: 10.1016/S1534-5807(03)00266-1
99. Nishioka M, Ueno K, Hazama S, Okada T, Sakai K, Suehiro Y, et al. Possible involvement of Wnt11 in colorectal cancer progression. Mol Carcinog (2013) 52(3):207–17. doi: 10.1002/mc.21845
100. He D, Yue Z, Liu L, Fang X, Chen L, Han H. Long noncoding RNA ABHD11-AS1 promote cells proliferation and invasion of colorectal cancer via regulating the miR-1254-WNT11 pathway. J Cell Physiol (2019) 234(7):12070–9. doi: 10.1002/jcp.27877
101. Sánchez-Vega F, Gotea V, Chen YC, Elnitski L. CpG island methylator phenotype in adenocarcinomas from the digestive tract: Methods, conclusions, and controversies. World J Gastrointest Oncol (2017) 9(3):105–20. doi: 10.4251/wjgo.v9.i3.105
102. Hu Z, Deng H, Qu H. Plasma SFRP5 levels are decreased in Chinese subjects with obesity and type 2 diabetes and negatively correlated with parameters of insulin resistance. Diabetes Res Clin Pract (2013) 99(3):391–5. doi: 10.1016/j.diabres.2012.11.026
103. Miyoshi T, Doi M, Usui S, Iwamoto M, Kajiya M, Takeda K, et al. Low serum level of secreted frizzled-related protein 5, an anti-inflammatory adipokine, is associated with coronary artery disease. Atherosclerosis (2014) 233(2):454–9. doi: 10.1016/j.atherosclerosis.2014.01.019
104. Veeck J, Geisler C, Noetzel E, Alkaya S, Hartmann A, Knüchel R, et al. Epigenetic inactivation of the secreted frizzled-related protein-5 (SFRP5) gene in human breast cancer is associated with unfavorable prognosis. Carcinogenesis (2008) 29(5):991–8. doi: 10.1093/carcin/bgn076
105. Su HY, Lai HC, Lin YW, Liu CY, Chen CK, Chou YC, et al. Epigenetic silencing of SFRP5 is related to malignant phenotype and chemoresistance of ovarian cancer through wnt signaling pathway. Int J Cancer (2010) 127(3):555–67. doi: 10.1002/ijc.25083
106. Zhao C, Bu X, Zhang N, Wang W. Downregulation of SFRP5 expression and its inverse correlation with those of MMP-7 and MT1-MMP in gastric cancer. BMC Cancer (2009) 9:224. doi: 10.1186/1471-2407-9-224
107. Zhou W, Ye C, Li L, Liu L, Wang F, Yu L, et al. Adipocyte-derived SFRP5 inhibits breast cancer cells migration and invasion through wnt and epithelial-mesenchymal transition signaling pathways. Chin J Cancer Res = Chung-kuo yen cheng yen chiu (2020) 32(3):347–60. doi: 10.21147/j.issn.1000-9604.2020.03.06
108. Evans DR, Venkitachalam S, Revoredo L, Dohey AT, Clarke E, Pennell JJ, et al. Evidence for GALNT12 as a moderate penetrance gene for colorectal cancer. Hum Mutation (2018) 39(8):1092–101. doi: 10.1002/humu.23549
109. Lee SH, Hatakeyama S, Yu SY, Bao X, Ohyama C, Khoo KH, et al. Core3 O-glycan synthase suppresses tumor formation and metastasis of prostate carcinoma PC3 and LNCaP cells through down-regulation of alpha2beta1 integrin complex. J Biol Chem (2009) 284(25):17157–69. doi: 10.1074/jbc.M109.010934
110. Iwai T, Kudo T, Kawamoto R, Kubota T, Togayachi A, Hiruma T, et al. Core 3 synthase is down-regulated in colon carcinoma and profoundly suppresses the metastatic potential of carcinoma cells. Proc Natl Acad Sci U States A (2005) 102(12):4572–7. doi: 10.1073/pnas.0407983102
111. Cao P, Wu Y, Sun D, Zhang W, Qiu J, Tang Z, et al. IGF2BP2 promotes pancreatic carcinoma progression by enhancing the stability of B3GNT6 mRNA via m6A methylation. Cancer Med (2022). doi: 10.1002/cam4.5096
112. Xiao S, Yang C, Zhang Y, Lai C. Downregulation of B3GNT6 is a predictor of poor outcomes in patients with colorectal cancer. World J Surg Oncol (2022) 20(1):110. doi: 10.1186/s12957-022-02561-x
113. Qi X, Zhang DH, Wu N, Xiao JH, Wang X, Ma W. ceRNA in cancer: possible functions and clinical implications. J Med Genet (2015) 52(10):710–8. doi: 10.1136/jmedgenet-2015-103334
114. Sun C, Li S, Zhang F, Xi Y, Wang L, Bi Y, et al. Long non-coding RNA NEAT1 promotes non-small cell lung cancer progression through regulation of miR-377-3p-E2F3 pathway. Oncotarget (2016) 7(32):51784–814. doi: 10.18632/oncotarget.10108
115. Shang A, Gu C, Wang W, Wang X, Sun J, Zeng B, et al. Exosomal circPACRGL promotes progression of colorectal cancer via the miR-142-3p/miR-506-3p- TGF-β1 axis. Mol Cancer (2020) 19(1):117. doi: 10.1186/s12943-020-01235-0
116. Wu H, Hu X, Li Y, Chen Q, Sun T, Qiao Y, et al. LNC473 regulating APAF1 IRES-dependent translation via competitive sponging miR574 and miR15b: Implications in colorectal cancer. Mol Ther Nucleic Acids (2020) 21:764–79. doi: 10.1016/j.omtn.2020.07.009
117. Lei X, Lei Y, Li JK, Du WX, Li RG, Yang J, et al. Immune cells within the tumor microenvironment: Biological functions and roles in cancer immunotherapy. Cancer Lett (2020) 470:126–33. doi: 10.1016/j.canlet.2019.11.009
118. Wang Y, Bhave MS, Yagita H, Cardell SL. Natural killer T-cell agonist α-galactosylceramide and PD-1 blockade synergize to reduce tumor development in a preclinical model of colon cancer. Front Immunol (2020) 11:581301. doi: 10.3389/fimmu.2020.581301
119. Wang M, Wang S, Desai J, Trapani JA, Neeson PJ. Therapeutic strategies to remodel immunologically cold tumors. Clin Trans Immunol (2020) 9(12):e1226. doi: 10.1002/cti2.1226
120. Overman MJ, McDermott R, Leach JL, Lonardi S, Lenz HJ, Morse MA, et al. Nivolumab in patients with metastatic DNA mismatch repair-deficient or microsatellite instability-high colorectal cancer (CheckMate 142): An open-label, multicentre, phase 2 study. Lancet Oncol (2017) 18(9):1182–91. doi: 10.1016/S1470-2045(17)30422-9
121. Le DT, Uram JN, Wang H, Bartlett BR, Kemberling H, Eyring AD, et al. PD-1 blockade in tumors with mismatch-repair deficiency. New Engl J Med (2015) 372(26):2509–20. doi: 10.1056/NEJMoa1500596
122. Chen EX, Jonker DJ, Loree JM, Kennecke HF, Berry SR, Couture F, et al. Effect of combined immune checkpoint inhibition vs best supportive care alone in patients with advanced colorectal cancer: The Canadian cancer trials group CO.26 study. JAMA Oncol (2020) 6(6):831–8. doi: 10.1001/jamaoncol.2020.0910
123. Yukami H, Kawazoe A, Lin YT, Koyama S, Fukuoka S, Hara H, et al. Updated efficacy outcomes of anti-PD-1 antibodies plus multikinase inhibitors for patients with advanced gastric cancer with or without liver metastases in clinical trials. Clin Cancer Res an Off J Am Assoc Cancer Res (2022) 28(16):3480–8. doi: 10.1158/1078-0432.CCR-22-0630
Keywords: weighted gene co-expression network analysis (WGCNA), tumor microenvironment, tumor immune infiltrating cells, copy number variation (CNV), nomogram, calibration curves, ceRNA networks
Citation: Yang W, Lu S, Peng L, Zhang Z, Zhang Y, Guo D, Ma F, Hua Y and Chen X (2022) Integrated analysis of necroptosis-related genes for evaluating immune infiltration and colon cancer prognosis. Front. Immunol. 13:1085038. doi: 10.3389/fimmu.2022.1085038
Received: 31 October 2022; Accepted: 08 December 2022;
Published: 22 December 2022.
Edited by:
Muhammad Suleman, University of Swat, PakistanReviewed by:
Farman Ullah, University of Swat, PakistanSusanna Ulahannan, University of Oklahoma, United States
Murad Ali Rahat, University of Swat, Pakistan
Copyright © 2022 Yang, Lu, Peng, Zhang, Zhang, Guo, Ma, Hua and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xiaobing Chen, zlyychenxb0807@zzu.edu.cn
†These authors have contributed equally to this work