NK cell-associated long non-coding RNAs reveal heterogeneity of colorectal cancer immune microenvironment

Abstract

Introduction: Individuals diagnosed with colorectal cancer (CRC) frequently confront a grave prognosis and exhibit poor responses to conventional treatment regimens. Immunotherapy, notably modalities centered on natural killer (NK) cells, represents a burgeoning frontier in the management of CRC. This study developed a validated prognostic model using NK-associated long non-coding RNAs (lncRNAs) to predict CRC outcomes.

Methods: Integrating single-cell RNA-seq (GSE146771_Smartseq2) and TCGA-COAD/READ bulk transcriptomic data, we identified NK-specific genes and correlated lncRNAs. A multi-step analytical approach-including univariate Cox regression for preliminary screening, LASSO regression to minimize overfitting, and multivariate Cox regression for final model optimization-yielded a robust 16-lncRNA prognostic signature with high predictive accuracy.

Results: This model demonstrated robust predictive performance across the training set, validation set, and 76 independent clinical samples. Mechanistic investigations revealed that AC010319.3 is highly expressed in NK cells, where it attenuates NK cell cytotoxicity by suppressing the expression of IFN-γ and granzyme B, thereby promoting the proliferation and invasion of CRC cells.

Discussion: This study systematically delineates the regulatory role of NK-associated lncRNAs within the CRC immune microenvironment, offering novel molecular targets and stratification strategies for CRC immunotherapy.

Keywords: NK cell-related lncRNAs; colorectal cancer; molecular subtyping; tumor immune microenvironment; tumor immune single-cell hub 2.

Figure 1

A detailed flowchart illustrates the construction, validation, and molecular subtyping of the NK cell-related lncRNA model in CRC. CRC: Colorectal Cancer; TISCH2: Tumor Immune Single-Cell Hub 2; TCGA: The Cancer Genome Atlas; TF: Transcription factor; LASSO: Least absolute shrinkage and selection operator; NK: Natural Killer cells.

Figure 2

Communication networks of NK cells in CRC. (A, B) UMAP plots show the distribution and abundance of different cell subpopulations in CRC. (C) A pie chart shows the percentage of NK cells. (D) Visualization of interaction probabilities between NK cells and other cells using CellChat. (E) A volcano plot shows differentially expressed genes within NK cells. Red indicates fold change > 1.5, FDR < 0.05; green indicates fold change < 1.5, FDR < 0.05.

Figure 3

Prognosis-associated NK cell-related lncRNAs. (A) A volcano plot shows 1,133 differentially expressed NK cell-related lncRNAs identified in CRC (red: logFC > 0.585, FDR-adjusted p < 0.05; green: logFC < 0.585, FDR-adjusted p < 0.05). (B) A heatmap visually displays the top 50 most significantly differentially expressed NK cell-related lncRNAs. (C) A forest plot shows the results of univariate Cox regression analysis, identifying 42 lncRNAs associated with CRC prognosis (green indicates hazard ratio < 1; red indicates hazard ratio > 1).

Figure 4

Prognosis-associated NK cell-related lncRNAs. (A) A heatmap shows the expression differences of identified lncRNAs between CRC and normal samples (*** indicates p < 0.001, ** indicates p < 0.01, * indicates p < 0.05). (B, C) Lasso regression analysis reveals overfitting phenomena in the model under different gene number settings and compares the severity of overfitting under these settings.

Figure 5

Construction and validation of the NK cell-related prognostic model. (A) A heatmap shows the expression of 16 lncRNAs in high- and low-risk groups in the training set. (B) Distribution and survival status of CRC patients in the training set. (C) Kaplan-Meier survival curve comparison between high- and low-risk groups in the training set. (D) ROC curve evaluation in the training set. (E) A heatmap shows the expression of 16 lncRNAs in high- and low-risk groups in the test set. (F) Distribution and survival status of CRC patients in the test set. (G) Kaplan-Meier survival curve comparison between high- and low-risk groups in the test set. (H) ROC curve evaluation in the test set. (I) A heatmap shows the expression of 16 lncRNAs in high- and low-risk groups in all patients. (J) Distribution and survival status of CRC patients in all patients. (K) Kaplan-Meier survival curve comparison between high- and low-risk groups in all patients. (L) ROC curve evaluation in all patients. (M) Heatmap displaying the expression of 16 lncRNAs in high- and low-risk groups within the external validation set(n=76). (N) Distribution and survival status of CRC patients in the external validation set(n=76). (O) Comparison of Kaplan-Meier survival curves between high- and low-risk groups in the external validation set set(n=76). (P) ROC curve evaluation in the external validation set(n=76).

Figure 6

Association between the prognostic model and clinical factors. (A) Univariate Cox regression shows factors affecting CRC prognosis. (B) Multivariate Cox regression shows independent factors affecting CRC prognosis. (C) ROC curve analysis evaluates the accuracy of various clinical variables and risk scores in predicting CRC prognosis. (D) Nomogram predicts CRC prognosis. (E) Calibration curve assesses the predictive ability of the nomogram. (F-O) Kaplan-Meier survival curves for various clinical subgroups based on risk scores.

Figure 7

Functional enrichment analysis in different risk groups. (A) A volcano plot shows differentially expressed genes between risk groups (red: logFC > 1, FDR-adjusted p < 0.05; green: logFC < 1, FDR-adjusted p < 0.05). (B) A heatmap shows the distribution of differentially expressed genes in different risk groups. (C) A Circos plot reveals changes in differentially expressed genes in the GO pathways. (D) Bubble plot reveals GO pathways enriched by significantly differentially expressed genes. (E) Bubble plot shows KEGG pathways enriched by significantly differentially expressed genes. (F) GSEA shows upregulated pathways in the high-risk group. (G) GSEA shows downregulated pathways in the low-risk group.

Figure 8

Immune and drug sensitivity analysis. (A-C) ESTIMATE scores, immune scores, and stromal scores in high- and low-risk groups. (D) Quantitative analysis of immune infiltration in the new CRC subtyping using various algorithms. (E) Box plots show immune function status in high- and low-risk groups. (F) Box plots show immune checkpoint status in high- and low-risk groups. (G) Box plots show drug sensitivity status in high- and low-risk groups. *p < 0.05, **p < 0.01, ***p < 0.001.

Figure 9

New CRC subtyping based on NK cell-related lncRNAs. (A) Sample distribution of different subtyping numbers. (B) CDF curves for different subtyping numbers. (C) Consensus CDF for different subtyping numbers. (D, E) Consensus matrices for three subtypes. (F) A Sankey diagram shows the relationship between different CRC subtypes and risk scores. (G) Survival curves for different CRC subtypes.

Figure 10

Immune and drug sensitivity analysis of the new CRC molecular subtyping. (A) Quantitative analysis of immune infiltration in the new CRC molecular subtyping using various algorithms. (B) ESTIMATE Score, ImmuneScore, and StromalScore for different CRC subtypes. (C) Immune checkpoint analysis for different CRC subtypes. (D) Drug sensitivity analysis for different CRC subtypes. *P < 0.05, **P < 0.01, ***P < 0.001.

Figure 11

AC010319.3 Promotes Colorectal Cancer Progression by Suppressing NK Cell-Related Functions. (A) Relative expression of lncRNA in NK cells within tumor tissues. (B) Comparison of lncRNA expression between NK cell lines and colorectal cancer cell lines. (C) Detection of intrinsic functions of NK cells through overexpression and knockdown of AC010319.3. (D) Transwell assay to validate invasive ability after overexpression and knockdown of AC010319.3. (E) Colony formation assay to validate proliferative ability after overexpression and knockdown of AC010319.3. (F) CCK-8 assay to validate proliferative ability after overexpression and knockdown of AC010319.3. *P < 0.05, **P < 0.01, ***P < 0.001.