GSE61635 and GSE50772 dataset of peripheral blood mononuclear cells (PBMC) from SLE patients, were downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database. The GSE61635 dataset, comprising 99 SLE samples and 30 control samples, was used as the training cohort. The GSE50772 dataset, including 61 SLE samples and 20 control samples, was used as the validation cohort. Two datasets of CRC patients, including GSE39582 and GSE17536, were retrieved from GEO database as well (20, 21). The dataset of colon adenocarcinoma (COAD) and rectum adenocarcinoma (READ) was obtained from The Cancer Genome Atlas (TCGA) database via UCSC Xena website (https://xena.ucsc.edu/) (22). The GSE39582 dataset, as a training set, consisted of 566 tumor samples and 19 control samples. GSE17536, as well as the COAD and READ datasets, were used as external validation cohorts. GSE17536 comprises 177 tumor samples. The COAD dataset includes 448 tumor and 39 control samples, while the READ dataset consists of 158 tumor and 9 control samples. The expression data from GEO database was generated using the GPL570 platform. Single-cell RNA sequencing data from 7 CRC patients in the E-MTAB-8107 dataset were downloaded from the TISCH2 database (http://tisch.compbio.cn) (23). Single-cell data from 3 SLE patients in the GSE142016 dataset were obtained from the GEO database (24).

DEGs in GSE39582 and GSE61635 datasets were screened out and visualized in the form of volcano plots using the “limma” and “ggplot2” packages in R software. |Log2 (fold change)| >1 and adjusted P value ≤ 0.05, with multiple testing correction performed using the Benjamini–Hochberg method were considered to indicate statistical significance. Subsequently, the common DEGs in both CRC and SLE datasets was shown in a Venn diagram. The expression patterns of common DEGs were visualized in the form of heatmaps with “pheatmap” package. Moreover, in order to analyze the biological functions and pathways involved in common genes, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using The Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://davidbioinformatics.nih.gov/), with terms/pathways with adjusted P value (Benjamini–Hochberg) < 0.05 were considered significantly enriched Additionally, a protein-protein interaction (PPI) network was constructed using STRING (https://string-db.org/) to evaluate potential functional associations among these shared genes.

Relationships between common DEGs and small molecule therapeutics were explored via the cMAP database (https://clue.io) (25). The 31 upregulated DEGs were incorporated into the cMAP platform, resulting in identification of the top 10 compounds exhibiting maximal enrichment scores. Chemical structure of small molecular compounds was illustrated via ChemSpider database (http://www.chemspider.com) and visualized with Chemdraw software (version 22.0.0.22).

Signature genes for CRC and SLE were identified using three machine learning methods: univariate Cox regression, LASSO, and SVM-RFE algorithm. First, prognostic genes were screened from common DEGs using univariate Cox regression implemented via the “survival” package. Subsequently, the LASSO regression model was applied using the “glmnet” package, where the hyperparameter lambda (λ) minimized regression coefficients toward zero to exclude weakly correlated features, with optimal parameters determined through 5-fold cross-validation (26). Concurrently, SVM-RFE was employed using the “mlbench” and “caret” packages to iteratively eliminate the least important feature vectors (27). Finally, characteristic genes identified by LASSO and SVM-RFE for CRC and SLE were intersected, yielding shared signature genes for downstream analysis.

NetworkAnalyst (http://www.networkanalyst.ca) is a comprehensive web-based tool designed to identify proteins of interest, construct networks, and perform network analysis and visualization (28). The characteristic genes for CRC and SLE were uploaded to the NetworkAnalyst website, mapped to the integrative interaction database to construct correlation networks. The interactions of significant genes with transcription factors (TFs) and microRNAs (miRNAs) were analyzed subsequently. Transcription factor targets were derived from the JASPAR TF binding site profile database. Comprehensive experimentally validated miRNA-gene interaction data were collected from TarBase database. Finally, coregulatory networks among genes, TFs, and miRNA were constructed by Cytoscape software (version 3.10.2).

The 22 types of immune cell composition in patients and controls were estimated with the cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) algorithm (29). The proportion and abundance of immune cells in each sample were shown as the stacked histogram with the “ggplot2” package. Subsequently, the differences between the proportions of immune cells in the both groups were shown in a box plot. Next, the association of the infiltrated immune cells was evaluated with the “corrplot”package and presented as a heatmap with the”pheatmap”package. Finally, the correlation between the four characteristic genes and the infiltrating immune cells was analyzed using the Spearman correlation coefficient and visualized as a heatmap.

The nomogram diagnostic model was constructed for four characteristic genes using the “rms”package. The risk score for each participant of SLE and controls was calculated based on the nomogram model. Receiver operating characteristic (ROC) curves via “pROC” package, violin plots via “ggbiplot” package, calibration curves via “calibration curve” package, and decision curve analysis (DCA) via “rmda” package evaluated the ability of the diagnostic model to discriminate between SLE and controls (30). GSE50772 was used as an external validation cohort to evaluate the robustness of the model.

Based on the normalized expression levels (Expi) and corresponding regression coefficients (Coei) of prognostic genes, the risk score of each patient was calculated according to the following formula:

First, the diagnostic performance of the four-gene-based risk score for CRC was evaluated using ROC analysis. Subsequently, its prognostic capability for predicting overall survival (OS) of CRC was assessed. CRC patients were divided into high-risk and low-risk groups based on the median cut-off value of the risk score. Subsequently, Kaplan-Meier analyses were performed using the “survival”package to compare the survival of patients in the both groups. Distribution in the proportion, survival time, and survival status of high-risk and low-risk groups in CRC patients were explored with the “ggrisk” package.

Furthermore, univariate and multivariate Cox regression analyses were also conducted to evaluate the impact of risk score and other clinical characteristics (including age, chemotherapy, clinical stage, T/N/M classification) on the prognosis of colorectal cancer patients. ROC, calibration, and DCA curves were evaluated the ability of the prognostic model to predict OS of CRC. The COAD and READ datasets served as external validation cohorts to assess the predictive capability of the nomogram.

Somatic mutations in the 15 most frequently mutated genes were analyzed using the “maftools” package to compare mutational profiles between high-risk and low-risk cohorts. Furthermore, copy number variation (CNV), specifically gene amplifications and deletions, were assessed in COAD and READ cohorts via “tidyverse” package.

Stemness gene sets were retrieved from StemChecker website (http://stemchecker.sysbiolab.eu/), which was based on the most comprehensive and updated collection of published stemness signatures defined by expression profiles, RNAi screens, literature curation, computationally derived, Nanog, and transcription factor target gene sets. Then, quantitative analysis of stemness enrichment score of each participant was elucidated with the single-sample gene set enrichment analysis (ssGSEA) method. Adjusted P value < 0.05 was used as the significance threshold.

The differential expression of immune checkpoint and human leukocyte antigen (HLA) were assessed across high-risk and low-risk populations. The immune score, stromal score, ESTMATE score, tumor purity of each sample was calculated with the “ESTMATE” package (31). In addition, tumor immune dysfunction and exclusion (TIDE) score of each sample was evaluated to predict the possibility of immune escape and immune response occurring in the tumor via the website (http://tide.dfci.harvard.edu) (32). Immunophenotypic score (IPS) data were obtained from The Cancer Immunome Atlas (TICA) database (https://tcia.at/home) (33). The responsiveness to anti-PD-1 or anti-CTLA-4 therapy was evaluated by comparing IPS between high-risk and low-risk groups via the “ggpubr” package.

The DEGs between the high-risk and low-risk groups in the GSE39582 dataset were defined and presented as the volcano plot. Potential biological functions were examined through GO enrichment analysis via Metascape database (https://metascape.org) (34). Enriched pathways were further integrated using the gene set enrichment analysis (GSEA) was performed with the “clusterProfiler” package, applying thresholds of |normalized enrichment score (NES)| > 1, NOM P value < 0.05, and q value < 0.25. Additionally, gene set variation analysis (GSVA) was conducted to identify pathways enriched in the two risk cohorts via the “GSVA” package.

Single-cell RNA sequencing data from 7 CRC patients in the E-MTAB-8107 dataset were processed using the “Seurat” package. Cells meeting quality control criteria (featuring more than 400 and fewer than 6000 expressed genes, with a mitochondrial gene percentage below 20%) were retained for subsequent analysis. Data normalization was performed using the “NormalizeData” function. Highly variable genes across cells were identified using the “FindVariableFeatures” function. Subsequently, data scaling was conducted via the “ScaleData” function, followed by principal component analysis (PCA). Cell clustering analysis was performed using the “FindNeighbors” and “FindClusters” functions, and the resulting clusters were visualized by applying uniform manifold approximation and projection (UMAP). Cell-cell communication was analyzed using the “CellChat” package.

We aimed to investigate the expression patterns of the four core genes across eight annotated cell types: epithelial, plasma, T, myeloid, mast, natural killer (NK), fibroblast, and endothelial cells. Furthermore, these cell types were categorized into three functional groups—epithelial, stromal (fibroblasts and endothelial cells), and immune cells (plasma, T, myeloid, mast, and NK cells)—to systematically compare and analyze the expression levels of the four hub genes within these broader functional compartments. Based on the positive or negative expression status of the deoxyribonuclease 1 like 3 (DNASE1L3) gene, all cells were stratified into a positive group and a negative group. The proportion of each cell subpopulation within these two groups was calculated to elucidate the expression pattern of DNASE1L3. Differential expression analysis was further performed between these two groups, and the results were visualized using a volcano plot. KEGG analysis was conducted on the DEGs to identify significantly enriched signaling pathways. Pseudotime analysis was analyzed using the “monocle” package to investigate developmental trajectories of cells (35). A similar analytical workflow was applied to single-cell RNA sequencing data from 3 SLE patients in the GSE142016 dataset.

Spatial transcriptome analysis of CRC was performed via the Sparkle website (https://grswsci.top/) (36). Briefly, inverse convolutional analysis was employed to accurately assess cellular composition on different spatial localizations. Strict quality control measures were implemented on the collected single-cell transcriptome data based on the number of expressed genes, the count of unique molecular identifiers, and the percentage of mitochondrial RNA in each cell. Subsequently, a signature score matrix was constructed by calculating the average expression of the top 25 specifically expressed genes in each locus for each cell type. Finally, an enrichment scoring matrix was generated using the “get enrichment matrix” and “enrichment analysis” functions in the “Cottrazm” package. The enrichment scores for each cell type were visualized using the “SpatialFeaturePlot” function in the “Seurat” package. The malignant group was defined if the score of malignant cells in the microregion was 1, it is defined as the non-malignant group otherwise. Spearman correlation between gene expression and cell content was visualized by the “linkET” package.

A tissue microarray containing tumor samples and adjacent tissues from 76 CRC patients was obtained from Aifang company (AF-CocSur2201, China). The pathological sections were first subjected to deparaffinization and rehydration. Antigen retrieval was then performed using EDTA buffer. Endogenous peroxidase activity was subsequently blocked with hydrogen peroxide, followed by blocking with 3% bovine serum albumin. Tissue sections were incubated overnight at 4 degrees with the primary antibody against DNASE1L3 (GTX114363, GeneTex, USA). A secondary antibody (GB23303, Servicebio, China) was applied and incubated at room temperature for 1 hour. Diaminobenzidine was used as the chromogen to visualize immunoreactivity, and hematoxylin was employed for nuclear counterstaining. Stained specimens were imaged and analyzed using a microscope (Model KF-FL-020, KFBIO, China). The percentage of positively stained cells was quantified using Visiopharm software.

SLE patients, CRC patients, and healthy volunteers were recruited from the First Affiliated Hospital of Zhengzhou University. Peripheral blood was collected and centrifuged to obtain plasma. Peripheral blood mononuclear cells (PBMC) were extracted by density gradient centrifugation method using Lymphoprep reagent (1858, Serumwerk Bernburg AG, Germany).

The human monocyte cell line THP-1, immortalized colon epithelial cell line NCM460, and the embryonic kidney cell line HEK-293T were purchased from HyCyte Company (China). THP-1, NCM460 and PBMC were cultured in RPMI 1640 medium (11875093, Gibco, USA) containing 10% fetal bovine serum (FBS, A5256701, Gibco, USA) and 1% penicillin/streptomycin (P/S, 15140122, Gibco, USA). HEK-293T cells were cultured in DMEM medium (11965092, Gibco, USA) containing FBS and P/S. All cells were cultured in 37 degrees incubator containing 5% CO2.

The DNASE1L3-targeting plasmids (sh-DNASE1L3), control plasmid (sh-NC), and helper plasmids (pSPAX2 and pMD2G) were synthesized by Tsingke company (China). These plasmids were transfected into HEK-293T cells using polyethylenimine (PEI; 40816ES01, Yeasen, China). After 48–72 hours of incubation, lentivirus-containing supernatants were harvested. THP-1 cells were then incubated with the lentiviral supernatant supplemented with 2 μg/mL polybrene (C0351, Beyotime, China) for 48 hours. Stable transductants were selected using 3 μg/mL puromycin (P8230, Solarbio, China). Successful knockdown of DNASE1L3 was confirmed by reverse transcription-quantitative polymerase chain reaction (RT-qPCR) and Western blotting analysis.

Total RNA of cells was extracted by Trizol reagent (15596026CN, Thermofisher Scientific, USA). Then, RNA was reverse transcribed to cDNA according to the instructions of a Reverse Tran scription Kit (R433, Vazyme, China). A total of 15 pairs of cDNA microarrays of CRC tumor and matched paraneoplastic tissues were obtained from Shanghai Outdo Biotech Company (MecDNA-HColA030CS03, China). Next, qPCR was performed according to the instructions of a SYBR Green PCR Kit (Q712, Vazyme, China). Finally, the mRNA expression was obtained by the 2^-ΔΔCt method and standardized β-actin as the internal control. Primer sequences were listed in Supplementary Table S1.

Plasma DNASE1L3 levels in patients with SLE (n = 19) or CRC (n = 19), and healthy controls (n = 19) were quantified using an ELISA kit (CSB-EL007052HU, CUSABIO, China) according to the manufacturer’s protocol. Optical density at 450 nm was immediately measured using a Synergy H1 microplate reader (BioTek Instruments, USA).

Cells were lysed on ice for 30 minutes in RIPA buffer (CW2333S, CWBIO, China) supplemented with 1% protease (CW2200S, CWBIO, China) and phosphatase inhibitor cocktail (CW2383S, CWBIO, China). Protein concentration was determined using a BCA assay kit (23227, Thermo Fisher Scientific, USA). Proteins were separated by SDS-PAGE and subsequently transferred onto PVDF membranes (IPVH00010, Millipore, USA). After blocking with 5% non-fat milk for 1 hour, membranes were incubated overnight at 4°C with primary antibodies against DNASE1L3 or β-tubulin (80762-1-RR, Proteintech, China). Next, membranes were probed with HRP-conjugated secondary antibodies (GAR1007, MultiSciences, China) for 1 hour at room temperature. Protein signals were detected using chemiluminescent reagent (GK10008, GLPBIO, USA). Relative expression levels of DNASE1L3 normalized to β-tubulin were quantified using ImageJ software.

The NCM460 cells were labeled with 1 μM CMFDA (HY-126561, MedChemExpress, USA) for 30 minutes and then washed with PBS before being resuspended in complete culture medium. Apoptosis was induced in NCM460 cells by stimulation with 10 μM staurosporine (HY-15141, MedChemExpress, USA) for 12 hours, followed by washing with PBS. Transfected THP-1 cells were differentiated into macrophages by incubation with 100 ng/mL phorbol myristate acetate (PMA; HY-18739, MedChemExpress, USA) for 48 hours, as previously reported (37, 38). THP-1-derived macrophages were stained with 2 μM DiI (HY-126561, MedChemExpress, USA) for 20 minutes and washed with PBS. Labeled NCM460 cells and macrophages were co-cultured at a ratio of 5:1 for 4 hours. Non-phagocytosed and unbound apoptotic cells were then removed by washing with PBS. Total RNA was extracted for RT-qPCR analysis. Fluorescence signals were acquired on a FACS Canto II flow cytometer (BD Biosciences, USA). Data analysis was performed using FlowJo software (10.8.1, USA). Efferocytosis was quantified as the proportion of macrophages that had phagocytosed NCM460 cells, represented by the percentage of FITC-positive cells.

All expression data from public databases was processed in R software (version 4.2.0). Unpaired t-tests were employed to compare two groups of variables exhibiting normal distributions. For non-normally distributed metrics, comparisons were made using the Mann-Whitney U test. Comparisons among multiple groups were performed using the Kruskal-Wallis test. Survival differences were assessed utilizing the log-rank test and visualized via Kaplan-Meier survival curves. The correlation was analyzed using the Spearman method.

The strategy of data processing was shown in Figure 1. Differential analysis showed 1967 DEGs between CRC and control tissues in the GSE39582 cohort, including 1087 up-regulated genes and 880 down-regulated genes. In the GSE61635 cohort, 721 DEGs, including 474 up-regulated genes and 247 down-regulated genes, were identified in PBMCs from SLE and controls. Volcano plots were applied to demonstrate the expression pattern of DEGs in GSE39582 and GSE61635 datasets (Figures 2A, B). To further explore the shared key genes of CRC and SLE, Venn diagrams were used to showed the 58 common DEGs of CRC and SLE (Figure 2C). The protein-to-protein interaction (PPI) network illustrated the interactions among proteins associated with the common DEGs (Figure 2D). Heatmaps were carried out to identify the expression pattern of common DEGs in GSE39582 and GSE61635 cohorts (Figures 2E, F). Biological process of GO analysis revealed that the common DEGs were enriched in regeneration, regulation of extrinsic apoptotic signaling pathway in absence of ligand, and cell aging (Figure 2G). Cellular component of GO analysis showed that the DEGs were enriched in tertiary granule lumen, RNA polymerase II transcription regulator complex, and anchored component of membrane (Figure 2H). Molecular function of GO analysis displayed that the DEGs were enriched in ATPase-coupled intramembrane lipid transporter activity, peroxidase activity, and oxidoreductase activity, acting on peroxide as acceptor (Figure 2I). KEGG analysis demonstrated that the DEGs were enriched in glutathione metabolism, chemokine signaling pathway, and thyroid hormone synthesis (Figure 2J).

The 31 up-regulated DEGs were uploaded to the cMap database to identify potential therapeutic drugs. After significance screening, the top 10 compounds—including treprostinil, IKK-2-inhibitor, PRL-3-inhibitor-I, idebenone, alrestatin, GW-843682X, BRD-K15107389, AT-9283, pyrimethamine, and MAZ-51, were identified as potential therapeutic agents for CRC and SLE (Figure 2K). These compounds primarily function as phosphatase inhibitors, VEGFR inhibitors, calcium channel modulators, etc. (Figure 2L). The chemical structures of these 10 compounds are shown in Supplementary Figure S1.

The univariate Cox, SVM-RFE, and LASSO were applied to explore the hub genes in CRC and SLE. First, the univariate Cox regression was carried out to identify 15 potential candidate genes from 58 common DEGs. The LASSO regression was used to identify 9 potential genes from CRC and 10 genes from SLE. The SVM-RFE algorithm was constructed to identified 14 potential genes from CRC and 13 potential genes from SLE (Figures 3A–F). The overlapping genes among all subsets were defined as hub genes, i.e., genes associated with SLE diagnosis as well as the diagnosis and prognosis of CRC. A Venn diagram identified the four hub genes—DNASE1L3, protein tyrosine phosphatase non-receptor type 14 (PTPN14), selenium-binding protein 1 (SELENBP1), and esophageal cancer related gene 4 (ECRG4)—which were subsequently included in further analyses (Figure 3G). Network analysis via the NetworkAnalyst database revealed that the four signature genes are regulated by transcription factors such as SRF, YY1, and IRF2, as well as miRNAs including hsa-mir-335-5p, hsa-mir-27a-3p, and hsa-mir-124-3p (Figures 3H, I).

The proportion of immune cells in each sample of CRC cohort was calculated by the “CIBERSORT”algorithm (Figure 4A). There were significant differences in the proportions of the 11 types of immune cells between the control and CRC groups. Compared to the control group, the CRC group had a higher proportion of M0 macrophages, M1 macrophages, activated mast cells, neutrophils, activated memory CD4 T cells, while a lower proportion of memory B cells, resting dendritic cells, resting mast cells, monocytes, plasma cells, and resting memory CD4 T cells (Figure 4B). In addition, correlation analysis of immune cells showed a significant positive correlation between M1 macrophages and follicular helper T cells (r = 0.449, P < 0.001), and a significant negative correlation between activated NK cells and resting NK cells (r = -0.566, P < 0.001) (Figure 4C). Moreover, the correlation revealed that the four genes were closely associated with immune cell infiltration. Plasma cells were positively correlated with SELENBP1, and DNASE1L3. Neutrophils were negatively correlated with SELENBP1 (Figure 4D).

Similarly, the proportion of immune cells in each sample of SLE cohort was displayed in Figure 4E. There were significant differences in the proportions of the 7 types of immune cells between the control and SLE groups. Compared to the control group, the SLE group had a higher proportion of M0 macrophages, neutrophils, yet a lower proportion of naive B cells, resting mast cells, resting memory CD4 T cells, and CD8 T cells (Figure 4F). There was a significant positive correlation between neutrophils and naive CD4 T cells (r = 0.481, P < 0.001), and a significant negative correlation between neutrophils and CD8 T cells (r = -0.695, P < 0.001) (Figure 4G). PTPN14 expression showed the highest positive correlation with neutrophil infiltration and highest negative correlation with eosinophil infiltration. SELENBP1 expression showed the highest negative correlation with neutrophils. DNASE1L3 expression demonstrated the strongest positive correlation with resting memory CD4T cells and the strongest negative correlation with plasma cells. ECRG4 expression shows the highest positive correlation with CD8 T cells, while it shows the highest negative correlation with neutrophils (Figure 4H).

In order to diagnose and predict SLE, four signature genes were used to construct a diagnostic nomogram model of SLE with the “rms” package (Figure 5A). ROC curves were applied to calculate the area under curve (AUC)values of the nomogram model to determine their sensitivity and specificity for diagnosing SLE. The results showed that the nomogram model had an AUC value of 0.993 for the diagnosis of SLE (Figure 5B). Based on the nomogram model, risk scores were calculated for each sample, and the violin plot demonstrated significantly higher risk scores for SLE compared to the control group (Figure 5C). Additionally, DCA calculated the clinical “net benefit” of the nomogram model at different threshold probabilities compared to a strategy of treating all patients or no patients (Figure 5D). Furthermore, the calibration plot revealed the predicted probabilities of the nomogram model were almost identical to the actual probabilities (Figure 5E). Similarly, the model showed good performance in discriminating SLE patients from controls in the GSE50772 dataset (Figures 5F–I).

The risk score for each sample was calculated based on the expression levels of the four genes. The risk score was defined as follows: Risk score = (0.2120) × PTPN14 + (-0.0354) × SELENBP1 + (-0.0788) × DNASE1L3 + (0.1093) × ECRG4. The histogram displayed the coefficients of the four genes in the risk score model (Figure 6A). The AUC values for the risk score in diagnosing CRC were 0.894 in GSE39582 (Figure 6B), 0.691 in COAD (Supplementary Figure S2A), and 0.668 in READ (Supplementary Figure S2B). The forest plot indicated that the four gene significantly influenced the survival prognosis of CRC (Figure 6C). CRC patients were stratified into high and low risk groups based on the median value of risk score. The risk of death was increased in the high-risk group compared to the low-risk group (Figure 6D). Notably, survival analysis uncovered that OS was significantly lower in the high-risk group compared with the low-risk group in the GSE39582 and GSE17536 cohorts. Patients in the high-risk group demonstrated a trend toward reduced OS compared to those in the in the low-risk group in the COAD and READ cohorts (Figure 6E).

Risk score, age, chemotherapy, clinical stage, T classification, N classification, and M classification were included in univariate Cox regression, and the results revealed that risk score, age, clinical stage, T classification, N classification, and M classification significantly affected prognosis of patients (Figure 7A). Multivariate Cox regression analysis verified that risk score, age, and M classification were significant prognostic factors for patients (Figure 7B). The above results showed that risk score, age, and M classification were independent prognostic factors for CRC patients. The associations between the four hub genes and key clinicopathological features of CRC were further systematically analyzed. The results indicated that the expression levels of PTPN14 and ECRG4 were higher in patients with T3–4 classification (Supplementary Figure S3A). Additionally, PTPN14 expression was elevated in patients with N1–2 classification (Supplementary Figure S3B). We also observed that DNASE1L3 expression was lower in M0 patients, whereas PTPN14 expression was higher in M1 patients (Supplementary Figure S3C). Furthermore, PTPN14 levels were increased in patients with stage III–IV (Supplementary Figure S3D), and ECRG4 expression was higher in patients who received chemotherapy (Supplementary Figure S3E).

Based on the significant factors, a nomogram model was constructed to predict the prognosis of CRC (Figure 7C). The nomogram demonstrated AUC values of 0.778, 0.752, and 0.728 for predicting 1-, 3-, and 5-year survival of CRC, respectively (Figure 7D). Patients in the high-risk group demonstrated a significantly poorer prognosis compared to those in the low-risk group (Figure 7E). DCA showed that the clinical “net benefit” of the predictive model at different threshold probabilities was higher than that of the strategy of treating all patients or not treating patients (Figure 7F). Moreover, the calibration curve revealed minimal deviation between the actual and predicted probabilities (Figure 7G).

Similar results were obtained in the validation cohort. The nomogram model demonstrated time-dependent predictive accuracy for COAD, with AUC values of 0.752, 0.728, and 0.679 at 1, 3, and 5 years, respectively. Furthermore, survival, DCA, and calibration curves collectively validated its robust predictive performance (Supplementary Figures S4A–D). Due to the lack of sufficient data on the 5-year survival of patients in the READ cohort, we analyzed the accuracy of the nomogram model in predicting 1-, 2-, and 3-year survival for READ. The AUC values were 0.842, 0.802, and 0.806, respectively (Supplementary Figure S4E). Furthermore, survival, DCA, and calibration curves further validated the model’s robust predictive capability (Supplementary Figures S4F–H). These results suggested that the nomogram model had good diagnostic and prognostic prediction efficacy for CRC.

Somatic mutation analysis revealed that APC, TP53, and TTN were the three most frequently mutated genes in both high-risk and low-risk groups across the COAD and READ cohorts (Figure 8A). In the COAD cohort, copy number amplifications were predominantly observed in PCDH17, TWIST1, and LY6E genes, whereas in the READ cohort, amplifications were more commonly detected in PCDH17, INHBA, and TWIST1. Copy number deletions were consistently prevalent in ZDHHC2, EPHX2, and IFI6 genes in both cohorts (Figure 8B).

A total of 204 DEGs were identified between the low- and high-risk group in the GSE39582 cohort. The distribution of these DEGs was shown in the volcano plot (Figure 8C). To gain a better understanding of potential molecular processes and functions, GO, KEGG, GSVA, and GSEA enrichment analysis was then conducted based on the DEGs. GO analysis illustrated that pathogenic genes were closely enriched in extracellular matrix, glycosaminoglycan binding, and collagen fibril organization (Figure 8D). KEGG analysis revealed that pathogenic genes were mostly enriched in protein digestion and absorption, ECM-receptor interaction, and cytoskeleton in muscle cells (Figure 8E). GSVA analysis indicated that DEGs were mostly clustered in immunoglobulin complex circulating, IgA immunoglobulin complex, and pH elevation (Figure 8F). GSEA analysis revealed that the antigen processing and presentation, complement and coagulation cascades, leukocyte transendothelial migration, ECM-receptor interaction, and phagosome were enriched in the high-risk group, while drug metabolism cytochrome P450, fatty acid degradation, folate biosynthesis, glycolysis/gluconeogenesis, and oxidative phosphorylation were enriched in the low-risk group (Figure 8G).

Cancer stem cells are a very small number of tumor cells in tumor tissue with biological properties including unlimited proliferation, self-renewal and multidirectional differentiation. There was a significant difference in stemness enrichment scores between the high- and low-risk groups by 26 stemness gene sets, suggesting that the pattern of tumor genesis, progression, infiltration and recurrence after treatment may be different between the two groups (Figure 9A).

Functional enrichment analysis of signature genes showed a close association with immune cell infiltration in CRC. Therefore, we explored the correlation between risk score and immune cell characteristics. Analyses of HLA typing showed significant differences in the 9 HLAs between the high- and low-risk groups. Eight HLA were upregulated in the high-risk group, whereas only one HLA was upregulated in the low-risk group (Figure 9B). The expression levels of CD40, CD86, CD200, CD274, CD276, HAVCR2, LAIR1, NRP1, PDCD1LG2, TNFRSF4, TNFRSF8, TNFRSF9, TNFSF4, TNFSF9, and TNFSF14 were significantly higher in the high-risk group, while BTNL2, CD40LG, CD160, CD244, CTLA4, HHLA2, LGALS9, PDCD1, TMIGD2, and TNFRSF25 were lower (Figure 9C). Differences in immune checkpoint gene expression between high- and low-risk groups may result in differential susceptibility to immune checkpoint inhibitors. Furthermore, analysis using the estimation of stromal and immune cells in malignant tumor tissues using expression data (ESTIMATE) algorithm revealed that the high-risk group exhibited significantly higher immune, stromal, and ESTIMATE scores, along with lower tumor purity (Figure 9D). In addition, TIDE score in high-risk patients was significantly higher than that of low-risk patients, suggesting a worse efficacy of immunotherapy (Figure 9E). Finally, we compared the IPS between the two risk groups in COAD and READ cohorts. The low-risk group exhibited a significantly elevated IPS, indicating a more promising potential for immune response. In contrast, the high-risk group showed a markedly reduced IPS, suggesting a potential resistance to immune checkpoint inhibitors (Figures 9F, G). In conclusion, there were significant differences in the microenvironment between high and low risk patients. These heterogeneities may contribute to more precise tumor immunotherapy.

Analysis of single-cell RNA sequencing data from 7 CRC patients in the E-MTAB-8107 dataset identified 23 clusters, which were further classified into 8 distinct cell types, including epithelial cells, plasma cells, T cells, myeloid cells, mast cells, NK cells, fibroblast cells, and endothelial cells (Figures 10A, B). The top three marker genes for each cell type are shown in Figure 10C. Cell communication analysis revealed that NK cells were the primary signal receivers, while fibroblast cells were the major signal senders (Figures 10D, E, Supplementary Figure S5). UMAP visualization illustrated the distribution and expression levels of the four feature genes across the eight cell types (Figure 10F). These findings indicate that DNASE1L3 was predominantly expressed in myeloid cells, PTPN14 in endothelial cells, SELENBP1 in epithelial and fibroblast cells, and ECRG4 mainly in fibroblastic cells (Figure 10G). The eight cell types were categorized into three broad classes: immune cells, epithelial cells, and stromal cells. We observed that DNASE1L3 was primarily distributed in immune cells, PTPN14 in stromal cells, SELENBP1 in both epithelial and stromal cells, and ECRG4 mainly in stromal cells (Figure 10H). Pseudotime analysis was subsequently performed (Figure 10I, J), revealing distinct developmental patterns for the four genes (Figure 10K). The expression of PTPN14, SELENBP1, and ECRG4 gradually increased along the pseudotime trajectory, whereas DNASE1L3 expression initially increased and then decreased (Figure 10L). Furthermore, all cells were stratified into two subgroups based on DNASE1L3 expression status (high vs. low) to evaluate the contribution of different cell types to gene expression. Myeloid cells were identified as the primary contributors to DNASE1L3-high expression (Figure 10M). A total of 360 DEGs were identified between DNASE1L3- high and - low cells, as shown in Figure 10N. KEGG enrichment analysis indicated that these DEGs were significantly enriched in pathways such as phagosome, Epstein–Barr virus infection, and staphylococcus aureus infection (Figure 10O).

Analysis of single-cell RNA sequencing data from 3 SLE patients in the GSE142016 dataset identified 11 clusters, which were annotated into 5 immune cell types, including CD4 T cells, CD8 T cells, NK cells, B cells, and myeloid cells (Figures 11A, B). The top three marker genes for each cell type are presented in Figure 11C. Cell communication analysis demonstrated that CD8 T cells were both the primary senders and receivers of signals (Figures 11D, E, Supplementary Figure S6). UMAP plots displayed the distribution and expression levels of the four feature genes across the five cell types (Figure 11F). The results showed that DNASE1L3 was mainly expressed in B cells, PTPN14 in NK cells, SELENBP1 in CD4 T cells and myeloid cells, and ECRG4 primarily in CD4 T cells (Figure 11G). Pseudotime analysis was conducted (Figures 11H, I); however, due to low expression levels of ECRG4 and SELENBP1, only the trajectories of PTPN14 and DNASE1L3 were explored (Figure 11J). PTPN14 expression gradually decreased along pseudotime, whereas DNASE1L3 expression progressively increased (Figure 11K). Further stratification of all cells into DNASE1L3-high and -low subgroups revealed that B cells were the predominant cell type in the high-expression group, while myeloid cells were most abundant in the low-expression group (Figure 11L). A total of 111 DEGs were identified between the two groups (Figure 11M). KEGG pathway analysis indicated significant enrichment in salmonella infection, shigellosis, diabetic cardiomyopathy, and chemical carcinogenesis-reactive oxygen species (Figure 11N).

Integrated spatial transcriptome sequencing analysis was performed to explore the spatial location and gene expression of four prognostic genes in CRC tissue. As shown in Figure 12A, Hematoxylin and eosin (H&E) staining of pathological sections was performed. Ten distinct cell types were identified, and their spatial localization was presented in Figure 12B. The demarcation between malignant and non-malignant regions was displayed in Figure 12C. The spatial distribution of the four signature genes was shown in Figure 12D. Further investigation into their expression revealed that DNASE1L3, PTPN14, and SELENBP1 were primarily localized to the malignant regions, whereas the expression of ECRG4 showed no significant difference between malignant and non-malignant areas (Figure 12E). Correlation analysis demonstrated that all four signature genes were closely associated with a variety of cell types (Figure 12F).

The mRNA expression levels of DNASE1L3, SELENBP1 and ECRG4 were significantly lower in the tumor group compared to the control group in the GSE39582 dataset, whereas the mRNA expression levels of PTPN14 were significantly higher in the tumor group (Figure 13A). Consistent with the aforementioned results, similar expression trends for DNASE1L3, ECRG4, and SELENBP1 were observed in the COAD and READ datasets, except that PTPN14 showed no significant difference in the COAD dataset but was markedly downregulated in the READ dataset (Supplementary Figure S7). Patients with low levels of DNASE1L3 exhibited a significantly poorer prognosis than those with high levels, whereas the opposite trend was observed for PTPN14. In contrast, the expression of SELENBP1 and ECRG4 showed no significant association with patient survival (Figure 13B). The qPCR results of the 15 pairs of tumor and pericarcinomatous tissues confirmed the above results (Figure 13C).

These findings suggest that DNASE1L3 and PTPN14 may represent potential key targets for improving prognosis in CRC. Given that immune dysregulation is a critical pathological feature shared by CRC and SLE, our study, along with existing evidence, indicates that DNASE1L3—primarily expressed in macrophages—plays a central role in mediating the clearance of apoptotic cells, maintaining self-antigen homeostasis, and preventing immune dysregulation. Additionally, IHC staining of tumor tissues and pericarcinomatous tissues from 76 CRC patients further confirmed the lower expression level of DNASE1L3 protein in tumor tissues (Figure 13D).

The expression levels of DNASE1L3, PTPN14, SELENBP1, and ECRG4 in the GSE61635 dataset were consistent with those in the GSE39582 dataset (Figure 14A). Similarly, RT-qPCR analysis of PBMCs from 25 SLE patients and 23 healthy controls confirmed the mRNA expression levels of these four genes (Figure 14B). Subsequently, a diagnostic nomogram was developed to predict the probability of SLE development in both the control and SLE group (Figure 14C). The ROC curve demonstrated that the nomogram achieved an AUC of 0.920 for diagnosing SLE (Figure 14D). Nomogram assessment revealed a significant discrimination in SLE risk between the control and SLE groups (Figure 14E). Both the calibration and DCA curves indicated that nomogram-based decisions could provide clinical utility for SLE prediction (Figures 14F, G).

Significantly reduced DNASE1L3 levels were observed in both the SLE and CRC groups compared to the control group via ELISA method (Figure 15A). To simulate the pathological microenvironment of low DNASE1L3 expression observed in both diseases, we knocked down DNASE1L3 expression in the human monocytic cell line THP-1 to further investigate the regulatory mechanisms of this gene in disease pathogenesis. RT-qPCR and Western blotting were used to verify the knockdown efficiency of DNASE1L3 in the constructed THP-1 cell line (Figures 15B, C). RT-qPCR results demonstrated that low levels of DNASE1L3 promoted the differentiation of THP-1 cells toward the M1 phenotype while suppressing their differentiation into M2 macrophages (Figure 15D). Concurrently, there was a significant increase in IL-1β and TNF-α secreted by M1 macrophages, along with a marked reduction in LOX protein expression by M2 macrophages (Figure 15E). Flow cytometric analysis confirmed that reduced DNASE1L3 impaired macrophage efferocytosis (Figure 15F).