Construction and validation of a prognostic model for colorectal cancer based on migrasome-related long non-coding RNAs

Data acquisition and MRLs screening

High clinical and pathological information as well as high-throughput RNA sequencing data for 39 normal colorectal tissues and 396 COAD tissues were acquired from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). The expression data for lncRNA and mRNA were then normalized, while clinical information, including age, sex, survival time, survival status, AJCC stage and stage (T, N, M), were collected. Previous research identified nine migrasome-related genes (MRGs) (Qin et al., 2022; Zhao et al., 2019; Zheng et al., 2023). Pearson correlation coefficients were used to identify lncRNAs associated with MRG expression, with p < 0.001 and correlation coefficient r > 0.4 selected as screening criteria. The “limma” package was eventually used to analyze the normalized data before visualizing the results with the “ggalluvial” package.

Prognostic model development and validation for MRLs

After randomly assigning clinical samples to training and testing cohorts, univariate Cox regression analysis was performed on the training group, using p < 0.001 as the significance threshold to identify seven potential candidate MRLs. In this case, overfitting was minimized by applying the Least Absolute Shrinkage and Selection Operator (LASSO), thereby selecting only the most relevant predictive factors. Multivariable Cox regression analysis identified four MRLs based on which a risk feature model was constructed for predicting the survival rate of COAD patients. The risk score was then calculated as follows:

$$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\sum _{\mathrm{N}=1}^{\mathrm{i}}\left(\mathrm{C}\mathrm{o}\mathrm{e}\mathrm{f}\left(\mathrm{l}\mathrm{n}\mathrm{c}\mathrm{R}\mathrm{N}\mathrm{A}\mathrm{i}\right)\times \mathrm{E}\mathrm{x}\mathrm{p}\mathrm{i}\left(\mathrm{l}\mathrm{n}\mathrm{c}\mathrm{R}\mathrm{N}\mathrm{A}\mathrm{i}\right)\right).$$

In the risk model, N represents the number of prognosis-associated lncRNAs, Coef refers to the regression coefficients calculated using multivariable Cox regression analysis, and Expi denotes each lncRNA’s expression level. Correlation studies were conducted using data from the four MRLs, with the relationships subsequently visualized on a heatmap. Using the prognostic model, risk scores were calculated for both the training and testing groups, with the median score subsequently selected for classifying samples into low-risk and high-risk groups. Risk curves were generated prior to visualization with the “pheatmap” tool in R. Kaplan-Meier (K-M) survival analysis as well as receiver operating characteristic (ROC) curve analysis were performed using the “survminer” and “timeROC” packages.

The “survival” and “RMS” packages were used to develop predictive plots which combined risk levels and clinical pathological features for estimating the 1-, 3-, and 5-year survival probabilities for COAD patients. The relationship between expected survival rates and observed outcomes were eventually compared using calibration curves.

Principal component analysis and functional enrichment analysis

Principal component analysis (PCA), a technique that is widely applied for dimensionality reduction, was performed for extracting primary feature components, with the “scatterplot3d” tool in R subsequently used to visualize differences between the low-risk and high-risk groups.

These two categories were then compared using the “limma” package to identify differentially expressed genes (DEGs), with a false discovery rate of <0.05 and a |log2(fold change)|of >1 selected as selection criteria (Ritchie et al., 2015). For functional enrichment analysis, enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) were assessed using the “clusterProfiler” package. GO analysis provides insights into cellular components (CC), molecular functions (MF), and biological processes (BP), while KEGG analysis integrates genomic, chemical, and functional data to analyze biological pathways. Additionally, the “limma” package was used for gene set enrichment analysis (GSEA), with GO and KEGG enrichment analysis for differentially expressed lncRNAs between the two risk groups subsequently conducted with the “enrichplot”, “org.Hs.eg.db”, “limma”, and “clusterProfiler” packages. Functional enrichment was considered significant at p < 0.05.

Tumor mutational burden and tumor immune analysis

Tumor mutation burden (TMB) is defined as the number of mutations per million bases and serves as a key biomarker in tumor progression and immunotherapy response. Mutation data for COAD were obtained from TCGA and analyzed using the “maftools” package to compare TMB differences between high- and low-risk groups and visualize the top 15 mutated genes. Survival analysis was conducted to evaluate its prognostic significance.

To explore immune characteristics, we used the CIBERSORT algorithm to estimate immune cell infiltration and assessed immune checkpoint expression differences between risk groups. Additionally, we performed Tumor Immune Dysfunction and Exclusion (TIDE) analysis (http://tide.dfci.harvard.edu/) to predict tumor immune escape. TIDE scores were compared between the high- and low-risk groups to evaluate potential differences in response to immune checkpoint blockade therapy. Statistical comparisons were performed using “limma” and “ggpubr”, with significance levels set at ^*p < 0.05, ^**p < 0.01, and ^**p < 0.001.

Prediction of potential drugs for the treatment of COAD

To predict potential antitumor drugs for treating colorectal cancer patients, we calculated the half-maximal inhibitory concentration (IC50) of antitumor drugs obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) website. The “limma,” “ggpubr,” and “ggplot2” packages were employed to identify antitumor drugs that may be effective in treating colorectal cancer tumors.

Cell culture and CRC patient samples

Human COAD cell lines HCT116, DLD-1, SW620, and RKO and the normal human intestinal epithelial cell line NCM460 were used in this study. The selection of these cell lines was based on their genetic background, invasiveness, and relevance to COAD progression. Specifically, HCT116 and RKO exhibit high invasiveness, SW620 originates from lymph node metastases, and DLD-1 represents a well-characterized COAD model. The normal cell line NCM460 was used as a control.

Cells were obtained from Wuhan Pricella Biotechnology Co., Ltd. (Wuhan, China). All cell lines were authenticated before use, and routine mycoplasma contamination tests were conducted to ensure experimental validity. HCT116 and DLD-1 cells were cultured in RPMI-1640 medium, while SW620 and RKO cells were maintained in high-glucose DMEM medium. All culture media were supplemented with 10% fetal bovine serum (FBS) and 100 U/mL penicillin-streptomycin. Cells were incubated at 37 °C in a humidified atmosphere with 5% CO₂.

For clinical sample validation, 20 paired colorectal adenocarcinoma (COAD) tumor tissues and adjacent normal tissues were collected from patients who underwent surgical resection at The Fourth Affiliated Hospital of Harbin Medical University. The inclusion and exclusion criteria for patient selection were as follows:

Inclusion criteria: (1) Patients with histologically confirmed primary COAD; (2) No history of preoperative chemotherapy, radiotherapy, or targeted therapy; (3) Availability of complete clinical and follow-up data; (4) Patients who provided written informed consent. Exclusion criteria: (1) Patients diagnosed with metastatic COAD at the time of surgery; (2) Patients with prior malignancies or concurrent cancers; (3) Presence of severe systemic diseases affecting study outcomes; (4) Poor RNA quality or insufficient tissue samples for further analysis. This study was approved by the Institutional Review Board (IRB) of The Fourth Affiliated Hospital of Harbin Medical University (Approval Reference No.: YXLLSC-2018-01). All procedures followed the International Ethical Guidelines for Biomedical Research Involving Human Subjects.

RNA extraction and cDNA synthesis

Total RNA was extracted from cells and tissue samples using TRIzol reagent (Takara, Shiga, Japan) following the manufacturer’s protocol. RNA concentration and purity were assessed with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), ensuring an A260/A280 ratio between 1.8 and 2.0. Subsequently, 1 µg of total RNA was converted to cDNA using the ReverTra Ace™ qPCR RT Kit (Toyobo Bio, Osaka, Japan). The reverse transcription reaction was performed at 37 °C for 15 min and terminated by heating at 85 °C for 5 s. The resulting cDNA was stored at −20 °C for future experiments.

Quantitative reverse transcription polymerase chain reaction

Quantitative reverse transcription polymerase chain reaction (two-step RT-qPCR) was conducted using the SYBR Green PCR kit (2× S6 Universal SYBR qPCR Mix, EnzyArtisan, Shanghai, China) on a QuantStudio™ 6 Flex real-time PCR system (Applied Biosystems, Waltham, MA, USA). Each reaction had a total volume of 10 µL, comprising 5 µL of SYBR Green mix, 0.5 µL of each primer, 1 µL of diluted cDNA template, and 3 µL of nuclease-free water. The PCR conditions included pre-denaturation at 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s and 60 °C for 30 s. Reaction specificity was confirmed by melt curve analysis.

Primer design and data analysis

All primers were designed using Primer-BLAST, and their specificity was verified using NCBI BLAST. The primers were synthesized by Synbio Technologies (Beijing, China). The sequences of primers used for PCR amplification in this study are listed below:

AC004846.1, forward: 5′-ATCCGTGTGGGTGCTCCTTCC-3′,

reverse: 5′- CGCTGGTACTTCGTTGCCTCTG-3′;

AC010789.2, forward: 5′-CCACCATGAAGAGAAAAGCAGG-3′,

reverse: 5′-TTTCTCTGACCAGTGCTTGTTCTG-3′;

ZEB1-AS1, forward: 5′-TGTCGGAGTTGGAAAGGGAC-3′,

reverse: 5′-CTACTAAGGAGGCTGCTGGC-3′;

LCMT1-AS1, forward: 5′-GTTTCTGCTTGCTGGTCGC-3′,

reverse: 5′-TGGATTCTCTGCCCTTTGCC-3′;

Vimentin, forward: 5′-AGGCAAAGCAGGAGTCCACTGA-3′,

reverse: 5′- ATCTGGCGTTCCAGGGACTCAT-3′;

E-cadherin, forward: 5′-GCCTCCTGAAAAGAGAGTGGAAG-3′,

reverse: 5′-TGGCAGTGTCTCTCCAAATCCG-3′.

Each sample was tested in triplicate, with β-actin serving as the internal reference control. Relative expression levels were determined using the 2^−ΔΔCt method. The expression differences of MRLs were assessed using t-tests. Graphs were generated using GraphPad Prism (version 9.5.1).

RNA transfection by lipofectamine

The LCMT1-AS1-set shRHA was synthesized by Applied Biological Materials company (Richmond, BC, Canada). Oligonucleotide transfection was conducted using MaxFect Lipo 3000 Transfection Reagent according to the manufacturer’s protocol. Forty-eight hours post-transfection, cells were harvested for further analysis. Knockdown efficiency was confirmed by RT-qPCR.

Proliferation assay

The CCK-8 assay was conducted to evaluate cell proliferation. A total of 2,000 transfected HCT116 and RKO cells were seeded into 96-well plates. CCK-8 reagent was added at specified time points and incubated at 37 °C for 1 h. Absorbance was measured at 450 nm.

A colony formation assay was conducted to evaluate cell growth. A total of 400 transfected HCT116 and RKO cells were seeded into six-well plates and maintained until colonies formed. Colonies were then fixed with paraformaldehyde for 30 min and stained with crystal violet. The number of colonies was counted for further analysis.

Transwell migration assay

A tumor cell migration assay was conducted following the manufacturer’s instructions. Briefly, cells were collected, resuspended in serum-free medium, and seeded onto Transwell inserts at a density of 100,000 cells per well. The inserts were then placed into lower chambers filled with 600 μL of medium containing 10% FBS and incubated at 37 °C for 24 h. Cells on the upper surface of the Transwell insert were removed with a cotton swab, while cells that had migrated to the underside were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet. Images were captured from five random areas, and cells were counted to determine the average number of migrated cells.

Statistical analysis

Cox regression, Pearson correlation coefficients, and K-M survival analysis were conducted using R (version 4.4.1; https://www.r-project.org/). Comparisons among groups were performed using GraphPad Prism software (version 9.5.1), with statistical tests conducted using analysis of variance (ANOVA). Correlation analysis was carried out using Spearman correlation coefficients. A p-value that is not significant (ns) is indicated as non-significant. A p-value less than 0.05 (^*) indicates statistical significance, while p-values less than 0.01 (^**) and 0.001 (^***) indicate higher degrees of significance.

Identification of MRLs associated with prognosis in COAD

Based on the TCGA-COAD data, 16,876 lncRNAs were identified. Using the expression profiles of the nine MRGs obtained from previous studies, the expression profiles of MRGs were then extracted from the TCGA-COAD dataset. Through correlation analysis between colorectal cancer lncRNAs and MRGs, a total of 2,317 MRLs were identified (|or| > 4, p < 0.001). The relationship between MRLs and MRGs in COAD patients was illustrated in a Sankey diagram (Fig. 1). Overall, 377 COAD patients were randomly divided into a training (N = 189) or a testing group (N = 188), and as shown in Table 2, the two groups were not significantly different in terms of their clinical characteristics. Analysis of the training set using univariate Cox regression models identified seven MRLs (p < 0.001) significantly associated with COAD patients’ prognosis (Fig. 2A). This was followed by LASSO regression analysis to screen for overfitting MRLs before investigating any potential association between four MRLs (AC004846.1, AC010789.2, ZEB1-AS1, LCMT1-AS1) and the clinical outcomes of COAD patients based on multivariate Cox regression analysis (Figs. 2B, 2C). Correlations between these four MRLs and MRGs were visualized in a heatmap (Fig. 2D).

Prognostic model development and assessment for MRLs

Using the following risk score formula, a novel prognostic model was developed for the MRLs:

$$\begin{array}{rl}\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=& {\displaystyle \mathrm{A}\mathrm{C}004846.1\times 1.903+\mathrm{A}\mathrm{C}010789.2\times 0.915+\mathrm{Z}\mathrm{E}\mathrm{B}1\text{-}\mathrm{A}\mathrm{S}1\times 0.870}\\ & +\mathrm{L}\mathrm{C}\mathrm{M}\mathrm{T}1\text{-}\mathrm{A}\mathrm{S}1\times \mathrm{1.176.}\end{array}$$

Risk scores were calculated for each patient in both the training and testing sets. Based on these scores, patients were categorized into high-risk and low-risk groups, using the median risk score as the cutoff value (Fig. 3A). Patients in the high-risk group exhibited a higher mortality rate compared to those in the low-risk group (Fig. 3B). Additionally, the heatmap demonstrated significant differences in MRL expression levels between the two patient groups (Fig. 3C), with the expression levels of the four MRLs being significantly higher in the high-risk group. K-M survival curve analysis revealed that the overall survival (OS) time for COAD patients in the high-risk group was significantly shorter than that of the low-risk group (Fig. 4A). To further validate the predictive ability of the MRL risk score, ROC analysis was conducted based on the risk model and clinical characteristics. The area under the ROC curve (AUC) effectively demonstrated the accuracy of this risk score. In the training set, we calculated AUC values for three different time points: 1 year (0.744), 3 years (0.760), and 5 years (0.861) (Fig. 4B). Similarly, the AUC values in the testing set were 0.665, 0.646, and 0.628, respectively. Furthermore, the AUC value of the MRL risk score (0.762) was higher than those for age (0.627), tumor stage (0.684), and gender (0.487) (Fig. 4C). These results indicate that the MRL risk score has good predictive ability for the prognostic outcomes of COAD patients.

Evaluation of the MRLs prognostic model in COAD patients

We conducted univariate and multivariate Cox regression analyses to assess the correlation between the MRLs risk score and the prognosis as well as clinical characteristics of COAD patients. The results indicated that the MRLs risk score (p < 0.001) could serve as an independent prognostic factor, similar to age (p < 0.001) and tumor stage (p = 0.041) (Figs. 5A, 5B). We constructed a nomogram model that integrates clinical characteristics and risk scores to predict the survival rates of COAD patients at 1 year, 3 years, and 5 years (Fig. 6). The calibration plot demonstrated that the nomogram model exhibited excellent predictive performance (Fig. 7A). Moreover, the MRLs risk score showed a higher concordance index (C-index) compared to other risk factors, further validating its superior predictive capability (Fig. 7B).

Clinical characteristic subgroups of the MRLs prognostic model and PCA

The significance of the MRLs prognostic model for clinical subgroups was further assessed by assigning patients to three subgroups based on sex (male and female), age (≤65 years and >65 years) and tumor stage (I–II and III–IV), with their survival outcomes subsequently analyzed. The results indicated that prognosis was not significantly affected by factors such as age, gender, or tumor stage (Figs. 8A–8F).

Four gene sets involving all genes, migrasome-related genes, model lncRNAs, and migrasome-related lncRNAs, were analyzed using PCA to determine their potential in distinguishing between risk groups of patients, with the results presented in Figs. 8G–8J. The analysis demonstrated that low-risk and high-risk patients could not be differentiated using three gene sets (all genes, migrasome-related genes, and migrasome-related lncRNAs), whereas the PCA plot for model lncRNAs revealed a clear separation between these groups. These findings indicate that model lncRNAs possess superior discriminative ability, effectively distinguishing high-risk from low-risk populations.

Enrichment analysis of the MRLs prognostic model

To investigate the functional relevance and signaling pathways associated with the prognostic model, we conducted GO functional analysis, KEGG enrichment analysis, and GSEA on the differentially expressed genes between the high-risk and low-risk groups. The GO enrichment analysis indicated significant biological activities of the model-related genes in the extracellular matrix, connective tissue development, and molecular interactions, which may be related to tissue repair, structural stability, and cellular development (Fig. 9A). The KEGG enrichment analysis revealed several pathways significantly associated with cell migration, adhesion, and signaling, particularly those closely related to the cytoskeleton, PI3K-Akt signaling pathway, focal adhesion, and extracellular matrix interactions (Fig. 9B). The GSEA demonstrated that functions such as extracellular matrix remodeling, plasma membrane protein complexes, and receptor complexes were significantly enriched in the high-risk group, indicating more active intercellular signaling, cell adhesion, and invasive biological processes. At the pathway level, cell adhesion molecules, cytokine-receptor interactions, and immune response-related pathways were significantly enriched. In contrast, the functional enrichment of the low-risk group focused on ribosomes and the respiratory chain, with metabolic pathways such as the tricarboxylic acid (TCA) cycle, oxidative phosphorylation, and ribosome biogenesis significantly enriched, suggesting that cellular metabolic functions and protein synthesis in the low-risk group were in a healthy state (Figs. 10A–10D).

TMB and immune correlation analysis

TMB quantifies the number of genetic mutations within a tumor. Elevated TMB levels are often linked to enhanced immune responses and can serve as predictors for immunotherapy efficacy. Utilizing somatic mutation data from TCGA, TMB scores were calculated, highlighting the 15 genes with the highest mutation frequencies. The waterfall plot results indicate that the mutation frequency of these 15 genes is higher in the high-risk group compared to the low-risk group. Approximately 97.73% of COAD patients in the high-risk group exhibited TMB, with the most common mutated genes being APC (70%), TP53 (62%), TTN (48%), KRAS (43%), and PIK3CA (26%) (Fig. 11A). In the low-risk group, about 93.18% of COAD patients were affected, with APC (73%), TP53 (52%), TTN (43%), KRAS (42%), and PIK3CA (33%) as the top five genes with the highest mutation frequency (Fig. 11B). Furthermore, additional analyses showed that TMB scores in the high-risk group were significantly greater than those in the low-risk group, indicating an increase in genetic mutation levels (Fig. 11C). To assess the predictive value of risk scores for the prognosis of COAD patients, we categorized patients into high and low mutation groups based on the median TMB of each sample. Kaplan-Meier analysis revealed that patients in the low mutation group had better survival outcomes compared to the high mutation group (Fig. 11D). When combining TMB with risk scores for survival analysis, significant differences in survival outcomes were observed among the four groups (p < 0.001), with patients exhibiting high TMB and high risk showing the poorest survival outcomes (Fig. 11E). These results suggest that risk scores possess good predictive value for the prognosis of COAD patients.

Analysis of tumor microenvironment differences revealed significant disparities between high-risk and low-risk groups. Compared to the low-risk group, samples from the high-risk group showed higher levels of stromal and immune cells, suggesting that tumor growth, invasion, and metastasis are more likely to occur in high-risk patients (Fig. 12A). TIDE analysis results showed that the TIDE score was significantly higher in the high-risk group than in the low-risk group, suggesting that high-risk patients may have stronger immune evasion capabilities and potential resistance to immune checkpoint inhibitor therapy (Fig. 12B). Immune infiltration analysis indicated a higher proportion of naïve B cells in the high-risk group, implying that these patients face a stronger immune challenge, while the low-risk group had a higher proportion of memory B cells, indicating a more effective immune response and potentially better prognosis (Fig. 12C). Analysis of immune-related functional differences revealed marked differentiation between high-risk and low-risk groups across several immune functions. The high-risk group scored higher on certain immunosuppressive mechanisms, such as co-inhibition of antigen-presenting cells (APCs) and pro-inflammatory cells, suggesting that this group may be more prone to immunosuppression and inflammation promotion. This state could enable tumor cells to evade immune surveillance and accelerate disease progression. In contrast, the low-risk group demonstrated relatively high scores for B cells and CD8+ T cells, indicating that these immune effector cells play a more active role. CD8+ T cells are critical for anti-tumor immunity, suggesting that the low-risk group may have stronger anti-tumor immune responses (Fig. 12D).

Correlations between the MRLs prognostic model and sensitivity to targeted therapy drugs

To further investigate the relationship between the MRLs prognostic model and drug sensitivity, we compared the IC50 of various drugs between the high-risk and low-risk groups using box plots. Among the eight drugs assessed, significant differences were found between the two groups (p < 0.05). The results indicated that patients in the high-risk group exhibited greater sensitivity to Cediranib, Dasatinib, Entospletinib, and Ruxolitinib (Figs. 13A–13D), while those in the low-risk group were more sensitive to Dabrafenib, Erlotinib, Oxaliplatin, and Sorafenib (Figs. 13E–13H). These findings suggest that our MRLs-associated prognostic model may serve as a valuable predictor of drug sensitivity.

Validation of genes in the MRLs prognostic model

We first evaluated the expression levels of four MRLs in colorectal cancer cell lines (HCT116, DLD-1, SW620, RKO) compared to the normal colorectal cell line NCM460. Notably, the expression levels of the four MRLs were significantly higher in the HCT116 and RKO cell lines (Fig. 14A). Subsequently, we analyzed the expression of these MRLs in 20 pairs of colorectal cancer tissues and adjacent normal tissues. As illustrated in Figs 14B–14E, the expression levels of the four MRLs in tumor tissues were elevated compared to those in adjacent normal tissues. These experimental results further confirm the reliability of the risk model.

LCMT1-AS1 enhances the proliferation and migration of COAD cells

Among the four lncRNAs in the model, LCMT1-AS1 has limited prior research. To investigate its effect on COAD, we knocked down LCMT1-AS1 in HCT116 and RKO cells, which exhibit high expression of this gene. The efficiency of the transient knockdown was confirmed by RT-qPCR (Fig. 15A). Following LCMT1-AS1 knockdown, Vimentin expression decreased while E-cadherin expression increased in HCT116 and RKO cells (Fig. 15B), suggesting that LCMT1-AS1 may play a role in promoting epithelial-mesenchymal transition (EMT). Knocking down LCMT1-AS1 may inhibit EMT, thereby suppressing the migratory and invasive abilities of colorectal cancer cells.

To assess the proliferation of HCT116 and RKO cells, we conducted CCK-8 and colony formation assays. As shown in results, LCMT1-AS1 knockdown significantly inhibited the growth of these cells (Figs. 15C, 15D). Additionally, Transwell assays were performed to evaluate the effect of LCMT1-AS1 on COAD cell migration. The results indicated that LCMT1-AS1 knockdown reduced the migratory ability of HCT116 and RKO cells (Fig. 15E).