Xinyun Li1, Lin Yang2, Wen Wang1, Xiangshu Rao2, Yu Lai2. 1. School of Traditional Chinese Medicine, Sichuan College of Traditional Chinese Medicine, China. 2. School of Basic Medicine, Chengdu University of Traditional Chinese Medicine, Chengdu, China.
Abstract
Colon cancer is a common digestive tract tumor. Although many gene prognostic indicators have been used to predict the prognosis of colon cancer patients, the accuracy of these prognostic indicators is still uncertain. Thus, it is necessary to construct a model for the prognostic analysis of colon cancer. We downloaded the original transcriptome data of colon cancer and performed a differential coexpression analysis of immune-related genes to obtain different immune-related long noncoding RNAs, which were paired as differentially expressed immune-related lncRNA pairs (DEirlncRNAPs). Then, the 1-year overall survival rate receiver operating characteristic curve was calculated, and the Akaike information criterion value was evaluated to determine the maximum inflection point, which was used as the cutoff point to identify groups of colon cancer patients at high and low risk for death. Subsequently, the optimal prediction model was established. Finally, we used the patients' survival times, clinicopathological features, tumor infiltrating immune cells, chemotherapy responses, and immunosuppressive biomarkers to verify the DEirlncRNAP model. Seventy-one DEirlncRNAPs were obtained to build the risk assessment model. The patients were divided into a high-risk group and a low-risk group according to the cutoff point. Then, the DEirlncRNAP model was verified using patient survival times, clinicopathological features, tumor-infiltrating immune cells, chemotherapy responses, and immunosuppressive biomarkers. A new DEirlncRNAP model for predicting the prognosis of colon cancer patients was established, which could reveal new insights into the relationships of colon cancer with tumor-infiltrating immune cells and antitumor immunotherapy.
Colon cancer is a common digestive tract tumor. Although many gene prognostic indicators have been used to predict the prognosis of colon cancer patients, the accuracy of these prognostic indicators is still uncertain. Thus, it is necessary to construct a model for the prognostic analysis of colon cancer. We downloaded the original transcriptome data of colon cancer and performed a differential coexpression analysis of immune-related genes to obtain different immune-related long noncoding RNAs, which were paired as differentially expressed immune-related lncRNA pairs (DEirlncRNAPs). Then, the 1-year overall survival rate receiver operating characteristic curve was calculated, and the Akaike information criterion value was evaluated to determine the maximum inflection point, which was used as the cutoff point to identify groups of colon cancer patients at high and low risk for death. Subsequently, the optimal prediction model was established. Finally, we used the patients' survival times, clinicopathological features, tumor infiltrating immune cells, chemotherapy responses, and immunosuppressive biomarkers to verify the DEirlncRNAP model. Seventy-one DEirlncRNAPs were obtained to build the risk assessment model. The patients were divided into a high-risk group and a low-risk group according to the cutoff point. Then, the DEirlncRNAP model was verified using patient survival times, clinicopathological features, tumor-infiltrating immune cells, chemotherapy responses, and immunosuppressive biomarkers. A new DEirlncRNAP model for predicting the prognosis of colon cancer patients was established, which could reveal new insights into the relationships of colon cancer with tumor-infiltrating immune cells and antitumor immunotherapy.
Colon cancer is a common intestinal epithelial malignant tumor with a high incidence and poor prognosis that seriously endangers human life. According to the statistics of the International Agency for Research on Cancer of the World Health Organization, the global incidence of colon cancer was approximately 1.93 million, and the death toll was approximately 940,000 in 2020.[ More than half of the cases of colon cancer are related to smoking, unhealthy eating habits, alcoholism, physical inactivity, and obesity.[ The early clinical symptoms of colon cancer are hidden and lack highly specific markers. Therefore, it is important to construct a reliable prognostic model for colon cancer patients.Recently, studies have suggested that the immune system plays a key role in the occurrence, development and treatment of colon cancer, especially immune checkpoint inhibitor (ICI) therapy, which is of great significance for prolonging the survival time of patients with colon cancer.[ Pembrolizumab, an inhibitor of PD-1, was used to treat colorectal cancer patients with type dMMR-MSI-H cancer in the KEYNOTE(KN)-016 trial. The results indicated that the 2-year overall survival (OS) rate was 72%.[ Navuzumab, another PD-1 inhibitor, was also shown to have a significant effect in patients with colorectal cancer in multicenter, open-label II clinical trials.[ In addition, the combination of navuzumab and the CTLA-4 inhibitor ipilimumab achieved 84% control of colon cancer and an 83% OS rate at 12 months, and the degree of side effects associated with the treatment was acceptable.[Long noncoding ribonucleic acids (lncRNAs) are transcripts over 200 nucleotides in length that have no protein coding function. Because lncRNAs are closely related to the occurrence and development of tumors and have obviously abnormal expression levels in tumor tissues and precancerous lesions,[ lncRNAs are ideal tumor markers with high specificity and sensitivity that are easy to repeatedly detect.[A recent study showed that the immune-related lncRNA (irlncRNA) signature possessed independent prognostic significance in colon cancer.[ According to a previous report by Hong et al,[ we screened differentially expressed immune-related lncRNA pairs (DEirlncRNAPs), established and verified a risk assessment model for colon cancer, and analyzed the relationships between the paired DEirlncRNAPs, tumor infiltrating immune cells and clinical treatment outcomes.
2. Methods
2.1 Data collection and preprocessing
The main process of this research is demonstrated in Figure 1. Transcriptome profiling (RNAseq) data and fragments per kilobase per million (FPKM) data from 449 patients with colon adenocarcinoma (COAD) and clinical data from 452 COAD patients were downloaded from The Cancer Genome Atlas (TCGA) database (https://tcga-data.nci.nih.gov/tcga/). Clinical data with no survival times available and duplicate data were deleted. The gene transfer format (GTF) annotation file was downloaded from Ensembl (http://asia.ensembl.org). Perl language was used to distinguish mRNA and lncRNA for subsequent analysis. A total of 2483 identified immune-related genes (ir-genes) were downloaded from ImmPort (http://www.immport.org). Then, we performed a coexpression analysis of the ir-genes and the previously obtained lncRNAs to identify irlncRNAs. We used the limma R package to analyze the differential expression of irlncRNAs and set the filter parameters as log fold-change (FC) > 1.5 and false discovery rate (FDR) < 0.05 to acquire differentially expressed irlncRNAs (DEirlncRNAs).
Figure 1.
The main flow chart of this study. AUCs = area under the curves, DEirlncRNAs = differentially expressed irlncRNAs, GTF = gene transfer format, irlncRNA = immune-related lncRNA, lncRNAs = Long noncoding ribonucleic acids, ROC = receiver operating characteristic, TCGA = The Cancer Genome Atlas.
The main flow chart of this study. AUCs = area under the curves, DEirlncRNAs = differentially expressed irlncRNAs, GTF = gene transfer format, irlncRNA = immune-related lncRNA, lncRNAs = Long noncoding ribonucleic acids, ROC = receiver operating characteristic, TCGA = The Cancer Genome Atlas.
2.2 Screening DEirlncRNAPs
We paired the DEirlncRNAs and called the paired lncRNA pairs DEirlncRNAP. Then, we used DEirlncRNAPs to build a 0 or 1 matrix. According to Hong et al’s algorithm, by comparing the expression of 2 lncRNAs in DEirlncRNAP, which can predict the prognosis of patients with COAD, was selected to minimize the error caused by the change in expression.
2.3 Establishment of a risk assessment model of colon cancer
First, we performed univariate analysis of the validated DEirlncRNAPs and used the glmnet R package for Lasso regression with 10-fold cross validation, and a P value < .05 was considered statistically significant. The number of iterations was set to 1000, and 1000 random stimuli were set for each cycle. Additionally, the frequency of each DEirlncRNAP in 1000 Lasso regression cycles was recorded, and DEirlncRNAPs with frequencies over 100 were subjected to Cox proportional hazard regression analysis. Then, we calculated the coefficient of DEirlncRNAPs and constructed the risk model. We drew receiver operating characteristic (ROC) curves of the 1-, 2-, and 3-year OS rates of the patients and calculated the area under the curves (AUCs). If the curve met the point pointing to the maximum AUC, the model was regarded as the best candidate, and we ended the calculation process. We used the following formula to determine the risk score of all clinical cases of COAD in the TCGA database: RiskScore = . Finally, Akaike information criterion (AIC) values were calculated to evaluate each point of the 1-year ROC curve to determine the maximum inflection point. We divided COAD patients into a high-risk group and a low-risk group according to the cutoff value.
2.4 Verification of the DEirlncRNAP risk assessment model
We used Kaplan–Meier analysis to show the difference in survival between patients in the high-risk and low-risk groups and R language to evaluate the prognostic value of the cutoff point. The R packages used in these operations are survival, glmnet, survivalROC, and survminer.Then, we analyzed the relationship between the DEirlncRNAP model and clinicopathological features by the chi-square test and drew a band diagram to show the results. The markers were as follows: *P < .05, **P < .01, ***P < .001. The Wilcoxon signed rank test was used to calculate the difference in risk scores in different groups with clinicopathological features, and box diagrams were drawn for visualization. Finally, univariate and multivariate Cox proportional hazards regressions were used to compare the DEirlncRNAP model with other clinicopathological features in predicting the prognosis of COAD patients and to evaluate the accuracy of the DEirlncRNAP model. We drew forest maps to show the results. The R packages used in these steps include limma, ComplexHeatmap, survivalROC, and ggpubr.
2.5 Study of the relationship between the DEirlncRNAP model and tumor-infiltrating immune cells
We analyzed the relationship between the model and immune cells by the Wilcoxon signed rank test and used 7 common methods to calculate the immune infiltrating cells of COAD samples in the TCGA database, including CIBERSORT,[ CIBERSORT-ABS,[ XCELL,[ QUANTISEQ,[ MCPCOUNTER,[ EPIC,[ and TIMER.[ The results are displayed as box plots. Then, a detailed Spearman correlation analysis was performed to study the relationship between immune cells and the risk score. The results are shown as a bubble chart. We set the threshold value as P < .05. These steps were completed in R language, and the R packages used included limma, scales, ggplot2, and ggtext.
2.6 Study of the relationship between the DEirlncRNAP model and clinical treatment
To improve the clinical relevance of the DEirlncRNAP model, we studied the relationship between the model and genes related to immunosuppression checkpoints and drew a violin diagram. This operation used the R packages gamma and ggpubr. In addition to ICI therapy, chemotherapy is also a common treatment for colon cancer. To study the correlation between the risk assessment model and the response to chemotherapy, we calculated the inhibitory concentration 50 (IC50) of some commonly used chemotherapy drugs for patients with colon cancer. These drugs include efitinib, thapsigargin, vinorelbine, 5-fluorouracil, and irinotecan. We used the Wilcoxon signed rank test to compare the IC50 values between the high- and low-risk groups, and the results are shown as box charts. In this step, a ridge regression model was constructed by using the R package pRRophetic to predict the IC50 value of chemotherapy drugs.
3. Results
3.1 Obtaining DEirlncRNAs
In this study, we collected transcriptome data of 514 patients with COAD from the TCGA database, consisting of 41 normal and 473 tumor samples. We used the GTF of Ensemble to annotate the data and analyzed the coexpression of known ir-genes and lncRNAs to obtain irlncRNAs. Then, we analyzed the differential expression of the irlncRNAs. We obtained 894 irlncRNAs (Table S1, Supplemental Digital Content 1, http://links.lww.com/MD/H261) and 227 DEirlncRNAs (Fig. 2A), of which the expression of 208 lncRNAs was upregulated and the expression of 19 lncRNAs was downregulated (Fig. 2B).
Figure 2.
Screening DEirlncRNAs and DEirlncRNAPs. The DEirlncRNAs selected from 473 COAD tumor and 41 normal samples from TCGA were analyzed by heatmap (A) and volcano plots (B). (C) The forest map shows 34 DEirlncRNAPs obtained by multivariate stepwise regression analysis, which were used to construct a Cox proportional hazards model. COAD = colon adenocarcinoma, DEirlncRNAPs = differentially expressed immune-related lncRNA pairs, DEirlncRNAs = differentially expressed irlncRNAs, TCGA = The Cancer Genome Atlas.
Screening DEirlncRNAs and DEirlncRNAPs. The DEirlncRNAs selected from 473 COAD tumor and 41 normal samples from TCGA were analyzed by heatmap (A) and volcano plots (B). (C) The forest map shows 34 DEirlncRNAPs obtained by multivariate stepwise regression analysis, which were used to construct a Cox proportional hazards model. COAD = colon adenocarcinoma, DEirlncRNAPs = differentially expressed immune-related lncRNA pairs, DEirlncRNAs = differentially expressed irlncRNAs, TCGA = The Cancer Genome Atlas.
3.2 Obtaining DEirlncRNAPs and constructing the risk assessment model
We conducted a pairwise comparison of the DEirlncRNAs and screened the 0-or-1 matrix and obtained 14,874 pairs of valid DEirlncRNAs. Then, we used modified Lasso regression analysis and single factor analysis and extracted 71 DEirlncRNAPs (Table 1). We used a stepwise method and obtained 34 DEirlncRNAPs to construct a Cox proportional hazards model (Fig. 2C). Next, we used the 71 DEirlncRNAPs to draw 1-, 2- and 3-year survival rate ROCs and calculated the AUCs. All AUCs were above 0.95 (Fig. 3A). The maximum AUC of the 1-year ROC curve was 0.955 (Fig. 3B), which was compared with other clinical features of COAD patients (Fig. 3C). Finally, we used AIC values to determine the cutoff point of the 1-year survival rate ROC curve (Fig. 3D). We evaluated the clinical data of 426 COAD patients collected from TCGA. Then, we calculated the risk scores of all patients and divided them into groups. Those with scores higher than the cutoff value were categorized as high-risk; otherwise, they were categorized as low-risk. As a result, we stratified 367 patients into the low-risk group and into the 59 high-risk group (Table S2, Supplemental Digital Content 2, http://links.lww.com/MD/H262).
Table 1
Seventy one DEirlncRNAPs for establishing risk assessment models.
Constructing the risk assessment model with DEirlncRNAPs. (A) For the 1-, 2-, and 3-year survival rate ROC curves, the AUCs were all >0.95. (B) The AUC value of the 1-year survival rate ROC curve was 0.955. (C) Comparison of the DEirlncRNAP model with ROC curves of other clinical features shows that the DEirlncRNAP model has the highest accuracy. (D) AIC was used to determine the maximum inflection point of the 1-year survival rate ROC curve, and the cutoff point was 46.974. AIC = Akaike information criterion, AUCs = area under the curves, DEirlncRNAPs = differentially expressed immune-related lncRNA pairs, ROC = receiver operating characteristic.
Seventy one DEirlncRNAPs for establishing risk assessment models.DEirlncRNAPs = differentially expressed immune-related lncRNA pairs.Constructing the risk assessment model with DEirlncRNAPs. (A) For the 1-, 2-, and 3-year survival rate ROC curves, the AUCs were all >0.95. (B) The AUC value of the 1-year survival rate ROC curve was 0.955. (C) Comparison of the DEirlncRNAP model with ROC curves of other clinical features shows that the DEirlncRNAP model has the highest accuracy. (D) AIC was used to determine the maximum inflection point of the 1-year survival rate ROC curve, and the cutoff point was 46.974. AIC = Akaike information criterion, AUCs = area under the curves, DEirlncRNAPs = differentially expressed immune-related lncRNA pairs, ROC = receiver operating characteristic.
3.3 Clinical evaluation based on the DEirlncRNAP risk assessment model
We combined the risk scores of all COAD patients with survival time and visualized the results (Fig. 4A and B) and found that the survival time and number of survivors in the low-risk group were significantly higher than those in the high-risk group. Kaplan–Meier analysis also indicated that the survival time of patients in the low-risk group was better than that of patients in the high-risk group, and the difference was statistically significant (P < .001; Fig. 4C).
Figure 4.
Prognostic prediction ability of the DEirlncRNAP model. (A) The risk score of each COAD patient; those with a risk score >46.974 were categorized as high-risk; otherwise, they were categorized as low-risk. (B) Comparison of survival time between the high- and low-risk groups. (C) Kaplan–Meier analysis indicated that the survival time of patients in the high-risk group was shorter than that of patients in the low-risk group (P < .001). COAD = colon adenocarcinoma, DEirlncRNAPs = differentially expressed immune-related lncRNA pairs.
Prognostic prediction ability of the DEirlncRNAP model. (A) The risk score of each COAD patient; those with a risk score >46.974 were categorized as high-risk; otherwise, they were categorized as low-risk. (B) Comparison of survival time between the high- and low-risk groups. (C) Kaplan–Meier analysis indicated that the survival time of patients in the high-risk group was shorter than that of patients in the low-risk group (P < .001). COAD = colon adenocarcinoma, DEirlncRNAPs = differentially expressed immune-related lncRNA pairs.According to the risk score of the DEirlncRNAP model, we conducted a sequence of chi-square tests to study the relationship between clinicopathological features and risk in COAD patients. We obtained a band diagram (Fig. 5A). We obtained a group of scatter plots using the Wilcoxon signed rank test. The results suggested that risk was related to clinical stage (Fig. 5B), T stage (Fig. 5C), M stage (Fig. 5D), N stage (Fig. 5E) and age (Fig. 5F). Univariate Cox regression analysis (Fig. 5G) showed that the differences in patient age (P = .04, HR = 1.021, 95% CI [1.001–1.042]), clinical stage (P < .001, HR = 2.222, 95% CI [1.714–2.881]) and risk score (P < .001, HR = 1.001, 95% CI [1.000–1.001]) were statistically significant. Multivariate Cox regression analysis (Fig. 5H) showed significant differences in age (P = .012, HR = 1.027, 95% CI [1.006–1.048]), clinical stage (P < .001, HR = 2.417, 95% CI [1.851–3.155]) and risk score (P < .001, HR = 1.001, 95% CI [1.000–1.001]). Univariate and multivariate Cox regression analyses indicated that the risk score could be used as an independent prognostic indicator.
Figure 5.
Clinical validation of the DEirlncRNAP risk assessment model. (A) Band diagram indicating that the risk score was significantly correlated with clinical stage, T stage, M stage and N stage. Scatter charts show that the risk score is related not only to clinical stage (B), T stage (C), M stage (D), and N stage (E) but also to patient age (F). (G) The forest map suggested that there were significant differences in age, clinical stage, and risk score in the univariate Cox hazard ratio analysis. (H) The multivariate Cox hazard ratio analysis indicated that there were significant differences in age, clinical stage, and risk score. DEirlncRNAPs = differentially expressed immune-related lncRNA pairs.
Clinical validation of the DEirlncRNAP risk assessment model. (A) Band diagram indicating that the risk score was significantly correlated with clinical stage, T stage, M stage and N stage. Scatter charts show that the risk score is related not only to clinical stage (B), T stage (C), M stage (D), and N stage (E) but also to patient age (F). (G) The forest map suggested that there were significant differences in age, clinical stage, and risk score in the univariate Cox hazard ratio analysis. (H) The multivariate Cox hazard ratio analysis indicated that there were significant differences in age, clinical stage, and risk score. DEirlncRNAPs = differentially expressed immune-related lncRNA pairs.
3.4 Study of tumor-infiltrating immune cells with the DEirlncRNAP model
The DEirlncRNAP model is related to ir-genes, so we used 7 different online analysis platforms to predict the enrichment of immune cells in different risk groups. The results of the Wilcoxon signed rank test showed that high risk negatively correlated with the infiltration of plasmacytoid dendritic cells and CD8+ T cells (Figure S1, Supplemental Digital Content 3, http://links.lww.com/MD/H263). Then, we carried out a Spearman correlation analysis and found that high risk positively correlated with hematopoietic stem cells, macrophages, monocytes, endothelial cells and cancer-associated fibroblast cells and negatively correlated with class switch memory B cells, plasmacytoid dendritic cells, CD8+ T cells and mast cells (Fig. 6) (Table S3, Supplemental Digital Content 4, http://links.lww.com/MD/H264).
Figure 6.
Analysis of tumor infiltrating cells in the DEirlncRNAP model. The results of Spearman correlation analysis suggested that high risk positively correlated with hematopoietic stem cells, macrophages, monocytes, endothelial cells, and cancer-associated fibroblast cells and negatively correlated with class switch memory B cells, plasmacytoid dendritic cells, CD8 + T cells and mast cells. DEirlncRNAPs = differentially expressed immune-related lncRNA pairs.
Analysis of tumor infiltrating cells in the DEirlncRNAP model. The results of Spearman correlation analysis suggested that high risk positively correlated with hematopoietic stem cells, macrophages, monocytes, endothelial cells, and cancer-associated fibroblast cells and negatively correlated with class switch memory B cells, plasmacytoid dendritic cells, CD8 + T cells and mast cells. DEirlncRNAPs = differentially expressed immune-related lncRNA pairs.
3.5 Study of the relationship between clinical treatment and the DEirlncRNAP model
Because ICIs have been shown to have a beneficial effect on the treatment of colon cancer in clinical trials, we studied whether the DEirlncRNAP model is related to immunosuppressive point biomarkers. The results showed that there was a positive correlation between high risk and high expression of CTLA-4 (P < .05; Fig. 7A). BTLA (P > .05; Fig. 7B), HAVCR-2 (P > .05; Fig. 7C), and LAG-3 (P > .05; Fig. 7D) expression levels were not significantly different. In addition to ICI therapy, chemotherapy is also a common treatment for colon cancer; thus, we studied the sensitivity of patients with high-risk colon cancer to common chemotherapeutic drugs. The results showed that the IC50 of gefitinib (P = .019), thapsigargin (P = .025) and vinorelbine (P = .012) negatively correlated with high risk, but there was no significant difference with 5-fluorouracil (P = .42) and irinotecan (P = .9). The results showed that the risk score model can also be used to predict drug sensitivity in COAD patients. (Fig. 7E) (Table S4, Supplemental Digital Content 5, http://links.lww.com/MD/H265).
Figure 7.
Clinical treatment analysis of the DEirlncRNAP model. (A) Violin plot showing that high risk negatively correlated with CTLA-4. (B-D) There was no significant difference in the expression levels of BTLA (B), HAVCR-2 (C), and LAG-3 (D) between the high- and low-risk groups. (E) Box plots showed that a high-risk score was associated with a lower IC50 of chemotherapy drugs such as gefitinib, thapsigargin and vinorelbine but not 5-fluorouracil or irinotecan. DEirlncRNAPs = differentially expressed immune-related lncRNA pairs.
Clinical treatment analysis of the DEirlncRNAP model. (A) Violin plot showing that high risk negatively correlated with CTLA-4. (B-D) There was no significant difference in the expression levels of BTLA (B), HAVCR-2 (C), and LAG-3 (D) between the high- and low-risk groups. (E) Box plots showed that a high-risk score was associated with a lower IC50 of chemotherapy drugs such as gefitinib, thapsigargin and vinorelbine but not 5-fluorouracil or irinotecan. DEirlncRNAPs = differentially expressed immune-related lncRNA pairs.
4. Discussion
As a common malignant tumor of the digestive tract, colon cancer has the characteristics of a high incidence rate, frequent postoperative recurrence and high mortality. At present, pathophysiological evaluation, treatment decision-making and prognosis evaluation of colon cancer mainly depend on cancer cell-centered evaluation factors, such as the TNM staging system, and serum molecular markers, such as CEA and CA19-9.[ Although many gene prognostic indicators have been used to predict the prognosis of patients with colon cancer, their accuracy is still uncertain.[ Therefore, it is urgent to explore reliable prognostic biomarkers to predict the survival of COAD patients.Interestingly, lncRNAs can regulate the expression of coding genes by affecting adjacent genes or distant genes on other chromosomes.[ Because of the obvious abnormal expression of lncRNAs in the tumors of patients with colon cancer, lncRNAs such as SNHG7,[ SNHG11,[ and POU6F2[ are highly specific and sensitive potential biomarkers. As a specific biomarker of colon cancer, lncRNAs have prominent effects on diagnosis, curative effect prediction and prognosis. Importantly, it was suggested that the accuracy of a cancer diagnosis model with a combination of multiple biomarkers was significantly better than that of a single biomarker prediction model.[ Due to the biological heterogeneity of tumors and the technical deviation caused by the use of different sequencing platforms, previous risk models of colon cancer need to be standardized to gene expression profiles. To achieve robustness of prediction, we used a new method without the need to consider the data deviation arising from the use of different platforms. Moreover, the new analysis method does not require the specific expression level of lncRNAs or the scaling and normalization of the expression levels of lncRNAs but rather uses relative ranking and paired comparisons of lncRNA expression levels. This method obtained reliable results in a previous study of human hepatocellular carcinoma.[In this study, we identified a prognostic DEirlncRNAP model to predict the prognosis of COAD patients. The prognostic model consisted of 71 DEirlncRNAPs, including 90 independent DEirlncRNAs. It has been shown that lncRNAs are key regulators of the immune response, participate in gene activation and regulate the immunophenotype.[ With the remarkable achievements of immunotherapy for cancer treatment, irlncRNAs have gradually become a new hot spot and have been found to be prognostic factors for breast cancer,[ renal cancer[ and pancreatic cancer.[A model for predicting the prognosis of colon cancer with 10 lncRNAs was established, which suggested that lncRNAs with high expression levels have important biological functions.[ Wu et al[ created a prognostic model composed of 12 irlncRNAs, and the internal verification C index of the nomogram was 0.807. We obtained DEirlncRNAPs by differential coexpression analysis and established a DEirlncRNAP model by modified Lasso penalized regression. This algorithm only requires analysis of the expression level of paired lncRNAs, and it is not necessary to calculate the specific expression level of each lncRNA, so it has higher clinical practical value than other models. Some DEirlncRNAs involved in our study have been found to play an important role in the occurrence, progression and metastasis of various cancers, such as CDKN2B-AS1 in lung cancer and cervical cancer,[ MNX1-AS1 in esophageal squamous cell carcinoma and gastric cancer,[ B4GALT1-AS1 in colon cancer,[ NKILA in hepatocellular carcinoma and breast cancer[ and MIR17HG in glioma and osteosarcoma.[ Among them, the expression of NKILA in colorectal cancer was lower than that in normal tissues and adenomas.[ MNX1-AS1 regulates the miR-218-5p/SEC61A1 axis to promote the development of colon cancer.[ Both B4GALTI-AS1 and MIR17HG promote the migration of colon cancer cells.[ The abnormal expression of these lncRNAs facilitates the development of colon cancer. We are the first to identify the other lncRNAs used to build this model that may be associated with cancer and could be novel biomarkers for further colon cancer research.We used 7 common methods to evaluate the difference in infiltrating immune cells between the high- and low-risk groups, including XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT. Because the algorithms of each platform are complex and different, we integrated the results. The results suggested that high risk positively correlated with hematopoietic stem cells, macrophages, monocytes, endothelial cells and cancer-associated fibroblast cells and negatively correlated with class-switch memory B cells, plasmacytoid dendritic cells, CD8+ T cells and mast cells. It has been reported that endothelial cells can form new tumor blood vessels and provide nutrition for tumors.[ Macrophages induce the proliferation and migration of colon cancer cells.[ Monocytes have also been found to be closely associated with liver metastasis of colon cancer.[ Plasmacytoid dendritic cells produce type I interferon, which has a strong immune response to tumors and contributes to the establishment of an immunosuppressive tumor microenvironment.[ Notably, tumor-infiltrating CD8+T cells often indicate a better immunotherapy response and prognosis.[ It was also found that the prognosis of high-risk patients with tumor-infiltrating immune cells was poor in the DEirlncRNAP model. Galon et al[ proposed that an immune score based on tumor immune cell density analysis can predict the prognosis of patients and that this method is more accurate than TNM analysis.This immune score also contributes to individualized immunotherapy and chemotherapy in the clinic[; thus, we analyzed the IC50 of chemotherapy drugs and immune checkpoint inhibitors in the high-risk group. The results showed that the high-risk group was more sensitive to gefitinib, thapsigargin and vinorelbine than to 5-fluorouracil and irinotecan. Liu et al[ found that vinorelbine can inhibit the metastasis of colon cancer cells by reducing metastasis, invasion and epithelial-mesenchymal transition. This model may provide a new direction for the further study of chemotherapeutic drugs for colon cancer. However, we believe that immunotherapy can eliminate cancer cells and produce new antigens, which may have an advantage over chemotherapy in inhibiting tumor development. Although the results of our model and immune checkpoint inhibitor analysis only showed that the difference in the expression of CTLA4 was statistically significant, the tumor-infiltrating immune cells and immunophenotype of the patients were different, and the expression of the verified biomarkers should also be different. Due to the heterogeneity of tumors, we need to effectively and accurately identify the tumor type, tumor invasion matrix, and molecular and functional characteristics of COAD, identify more biomarkers, develop different immunotherapy regimens, and carry out personalized precision treatment.However, there are still some deficiencies in our research. The TCGA data we adopted have some defects, mainly because the data are only from a small number of countries, the sample size is limited, the data are not updated in a timely manner, and the sequencing technologies and quality control methods are different across the samples, such that all of which have a certain impact on the accuracy of the data. Moreover, we only searched the RNA-seq dataset, which does not fully represent the status of the tumor genes, and we were not able to collect the lncRNA expression level data of the colon cancer patients. The same dataset was used to build and verify the model. The prognostic efficacy and potential mechanism of this model still need to be externally verified through the analysis of real clinical data and basic research. Therefore, we used a 0-or-1 matrix to screen the qualified DEirlncRNAPs without calculating the expression level of lncRNAs to minimize sample error. In addition, the algorithm was verified through several parameters, such as the survival time of the patients, clinicopathological features, tumor infiltrating immune cells and clinical therapeutics. Nevertheless, external validation with clinical data is still needed to confirm the reliability of the model. Therefore, more functional studies and in vitro and in vivo experiments should be carried out to test the accuracy of the DEirlncRNAP risk assessment model to improve its clinical applicability.
5. Conclusion
In summary, the findings of this study provide an immune-related lncRNA model that does not require the prediction of lncRNA expression levels to predict prognosis for patients with colon cancer and increase the understanding of the relationships of colon cancer with tumor infiltrating immune cells and antitumor immunotherapy.
Author contributions
Conceptualization: Lin Yang.Data curation: Xinyun Li.Formal analysis: Xinyun Li, Lin Yang.Methodology: Xiangshu Rao.Validation: Wen Wang.Visualization: Xinyun Li, Lin Yang, Yu Lai.Writing – original draft: Xinyun Li, Lin Yang.Writing – review & editing: Yu Lai.
Authors: Ewan A Gibb; Katey S S Enfield; Greg L Stewart; Kim M Lonergan; Raj Chari; Raymond T Ng; Lewei Zhang; Calum E MacAulay; Miriam P Rosin; Wan L Lam Journal: Oral Oncol Date: 2011-08-10 Impact factor: 5.337
Figure 1.
Figure 2.
Figure 3.
Figure 4.
Figure 5.
Figure 6.
Figure 7.