To comprehensively identify potential therapeutic targets for CRC, we compiled a broad list of druggable genes, which encode proteins amenable to modulation by small molecules or biologics (21). Two primary resources were used, including the Drug–Gene Interaction Database (DGIdb v4.2.0) (22) and a comprehensive review of human druggable genes (21). These sources integrate data from drug-target databases, clinical studies, and chemical biology literature to define the druggable genome. From these databases, we curated a total of 4,479 unique druggable genes ( Supplementary Table S1 ). This wide-scope inclusion strategy ensured that we captured both clinically validated targets and exploratory candidates relevant to complex diseases, including cancer. These genes were then used as input for subsequent eQTL and MR analyses to evaluate their causal relevance to CRC risk.

Cis-eQTLs were defined as genetic variants located within 1 megabase (Mb) of the transcription start site of a gene and significantly associated with its expression (minor allele frequency > 0.01). We obtained cis-expression quantitative trait loci (cis-eQTLs) data from the eQTLGen Consortium (https://eqtlgen.org/), which provides large-scale eQTL information derived from blood samples of 31,684 individuals of European ancestry (23). This dataset includes expression data for 16,987 genes and over 30,000 significant cis-eQTLs. Among the 4,479 druggable genes initially identified, 2,525 had corresponding cis-eQTLs available in the eQTLGen dataset and were included in downstream MR analysis. This dataset served as the basis for selecting genetic instruments to evaluate causal relationships between gene expression and CRC risk.

Summary-level GWAS data for CRC were obtained from previously published meta-analyses and biobank studies, all based on individuals of European ancestry (24). The GWAS statistics could be further obtained from the GWAS catalog (https://www.ebi.ac.uk/gwas/home). The discovery analysis included 78,473 CRC cases and 107,143 controls (GCST90255675). Two replication datasets were used: one from Sakaue et al. (6,581 cases and 463,421 controls; GCST90018808) (25) and one from the FinnGen project (11,790 cases and 378,749 controls) (26). To explore site-specific effects, we also conducted stratified analyses using GWAS data for colon cancer (3,793 cases and 410,350 controls; GCST90011811) and rectal cancer (2,091 cases and 410,350 controls; GCST90011810) (27). All participants provided informed consent, and ethical approvals were obtained from the relevant institutional review boards. Detailed cohort characteristics are provided in Table 1 .

To ensure robust causal inference in MR, we applied strict criteria for selecting single nucleotide polymorphisms (SNPs) as instrumental variables (IVs) for each druggable gene. First, we extracted cis-eQTL SNPs with genome-wide significance (P < 5.0 × 10^-8) (28) from the eQTLGen dataset, ensuring strong association with gene expression. To eliminate linkage disequilibrium (LD), SNPs were clumped using the 1000 Genomes Project (European population) reference panel, with an LD threshold of r² < 0.1 and a clumping window of 10,000 kb (29). SNPs incompatible across exposure and outcome datasets were excluded. For palindromic SNPs, we inferred strand orientation based on allele frequencies or removed them if this information was unavailable. The F-statistic was calculated for each IV using the formula F = (N − k − 1)/k × [R²/(1 − R²)] to assess instrument strength (30). Only SNPs with F-statistics > 20 were retained to avoid weak instrument bias. After filtering, a total of 40,356 SNPs were selected as IVs corresponding to 2,525 druggable genes, forming the basis for downstream MR analysis (31). Details about the IVs are shown in Supplementary Table S2 .

We conducted two sample MR analysis using the TwoSampleMR R package (version 0.6.6) to estimate the causal effects of gene expression on CRC, colon cancer, and rectal cancer risk. For genes with only a single valid instrumental SNP, we applied the Wald ratio method. When two or more independent SNPs were available, the inverse-variance weighted (IVW) method was used to derive overall causal estimates (32). To control for multiple hypothesis testing, we applied Bonferroni correction based on the total number of gene tested (33). Statistical significance thresholds were set at P < 2.02 × 10^-5. This corrected threshold ensured robust identification of causal genes while minimizing the false discovery rate.

To assess the robustness of the MR results and detect potential violations of MR assumptions, we conducted a series of sensitivity analyses for genes with significant causal associations. For genes with multiple IVs, MR-Egger intercept tests were performed to evaluate the presence of horizontal pleiotropy, with P < 0.05 indicating potential bias (34). We also calculated Cochran’s Q statistic under both the IVW and MR-Egger models to assess heterogeneity among the SNP-specific estimates (35). Additionally, a leave-one-out analysis was carried out by iteratively excluding one SNP at a time to identify any disproportionately influential variants. Only genes that met all quality control criteria, including consistent effect directions across MR methods, non-significant MR-Egger intercepts, and stable leave-one-out results, were retained as high-confidence targets for further investigation.

To explore the functional interactions among genes identified by MR, we constructed a protein-protein interaction (PPI) network using the STRING database (https://string-db.org), limiting results to “Homo sapiens” and applying a minimum confidence score of 0.15. The network included only interactions supported by experimental or database evidence. To further investigate the biological relevance of these genes, we performed functional enrichment analysis using the clusterProfiler R package (36). Gene Ontology (GO) enrichment was used to categorize genes based on biological process (BP), cellular component (CC), and molecular function (MF) (37). In parallel, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was conducted to identify significantly enriched signaling pathways (38). Only terms and pathways with P < 0.05 were considered statistically significant. These analyses provided insights into the potential mechanisms through which the causal genes may influence CRC development and progression.

Colocalization analysis was performed using the coloc R package (39) to determine whether gene expression and CRC risk share the same causal variant. SNPs were harmonized using the TwoSampleMR pipeline to ensure alignment between exposure and outcome datasets. Default prior probabilities were set to P1 = 1 × 10^-4 (association with gene expression), P2 = 1 × 10^-4 (association with CRC), and P12 = 1 × 10^-5 (association with both). The posterior probabilities (PP) were used to evaluate five hypotheses, with strong colocalization defined as PP.H4/(PP.H3 + PP.H4) > 0.7 (40, 41), indicating a high likelihood that both traits are influenced by the same causal variant. Genes meeting this threshold were prioritized as high-confidence druggable targets.

To assess potential off-target effects and horizontal pleiotropy of the prioritized druggable genes, we conducted a phenome-wide association study (PheWAS) using the AstraZeneca PheWAS Portal (https://azphewas.com) (42). This platform includes genotype-phenotype associations derived from approximately 450,000 UK Biobank participants, covering ~15,500 binary and ~1,500 continuous phenotypes (43). For each candidate gene, we queried associated variants to identify significant relationships with traits unrelated to CRC. A significance threshold of P < 1 × 10–⁶ was applied, in accordance with the portal’s recommended cutoff for suggestive associations, to control for multiple testing and reduce false positives. This analysis provided insight into the broader phenotypic impact of each gene, helping to evaluate the safety and specificity of targeting these genes in CRC therapy.

Given that most therapeutic agents target small-molecule proteins, we aimed to identify druggable targets among the 6 causal genes showing strong colocalization signals with CRC. To explore potential repurposing opportunities, we integrated our results with evidence from OpenTargets (44), DrugBank (45), and DGIdb (22). This integrative approach allowed us to prioritize candidate compounds, both approved and investigational, with favorable safety profiles that could be repositioned for CRC therapeutic applications.

The scRNA-seq data, comprising 33 samples from pre- and post-treatment tumor tissues and adjacent normal tissues of CRC patients, were retrieved from the GEO database (accession ID: GSE205506). We applied a strict quality control criteria to ensure data reliability. Cells were retained if they contained fewer than 20% mitochondrial gene content, expressed more than 200 genes, and had between 200 and 6000 detected genes, appearing in at least three cells. A total of 175,566 high-quality cells were selected for downstream analysis. To correct batch effects and improve clustering accuracy, we utilized the Harmony algorithm for data integration (46). The data was normalized using log-normalization, and the FindVariableFeatures function identified the top 2000 highly variable genes. Principal Component Analysis (PCA) was performed for dimensionality reduction, followed by soft k-means clustering using the Harmony package. Cells were then grouped into distinct clusters using the FindClusters function with a resolution of 0.6. Cell type annotation was conducted based on canonical marker genes, differential expression patterns, and known cellular profiles. Additionally, the expression patterns of potential targets were identified in bulk RNA-seq and scRNA-seq data. Bulk RNA-seq data, the TCGA-COAD cohort, was downloaded from UCSC Xena data portal (https://xena.ucsc.edu/). We first compared the expression levels of the targets between tumor and adjacent normal tissues. Subsequently, featureplots were utilized to visualize the expression patterns of target genes in CRC microenvironment. Additionally, we analyzed the association between the expression levels of target genes and key T cell exhaustion markers (PD-1, PD-L1, TIM-3) as well as the immune-suppressive cytokine IL-10 via GEPIA2 (http://gepia2.cancer-pku.cn/).

Paired tumor and adjacent normal tissues were collected from 31 patients with histologically confirmed CRC at Renji Hospital, Shanghai Jiao Tong University School of Medicine (Approval No.: KY2024-188-C). Tissue samples were immediately snap-frozen in liquid nitrogen and stored at -80°C until further processing. Total RNA was extracted using TRIzol reagent (Invitrogen) in accordance with the manufacturer’s instructions. Reverse transcription was performed using Superscript II reverse transcriptase (Invitrogen) and gene-specific primers. Real-time quantitative PCR (RT-qPCR) was conducted using the ABI Prism 7300 Sequence Detection System (Applied Biosystems) to assess the expression of LAMC1, TFRC, TNFSF14, PLK1, TYMS, TSSK6, and the reference gene GAPDH. Relative transcript levels were normalized to GAPDH expression using the relative standard curve method, following the manufacturer’s technical guidelines (Applied Biosystems). All data analysis adhered to the MIQE guidelines (47), ensuring experimental transparency and reproducibility. Detailed primer sequences were provided in Supplementary Table S3 . The statistical differences were analyzed by one-way ANOVA and t-test using GraphPad Prism 10.1.2 (San Diego, CA, USA). P < 0.05 were considered significantly different.

Immunohistochemical (IHC) images of both tumor and normal tissues were obtained from the Human Protein Atlas (HPA, https://www.proteinatlas.org/). However, IHC data of both normal and tumor tissues for TSSK6 were not available in the database and were therefore not included in the analysis.

Using a two-sample MR framework, we systematically assessed the causal relationship between the expression of 2,525 druggable genes and CRC risk. After Bonferroni correction (P < 2.02 × 10^-5), 47 genes demonstrated significant causal associations in the discovery cohort, including 45 via the IVW method and 2 via the Wald ratio method ( Figure 2 ; Supplementary Table S4 ). To ensure robustness, we conducted sensitivity analyses, including MR-Egger intercept, Cochran’s Q test, and leave-one-out tests. Four genes (HLA-DRB5, CTSF, TSPO, and TRPC6) were excluded due to evidence of horizontal pleiotropy (P < 0.05). The remaining 43 genes passed all quality control criteria and were retained for further validation and downstream analyses. These results provide a high-confidence set of candidate genes potentially involved in CRC pathogenesis and highlight their potential as therapeutic targets.

To elucidate the functional relationships and interaction landscape among the genes causally linked to CRC, we constructed a PPI network using the STRING database. The network was restricted to interactions in “Homo sapiens” and filtered with a minimum interaction score of 0.15 to ensure confidence in the associations ( Figure 3 ). The resulting network revealed a highly interconnected core, with genes such as PTPN11, CDC42, TFRC, HSP90AA1, and PLK1 serving as central hubs, suggesting key regulatory roles in CRC-related pathways. To further explore the biological relevance of these candidate genes, GO and KEGG pathway enrichment analyses were conducted using the R package clusterProfiler. GO analysis identified significant enrichment in terms associated with immune function and cell signaling, particularly lymphocyte-mediated immunity, dendritic cell migration, and protein localization to the nucleus ( Figure 4A ). Cellular component analysis highlighted localization to secretory granules, cytoplasmic vesicle lumens, and the external side of the plasma membrane, while molecular function analysis emphasized roles in protein tyrosine kinase activity and cytokine receptor activity. KEGG pathway enrichment revealed that the candidate genes were significantly associated with cancer-relevant signaling pathways ( Figure 4B ). The most enriched pathways included cytokine-cytokine receptor interaction, epithelial cell signaling in Helicobacter pylori infection, and leukocyte transendothelial migration. These findings underscore the potential involvement of the identified genes in inflammation, immune regulation, and tumor microenvironment remodeling-hallmarks of CRC development and progression.

To further refine the causal gene list and minimize potential false positives due to LD, we performed colocalization analysis on the 43 genes identified from MR. This approach assessed whether CRC-associated variants and gene expression signals shared the same causal variant. Six genes (TFRC, TNFSF14, LAMC1, PLK1, TYMS, and TSSK6) met this stringent criterion, indicating a high likelihood that gene expression and CRC risk are driven by the same genetic variants ( Table 2 ). These genes were subsequently prioritized as high-confidence druggable targets. The remaining 37 genes, while suggestive, did not reach the colocalization threshold and were retained as secondary candidates for further evaluation. Sensitivity analyses supported the robustness of these six genes, showing no evidence of horizontal pleiotropy or significant heterogeneity ( Table 3 ). Consistent effects across alternative MR models and leave-one-out analyses further confirmed the stability of these associations ( Supplementary Figures S1 , S2 ). These results highlight a subset of genes with both strong causal and colocalization evidence, warranting focused investigation in downstream validation and translational studies.

In the replication phase, among the 6 genes with significant colocalization evidence, 2 genes (LAMC1 and TNFSF14) were validated to be significant in all three replication cohorts (P < 0.05), and another 2 genes (TFRC and PLK1) were validated to be significant in at least one replication cohort (P < 0.05) ( Figure 5 ). Additionally, the 6 genes with significant colocalization evidence were utilized for further stratified analysis in rectal cancer and colon cancer, respectively. In rectal cancer, two of the 6 genes (LAMC1 and PLK1) were identified with significant casual effect (P < 0.05) ( Figure 6 ). In colon cancer, 2 of the 6 genes (LAMC1, and TSSK6) were identified with significant casual effect (P < 0.05) ( Figure 6 ). Sensitivity analyses supported these findings, as no evidence of horizontal pleiotropy was detected for these genes ( Supplementary Tables S5 , S6 ). Consistent trends across four additional MR models and stability in the leave-one-out analysis confirmed the robustness of these findings ( Supplementary Figures S3 – S6 ).

The phenome-wide scan explored the relationships between potential drug targets and a variety of traits, providing valuable insights into possible side effects during drug development. This analysis revealed that none of the thirteen genes with strong colocalization evidence were linked to multiple phenotypes (P < 5 × 10^-8) ( Figures 7A–F ).

To assess therapeutic potential, we investigated drugs targeting the six key genes (TFRC, TNFSF14, LAMC1, PLK1, TYMS, and TSSK6) using DGIdb, OpenTargets, and DrugBank. Several approved or investigational drugs were identified ( Table 4 ). LAMC1 is targeted by agents like Ocriplasmin and Lanoteplase, mainly used for non-cancer conditions, but may be repurposed due to its role in the tumor microenvironment. TFRC, involved in iron metabolism, is already targeted by various approved iron supplements, and Pabinafusp alfa is under investigation for other diseases, showing potential for CRC applications. TNFSF14, a key immune regulator, is targeted by Baminercept, suggesting possible synergy with immunotherapies. PLK1 has multiple targeted drugs, including investigational anticancer agents like Volasertib and Onvansertib, with several approved drugs also showing potential off-target interactions. TYMS is already a standard chemotherapy target in CRC, with approved drugs such as Fluorouracil, Capecitabine, and Raltitrexed confirming its therapeutic relevance. These findings support drug repurposing strategies for CRC, potentially accelerating the development of targeted therapies.

After applying the Harmony algorithm, the cellular distribution within each sample remained largely consistent, indicating the absence of significant batch effects among the samples, making them suitable for downstream analysis ( Figure 8A ; Supplementary Figure S7 ). A total of 18 distinct cell clusters were identified using a resolution parameter of 0.4 ( Figure 8B ). We reanalyzed 175,566 scRNA-seq cells spanning 26 samples, including primary tumor, adjacent normal tissues and PBMC samples. Based on classical marker genes, we successfully classified various cell populations, including epithelial cells, B cells, plasma cells, T/NK cells, myeloid cells, fibroblasts, myofibroblasts and endothelial cells ( Figure 8C ). The marker genes corresponding to each cell type displayed distinct expression patterns, reinforcing the accuracy of our cell annotation ( Figures 8D, E ).

Detailed information on the expression levels of the 6 candidate genes is presented, along with information regarding the different expression patterns within tumor and normal samples ( Figures 9A–F ). The expression levels of TFRC, LAMC1, PLK1, TYMS, and TSSK6 were significantly higher in tumor tissues while TNFSF14 exhibited higher expression levels in normal tissues. Notably, LAMC1 showed high expression levels in endothelial cells, fibroblasts and myofibroblasts, and TFRC in epithelial cells. The correlation analysis of the expression levels of target genes with key T cell exhaustion markers (PD-1, PD-L1, TIM-3) and the immunosuppressive cytokine IL-10 demonstrated that TNFSF14, LAMC1, and TYMS were significantly and positively associated with these immune markers ( Supplementary Figure S8 ).

We performed RT-qPCR to validate the results of our previous analyses, where consistent expression trends were observed ( Figure 10A ). Specifically, TFRC, PLK1, TYMS, and TSSK6 were remarkably upregulated in tumor tissues compared to normal controls (P < 0.0001). Similarly, LAMC1 showed a significant increase (P < 0.05), while TNFSF14 was significantly downregulated in tumor tissues (P < 0.0001). Consistent with these findings, the IHC staining results of the five targets demonstrated similar trends at the protein level ( Figure 10B ). Strong immunoreactivity of TFRC, LAMC1, PLK1, and TYMS was observed in tumor tissues, whereas weak or no staining was observed in normal tissues. TNFSF14 showed similar staining intensities in tumor and normal tissues.