Mixed-model and transcriptome-wide association analyses identify transcription factors and genes associated with colorectal cancer susceptibility

Domain: biology
Published: 2026-01-15
Editorial outcome: major_revision

Abstract
Susceptibility transcription factors (TF) whose DNA bindings are altered by genetic variants regulating colorectal cancer (CRC) risk genes remain poorly defined. Using generalized linear mixed models, we analyze 218 TF ChIP-Seq datasets alongside GWAS data from 100,204 CRC cases and 154,587 controls of East Asian and European ancestries. We identify 51 TFs and TF-cofactor interactions, including VDR-cofactors, as key regulators of CRC risk. Integrating these TF insights with transcriptome-wide association studies (TWAS), we further evaluate associations between genetically predicted gene expression, alternative splicing, and alternative polyadenylation with CRC risk, using RNA-seq data from 364 Asian-ancestry and 707 European-ancestry individuals. Multi-ancestry TWAS identify 222 risk genes, including 95 novel genes and 48 potentially druggable targets. Single-cell analysis provides additional functional evidence supporting ~45% of these genes, and experimental validation confirms oncogenic roles for RHPN2, IRS2, and TXN. Our findings elucidate key TF–gene regulatory networks and uncover novel CRC risk genes.

Introduction
Transcription factors (TFs) play a pivotal role in mediating the effects of genetic risk variants by regulating gene expression in a context-specific manner, thereby influencing disease risk1,2. A better understanding of the regulatory mechanisms by which TFs mediate these genetic variants is critical for elucidating the molecular basis of complex traits, including cancer. Previous studies have identified susceptible TFs through enrichment analyses of genome-wide association studies (GWAS)-identified risk variants in TF binding sites through chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq)3–6. However, these approaches often suffer from suboptimal statistical power, as they focus primarily on a narrow range of GWAS-identified risk variants. Recently, we developed a statistical method using a generalized linear mixed model (GLMM) to enhance power by analyzing TF ChIP-seq data in conjunction with GWAS summary statistics of variants across the genome7. This method evaluated associations between variants’ risk (e.g., Chi-squared values from GWAS) and their binding status by a TF, and has identified known susceptible TFs, such as FOXA1, ESR1 and GATA3 for breast cancer, by showing significantly higher risk for variants occupied by these TFs compared to genome-wide controls7. However, the specific TFs involved in colorectal cancer (CRC) susceptibility remain unclear. Thus, identifying TFs that interact with CRC risk variants could provide insights into TF-mediated gene regulation and enhance our understanding of CRC susceptibility.

GWASs have successfully identified more than 250 risk loci for CRC1,2,8, and integrative analyses of expression quantitative trait loci (eQTL) and these GWAS loci have revealed numerous candidate susceptibility genes9–12. Further transcriptome-wide association studies (TWAS) have linked predicted gene expression levels with disease risk, offering higher statistical power and potentially revealing novel associations overlooked in GWAS13,14. Prior research, including our own, has conducted TWAS in CRC and identified more than 200 putative susceptibility genes8,15,16. However, conventional TWAS methods are prone to false positives due to the inclusion of non-regulatory variants17–19. To address this, we developed the sTF-TWAS approach, which refines gene expression prediction by focusing only on regulatory variants within TF-occupied cis-regulatory elements (STFCREs)7,20. This approach significantly enhances disease gene detection, outperforming conventional TWAS methods20. In addition, our recent alternative polyadenylation (APA)-WAS framework has uncovered new cancer risk genes21, while its application to CRC has yet to be explored. Moreover, while TWAS has primarily been conducted in European populations, CRC studies in East Asian populations are limited due to a lack of relevant RNA-seq data from colon tissue.

In this study, we integrate 218 ChIP-seq datasets from CRC-related cells with GWAS data from 100,204 CRC cases and 154,587 controls of both European and East Asian ancestries to identify susceptible TFs for CRC. We further utilize RNA-seq data from 364 Asian-ancestry and 707 European-ancestry individuals to build ancestry-specific models for gene expression, alternative splicing (AS), and APA) under the sTF-TWAS framework. Our goal is to identify novel CRC susceptibility genes through multi-ancestry TWAS analyses, providing a deeper understanding of the genetic basis of CRC.

Results
Identification of susceptible CRC-associated TFs and key TF-cofactor interactions
We integrated 218 ChIP-seq data for 84 TFs (Supplementary Data 1) with CRC GWAS summary statistics from European and East Asian ancestry to identify CRC susceptible TFs (Fig. 1A). At Bonferroni-corrected P < 0.05, we identified 51 TFs associated with CRC susceptibility. Among the top TFs were VDR, MYC, JUN, CDX2, NIPBL, MED12, and ETV5 (Fig. 1B and Supplementary Data 2, “Methods”). We conducted motif analysis for each TF-binding region derived from ChIP-seq data and observed significant enrichment for host motifs, including those for VDR (24.2%), JUN (11.1%), KLF5 (49.7%), and HNF4A (36.4%) (Fig. 1C). We further explored the impact of genetic variations in co-occupied cis-regulatory regions of these 51 TFs and observed a substantial proportion of genetic variants located in co-occupied binding sites (Fig. 1D and Supplementary Data 3). Pair-wise analyses identified 154 significant interactions at Bonferroni-corrected P < 0.05, showing that genetic variants in regions co-occupied by multiple TFs exhibited stronger associations with CRC risk compared to variants bound by a single TF (Fig. 1D and Supplementary Data 3). Notably, VDR emerged as a key player in these interactions, highlighting its central role in CRC susceptibility (Fig. 1D and Supplementary Data 3). In addition, we applied stratified LD Score Regression (S-LDSC)22 as an alternative approach and found that 40 of the 51 risk TFs identified by our method were also detected by S-LDSC at Bonferroni-corrected P < 0.05. Most of the remaining 11 TFs showed nominal significance (P < 0.05) (Supplementary Fig. 1A), providing further support for the validity of our findings. We also conducted simulation analyses to compare the statistical power of our GMM with S-LDSC. Our results demonstrate that GLMM consistently outperformed S-LDSC in detecting causal TFs across varying levels of heritability (Supplementary Fig. 1B). Quantile-quantile (Q-Q) plots demonstrate that Type I error rates are under control for both approaches, with overall comparable performance (Supplementary Fig. 1C–F). These simulation results are consistent with our observations in real data, supporting GLMM as a robust approach for identifying risk TFs in colorectal cancer.

Associations of genetically predicted gene expression with CRC risk
For the multi-ancestry TWAS under the sTF-TWAS framework, we built gene expression prediction models using RNA-seq data from 364 participants of East Asian ancestry (from the ACCC) and 707 participants of European ancestry (N = 423 from the BarcUVa-seq and N = 284 from the GTEx projects) (Methods)16. A total of 7096 (for ACCC), 14419 (for BarcUVa-seq), 11361 (for GTEx) gene expression prediction models were built with a coefficient of determination R2 > 0.01 (10% correlation) using the elastic net approach (Supplementary Data 4). We then applied these ancestry-specific prediction models to summary statistics of CRC GWAS from 80,774 cases and 105,298 controls of European ancestry and from 23,572 cases and 48,700 control of East Asian ancestry. At a Bonferroni-corrected P < 0.05, we identified 138 genes whose genetically predicted expression was associated with CRC risk after applying a meta-analysis of TWAS results based on models built in both European and East Asian populations (Supplementary Fig. 2, “Methods”). In addition, we evaluated the associations separately in European and Asian populations. At the Bonferroni-corrected threshold, we identified two additional genes (SENP8 and PCNT) from the European population, while no additional significance was identified in the East Asian population (Supplementary Data 5).

Together, 140 genes showed an association between their genetically predicted expression and CRC risk (Fig. 2A). Of them, 41 (29% of 140) genes were not previously reported to associate with CRC risk1,2,8,10,15,16 (Table 1). Of note, seven putative susceptibility genes (PPAN, DERL3, PROX1, TMEM87B, SV2A, FBXL20, and MED1) are located at loci, with >2 Mb away from the closest GWAS-identified risk variants (Table 1). To assess whether the genes identified were independent of the established GWAS association signals, we conducted conditional analyses for their associations with CRC risk, adjusting for the associations with the closest lead SNP for each locus (“Methods”). Of the 99 previously reported genes located in GWAS-identified risk loci, conditional analysis showed that 39 genes remained statistically significant at a nominal P < 0.05 (Supplementary Data 6).

Additional risk genes identified by AS-WAS and APA-WAS
We built prediction models for AS and APA using the aforementioned RNA-seq data in colon tissues, generating 9680 AS models for ACCC, 10,025 for BarcUVa-seq, and 20,240 for GTEx, at prediction performance R2 > 0.01 (Methods, Supplementary Data 7). A meta-analysis of AS-WAS results from these models, we identified 82 risk genes at a Bonferroni-corrected threshold P < 0.05 (Fig. 2B). We identified three additional genes, POLD3, PARP11, and CEP89, specific in European populations (Supplementary Data 8). In total, we identified 85 genes after combining the results from the meta-analysis and population-specific analysis. For APA, a total of 2569 (for BarcUVa-seq), 4399 (for GTEx), and 2744 (for ACCC) prediction models were built successfully with R2 > 0.01 (Supplementary Data 9). After meta-analysis of APA-WAS results from these models, we identified 29 risk genes, at a Bonferroni-corrected P < 0.05 (Fig. 2C). Similarly, we identified four additional genes (ARPC5, PHLDB2, COA4, NEIL1) specifically in the European population (Supplementary Data 10). In total, we identified 33 genes after combining the results from the meta-analysis and population-specific analysis.

In total, 112 genes were identified for CRC risk combining results of SP-TWAS and APA-WAS. Of them, 58 were not previously reported for CRC risk (Table 2, 3). Among the remaining 54 previously reported genes, conditional analysis showed that 24 genes remained statistically significant at a nominal P < 0.05 (Supplementary Data 11).

Risk genes identified by combining TWAS, AS-WAS, and APA-WAS analyses
By integrating results from TWAS, AS-WAS, and APA-WAS, we identified 222 unique CRC risk genes, 95 of which had not been previously associated with CRC (Fig. 2). Among them, four genes (COLCA2, METRNL, ARPC5, COX14) were commonly identified by all three TWAS analyses and 28 genes commonly detected by any two of the TWAS analyses (Supplementary Fig. 3). By comparing genes reported in previous CRC TWAS and other genetic studies1,2,8,10,15,16, 95 of them were not previously identified as CRC susceptibility genes, particularly including 13 genes reported in new loci (Tables 1–3).

Using scRNA-seq data to explore risk genes potentially contributing to CRC development
We next analyzed scRNA-seq datasets spanning CRC developmental stages, including colon normal, serrated polyps, conventional adenomas, and carcinomas tissues from the COLON MAP23. Our scRNA-seq data analysis identified several major cell types (e.g., ABS, GOB, STM) in normal tissues, SSC in serrated polyps, and ASC in conventional adenomas, consistent with previous studies23 (Supplementary Fig. 4). Differential analysis of these major cell types was performed using DESeq224 to evaluate their contribution to two CRC etiology pathways: the normal-serrated polyp-carcinoma pathway and the normal-adenoma-carcinoma pathway (“Methods”). At an FDR-adjusted P < 0.05, 48 genes exhibited significantly differential expression with consistent associations from TWAS analysis across one of these major cell types in the normal-serrated polyp-carcinoma pathway (Supplementary Data 12), while 92 genes showed significant differential expression in the normal-adenoma-carcinoma pathway (Supplementary Data 13). Combined with these results, we provide additional evidence that 99 genes (44.6% of the identified 222 risk genes) are potentially involved in CRC development.

Druggable genes for CRC prevention and intervention
We comprehensively annotated 222 identified risk genes as therapeutic targets for existing drugs, including those either approved or in clinical development (“Methods). Our analysis identified nine CRC susceptibility genes that are druggable by 125 drugs approved or in clinical trials for CRC treatment (Fig. 3 and Supplementary Data 14). In addition, we identified 26 genes targeted by 519 drugs approved or in clinical trials for the treatment of other cancers (Supplementary Data 15). Furthermore, 40 genes were identified as targets of 585 drugs approved or in clinical trials for non-cancer indications (Supplementary Data 16). In total, 48 of the 222 CRC susceptibility genes (21.6%) were found to be druggable, with potential links to 1216 drugs, either approved or under investigation.

Functional assays for three putative oncogenes, RHPN2, IRS2, and TXN
To functionally examine the CRC susceptibility genes, we performed CRISPR/Cas9-mediated gene knockout. Three risk genes, including one previously reported gene, RHPN2, and two newly identified risk genes, IRS2 and TXN, were selected for downstream experiments based on evidence of their potential oncogenic functions and literature review (Methods). This was supported by their predicted elevated expressions being associated with an increased risk of CRC, as well as differential analysis of single-cell data highlighting their oncogenic roles. For each gene (RHPN2, IRS2, or TXN), we employed two independent sgRNAs for protein coding region targeting (Supplementary Data 17). We used western blotting for confirming target protein reduction and CCK-8 proliferation and colony formation (Fig. 4) assays for gene’s functions. TXN proteins displayed 65-98% reductions between knockout cells and control cells (sgRNA targeting LacZ) in the CRC cell lines (RKO, HCT116, and SW480) (Supplementary Fig. 5). We showed that TXN knockout significantly inhibited cell proliferation and colony formation in all CRC cell lines (Fig. 4A–C). IRS2 proteins had 66–86% reductions between knockout cells and control cells in three CRC cell lines (Supplementary Fig. 5). There was only significant inhibition of colony formation (Fig. 4E, F), indicating a long-term effect of knocking out the IRS2 protein. RHPN2 proteins had 34–61% reductions between knockout cells and control cells (Supplementary Fig. 5). However, the data showed a significant inhibition of cell proliferation and colony formation only in RKO cells, together with a trend of inhibition in HCT116 cells (Fig. 4G–I). There was a minimal effect in SW480 cells by knocking out RHPN2. The variable extents of RHPN2 dependency in cell proliferation and growth were matched to the consistent protein expression levels in these cells. Our experimental verification demonstrated the potential oncogenic functions of TXN, IRS2 and RHPN2, aligning with our TWAS analysis.

Methods
Study populations
This study utilized a meta-analysis summary statistics dataset of 31 GWAS comprised of 186,072 individuals of European ancestry (80,774 CRC cases and 105,298 controls) and 72,272 individuals of Asian ancestry (23,572 CRC cases and 48,700 controls). The meta-analysis was conducted using the fixed-effects inverse variance weighted model implemented in METAL58. Details on sample selection and matching, sample numbers, and demographic characteristics of study participants have been described previously59–62. All participants provided written informed consent, and all study protocols were approved by the relevant research ethics committee or institutional review board.

We included RNA-seq data from tumor-adjacent normal colon tissues and genotyping data of 364 individuals of East Asian ancestry who participated in the Asia Colorectal Cancer Consortium (ACCC). We also included data from individuals of European ancestry, including RNA-seq in normal colon tissues and genotyping data of 423 individuals in the BarcUVa-Seq project and 284 individuals in the GTEx project. The sample collection for these data has been described in detail in previous studies1,63,64.

Generalized linear mixed models
We systematically analyzed TF ChIP-seq data generated in normal colorectal epithelium and CRC cell lines from the Cistrome database65. Only ChIP-seq datasets meeting Cistrome’s recommended quality control standards were included. After evaluation, we curated 218 ChIP-seq datasets corresponding to 84 TFs. To ensure robust analysis, we focused on datasets with more than 1000 peaks and a fold change greater than 10. For TFs with multiple ChIP-seq datasets, only the dataset with the highest number of peaks was retained. This process resulted in a high-quality set of ChIP-seq data for 69 TFs, which were used in downstream analyses.

We used our recently developed generalized linear mixed model approach7 to investigate associations of CRC risk with variations of TF-DNA binding by a single TF (Eq. 1). In brief, a N x M matrix for all genetic variants from the GWAS summary statistics and the annotation from all available TF-DNA binding regions (“N” refers to the number of genetic variants and “M” refers to the number of TF ChIP-seq datasets; the variant in ChIP-seq peak is coded as “1” and the one not in ChIP-seq peak is coded as “0”). We used the Chi-squared value for each genetic variant reported in the GWAS summary data to measure its association with CRC risk. We then used generalized linear mixed models to estimate the associations between the Chi-squared values (Y) and TF binding status of genetic variants located in binding sites of each TF, given LD blocks (non-overlapping segments of 100kb) of genetic variants to handle the dependence between genetic variants (Eq. (1)).

Yij = β0 + β1 TFij + Vi + εij (1)

In Eq. (1), Yij is the Chi-squared value for the j-th variant in the i-th LD block; β0 is the fixed effect of the Chi-squared values (ΔX2), β1 is the fixed effect of the TF status; TFij is the j-th TF value (i.e., 1 for a variant located in a TF binding site, 0 otherwise) in the i-th LD block; Vi is the random intercept for the i-th LD block; and εij is the error term. Based on this statistical model, we identified cis-regulatory elements occupied by TFs whose genetic variations of TF-DNA bindings are associated with target cancer risk at Bonferroni-corrected P < 0.05. We additionally used generalized linear mixed models to estimate the associations of the Y values of variants with the TF-pair occupancy if they showed a significant interaction (Eq. (2)). Of note, we used genetic variants non-occupied by any of the TFs investigated as the referenced control group for all the analyses.

Yij = β0 + β1 TF1ij + β2 TF2ij + β3 TF1ij × TF2ij + Vi + εij (2)

Simulation studies
To benchmark the statistical power of GLMM against S-LDSC, we conducted simulations under a polygenic causality scenario where genotypes influence phenotypes through TF binding sites. Specifically, we assumed that only SNPs located within TF binding sites contribute to disease risk, while SNPs outside these regions have no effect.

For each replicate, we randomly selected K causal variants (K=500, 1000, 5000, or 10,000) from a pool of 600,000 SNPs annotated to 84 TFs in CRC cell lines (Supplementary Data 1). Under an additive architecture, phenotypes were simulated using Eq. 3:

Y*i = K j=1 β1 × XTFij + ε (3)

where, XTF is the genotype (0, 1, 2) of j-th variant for subject i use data from 489 European ancestry from 1000G project phase 366. β1 was set to 1 and ε is sampling from a normal distribution with heritability h2 = 0.5, 0.25, 0.1.

To facilitate downstream GWAS, we binarized the continuous trait into a case–control phenotype C* based on the mean value of Y* across individuals using Eq. (4).

C*i = 1, if Y*i ≥ median(Y*i); 0, if Y*i < median(Y*i) (4)

We performed GWAS using PLINK 2.067 and applied both GMM and S-LDSC to the resulting summary statistics to estimate TF–disease associations. For each K, we conducted 50 replicates, each with a different random subset of causal SNPs, ensuring reproducibility while keeping K fixed. Results were summarized as the number of true-positive TFs detected at a nominal threshold (P < 0.05) by each method.

To evaluate the type I error rate of GLMM, we conducted simulation analyses under the null. Specifically, we randomly assigned phenotype values (0 or 1) to 489 individuals of European ancestry and performed logistic regression to generate GWAS summary statistics under the null. We then applied GLMM and S-LDSC using different numbers of TF binding variants (500, 1000, 5000, and 10000) and repeated the entire procedure 1000 times. Finally, we compared the Type I error rates of GLMM and S-LDSC using quantile-quantile (QQ) plots of p-values across the 1000 replicates for each binding variant set size.

Gene expression profiling and data processing
The RNA-seq library preparation, sequencing, data processing, and quality control for the ACCC1 and the BarcUVa-Seq project63 have been described previously. In brief, paired-end sequencing was conducted for mRNA obtained from the BarcUVa-Seq project using either an Illumina HiSeq 2500 or a NovaSeq 6000 instrument and for mRNA obtained from the ACCC project using a HiSeq 3000 Sequencer. The raw RNA-seq data were processed following the GTEx Consortium pipeline64. BAM files containing aligned reads to the reference genome (GRCh37 for ACCC and GRCh38 for BarcUVa-Seq and GTEx project), were generated using STAR68 (v2.5.4). These BAM files, along with GENCODE release 19 (for ACCC) and 26 (for BarcUVa-Seq and GTEx) were used for quantification of gene expression, AS and APA.

Building on our previous work in processing gene expression and AS16, gene expression quantification was performed using RNA-SeQC69. Read counts and Transcripts Per Kilobase Million (TPM) values were calculated from aligned reads based on the following criteria: (1) reads were uniquely mapped; (2) aligned reads were properly paired; (3) the read alignment distance was <6. Genes with expression thresholds of ≥0.1 TPM in ≥20% of samples and ≥6 reads (unnormalized) in ≥20% of samples were selected. Quantile normalization was applied to standardize the gene expression data within each sample to the same distribution, followed by rank-based inverse normal transformation across all samples.

AS was quantified using LeafCutter70 and RegTools71. In brief, the intron usage was quantified from aligned BAM files using the junction extract from the RegTools with the default options. Next, the introns were clustered together using the leafcutter_cluster_regtools.py script from LeafCutter, with the following options: -m 50 -l 500000. The clusters were mapped to each gene using the map_clusters_to_genes.R script from the cluster_prepare_fastqtl.py provided by the GTEx. Introns with low counts or low complexity were filtered based on the criteria in the in the GTEx pipeline (e.g., cluster_prepare_fastqtl.py). The filtered counts were normalized using the prepare_phenotype_table.py script from LeafCutter.

For APA analysis, following our recent APA-WAS study21, APA was qualified using the percentage of PDUI estimated from DaParsv2.021. We performed quantile normalization to transform the quantified PDUI values of APA for each sample to the same distribution. We used 60 probabilistic estimation of expression residuals (PEER) factors72 to correct batch effects and experimental confounders in downstream prediction model building for gene expression, AS and APA.

Genotyping data processing
The generation and processing of genotyping data have been described in our previous work1,16. Briefly, CRC patients of East Asian ancestry from the ACCC were genotyped using the Illumina OncoArray and the expanded Illumina MEGA Array. For genotype quality control, genetic variants were excluded based on the following criteria: (1) genotype call rate <95%; (2) ambiguous variants; (3) duplicated variants; (4) P for Hardy-Weinberg equilibrium (HWE) <1.0 × 10−6. For quality control of samples, we excluded those with: (1) genotype call rate <95%; (2) genetically identical or duplicated samples; (3) first- or second-degree relatives; (4) racial outliers. To estimate the genetic relatives among samples, the quality-controlled common SNPs (minor allele frequency [MAF] >0.01) were first pruned based on the window size=1500 SNPs and pairwise r2 >0.2. Identity-by-descent (IBD) was calculated using the 1000 Genomes Project Phase 3 mixed reference haplotypes with the Michigan Imputation Server (Minimac4 for imputation and Eagle v2.4 for phasing). Variants with MAF <0.01 were excluded before imputation. Approximately 4.5 million variants with MAF >0.05 and R2 >0.3 were included for prediction model building.

For the BarcUVa-Seq project, approximately 400,000 SNPs were genotyped from DNA isolated from peripheral blood samples, using the OncoArray BeadChip. Quality controls and genotyping data processing details were provided in a previous study63. Genotypes were aligned to the reference genome GRCh38 and imputed using the TOPMed (Version R2) reference panel63. Genetic variants with an imputation quality of R2 <0.3 were excluded.

In the GTEx project, whole genome sequencing (WGS) was performed on DNA samples from 284 individuals of European ancestry. WGS-based genotype data for these individuals, in VCF format, were accessed from the dbGap (phs000424.v8.p2). Detailed information on WGS variant quality control is provided in GTEx documentation64.

Conditional analyses based on adjusted GWAS association signals
To investigate whether the identified associations between genetically predicted gene expression (or AS, APA) and CRC risk are independent of the association signals identified in previous GWAS, we performed conditional analyses by adjusting for the lead variant using the GCTA-COJO method73. The lead variant in each risk locus was determined based on the strongest association from GWAS summary statistics derived from the meta-analysis of European and East Asian populations or from Europeans alone. Genotype data from either European or East Asian populations in the 1000 Genomes Project Phase 3 were used to estimate linkage disequilibrium (LD). We then conducted TWAS analyses (or AS, APA) using the adjusted GWAS summary statistics to explore the association between genetically predicted gene expression and CRC risk. In addition, joint analyses were performed to combine p-values for the European population, as well as datasets from both European and East Asian populations, for gene expression (or AS, APA).

Building genetically predicted models of gene expression, AS, and APA
To construct genetically predicted models for gene expression, AS, and APA, based on our sTF-TWAS framework20, we selected only genetic variants that located in the binding regions of any susceptible TFs for variants with MAF >0.05, and HWE P < 10−6. We also excluded variants missing in CRC GWAS summary statistics from either European or Asian populations, which resulted in approximately 378,000 genetic variants. Prediction models were trained using the elastic net approach based on their flanking variants within a ±1 Mb region, accounting for potential confounding factors, such as the top five principal components (PCs), gender, potential batch effects, and other factors derived from the PEERs method. A tenfold cross-validation approach was used to fine-tune the model parameters. The performance of the prediction models was evaluated by calculating the squared correlation (R2) between the predicted values and the observed levels for gene expression, AS and APA, respectively.

Association analyses of genetically predicted gene expression, AS, and APA with CRC risk
Using the weight matrix and summary statistics from CRC GWAS datasets consisting of 186,072 individuals of European ancestry and 72,272 individuals of East Asian ancestry, we evaluated the association between gene expression (or AS, APA) and CRC risk under the sTF-TWAS framework20. The statistical method was outlined in Eq. (3)

Z = ωl β^l / se(β^l) σ^g / σ^l (5)

Briefly, the Z-score was used to estimate the association between predicted gene expression (or AS, APA) and CRC risk, ωl is the weight of variant l for predicting the expression of gene g (or ASg, APA), β^l and se(β^l) are the association regression coefficient and its standard error, respectively, for variant l in GWAS, and σ^l and σ^g are the estimated variances of variant l and the predicted expression of gene g (or ASg, APA). We further conducted a joint analysis to combine TWAS p-values for each gene from the BarcUVa-Seq and GTEx datasets for the European population, as well as from datasets across European and East Asian populations (ACCC), using the aggregated Cauchy association test (ACAT) method74. Similarly, joint analyses were also performed to combine p-values for AS and APA, respectively. Statistically significant genes were determined based on Bonferroni corrected P < 0.05, considering the total number of tests for the models of gene expression, AS, and APA, respectively.

Differential gene expression analysis from single cell data
We downloaded and analyzed single-cell RNA sequencing (scRNA-seq) data spanning CRC developmental stages, including colon normal tissues (n=31), serrated polyps (n=19), conventional adenomas (n=29), and carcinomas (microsatellite stability [MSS], n=17; and microsatellite instability-high [MSI-H], n=15) from the Colorectal Molecular Atlas Project (COLON MAP)23. We filtered the cell count matrix for each sample by removing low-expressed genes (expressed in ≤3 cells) and excluded cells with high mitochondrial activity (≥20% of transcripts attributed to mitochondrial genes) using Seurat (version 4.2.1)75,76. We next identified the top 2000 highly variable genes for principal component analysis (PCA) and clustered cells based on the first 30 principal components using the Shared-Nearest-Neighbors (SNN) graph and the Louvain community detection algorithm. These clusters were visualized using the dimensionality reduction technique Manifold Approximation and Projection (UMAP).

For differential expression analysis across CRC developmental stages, we examined risk genes associated with the normal-serrated polyp-carcinoma pathway and the normal-adenoma-carcinoma pathway. We performed integrative analysis of the cell count matrix across stages, generating a pseudobulk gene count matrix for each cell type within each tissue sample. Using these pseudobulk count data, we applied DESeq224 to identify differentially expressed genes between developmental stages. DESeq2 accounts for heterogeneous variability (heteroscedasticity) in gene expression, where variance depends on mean values, by modeling this relationship. This approach effectively adjusts for differences in variability between low- and high-expressing genes, leading to more accurate differential expression analysis. Specifically, we focused on potential CRC cell origins in major cell types, including absorptive cells (ABS), goblet cells (GOB), and stem cells (STM), comparing normal tissues with serrated polyps and normal tissues with conventional adenomas. Furthermore, we analyzed differential gene expressions between adenoma-specific cells (ASC) and MSS, as well as ASC and MSI-H. Similarly, differential gene expression analyses were performed between serrated-specific cells (SSC) and MSS, and SSC and MSI-H. Statistically significant genes were considered as a False Discovery Rate (FDR)-adjusted p-value < 0.05.

Identifying druggable targets for the identified risk genes
We identified potentially druggable risk genes and associated drugs by integrating drug-protein interaction data from four well-established drug databases: DrugBank52, ChEMBL53, the Therapeutic Target Database54 (TTD) and Open Targets55. Specifically, we extracted key information for each drug, including commercial names, maximum clinical trial phases, and indications. Using this data, we annotated the identified risk genes as druggable targets and linked their corresponding drugs that are either approved or in clinical trials. In addition, we categorized therapeutic drugs based on their indications, grouping them into those targeting CRC, other types of cancer, and non-cancer diseases.

Cell culture and cell proliferation and colony formation assays
Human colorectal cancer cell lines (HCT116, RKO, and SW480) from ATCC were cultured in RPMI 1640 supplemented with 10% FBS and penicillin/streptomycin. Cells were kept in a humidified incubator at 37˚ with 5% CO2. CCK-8 cell proliferation was performed for four days. In detail, 5×103 cells were seeded into a 96-well plate and cultured for four days. The CCK-8 reagent (APExBIO) was then added into each well, and the plate was read at OD 450 nm using a microplate reader (BioTek). Relative cell growth was calculated by normalizing the OD value of knockout cells to the control cells. For long-term effect, colony formation assay was performed for 10–15 days. In detail, 200 cells (HCT116 and RKO) or 500 cells (SW480) were seeded into a 12-well plate. Cell colonies were visualized by crystal violet staining. The colonies were photographed, and numbers were counted.

CRISPR/Cas9-mediated gene knockout
For gene knockout, the lentiviral vector LentiCRISPRv2 (Addgene) was used to clone sgRNA of the individual target gene. Lentiviruses were amplified in HEK293T cells. Briefly, lentiviral vector and packaging vectors, psPAX2 and pVSV-G, were transfected into 293T cells with lipofectamine 2000 (Thermo Fisher Scientific). Viral supernatants were collected by centrifugation at 48 and 96 hours after transfection. CRC lines were infected with the corresponding lentiviruses for 24 h. Cells were changed to new growth medium for another 24 h and then selected with puromycin (2μg/ml, Sigma-Aldrich) for 5 days. The CRISPR/Cas9 edited cells were verified by western blotting of the target protein before using for the functional assays.

Western blotting
Total proteins were extracted from the corresponding cells lysed with RIPA buffer plus protease and phosphatase inhibitors (Roche). Twenty ug of total proteins were separated by SDS-PAGE, transferred to PVDF membrane and probed with the indicated primary antibodies. Anti-rabbit or anti-mouse secondary antibody was subsequently applied, and western blot signals were detected using enhanced chemiluminescence substrate (West Femto, Thermo Fisher Scientific). β-Tubulin was used as the loading control. Band intensity was calculated with NIH ImageJ software.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript integrates several ChIP-seq datasets with colorectal cancer (CRC) GWAS summary statistics to nominate transcription factors (TFs) of relevance to CRC by broadly searching for enrichment of TF binding in GWAS loci. The authors also perform extensive TWAS and TWAS-adjacent analyses, using an approach termed sTF-TWAS. Through these analyses, the authors nominate 222 genes as associated with CRC. The authors also perform functional validation for three genes (TXN, IRS2, and RHPN2) in CRC cell lines, reporting effects on cell viability and formation. Overall, a strength of the manuscript is the creativity involved in some of the analyses; for example, TWAS approaches using alternative polyadenylation are somewhat unusual, although interesting. The manuscript reports findings of potential value to the field.\nThe manuscript has a few major limitations:\n1. Among the key claims is the list of CRC-associated TFs. To perform this analysis, the authors use a non-standard approach introduced in a previous manuscript by a similar group of authors (Wen et al., 2021). This approach models GWAS chi-squared statistics in a linear mixed model framework, including as a fixed effect TF binding, and a random intercept term for each LD block. This approach is conceptually somewhat similar to LD score regression (Bulik-Sullivan et al., 2015) with some important differences; the approach in this manuscript does not include the LD score of a given SNP but instead uses a random intercept term for each LD block. Although a reasonable model, the results are clearly quite different from LD score regression (Supplementary Figure 1), which is a much more common approach in the field. The text refers to this in contradictory terms, first stating that the results are “consistent” but also that their “method outperform[ed] [LD score regression].” Without a much more careful analysis here, along with detailed simulations comparing their method to LD score regression, I find it quite difficult to accept this claim; there is no real justification provided for why the results are so different. This is important because it substantially changes the list of reported CRC relevant TFs.\n2. It is slightly difficult to parse the use of scRNA-seq data here to “evaluate” some of the novel genes. The authors report 99 out of 222 genes were differentially expressed in CRC relevant pathways, though this is somewhat difficult to interpret, because presumably much of this is association rather than causal relationships. Presumably there are large differences in power as well for highly vs lowly expressed genes.\n3. The third main figure contributes little to the manuscript as it contains many nodes in the graph, with leaf nodes that are very small. Most clear is a set of a small number of central nodes, though one wonders what the graph adds beyond the reported tables.\n4. The discussion contains text referencing the important of using “sTF-TWAS” instead of more common TWAS approaches, yet this reviewer could find no evidence of a direct comparison in this manuscript; only a reference to another manuscript that performed a direct comparison.\nMinor limitations:\n1. Supplementary Figure 6 should be removed or substantially revised. This reviewer can only read a central node “colorectal cancer.”\n2. The discussion contains musings about TWAS (“methods like transTF-TWAS and deep-learning method”) that seem slightly out of place for this manuscript.\n3. The text in the discussion about lookups in a drug-protein interaction database is quite long and unclear what the reader should take from this.", "reviewer_2": "The study conducted integrative analysis of multiple types of omics data, including ChIP-seq, gene expression, GWAS, and knockout experiments, to identify putative transcription factors and genes for colorectal cancer. Overall, this is a well-conducted study and well-written paper. I have the following comments.\n1. The first two subsections of Results, i.e., CRC-associated TFs (page 6, lines 172-189) and genetically predicted expression vs CRC risk (line 191-218), both presented interesting findings, but there was little effort to connect the two of them. Since the TWAS analysis used TF information, i.e., the underlying hypothesis is that cis-SNPs could affect gene expression by disrupting TF binding, are the TWAS significant genes likely to be the target genes of CRC-associated TFs? Are the CRC-associated TFs among the TWAS significant genes?\n2. It would be great to see more validation analysis for the TFs and genes that was identified. Although the knockout experiments serve as validation, they were only applicable to a few selected genes. The authors may consider using independent datasets to validate their findings or conduct gene-set enrichment analysis.\n3. Page 7, lines 240-247 (combining TWAS, AS-WAS, and APA-WAS) It would be great if the authors could elaborate on genes commonly identified by all three analyses. Could this be due to any interconnection between gene expression, alternative splicing and alternative polyadenylation? Or are they mostly independent signals?\n4. Figure 3 and Supplementary Figure 6. Circular presentation of the gene-drug relationship is a bit unnecessary. Especially, in Supplementary Figure 6, the texts and lines were so small that they are difficult to read. A table could represent his relationship, though the table will be large.\n5. Supplementary Figure 1. The x axis label should be “transcription” factors, not “transcriptome”"}, "author_response": "Response: We thank the reviewer for the positive feedback and insightful comments that helped improve this work. In the revised manuscript, we have included new simulation analyses to compare the statistical power of our proposed generalized mixed model (GMM) with stratified LD Score Regression (S-LDSC). Our simulations demonstrate that GMM generally outperforms S-LDSC in detecting significant TFs (Supplementary Fig. 1B), and these findings are consistent with our real data results.\nAs the primary goal of our study is to uncover CRC-associated TFs, we also provide the evidence that 40 of our reported 51 TFs (78.4%) were also detected by S-LDSC at Bonferroni-corrected P < 0.05. Additionally, most of the remaining 11 TFs had nominal significance (P < 0.05) (Supplementary Fig. 1A). These results provide further support the validity of our findings. We have included the new results in the revised Results section and added a detailed description of the simulation framework in the Supplementary Note. We have also uploaded our codes and the detailed information of our input and simulation data to our GitHub repository (https://github.com/XingyiGuo/sTF-TWAS-CRC).\nSimulation framework:\n“To benchmark the statistical power of GMM against S-LDSC, we conducted simulations under a polygenic causality scenario where genotypes influence phenotypes through TF binding sites. Specifically, we assumed that only SNPs located within TF binding sites contribute to disease risk, while SNPs outside these regions have no effect.\nFor each replicate, we randomly selected K causal variants (K = 500, 1,000, 5,000, or 10,000) from a pool of 600,000 SNPs annotated to 84 TFs in CRC cell lines (Supplementary Table 1).\nUnder an additive architecture, phenotypes were simulated as:\n𝐾\n𝑌∗ = ∑𝛽 ×𝑋𝑇𝐹 + 𝜀 (1)\n𝑖 1 𝑖𝐗 𝑖\n𝑗=1\nwhere, 𝑋𝑇𝐹 is the genotype (0,1,2) of j-th variant for subject i using data from 489 European ancestry from1000G project phase 3 [PMID: 26432245]. β was set to 1 and 𝜀 is sampling from a normal distribution 𝑁(0, 1−ℎ2 / ℎ2 * 𝑣𝑎𝑟(𝑋)), with heritability ℎ2= 0.5 and 𝑣𝑎𝑟(𝑋) is the variance from genotype part.\nTo facilitate downstream GWAS, we binarized the continuous trait into a case–control phenotype 𝐶∗ based on the mean value of 𝑌∗ across individuals using equation (2).\n𝑁\n1, 𝑌∗ ≥ 𝑚𝑒𝑑𝑖𝑎𝑛(∑𝑌∗)\n𝑖 𝑖\n𝐶∗ = 𝑖=1 (2)\n𝑖 𝑁\n0, 𝑌∗ < 𝑚𝑒𝑖𝑎𝑛(∑𝑌∗)\n𝑖 𝑖\n{\n𝑖=1\nWe performed GWAS using PLINK 2.0 [PMID: 25722852] and applied both GMM and S-LDSC to the resulting summary statistics to estimate TF–disease associations.\nFor each K, we conducted 50 replicates, each with a different random subset of causal SNPs, ensuring reproducibility while keeping K fixed. Results were summarized as the number of true-positive TFs detected at a nominal threshold (P < 0.05) by each method.” ( page 7, Supplementary Note)\nResult section:\n“In addition, we applied stratified LD Score Regression (LDSC) as an alternative approach and found that 40 of the 51 risk TFs identified by our method were also detected by S-LDSC at a Bonferroni-corrected P < 0.05. Most of the remaining 11 TFs showed nominal significance (P < 0.05) (Supplementary Fig. 1A), providing further support for the validity of our findings. We also conducted simulation analyses to compare the statistical power of our GMM with S-LDSC. The results demonstrate that GMM generally outperforms S-LDSC in detecting significant TFs (Supplementary Fig. 1B), consistent with the observations from our real data analysis (Supplementary Fig. 1A).” (page 6)\n2. It is slightly difficult to parse the use of scRNA-seq data here to “evaluate” some of the novel genes. The authors report 99 out of 222 genes were differentially expressed in CRC relevant pathways, though this is somewhat difficult to interpret, because presumably much of this is association rather than causal relationships. Presumably there are large differences in power as well for highly vs lowly expressed genes.\nResponse: We thank the reviewer for this thoughtful comment. To address concerns regarding differences in power between highly and lowly expressed genes, we have added clarification on how differential expression analysis is performed using DESeq2. We have included the clarification in our Supplementary Note. “DESeq2 accounts for heterogeneous variability (heteroscedasticity) in gene expression, where variance depends on mean values, by modeling this relationship. This approach effectively adjusts for differences in variability between low- and high-expressing genes, leading to more accurate differential expression analysis.” (Supplementary Note, page 6)\nWe have also revised the Discussion and included the limitation of differential expression analyses to reflect causal associations. “We also acknowledge that some genes identified through differential expression analysis across developmental stages may not reflect true causal associations, particularly given the high correlation among many gene expression levels.” (page 10). In addition, we have revised the paragraph subtitle to: “Using scRNA-seq data to explore risk genes potentially contributing to CRC development” (page 7).\n3. The third main figure contributes little to the manuscript as it contains many nodes in the graph, with leaf nodes that are very small. Most clear is a set of a small number of central nodes, though one wonders what the graph adds beyond the reported tables.\nResponse: We thank the reviewer for this comment and apologize for any difficulty in interpreting Figure 3. We have replaced it with a high-resolution version for improved clarity. Figure 3 summarizes a subset of data from Supplementary Table S14, focusing on drugs that are either approved or in phase II/III clinical trials for colorectal cancer. Our goal is to provide a visual overview linking actionable genes to their corresponding drugs, offering readers an intuitive understanding of the therapeutic landscape and the translational potential of our findings.\n4. The discussion contains text referencing the important of using “sTF-TWAS” instead of more common TWAS approaches, yet this reviewer could find no evidence of a direct comparison in this manuscript; only a reference to another manuscript that performed a direct comparison.\nResponse: We previously demonstrated that the sTF-TWAS framework identifies more cancer risk genes than conventional TWAS approaches in studies of breast, lung, and prostate cancer, as well as in simulation analyses [PMID: 36402776]. In the revised manuscript, we have included an additional comparison of our sTF-TWAS approach with conventional TWAS and AS-WAS [PMID: 37632791] in the Discussion section.\n“Compared to traditional TWAS analysis in our previous study16, the sTF-TWAS method improved prediction of gene expression and AS, increasing the number of predicted genes by approximately 1,500 and 2,000, respectively. Moreover, sTF-TWAS identified 94 additional genes beyond those detected by the traditional TWAS and AS-WAS.” (page 9).\nMinor limitations:\n1. Supplementary Figure 6 should be removed or substantially revised. This reviewer can only read a central node “colorectal cancer.”\nResponse: Thank you for pointing this out and we apologize for the difficulty in interpreting Supplementary Figure 6. As the underlying information is already presented in Supplementary Table S19, we have removed original Supplementary Figure 6 from the revised manuscript.\n2. The discussion contains musings about TWAS (“methods like transTF-TWAS and deep-learning method”) that seem slightly out of place for this manuscript.\nResponse: We have removed the sentence.\n3. The text in the discussion about lookups in a drug-protein interaction database is quite long and unclear what the reader should take from this.\nResponse: We have revised the description of the drug–protein interaction database and highlighted several novel risk genes that are druggable with existing compounds. Furthermore, we have incorporated additional biological evidence to facilitate the translation of our genetic findings into precision medicine strategies for CRC prevention and intervention.\n“Among these six novel genes, three of them (DDR1, SIK2, MAPK13) encode kinases involved in, or interacting with, the MAPK signaling pathway and are potential targets for approved multi-kinase inhibitors such as Regorafenib1, as well as other drugs including Axitinib2, Cabozantinib3, Pazopanib4, and Sorafenib5. The metabolic enzyme PKM is targetable by compounds such as Resveratrol, which inhibits tumor metabolism via PKM family proteins6, and Metformin, which may reduce CRC risk by suppressing proliferative signaling and enhancing chemotherapy efficacy7. PARP11, implicated in DNA repair, can be inhibited by approved PARP inhibitors including Olaparib, Niraparib, and Talazoparib8. Furthermore, numerous drugs approved for indications in other cancers or non-cancer diseases were found to target CRC susceptibility genes, highlighting opportunities for drug repurposing in CRC treatment. Overall, our protein–drug mapping expands the catalogue of druggable genes, candidate therapeutics and supports the development of precision medicine strategies for CRC prevention and intervention.” (page 11)\nReviewer #2 (Remarks to the Author):\nThe study conducted integrative analysis of multiple types of omics data, including ChIP-seq, gene expression, GWAS, and knockout experiments, to identify putative transcription factors and genes for colorectal cancer. Overall, this is a well-conducted study and well-written paper. I have the following comments.\n1. The first two subsections of Results, i.e., CRC-associated TFs (page 6, lines 172-189) and genetically predicted expression vs CRC risk (line 191-218), both presented interesting findings, but there was little effort to connect the two of them. Since the TWAS analysis used TF information, i.e., the underlying hypothesis is that cis-SNPs could affect gene expression by disrupting TF binding, are the TWAS significant genes likely to be the target genes of CRC-associated TFs? Are the CRC-associated TFs among the TWAS significant genes?\nResponse: We thank the reviewer for the positive feedback and insightful comments. Under the sTF-TWAS framework, TWAS-significant genes are linked to target genes of CRC-associated TFs, as the prediction model is constructed using only cis-SNPs bound by these TFs. By examining TF occupancy information for the genetic variants included in the prediction model, we found that 117 of the 222 putative susceptibility genes (52.7%) are commonly regulated by the top five TFs (RAD21, CBX3, MAX, SP1, and POLR2A). In addition, we conducted a new TF enrichment analysis on these 222 genes using publicly available TF ChIP-seq data from the Enrichr [PMID: 23586463, PMID: 27141961]. This analysis identified four CRC-associated TFs (VDR, EP300, MYC, and SOX2) significantly enriched at FDR < 0.05, and an additional set of 17 TFs enriched at nominal P < 0.05 (Supplementary Table 18). Of note, among the 222 risk genes, several TFs themselves, such as MYC, PROX1, SOX9, KLF5, and RAD21, were identified as significant genes. It is worth noting that most CRC-associated TFs were not detected by TWAS, likely because TWAS captures genetic heritability primarily from cis-variants. However, our generalized mixed model (GMM) analysis using genome-wide variants provided evidence of their susceptibility. Taken together, these findings underscore the role of CRC-associated TFs and TF-based transcriptional networks in CRC genetic susceptibility. We have added this information to the Discussion section to highlight these insights.\n“By examining TF occupancy information for the genetic variants included in the prediction model, we found that 117 of the 222 putative susceptibility genes (52.7%) are commonly regulated by the top five TFs (RAD21, CBX3, MAX, SP1, and POLR2A). In addition, we conducted a new TF enrichment analysis on these 222 genes using publicly available TF ChIP-seq data from the Enrichr9, 10. This analysis identified four CRC-associated TFs (VDR, EP300, MYC, and SOX2) significantly enriched at FDR < 0.05, and an additional set of 17 TFs (Supplementary Table 18) enriched at nominal P < 0.05. Of note, among the 222 risk genes, several TFs themselves, such as MYC, PROX1, SOX9, KLF5, and RAD21, were identified as significant genes. It is worth noting that most CRC-associated TFs were not detected by TWAS, likely because TWAS captures genetic heritability primarily from cis-variants. However, our GMM analysis using genome-wide variants provided evidence of their susceptibility. Taken together, these findings underscore the role of CRC-associated TFs and TF-based transcriptional networks in CRC genetic susceptibility.” (page 9)\n2. It would be great to see more validation analysis for the TFs and genes that was identified. Although the knockout experiments serve as validation, they were only applicable to a few selected genes. The authors may consider using independent datasets to validate their findings or conduct gene-set enrichment analysis.\nResponse: We thank the reviewer for the insightful comments. As noted above, we conducted TF enrichment analysis based on the 222 risk genes using Enrichr [PMID: 23586463; PMID: 27141961]. This analysis revealed four CRC-associated TFs (VDR, EP300, MYC and SOX2) with significant regulatory enrichment at FDR < 0.05, highlighting their potential roles in the regulation of these risk genes. In addition, 17 TFs showed enrichment at a nominal P < 0.05.\nWe acknowledge that future fine-mapping efforts are required to experimentally validate risk variant–TF–gene regulatory relationships using approaches such as luciferase reporter assays and/or enhancer knockout experiments. In our prediction models, we identified 144 lead variants, representing the strongest associations with CRC risk for each model, at P < 5 × 10⁻⁸. These variants are predicted to influence the binding affinity of CRC-associated TFs and are linked to 162 risk genes (Supplementary Table 19). These data provide valuable insights for further fine-mapping studies through in vitro/in vivo functional assays to pinpoint causal variants and genes for CRC. We have incorporated these descriptions and outlined the prospective directions in the Discussion section.\n“Moving forward, in our prediction models, we identified 144 lead variants—representing the strongest associations with CRC risk for each model, at P < 5 × 10⁻⁸. These variants are predicted to influence the binding affinity of CRC-associated TFs and are linked to 162 risk genes (Supplementary Table 19). Taken together, these findings offer valuable insights for prioritizing functional studies aimed at fine-mapping causal variants and genes for CRC through in vitro and in vivo assays.” (page 9)\n3. Page 7, lines 240-247 (combining TWAS, AS-WAS, and APA-WAS) It would be great if the authors could elaborate on genes commonly identified by all three analyses. Could this be due to any interconnection between gene expression, alternative splicing and alternative polyadenylation? Or are they mostly independent signals?\nResponse: We thank the reviewer for the insightful comments. Gene expression reflects overall transcript abundance, alternative splicing reflects exon usage, and APA reflects 3′ end site usage; therefore, these molecular traits are measured differently in RNA-seq data analysis. Among the 222 identified risk genes, only four genes (COLCA2, METRNL, ARPC5, COX14) were detected by TWAS, AS-WAS, and APA-WAS. This observation aligns with prior research showing that different molecular traits within the same gene can represent largely independent signals (PMID: 37601976; PMID: 38759092). For the commonly detected genes, we evaluated the genetic variants included in prediction model of each molecular. We found that very few of shared variants were observed across them. For example, for gene METRNL, the gene expression, APA, and AS models included 32, 32, and 13 variants, respectively. Three variants were shared by all models; eight overlapped between expression and APA, seven between expression and AS, and three between AS and APA. This indicates genetic heritability may independently contribute to each molecular trait, while a small proportion of risk variants may influence expression, splicing, and APA through interconnected mechanisms under multi-layered regulatory control. We have included additional discussion in the Discussion section.\n“Among the 222 identified risk genes, only four genes (COLCA2, METRNL, ARPC5, COX14) were detected by TWAS, AS-WAS, and APA-WAS. This observation aligns with prior research showing that different molecular traits within the same gene can represent largely independent signals11, 12. For these four commonly detected genes, we found that very few shared variants were observed among them, indicating genetic heritability may independently contribute to each molecular trait, while a small proportion of shared risk variants may influence expression, splicing, and APA through interconnected mechanisms under multi-layered regulatory control.” (page 9).\n4. Figure 3 and Supplementary Figure 6. Circular presentation of the gene-drug relationship is a bit unnecessary. Especially, in Supplementary Figure 6, the texts and lines were so small that they are difficult to read. A table could represent his relationship, though the table will be large.\nResponse: We thank the reviewer for this comment and apologize for the difficulty in interpreting Figure 3. We have provided a high-resolution version for better clarity. Figure 3 shows a subset of Supplementary Table S14, highlighting drugs approved or in phase II/III trials for colorectal cancer, to give readers a visual overview linking actionable genes to their target drugs. For Supplementary Figure 6, we have removed it contains substantial information already presented in Supplementary Table S19.\n5. Supplementary Figure 1. The x axis label should be “transcription” factors, not “transcriptome”\nResponse: Fixed."}, {"round": 2, "reviewer_comments": {"reviewer_1": "In this revised text, the authors substantially revised the discussion, and conducted a simulation to compare their \"GMM\" approach to a sLDSC. Although these revisions are helpful, the text would still benefit from consideration of a few major points.\n1. Regarding the simulation the authors added to compare sLDSC to their approach; overall this simulation is unconvincing for two reasons. This reviewer interprets equation (1) as indicating the variance explained by genotype is 50% on the liability scale for the case-control. This frankly is very large effect size. The variance of the noise term is also unclear, given that is specified as (0.5 / 0.5) * Var(X) where X is a vector of binomial random variables with distinct variances. This reviewer checked the provided code, and found that the code used to perform the simulation (https://github.com/XingyiGuo/sTF-TWAS-CRC/blob/e905b7c44c508265fad081c7676248902ee65184/src/simultation/1.simulate_Yvalues_basedonGeno.R#L73) does not in fact correspond to the equation presented in the rebuttal; in fact the code defines the the variance of the error term as the variance of the genetic liability (genotype * beta) for each individual, though that happens to be equivalent in this case to what the authors intend since beta is fixed to 1. Secondly, and more importantly, the authors do not assess false positives in their simulation, only statistical power. Finally - the authors should note the standard nomenclature is to refer to such models as \"generalized linear mixed models\" or GLMMs, rather than GMMs.\n2. The authors appears to have replaced Figure 3, yet the figure still has little value to the reader as the leaf nodes are still far too small to be read. The reader cannot be expected to zoom in 250% to read the labels. The text in Figure 1D is also much too small.\n3. The discussion now contains claims that sTF-TWAS increases power by 94 additional genes over 'traditional TWAS'. Have the authors interrogated the false positive rate of their approach ? In the absence of careful false positive estimation, it is difficult to accept the claims of power increase. It is also worth clarifying explicitly what the authors refer to as \"traditional TWAS\"; the authors appear to use an in house elastic pipeline, rather than using published weights from a published study (e.g., PredictDB).\n4. The revised discussion still contains very specific speculative callouts to multi-kinase inhibitors along with Metformin + others, which strikes this reviewer as overly speculative.\n(Remarks on code availability)\nI reviewed a subset of the code in this repo: https://github.com/XingyiGuo/sTF-TWAS-CRC, and did not attempt to execute the code myself. I encourage the authors to use renv to document their dependencies.", "reviewer_2": "The authors have adequately addressed my comments."}, "author_response": "Reponses: We appreciate reviewer’s insightful comments. We acknowledge 50% heritability is high and we conducted new simulations using additional three ranges of heritability (5%, 10%, 25%). As the reviewer noted, in our simulation based on equation (1), we sampled 𝜀 from a normal distribution 𝑁(0, 1−ℎ2 / ℎ2 * 𝑣𝑎𝑟(∑𝐾 𝛽 ×𝑋𝑇𝐹)), which corresponds to defining the variance of the genetic liability.\n𝐾\n𝑌∗ = ∑𝛽 ×𝑋𝑇𝐹 + 𝜀 (1)\n𝑖 1 𝑖𝐗 𝑖\n𝑗=1\nWe carefully reviewed our code, and the implementation of Equation (1) is carried out in lines 82–84 of the R script titled “1.simulate_Yvalues.R”. We apologize for the nomenclature, and we corrected all “generalized mixed models” to “generalized linear mixed models (GLMM)”. In addition, we also conducted analyses to assess the type I error rate for GLMM under null hypothesis that genotype is not associated with phenotype.\nSimulation studies:\n“To evaluate the type I error rate of GLMM, we conducted simulation analyses under the null hypothesis of no association between genotype and phenotype. Specifically, we randomly assigned binary phenotype values (0 or 1) to 489 individuals of European ancestry and performed logistic regression to generate GWAS summary statistics under the null. We then applied GLMM and S-LDSC using different numbers of TF binding variants (500, 1,000, 5,000, and 10,000) and repeated the entire procedure 1,000 times. Finally, we compared the Type I error rates of GLMM and S-LDSC using quantile-quantile (QQ) plots of p-values across the 1,000 replicates for each binding variant set size.” (page 8, Supplementary Note)\nResults:\n“Our results demonstrate that GLMM consistently outperformed S-LDSC in detecting causal TFs across varying levels of heritability (Supplementary Fig. 1B, Supplementary Table 21). Quantile-quantile (Q-Q) plots demonstrate that Type I error rates are under control for both approaches, with overall comparable performance (Supplementary Fig. 1C-F).These simulation results are consistent with our observations in real data, supporting GLMM as a robust approach for identifying risk TFs in colorectal cancer.”\nSupplementary Figure 1. B) Bar plots showing the number of transcription factors (TFs) detected by the GLMM and S-LDSC models across four sets of TF-binding variants (500, 1,000, 5,000, and 10,000) and four heritability levels (0.05, 0.1, 0.25, and 0.5)\nSupplementary Figure 1. C) Quantile–quantile (QQ) plots of p-values under the null hypothesis for the GLMM and S-LDSC models. Orange dots represent GLMM results, green dots represent S-LDSC results, and the gray dashed line indicates the expected distribution under the uniform null.\n2. The authors appears to have replaced Figure 3, yet the figure still has little value to the reader as the leaf nodes are still far too small to be read. The reader cannot be expected to zoom in 250% to read the labels. The text in Figure 1D is also much too small.\nReponses: We apologize for the lack of clarity in Figure 3. On our end, the figure appeared clearly in the PDF or PNG format without requiring zoom (see Figure below). To improve readability, we have increased the font sizes in both revised Figure 3 and Figure 1D\nSupplementary Figure 3. A summary of novel genes identified through TWAS, sp-TWAS, and APA-WAS. Numbers in Venn Diagram plots refer to the total number of genes; novel genes name are displayed in the boxes.\n3. The discussion now contains claims that sTF-TWAS increases power by 94 additional genes over 'traditional TWAS'. Have the authors interrogated the false positive rate of their approach ? In the absence of careful false positive estimation, it is difficult to accept the claims of power increase. It is also worth clarifying explicitly what the authors refer to as \"traditional TWAS\"; the authors appear to use an in house elastic pipeline, rather than using published weights from a published study (e.g., PredictDB).\nReponses: “Traditional TWAS” refers to transcriptome-wide association studies through S-prediXcan in our previous work (PMID: 37632791). Of note, our prior methodological study demonstrated that the sTF-TWAS approach maintains proper type-I error control under the null hypothesis in the extensive simulation analysis (PMID: 36402776). The elastic net model used here was implemented in-house to streamline computation and has been thoroughly validated for robustness. Notably, this model has been used in two previous publications in Nature Communications (PMID: 34518541, PMID: 36402776), further validating its utility and reliability.\n4. The revised discussion still contains very specific speculative callouts to multi-kinase inhibitors along with Metformin + others, which strikes this reviewer as overly speculative.\nReponses: We thank the reviewer for the thoughtful comments. We agree that speculative interpretations should be presented cautiously, especially when discussing potential therapeutic implications. In the revised discussion, we have carefully revised the relevant text to directly connections between our risk genes and drugs. We also emphasized these are exploratory findings rather than definitive therapeutic recommendations.\nDiscussion:\n“Among these six novel genes, PKM2 from PKM family, a member of the PKM family, has been reported to mediate VEGF-induced angiogenesis in human endothelial cells. Resveratrol inhibits this process by suppressing aerobic glycolysis, potentially through modulation of PKM2 nuclear translocation. PARP11, implicated in DNA repair, can be inhibited by approved PARP inhibitors including Olaparib, Niraparib, and Talazoparib. However, these results are observations are exploratory and require further validation.” (page 11)\nReviewer #1 (Remarks on code availability):\nI reviewed a subset of the code in this repo: https://github.com/XingyiGuo/sTF-TWAS-CRC, and did not attempt to execute the code myself. I encourage the authors to use renv to document their dependencies.\nReponses: We appreciate reviewer’s suggestion. We documented our R environment using the R package named “renv”. Moreover, we uploaded all code and example data. Genotype data cannot be redistributed due to policy restrictions but can be freely downloaded by reviewers from the official 1000 Genomes website (https://www.internationalgenome.org/data-portal/data-collection/30x-grch38) or the processed folder named “Variant_effects/1000G_EUR_Phase3_plink_hg38”(https://console.cloud.google.com/storage/browser/broad-alkesgroup-public-requester-pays) for verification."}, {"round": 3, "reviewer_comments": {"reviewer_1": "I have no further comments for the authors. They have addressed my previous comments."}, "author_response": ""}]

Decision letter