We downloaded RNA-seq data including lncRNAs and mRNAs expression data and clinical information of colon cancer patients from The Cancer Genome Atlas (TCGA) dataset (https://cancergenome.nih.gov/). In this study, inclusion criteria were as follows: (1) patients with a diagnosis of COAD; (2) patients with complete clinical information on age, gender, vital status, pathologic T stage, pathologic N stage, pathologic M stage, and tumor stage; (3) patients with specific follow‐up time. Then, we treated the RNA matrix in the following way: (1) same RNA with more than one row were averaged; (2) RNAs with low value (expression <1) were excluded; (3) edgeR package of R software was used to normalized the RNA-seq data of ‘Level‐3’ read counts.

WGCNA is a systems biology method for finding gene modules with highly correlated expression levels and for relating them to clinical traits. Therefore, WGCNA is widely used to identify and screen biomarkers or therapeutic targets (13). Top 5000 mad lncRNAs and mRNAs were selected for further analysis with WGCNA package. First, all samples were clustered to construct a sampleTree and detect outliers based on cut height. Sample dendrogram and trait heatmap were used to develop networks to investigate the relationship between the corresponding sample gene expression data and clinical phenotypes. Then we selected β = 4 for lncRNAs and β = 6 for mRNAs as the soft‐thresholding parameters to construct the adjacency matrix. Next, we transformed the adjacency matrix into the topological overlap matrix. According to the topological overlap matrix (TOM)‐based dissimilarity measure, genes with highly absolute correlations were clustered into the same module to generate a cluster dendrogram (deep-split = 2, minimum cluster size = 30, cut height = 0.25). To visually represent the relationships between modules and clinical features of COAD, we calculated the Pearson correlation coefficient and plotted heatmap. Gene modules were significantly correlated with traits when P-value < 0.05. High correlation modules were selected as key modules for further analysis.

We constructed the ceRNA network based on lncRNAs and mRNAs in the key modules. The MiRcode (http://www.mircode.org/) database was used to predict lncRNA-miRNA interactions according to lncRNAs in the key module. StarBase (http://starbase.sysu.edu.cn/index.php) database was used to transform candidate miRNAs into human mature miRNA names. TargetScan (http://www.targetscan.org/vert_72/), miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/search.php), and miRDB (http://www.mirdb.org/) databases were used to predict miRNA-mRNA interactions, which must exist in three databases simultaneously. Finally, we obtained the overlaps between the candidate mRNAs and mRNAs in the key module. Cytoscape 3.7.1 software was used to construct and visualize the lncRNA‐miRNA ‐mRNA network.

We used Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analyses to deeply understand the biological function of target mRNAs. Then, the clusterProfiler package was used to conduct functional enrichment analyses of target mRNAs, and the enrichplot package and ggplot2 package were used to visualize the results of the enrichment analyses. Top enriched terms with P-value <0.05 were selected.

In order to construct a lncRNA-based risk model to predict OS of COAD patients, the following steps were performed with survival package, glmnet package, survminer package, and time receiver operating characteristic (ROC) package. Firstly, univariate COX regression analysis was used to select lncRNAs from the original ceRNA network. Next, lncRNAs with P-value <0.1 were included in the least absolute shrinkage and selection operator (LASSO) regression analysis. Lambda.min of LASSO analysis was used to construct lncRNA prognostic signature. The following formula was used to calculate the risk score of patients: Risk score = β_RNA1 * exp_RNA1 + β_RNA2 * exp_RNA2 + … + β_RNAn * exp_RNAn, where β _RNA indicated the prognostic coefficient of lncRNA, and exp_RNA stood for the expression level of lncRNA. At last, the risk score of every COAD patient was calculated based on the prognostic signature. All patients were divided into high-risk and low-risk sets according to the median risk score value. The Kaplan–Meier curve was used to evaluate OS, and ROC curve and areas under the curve (AUC) were used to evaluate the predictive accuracy of the prognostic risk model. Heatmap of clinical traits and the expression levels of lncRNAs was constructed by pheatmap package.

Nomogram was used to demonstrate a risk model based on lncRNAs and clinical characteristics for assessing OS in COAD patients. Univariate and multivariate COX regression analyses were performed for clinical characteristics using survival package. The P-value equal to 0.05 was selected to be the significant threshold. Clinical features with a significant threshold in multivariate analysis were identified as independent prognostic factors and included in nomogram. Finally, rms package was used to establish a nomogram. The C-index, calibration curve, and ROC analysis were used to describe the prediction value of the nomogram.

We integrated mRNA expression and clinical prognostic information (vital status and follow‐up time) of COAD patients. Survival package in R software was used to perform survival analysis. The mRNAs with P-value <0.05 were considered as crucial prognostic factors.

Colon cancer and adjacent normal tissues were collected from patients undergoing radical resection and with a postoperative pathology diagnosis of colon cancer at Zhongnan Hospital of Wuhan University from October 2019 to September 2020. Written informed consent was obtained from each patient before surgery and the patient protocols were approved by the hospital’s ethics committee. The samples were stored at −80°C immediately after resection until use.

Colon cancer cell lines (SW480, HCT 116, SW620, HT-29, and DLD-1) and normal human colon mucosal epithelial cell line (NCM460) were purchased from the Cell Bank of the Chinese Academy of Sciences. The cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM; BI, Israel) supplemented with 10% fetal bovine serum (FBS; BI, Israel) at 37°C with 5% CO₂.

The EASYspin RNA Mini Kit (Aidlab Biotechnologies Co., Ltd, China) was used for total RNA extraction according to the manufacturer’s instructions. Reverse transcription was performed using the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, USA) with 2 µg total RNA. PCR was conducted using the SYBR Green qPCR Mix (High ROX) (Monad, China) on a CFX96 Touch Real-Time PCR Detection System (Bio-Rad, USA). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) and U6 were used as internal controls and relative RNA expression was quantified using the 2^-ΔΔCt method (14). The primers (Tsingke, China; RiboBio, China) were listed in Supplementary Table 1 .

To over-expressed DNMBP-AS1 and FAM87A, lentivirus was designed and constructed by the Genechem Company (China). Stable transfected SW480 cells were obtained by selection with puromycin after lentivirus infection according to the manufacturer’s protocol. The miR-93-5p/17-5p mimics and mimics negative control (GenePharma, China) were transfected into SW480 cells using Lipofectamine 3000 (Invitrogen, USA) following the manufacturer’s instructions.

In the Cell Counting Kit-8 (CCK-8) assay, stable transfected SW480 cells were seeded into 96-well plates at a concentration of 2x10³ cells/well. 10 µL of CCK-8 reagent were added to each well and the cells were incubated for a further 2 h. To estimate cell viability, the absorbance of each well was measured at 450 nm every 24 h for 4 days.

In clone-formation assay, 1x10³ cells were plated in each well of a six-well plate. After culturing for 10 days, the cells were fixed with 4% paraformaldehyde for 30 minutes, stained with 0.1% crystal violet for 30 minutes, photographed, and counted the number of clones in each plate.

In the wound-healing assay, stable transfected SW480 cells were seeded into six-well plates. After culturing for 24 h, the cells were scratched with a linear wound with a 200 µL pipette tip, washed with phosphate buffered saline (PBS), imaged with an inverted microscope (Nikon, Japan) at 0 and 24 h, and measured the wound width using Image J software (version 1.53).

In the transwell assay, cellular invasion and migration assays were performed using transwell inserts (Corning, USA) coated with or without Matrigel, respectively. Stable transfected SW480 cells suspension (5×10⁴ cells) was added to the upper chamber, and DMEM containing 20% FBS was added to the lower chamber. After culturing at 37°C for 24 h, cells that had migrated or invaded to the lower chambers were fixed with 4% paraformaldehyde for 30 minutes and stained with 0.1% crystal violet for 30 minutes. Afterward, the stained cells were photographed under an inverted microscope (Leica, Germany) and counted using Image J software (version 1.53).

Bioinformatics analysis was applied to predict the potential binding sites of DNMBP-AS1/NHLRC3 and miR-93-5p/17-5p. Wild-type or mutant fragment of DNMBP-AS1 and NHLRC3 were constructed and inserted into the pmirGLO plasmid. Then the report plasmids and negative control or miR-93-5p/17-5p mimics were co-transfected into SW480 cells using Lipofectamine 3000 (Invitrogen, USA). After 24 h, the Dual-Luciferase Reporter Assay System (Promega, USA) was applied to detect Firefly and Renilla luciferase activity.

SW480 cells transfected with DNMBP-AS1 or control vector were harvested and suspended in PBS. Five-week-old male BALB/c nude mice (SPF Biotechnology Co., Ltd, China) were injected subcutaneously with stably transfected SW480 cells (1×10⁶), termed as vector or over-DNMBP-AS1 group (n = 3 per group). Tumor volume was monitored weekly and calculated using the formula: volume (mm³) = width²×length/2. After 5 weeks following the inoculation, the mice were euthanized and tumor samples were excised and weighted. The experiment was approved by the Animal Ethical and Welfare Committee of Zhongnan Hospital of Wuhan University and conducted following the ethics guidelines of the use of laboratory animals.

The data were analyzed using R (version 3.6.3) and GraphPad Prism 9 software. Differences between the two groups were assessed by unpaired Student’s t-test and those among multiple groups were assessed by one-way analysis of variance (ANOVA). The correlations were analyzed by Pearson’s test (r). The data are expressed as the mean ± standard deviation. P-value <0.05 indicated that the difference is statistically significant.

The workflow of this study is shown in Figure 1 . Expression matrices of lncRNAs and mRNAs were collected from 320 patients, who were pathologically diagnosed as COAD. Specific clinical and pathological characteristics of all COAD patients are shown in Table 1 , including age, gender, vital status, pathologic T stage, pathologic N stage, pathologic M stage, and tumor stage.

We selected the top 5000 mad lncRNAs and mRNAs for co‐expression analysis with WGCNA package. The cut heights were set as 55000 (lncRNA) and 2000000 (mRNA) to remove the outliers in sampleTree ( Figures 2A, B ). The sample dendrogram and trait heatmap split the selected samples into the different clusters and provided the distribution map of clinical trait data ( Figures 2C, D ). We found that when the power was equal to 4 in lncRNAs, the scale R² = 0.91 ( Figures 3A, B ); when the power was equal to 6 in mRNAs, the scale R² = 0.91 ( Figures 3C, D ). Therefore, we selected β= 4 as soft-thresholding in lncRNAs, and β= 6 as soft-thresholding in mRNAs co‐expression analysis. Dendrogram of the gene modules was based on a TOM-based dissimilarity measure. Sixteen modules were screened out in lncRNAs and twenty-three modules in mRNAs were shown in different colors ( Figure 4 ). The heatmap of module-trait relationships was used to select the key modules for further research. The red and purple modules were the key modules for lncRNAs and mRNAs, respectively. Both of them showed a significantly positive correlation with stage and M stage ( Figures 5A, B ). Subsequently, scatter plots of Gene Significance (GS) vs. module membership (MM) were shown in the key module. The plots revealed that MM in red module was significantly correlated with the tumor stage (cor=0.42,p=7.5e-18), and M stage (cor=0.35,p=4.1e-13). We found that MM in the purple module was significantly correlated with the tumor stage (cor=0.53,p=7.1e-12), and M stage (cor=0.54,p=2.4e-12) ( Figures 5C-F ).

The ceRNA network was built using 405 lncRNAs in red module and 145 mRNAs in purple module. We used miRcode database to match miRNAs with 405 lncRNAs. We got 687 lncRNA‐miRNA pairs containing 57 candidate lncRNAs and 87 candidate miRNAs. Then, based on three databases of TargetScan, miRDB, and miRTarBase, we obtained 5813 miRNA‐mRNA pairs including 87 miRNAs and 2556 candidate mRNAs. Finally, we picked out 34 target mRNAs overall after matching the 2556 candidate mRNAs with the 145 mRNAs in the purple module. Based on the above results, we constructed a lncRNA‐miRNA‐mRNA network with 50 target lncRNAs, 41 target miRNAs, and 34 target mRNAs and visualized the network using Cytoscape ( Supplementary Figure 1A ). The detailed correspondences of the original ceRNA network are shown in Supplementary Table 2 . The GO analysis results showed that target mRNAs were chiefly enriched in steroid hormone receptor signaling pathway, negative regulation of GTPase activity, negative regulation of B cell apoptotic process, transcription factor complex, and steroid hormone receptor binding. KEGG analysis results were significantly involved in hippo signaling pathway, longevity regulating pathway - multiple species, and microRNAs in cancer ( Supplementary Figures 1B, C ).

Univariate COX regression analyses for all 50 target lncRNAs were conducted to evaluate the associations between lncRNAs expression and OS. Eleven lncRNAs with a P-value < 0.1 in the univariate analyses were included in the LASSO regression ( Figure 6A ). Then, 9 lncRNAs (SNHG11, STK24-AS1, AL590483.1, MIR210HG, DNMBP-AS1, AL928654.1, AC019330.1, FAM87A, and HAR1A) were selected to construct a prognostic risk model for COAD after LASSO regression analysis ( Figure 6B ). The risk score of the model for OS was calculated in the following way: risk score = (-0.117 × expression level of SNHG11) + (-0.103 × expression level of STK24-AS1) + (-0.153 × expression level of AL590483.1) + (0.090 × expression level of MIR210HG) + (-0.024 × expression level of DNMBP-AS1) + (0.122 × expression level of AL928654.1) + (-0.048 × expression level of AC019330.1) + (-0.110 × expression level of FAM87A) + (0.064 × expression level of HAR1A). All 320 COAD patients were separated into high-risk and low-risk sets based on the median risk score. Vital status, survival time, and risk score levels were presented in Figure 6C . Survival analysis of OS showed that COAD patients with high-risk scores had a worse prognosis (P-value < 0.0001) ( Figure 6D ). The ROC curves of the risk model indicated that the values of AUC for 3-, 5-, and 10- year were 0.7, 0.696, 0.774, respectively ( Figure 6E ). From the above, the risk model based on the 9 lncRNAs revealed good performance in estimating OS. To explore the relationship between the 9 lncRNAs expression and clinical factors, we made a hierarchical cluster diagram ( Supplementary Figure 2 ).

We identified age, stage, and risk score as independent clinical prognostic factors by univariable and multivariable Cox regression analyses. All of the factors showed that high expression had a worse prognosis ( Table 2 ). Moreover, a nomogram was built based on the three independent prognostic factors ( Figure 7A ). The C-index for the nomogram was 0.823 and the calibration curves for predicting 3-, 5-, and 10-year OS were closely matched with actual survival outcomes ( Figures 7B–D ). The values of age AUC for 3-year survival, 5-year survival, and 10-year survival were 0.569, 0.619, and 0.517, respectively. The values of stage AUC for 3-year survival, 5-year survival, and 10-year survival were 0.779, 0.71, and 0.681, respectively. The values of whole nomogram AUC for 3-year survival, 5-year survival, and 10-year survival were 0.862, 0.808, and 0.852, respectively ( Figures 7E–G ). These validate results indicated that the nomogram combined with the lncRNA prognostic signature and clinical factors provided a greater performance for OS prediction than a single factor.

Based on the lncRNAs in the prognostic risk model, the ceRNA network further shrank with 9 lncRNAs, 35 miRNAs, and 29 mRNAs ( Figure 8A ). The specific correspondences of the ceRNA network are shown in Supplementary Table 3 . Survival analyses for 29 mRNAs was then applied and Kaplan ‐Meier curves show that HOXA10 and NHLRC3 with low expression had a worse prognosis in OS ( Figure 8B ). Therefore, we finally constructed a survival-related ceRNA network including 5 lncRNAs (SNHG11, FAM87A, DNMBP-AS1, MIR210HG, and AL928654.1), 8 miRNAs (hsa-mir-424-5p, hsa-mir-195-5p, hsa-mir-497-5p, hsa-mir-16-5p, hsa-mir-93-5p, hsa-mir-519d-3p, hsa-mir-96-5p, and hsa-mir-17-5p), and 2 mRNAs (HOXA10 and NHLRC3) ( Figure 8C ). FAM87A could form a ceRNA network with any of the miRNAs and mRNAs.

Based on qRT-PCR results, DNMBP-AS1 and FAM87A were indicated to be down-regulated in colon cancer cells and tissues. Compared with that in NCM460 cells, DNMBP-AS1 expression was significantly down-regulated in SW480 and DLD-1, and FAM87A expression was notably down-regulated in SW480, HT-29, and DLD-1 ( Figure 9A ). Therefore, the SW480 cell line was selected for further analysis. The results of 30 pairs of human colon cancer tissues and their adjacent normal tissues showed that DNMBP-AS1 and FAM87A expression were elevated in normal tissues ( Figure 9B ).

To determine the function of DNMBP-AS1 and FAM87A in colon cancer, a series of assays was performed using stable transfected SW480 cells. Firstly, qRT-PCR was used to detect the expression levels of DNMBP-AS1 and FAM87A in SW480 cells transfected with control vector (vector), DNMBP-AS1 (over-DNMBP-AS1), or FAM87A (over-FAM87A). The results demonstrated that both DNMBP-AS1 and FAM87A expression significantly increased ( Figures 10A, B ). Then, CCK-8 assay and clone-formation assay were used to reveal that cell proliferative potential. Up-regulation of DNMBP-AS1 and FAM87A both significantly decreased the cell viability of SW480 cells ( Figure 10C ). Stable ectopic expression of DNMBP-AS1 and FAM87A reduced the number of clone spot than those stably transfected with the control vector ( Figure 10D ). Subsequently, wound-healing assay and transwell assay were carried out to evaluate the effects of DNMBP-AS1 and FAM87A on cell metastasis. After up-regulating lncRNA expression, both DNMBP-AS1 and FAM87A, the wound healing rates were remarkably decreased ( Figure 10E ). In transwell assay without Matrigel, migrated SW480 cells in over-lncRNA groups were less than those in the vector group ( Figure 10F ). In transwell assay with Matrigel, invaded cells were evidently decreased in over-lncRNA groups compared with that in vector group ( Figure 10G ). Finally, the animal experiment was performed to reveal the effect of stable transfected DNMBP-AS1 on cellular proliferation in vivo. Five weeks after subcutaneous injection, both xenograft tumors size and weight in over-DNMBP-AS1 group were obviously reduced ( Figure 10H ). On the contrary, these results suggested that over-expression of DNMBP-AS1 evidently suppressed the growth of colon cancer in vitro and in vivo.

Through the previous bioinformatics analysis, we predicted potential miRNAs and mRNAs that bind to DNMBP-AS1 and FAM87A. In the final ceRNA network, three miRNAs (hsa-mir-17-5p, hsa-mir-519d-3p, and hsa-mir-93-5p) and one mRNA (NHLRC3) were related to DNMBP-AS1; eight miRNAs (hsa-mir-424-5p, hsa-mir-195-5p, hsa-mir-497-5p, hsa-mir-16-5p, hsa-mir-93-5p, hsa-mir-519d-3p, hsa-mir-96-5p, and hsa-mir-17-5p), and two mRNAs (HOXA10 and NHLRC3) were related to FAM87A ( Figure 8C ). Afterward, those candidate miRNAs and mRNAs were evaluated in validation assays in SW480 cells. We found that over-expression of DNMBP-AS1 decreased the expression levels of miR-93-5p/17-5p and increased the expression level of NHLRC3, while the expression of miR-519d-3p was not affected ( Figure 11A ). As for FAM87A, the qRT-PCR results showed that the expression levels of most miRNAs were elevated and all mRNAs were significantly decreased in over-FAM87A groups ( Figure 11B ), which were not consistent with our prediction and the literature reported. Therefore, further analyses were performed to explore the effects of DNMBP-AS1 acting as ceRNA on miR-93-5p/17-5p and NHLRC3. In addition, to further excavate more miRNAs that were sponged by DNMBP-AS1 and more downstream mRNAs, we also verified the expression of miRNAs and mRNA in the secondary ceRNA network. Compared with that in the vector group, the expression of miR-106a-5p, miR-192-5p, miR-214-3p, and miR-217 were reduced, and the expression of AKAP11, RB1, NFIC, and DACH1 were increased in the over-lncRNA group ( Supplementary Figures 3A, B ). All of these results are consistent with the theoretical results of DNMBP-AS1 acting as ceRNA, so those candidate miRNAs and mRNAs are worth further vitro and vivo study.

The diagrammatic sketch of the binding sites for miR-93-5p/17-5p in DNMBP-AS1 is shown in Figure 11C . To find evidence of direct binding between DNMBP-AS1 and miR-93-5p/17-5p based on their complementary sequences, miRNA negative control (miR-NC), miR-93-5p mimics (miR-93-5p), or miR-17-5p mimics (miR-17-5p) were co-transfected with the wild-type or mutant DNMBP-AS1 dual-luciferase reporter, respectively. Obvious reductions in luciferase activity were observed when wild-type DNMBP-AS1 and miR-93-5p/17-5p mimics were co-transfected, but no significant change was observed in the mutant DNMBP-AS1 group ( Figure 11D ). Next, we identified miR-93-5p/17-5p were up-regulated both in colon cells (SW480) and colon tumorous tissue samples ( Figures 11E, F ). In colon cancer tissues, a negative correlation was observed between the expression of DNMBP-AS1 and miR-93-5p, and also of DNMBP-AS1 and miR-93-5p ( Figure 11G ). Taken together, this series of evidence strongly indicated that DNMBP-AS1 could act as miR-93-5p/17-5p sponges.

qRT-PCR analysis was performed to further confirm that over-expression of DNMBP-AS1 could decrease the miR-93-5p/17-5p expression levels ( Figure 12A ) while miR-93-5p/17-5p expression levels could markedly rise when transfected with miR-93-5p/17-5p mimics ( Figure 12A ). To confirm whether DNMBP-AS1 exerts its biological function by sponging miR-93-5p/17-5p, the CCK-8 assay and transwell assay were conducted for rescue experiments. We observed that when DNMBP-AS1 and miR-93-5p/17-5p mimics were co-transfected, the reduced cell viability mediated by DNMBP-AS1 over-expression could be reversed ( Figure 12B ). The number of migrated cells and invaded cells were elevated in DNMBP-AS1 and miR-93-5p/17-5p co-transfected groups compared to those in DNMBP-AS1 and miR-NC co-transfected group ( Figures 12C, D ). This evidence suggested that DNMBP-AS1 affected colon cancer cell proliferation and metastasis by sponging miR-93-5p/17-5p.

The bioinformatics prediction results demonstrated that there were binding sites between miR-93-5p/17-5p and 3′UTR of NHLRC3 ( Figure 13A ). The results of the dual luciferase reporter assay revealed that luciferase activity was reduced after co-transfection of wild-type NHLRC3 and miR-93-5p/17-5p mimics, while no reduction in luciferase activity was observed in the mutant group ( Figure 13B ). In order to further explore the interactions among DNMBP-AS1, miR-93-5p/17-5p, and NHLRC3, we detected NHLRC3 expression in co-transfected miR-NC, miR-93-5p, or miR-17-5p with vector or over-lncRNA. We found that in the cells transfected with DNMBP-AS1 and miRNA negative control (Over-lncRNA + miR-NC), the expression of NHLRC3 was up-regulated in SW480 cells ( Figure 13C ). On the contrary, in the cells transfected with DNMBP-AS1 and miR-93-5p/17-5p mimic (Over-lncRNA + miR-93-5p/17-5p), the expression of NHLRC3 was down-regulated compared with Over-lncRNA + miR-NC group ( Figure 13C ). The mRNA expression of NHLRC3 was examined by qRT-PCR, NHLRC3 were down-expressed both in colon cells and colon cancer tissues ( Figures 13D, E ). Pearson correlation results manifested that NHLRC3 expression were negatively correlated with miR-93-5p/17-5p and that was positively correlated with DNMBP-AS1 in 30 pairs of human colon samples ( Figure 13F ). All of the results illustrated that DNMBP-AS1 competitively bound to miR-93-5p/17-5p and acted as ceRNA to boost NHLRC3 expression.