A novel m7G-related lncRNA risk model for predicting prognosis and evaluating the tumor immune microenvironment in colon carcinoma

Abstract

N7-Methylguanosine (m7G) modifications are a common type of posttranscriptional RNA modifications. Its function in the tumor microenvironment (TME) has garnered widespread focus in the past few years. Long non-coding RNAs (lncRNAs) played an essential part in tumor development and are closely associated with the tumor immune microenvironment. In this study, we employed a comprehensive bioinformatics approach to develop an m7G-associated lncRNA prognostic model based on the colon adenocarcinoma (COAD) database from The Cancer Genome Atlas (TCGA) database. Pearson's correlation analysis was performed to identify m7G-related lncRNAs. Differential gene expression analysis was used to screen lncRNAs. Then, we gained 88 differentially expressed m7G-related lncRNAs. Univariate Cox analysis and Lasso regression analysis were performed to build an eight-m7G-related-lncRNA (ELFN1-AS1, GABPB1-AS1, SNHG7, GS1-124K5.4, ZEB1-AS1, PCAT6, C1RL-AS1, MCM3AP-AS1) risk model. Consensus clustering analysis was applied to identify the m7G-related lncRNA subtypes. We also verified the risk prediction effect of a gene signature in the GSE17536 test set (177 patients). A nomogram was constructed to predict overall survival rates. Furthermore, we analyzed differentially expressed genes (DEGs) between high-risk and low-risk groups. Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were conducted with the analyzed DEGs. At last, single-sample gene set enrichment analysis (ssGSEA), CIBERSORT, MCP-COUNTER, and Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithms were utilized to discover the relationship between the risk model and the TME. Consequently, the m7G-related lncRNA risk model for COAD patients could be a viable prognostic tool and treatment target.

Keywords: N7-methylguanosine(m7G); colon carcinoma; long noncoding RNA(LncRNA); risk model; tumor immune microenvironment.

Figure 1

The flowchart of the overall study design.

Figure 2

Differentially expressed m7G-related genes. (A) Violin plot showing the differential expression of 27 m7G-related genes between tumor and normal tissues from the COAD. (B) Heat map of 21 differentially expressed m7G-related genes between tumor and normal tissues (P < 0.05).The PPI network (C) and number of interaction nodes (D) of 21 differentially expressed m7G-related genes. (E) Pearson correlation analysis of 21 differentially expressed m7G-related genes. The red color represents a positive correlation; the blue color represents a negative correlation. *P < 0.05, **P < 0.01, and ***P < 0.001.

Figure 3

Construction of the prognostic m7G-related lncRNAs risk model. (A) Venn diagram of cor-lncRNA, diff-lncRNA, and GSE17536lncRNA. (B) Forest map of 11 prognostic m7G-related lncRNAs by univariate Cox analysis (P < 0.05). (C) The correlation between 21 differentially expressed m7G-related genes and 11 prognostic m7G-related lncRNAs. The red color represents a positive correlation; the blue color represents a negative correlation. *P < 0.05, **P < 0.01, and ***P < 0.001. (D) 1,000 cross-validation to determine the optimal penalty parameter lambda (λ). (E) Lasso regression of the 11 m7G-related lncRNAs. (F) The Sankey diagram displayed the relationship between the m7G regulators mRNA expression and the m7G-related lncRNAs.

Figure 4

Identification of m7G-associated clusters and prognostic analysis between clusters. (A) The consensus matrix (k = 3) of 452 COAD samples by Consensus Cluster analysis. (B) The relative change in area under the CDF curve for k = 2–10. The KM plot showing overall survival in three clusters (C), training set (D), and test set (F).The ROC curve of the training set (D) and test set (F). The 1-, 3-, and 5-year ROC analyses of risk score in the training set (E) and test set (G).

Figure 5

Validation of the prognostic risk model. Scatter plot revealing the risk score distribution of high risk and low risk and the relationship between survival time and risk score based on the training set (A) and test set (B). Heat map displaying the differential expression of the eight prognostic m7G-related lncRNAs in the high- or low-risk group. (C) Heat map showing clinicopathological features (TMN stage, stage, age, gender) and differences in the expression of eight lncRNAs in the high- and low-risk groups. *P < 0.05, **P < 0.01, and ***P < 0.001. The 3D scatter plot of PCA results of the training set (D) and test set (F).The t-SNE analysis of the training set (E) and test set (G).

Figure 6

Survival analysis and construction of a nomogram. (A) Survival analysis in subgroups including gender, age, and tumor stages. Univariate Cox regression analysis revealing the association between patients’ overall survival and clinicopathological parameters along with m7G-related lncRNA risk scores in the training set (B) and test set (D). Multivariate Cox regression analysis uncovering independent prognostic factors in the training set (C) and test set (E). (F) Nomogram depending on the m7G-related lncRNA risk score and other clinicopathologic feature predicting the 1-, 3-, and 5-year overall survival for COAD patients. (G) Calibration curves illustrating the consistency between predicted and observed 1-, 3-, and 5-year overall survival rates in COAD patients based on the nomogram.

Figure 7

Gene set enrichment analysis. DCA curve (A) and ROC curve (B) of the nomogram, risk, and other clinicopathologic feature in COAD. (C) GSEA results illustrating 10 significant enrichment of GO in low-risk and high-risk groups. The results of GO enrichment analysis of the differentially expressed genes shown by barplot (D), bubble chart (E), and chord diagram (F).

Figure 8

Immune infiltration analysis of the prognostic m7G-related lncRNA risk model. The infiltrating levels of 16 immune cell types (A) and 13 immune functions (B) in high-risk and low-risk groups estimated by ssGSEA. (C) The correlation of immune score and risk score calculated by CIBERSORT. (D) The violin diagram revealing the abundance of 10 types of immune and stromal cells between two groups via MCPcounter. (E) Heat map of 29 immune cells and functions displaying the difference of the immune score, stromal score, estimated score, and tumor purity in two groups through ESTIMATE. (F) Box plot of 24 MHC molecules’ expression level in two groups. ns, not significant, *P < 0.05, **P < 0.01, and ***P < 0.001.