Codebook Transcription Factor Plasmids, DNA-Binding Assays and Motif Analysis

Plasmids and Inserts

Sequences and supporting information are provided in Supplementary Tables 2–4. Briefly, Codebook transcription factors (TFs) and their DNA-binding domains (DBDs) were selected from data reported in a previous study¹ and available through the Human TFs database.

Inserts ending in “-FL” represent the full-length open reading frame (ORF) of a representative protein isoform. Inserts ending in “-DBD” contain all predicted DBDs, flanked by either 50 amino acids or the amino- or carboxy-terminal end of the protein. Inserts designated “-DBD1,” “-DBD2” or “-DBD3” contain subsets of the protein’s DBDs. These constructs were manually designed, particularly for large C2H2 zinc-finger arrays.

All inserts were produced as recoded synthetic ORFs by BioBasic. Each ORF was flanked by AscI and SbfI restriction sites and subcloned into up to three expression plasmids:

1. pTH13195: a tetracycline-inducible, N-terminal eGFP-tagged mammalian expression vector containing FLiP-in recombinase sites⁸;
2. pTH6838: a T7 promoter-driven, N-terminal GST-tagged bacterial expression vector⁵³; and
3. pTH16500 (pF3A-ResEnz-egfp): an SP6 promoter-driven, N-terminal eGFP-tagged bacterial expression vector modified from pF3A–eGFP⁷ to include the two restriction sites downstream of eGFP.

Protein Production

Each experiment used protein produced with one of the following expression systems:

1. FLiP-in HEK293 cells (Thermo Fisher Scientific, R78007), induced with doxycycline for 24 h. This system was used for inserts cloned into pTH13195;
2. PURExpress T7 recombinant in vitro transcription–translation system (NEB, E6800L), used for inserts in pTH6838; or
3. SP6-driven wheat germ extract-based in vitro transcription–translation system (Promega, L3260), used for inserts in pTH16500.

DNA-Binding Assays

We used previously established protocols for ChIP–seq⁸, protein-binding microarrays (PBMs)⁵⁰ and SMiLE-seq⁷. Detailed procedures for GHT-SELEX, HT-SELEX, ChIP–seq and SMiLE-seq data generation and preliminary analysis are described in the accompanying publications^10–13.

For PBM experiments, proteins were tested on two universal arrays, HK and ME, which contain different probe sequences⁵⁴. Most C2H2 zinc-finger proteins were excluded from PBM analysis because of the low success rate of this assay, presumably because these proteins recognize long DNA-binding sites.

Unless otherwise stated, each full-length protein construct was analysed twice by ChIP–seq¹¹. Each construct was analysed once by HT-SELEX and GHT-SELEX, including both full-length proteins and DBD constructs. Several constructs were tested multiple times to assess the effect of experimental variables during development of the GHT-SELEX method¹⁰.

SMiLE-seq assays were performed for 278 TFs, representing 388 constructs. Controls were not included because most selected TFs already had published SMiLE-seq data. This resulted in 299 TFs with motifs derived from SMiLE-seq. A subset of Codebook proteins, primarily those with unknown DBDs, was excluded because the proteins were unsuccessful in all other assays. To evaluate reproducibility, 82 randomly selected constructs were tested multiple times across independent SMiLE-seq experiments.

An anti-eGFP antibody (Ab290, Abcam) was used as the primary antibody in ChIP–seq, GHT-SELEX and western blot experiments, with assay-specific quantities as described previously^10,11.

For ChIP, 2 µl of undiluted polyclonal antiserum, equivalent to 10 µg total IgG, was immobilized on 60 µl of protein G magnetic bead suspension (Dynabead 10004D, Thermo Fisher Scientific). For HT-SELEX and GHT-SELEX, antibody–bead master mixes were prepared by immobilizing 6 µl of antiserum on 100 µl of protein G Sepharose bead slurry (Cytiva, 28-9670-70). One-microlitre aliquots were used for each selection, corresponding to approximately 0.06 µl antiserum, or 300 ng total IgG, per reaction. For western blotting, membranes were incubated in 15 ml of antibody solution diluted 1:5,000, corresponding to approximately 15 µg total IgG per membrane.

Data Processing and Motif Derivation

Motif derivation and evaluation are described in detail in the accompanying study¹². Following initial preprocessing, we identified likely bound, or “true positive,” sequences for each experiment. Of 4,873 experiments, 721 were excluded because they contained too few likely bound sequences or had other technical problems, as documented in Supplementary Table 4.

We applied multiple motif-discovery tools, listed in Supplementary Table 15, to training subsets from each experiment. The resulting motifs were evaluated using test sequences from the same experiment and independent datasets for the same TF, including test sets from other experiments involving that TF.

All experiments and motifs were assessed using binary classification. Motif performance was measured using several metrics, including the area under the receiver operating characteristic curve (AUROC) and the area under the precision–recall curve. The complete motif collection, provided as position-weight matrices (PWMs), is available through Zenodo⁵⁵. An interactive motif browser is available at MEX.

Systematic Filtering of Artefactual Motifs

During dataset curation, we specifically addressed recurrent artefactual motifs caused by systematic experimental noise or by characteristics of individual motif-discovery tools. For example, the ACGACG motif was frequently enriched in HT-SELEX experiments. This motif was considered an artefact because it matches the constant flanking sequence used in the assay.

Experiments using cell lysates also occasionally enriched motifs corresponding to abundant endogenous HEK293 proteins, including NFI, YY1 and ETS-family TFs. To reduce the impact of such artefacts, we:

1. manually assembled a catalogue of recurrent artefact motifs (Supplementary Table 16); and
2. scanned the complete motif collection with MACRO-APE⁵⁶, removing motifs highly similar to entries in the artefact catalogue.

We also removed motifs matching constant, non-variable DNA regions used in HT-SELEX, GHT-SELEX and SMiLE-seq experiments. ETS-related motifs were retained for ETS-family positive controls, including ELF3, FLI1 and GABPA. After expert curation, we confirmed that enriched k-mers in HT-SELEX experiments were not associated with potential artefacts from individual expression systems¹⁰.

Evaluation of Motif-Discovery Success Rates

To compare motif-discovery tools across experimental platforms, we began with successful experiments and TFs. For each platform X, such as ChIP–seq, and motif-discovery tool Y, such as Autoseed, we counted the number of experiments or TFs that produced a motif highly similar to the manually curated reference motif.

For each experiment or TF, we evaluated the complete set of candidate motifs generated by the relevant platform–tool combination. A one-pass MACRO-APE (v.3.06) ScanCollection analysis was used to determine whether any motif passed the similarity threshold, with the following parameters: -c 0.05 --rough-discretization 10 --precise 100500⁵⁶. Before comparison, motifs were converted to log-odds PWMs as described previously¹².

Expert Evaluation of Individual Experiments

We developed an expert-curation workflow to determine whether individual experiments were sufficiently reliable for inclusion in downstream analyses. Initially, a committee of annotators voted on the success or failure of each experiment. Every experiment was reviewed by at least three annotators.

All disagreements, involving approximately 300 experiments, were resolved jointly by A.J., I.V.K. and T.R.H. This subcommittee also reviewed every experiment classified as successful.

Annotators used an early version of the MEX portal (https://mex.autosome.org) to review PWMs scored against all experiments. They evaluated whether experiments produced consistent motifs across replicates or whether the motifs scored highly across datasets.

The review also considered:

• whether motifs were consistent with the expected specificity of the TF family, such as the E-box-like CAnCTG motif produced by BHLHA9;
• whether motifs were similar between closely related paralogues, such as ZXDA, ZXDB and ZXDC;
• the number and quality of peaks identified by ChIP–seq or GHT-SELEX;
• whether peaks were reproduced across independent assays, such as ChIP–seq and GHT-SELEX;
• similarity between Codebook PWMs and publicly available PWMs; and
• the enrichment of known or suspected contaminant motifs.

Post-Evaluation Peak Processing

After successful experiments were identified, we re-derived peak sets for ChIP–seq and GHT-SELEX to generate one peak set per TF, following procedures described in the accompanying publications^10,11.

For ChIP–seq, peak calling was repeated with MACS2 (v.2.2.9.1)⁵⁷ using experiment-specific background sets generated with a previously described method⁸. Peak sets from replicates of the same TF were then merged with BEDTools (v.2.30.0) merge⁵⁸, as described in the accompanying manuscript¹¹.

GHT-SELEX peaks were generated using MAGIX, a new method that calculates read enrichment at each selection cycle and treats separate experiments as independent statistical samples. This approach produces one enrichment coefficient for each peak¹⁰.

Expert Motif Curation

To identify one representative PWM for each TF, we assembled the highest-scoring candidate PWMs, performed additional testing with reprocessed peak data and manually reviewed the results.

For each TF, we began with the union of three sets of 20 PWMs:

• the 20 PWMs with the highest AUROC on successful ChIP–seq experiments;
• the 20 PWMs with the highest AUROC on successful GHT-SELEX experiments; and
• the 20 PWMs with the highest AUROC on successful HT-SELEX experiments.

PWMs were selected independently of the dataset from which they originated. We reassessed the candidates against ChIP–seq and GHT-SELEX data using two complementary approaches.

First, we recalculated AUROC values for the candidate PWMs using merged, thresholded ChIP–seq peak sets (P < 10⁻¹⁰)¹¹. Peaks were scored with AffiMX (v.1)²⁵. Negative sets were generated with BEDTools (v.2.30.0) shuffle⁵⁸ using the -noOverlapping option. These sets contained the same number of regions and the same peak-width distribution as the corresponding ChIP–seq sets.

The same procedure was used for GHT-SELEX, using thresholded peak sets defined with a “kneedle”⁵⁹ specificity value of 30 in the sorted enrichment values¹¹.

Second, we calculated the Jaccard index to quantify the overlap between PWM matches and ChIP–seq or GHT-SELEX peaks. PWM matches were identified with MOODS (v.1.9.4)⁵¹ using -p 0.001. For each assay, overlap was optimized by applying different peak thresholds and selecting the cutoff that produced the highest Jaccard index¹⁰.

A committee comprising A.J., T.R.H., K.U.L., A.F., R.R., M.A. and I.Y. then selected one representative PWM per TF. Selection prioritized high performance across all scores, consistency with the expected DBD class and, where appropriate, recognition-code-predicted motifs. PWMs also needed to have sufficient information content (IC).

Motif Similarity Analysis and Clustering

We used two complementary approaches to estimate the number of distinct motifs represented by Codebook TFs and to identify previously unknown motifs added to the human TF repertoire.

In the first analysis, we compared 1,582 PWMs representing Codebook TFs and TFs with previously characterized specificities. The Codebook set included 177 representative PWMs from 177 TFs. The comparison set contained 1,405 PWMs from 1,211 TFs. For previously characterized TFs, we used the “best” PWMs reported in an earlier study¹, except for MTF2 and PHF1, which lacked assigned PWMs. Their PWMs were retrieved from CisBP⁵³.

PWM similarity was measured using the correlation between pairwise affinities to 150,000 random sequences of length 100, calculated with MoSBAT (v.1)²⁵. The 177 Codebook TFs were clustered using Pearson correlation coefficient (PCC) distance and average linkage. The optimal number of clusters was 129, based on the highest silhouette value⁵² calculated with the silhouette function in the R cluster package (v.2.1.8.1), using R (v.4.3.2). This corresponded to a PWM distance threshold of 0.76.

We then clustered all 1,582 motifs using the same distance threshold. This produced 613 clusters, including 92 clusters containing only Codebook TFs. The motifs and cluster assignments are provided in Supplementary Table 10.

Independently, we generated a non-redundant motif set for MARA analysis using a previously described procedure⁶⁰, expanded to include Codebook motifs. HOCOMOCO (v.12)²⁹ motifs were merged with evaluated Codebook motifs. Where possible, we selected motifs generated with ChIPMunk and highly ranked during benchmarking¹² to maintain consistency with HOCOMOCO (v.12).

Motif similarity was estimated with MACRO-APE (v.3.0.6)⁵⁶ using a motif P-value cutoff of 0.0005 and the default matrix-discretization parameter (-d 1). The parameter was increased to 10 for improved precision only for motif pairs with a Jaccard similarity above 0.01 at -d 1.

Using the pairwise similarity matrix, we performed agglomerative clustering with average linkage in sklearn (v.1.8.0). The number of clusters was selected by maximizing the silhouette score. Low-quality HOCOMOCO “D” motifs were excluded. One representative motif was selected for each cluster based on the highest average similarity to all other motifs in that cluster. The resulting non-redundant motif set is available from Zenodo⁶¹, and cluster annotations are provided in Supplementary Table 11.

Motif Degeneracy Analysis

To determine whether low information content is an intrinsic characteristic of certain DNA-binding motifs, we adjusted PWM information content and evaluated binding-site prediction accuracy using ChIP–seq and GHT-SELEX data.

Information content was adjusted independently at each base pair by iteratively scaling base probabilities until the PWM reached an average information content of 1 bit per base pair. The logo_rescale.pl script is available from GitLab at https://gitlab.sib.swiss/EPD/pwmscan.

Comparison with External Peak Sets and PWMs

For comparison with external datasets, we downloaded ChIP–seq peak sets for all Codebook TFs from GTRD⁶² and ENCODE (4 December 2023)⁶³. Datasets were grouped into four cell-type categories: HEK293/HEK293T, HepG2, K562 and other cells.

GTRD was prioritized because it processed most ENCODE consortium experiments and included additional non-ENCODE studies. When multiple experiments were available for a TF within a cell-type category, we selected the experiment with the greatest number of peaks. If multiple peak callers were available, we used the following preference order: MACS, GEM, SISSRS, PICS and PEAKZILLA. Dataset identifiers and metadata are listed in Supplementary Table 7.

External peak sets were used as downloaded, except for GEM-derived sets, which contain summit-only peaks with a width of one base. These peaks were expanded by 250 bases in both directions.

For Codebook analyses, we used merged and thresholded Codebook ChIP–seq peak sets as described in the expert motif-curation section. Negative sets were generated with BEDTools (v.2.30.0) shuffle⁵⁸ using -noOverlapping, preserving the number and width distribution of the original peaks.

PWMs for Codebook TFs were downloaded from JASPAR (2024)⁶⁴, HOCOMOCO (v.12)²⁹ and Factorbook⁶⁵ on 15 December 2023. These data are listed in Supplementary Table 13. Codebook and external peak sets, together with corresponding negative sets, were scanned using representative expert-curated Codebook PWMs and AffiMX (v.1)²⁵. AUROC values were then calculated.

For 20 TFs with successful Codebook ChIP–seq experiments, Codebook PWMs, external ChIP–seq datasets and external PWMs, we compared PWM performance across datasets. For each TF, we selected the external PWM that produced the highest AUROC on the corresponding external peak set. These best-performing PWMs were then used to scan the Codebook peak sets and calculate AUROC values.

Curation of External Motifs for Previously Uncharacterized TFs

We examined proteins previously annotated as putative TFs¹ that lacked a credible motif at the time of the original annotation. All available motifs for these TFs were downloaded from JASPAR (2024)⁶⁴, HOCOMOCO (v.12)²⁹ and Factorbook⁶⁵, producing 484 PWMs.

Each PWM was manually evaluated according to the following criteria:

1. whether similar motifs were observed for the protein across multiple independent datasets;
2. whether the motif was consistent with the structural class of the TF; and
3. whether the motif reflected the protein’s intrinsic DNA-binding specificity rather than the target site of another TF or a likely technical artefact.

The curated motifs are listed in Supplementary Table 13, with additional information in Supplementary Table 10. The PWMs are available in Supplementary Data 1. Examples of artefactual and correctly curated external motifs are shown in Extended Data Fig. 10.

Identification and Scoring of Promoter Sequences for MARA

MARA requires promoter activity data across samples and motif scores across promoters. Promoter activity data were downloaded from the FANTOM5 resource at https://fantom.gsc.riken.jp/5/ using hg38_fair+new_CAGE_peaks_phase1and2_tpm_ann.osc.txt.

Expression values were log₂-transformed using a pseudocount of 0.05. We retained promoters associated with genes encoded in the nuclear genome and excluded time courses, perturbation experiments and human total RNA samples. The final dataset contained 209,374 promoters and 1,020 samples, including replicates. These represented 583 unique samples: 142 tissues, 187 primary cell types and 254 cell lines. Metadata are available in Supplementary Table 12 and Zenodo⁶¹.

Motif scanning was performed on regions from FANTOM5 hg38_fair+new_CAGE_peaks_phase1and2.bed, extending from 250 bp upstream to 10 bp downstream of the representative transcription start site (TSS).

SPRY-SARUS (v.2.2.3; https://github.com/autosome-ru/sarus) was used to calculate sum-occupancy scores⁶⁶ for each representative motif in the motif clusters described above. The final analysis included 632 motif clusters: 130 containing only Codebook motifs, 471 containing known motifs and 31 mixed clusters. Only TFs jointly expressed at levels greater than zero in at least one FANTOM5 sample were included.

We used MARADONER (v.0.13)⁶⁷, a Python command-line tool available through PyPI. Analyses were performed with maradoner create, maradoner fit and maradoner export using default parameters (https://github.com/autosome-ru/MARADONER).

Conceptually, MARADONER extends the original MARA³⁶ and isMARA⁶⁸ frameworks. The model assumes that promoter activity in each sample is a linear function of sample-specific motif activities, where each motif represents TFs with shared DNA-binding specificity. We used the following matrix-variate linear mixed model:

$$Y={{{\boldsymbol{\mu }}}_{{p}}{{\bf{1}}}_{{s}}}^{{\rm{T}}}+{{\bf{1}}}_{{p}}{{\boldsymbol{\mu }}}_{{s}}+B\,U+E,\,E \sim {\rm{M}}{\rm{N}}(0,{I}_{p},D),\,U \sim {\rm{M}}{\rm{N}}({{{\boldsymbol{\mu }}}_{{\rm{m}}}{{\bf{1}}}_{{s}}}^{{\rm{T}}},\varSigma ,{G}),$$

Here, Y is a p × s matrix of log-scale promoter activity, where p is the number of promoters and s is the total number of samples. The vectors μ_p and μ_s represent promoter-specific and sample-specific means, respectively, and 1_n is a vector of ones of length n.

B is a p × m matrix of promoter-level motif scores, and U is an m × s random matrix of motif activities. E is a p × s noise matrix. MN denotes a matrix-variate normal distribution, I_p is the p × p identity matrix, D is a diagonal matrix of noise variances, Σ is an m × m matrix of motif variances, G is a diagonal sample-scaling matrix and μ_m is a vector of motif-specific mean activities.

Unlike classical MARA, modelling μ_m enables explicit separation of activator and repressor activity. MARADONER estimates parameters using a four-stage restricted maximum-likelihood procedure. First, parameters in E are isolated through a transformation orthogonal to the 1_p, 1_s and B matrices, enabling estimation of D. Second, a transformation orthogonal only to the 1_p and 1_s vectors is used to estimate Σ and G. Third, a transformation orthogonal to B is used to estimate the total mean effect of the μ_p and μ_s terms. Finally, μ_m is estimated using the remaining known parameters.

After separating U from \(B{{{\boldsymbol{\mu }}}_{{\rm{m}}}{{\bf{1}}}_{{s}}}^{{\rm{T}}}\), the deviation-from-mean matrix \(\hat{U}=U-B{{{\boldsymbol{\mu }}}_{{\rm{m}}}{{\bf{1}}}_{{s}}}^{{\rm{T}}}\) represents sample-specific motif activity. \(\hat{U}\) is estimated using a maximum a posteriori procedure. MARADONER reports both raw maximum a posteriori estimates and standardized motif activities, calculated by dividing each activity by the square root of its posterior variance.

The availability of maximum-likelihood estimates for motif-specific variances enables an ANOVA-like test based on the asymptotic properties of maximum-likelihood estimates. MARADONER performs a Wald test using standard errors derived from the square roots of the diagonal elements of the asymptotic covariance matrix corresponding to Σ.

MARADONER assumes that gene-expression or promoter-activity data are supplied on a log scale. For the motif-score matrix B, each column is normalized by default using the negative logarithm of its empirical survival function.

TOP and CTOP Peak-Set Analyses

To identify TOP sites, we first determined thresholds for ChIP–seq peaks, GHT-SELEX peaks and PWM-derived peaks that maximized the three-way Jaccard metric. Thresholds were calculated independently for each TF.

PWM matches were generated with MOODS (v.1.9.4)⁵¹ using a P-value cutoff of 0.001. Neighbouring matches separated by less than 200 bp were merged and rescored using sum-occupancy scores. TOPs were defined as peaks that exceeded the relevant thresholds in all three datasets and overlapped across all three peak sets.

To identify CTOP sites, we extracted phyloP scores from the Zoonomia consortium⁶⁹ for each base within TOP sites and 100 flanking bases. Sites overlapping the ENCODE Blacklist⁷⁰ or protein-coding sequences were removed because codons can cause biased phyloP scores.

We used three statistical tests to evaluate phyloP enrichment across PWM matches:

• two tests of the association between PWM information content and phyloP scores at each PWM position, using either PCC or a likelihood-ratio test; and
• a Wilcoxon test for increased phyloP scores across the PWM match.

Additional methodological details are provided in the accompanying manuscripts^10,11. The custom R (v.4) script is available at GitHub.

Intersection of TOPs, CTOPs and Genomic Features

All CTOPs were first clustered with BEDTools (v.2.30.0) merge⁵⁸ using a maximum distance of 100 bp. The resulting clusters were intersected with the following genomic feature sets:

• basic canonical protein-coding promoters from GENCODE (v.44)⁷¹, defined as 1,000 bp upstream and 500 bp downstream of the canonical TSS;
• the unmasked CpG island track;
• PhastCons conserved elements from the Multiz 470 Mammalian alignment;
• the RepeatMasker track from UCSC⁷²; and
• ChromHMM HEK293 enhancers¹¹.

Promoters were classified as CpG-island or non-CpG-island promoters according to whether the GENCODE basic TSS was located within ±50 bp of a CpG island in the unmasked track.

CTOP clusters were assigned to one genomic-feature category according to the following priority order:

1. CpG-island-associated protein-coding promoter;
2. other CpG islands;
3. non-CpG-island-associated protein-coding promoter;
4. enhancer;
5. clusters containing a CTCF-binding site but not overlapping a CpG island, promoter or enhancer;
6. transposable elements when none of the previous categories applied;
7. non-transposable-element repeats when none of the previous categories applied; and
8. other, for clusters that did not overlap any examined feature.

Analysis of Allele-Specific Binding

GHT-SELEX and ChIP–seq data enabled direct assessment of allele-specific binding (ASB) by measuring allelic differences in read counts at single-nucleotide variants (SNVs). These datasets were not originally generated for ASB analysis, so several limitations apply, including relatively low read counts, linked SNVs and the abnormal karyotype of HEK293 cells, which originated from a single individual.

SNV calling identified 924,997 variant calls overlapping common dbSNP variants. These included 889,814 calls from 361 ChIP–seq experiments and 35,183 calls from 370 multicycle GHT-SELEX experiments, representing 122,364 unique genomic locations and distinct rsSNP IDs.

Among these variants, 10,009 SNPs corresponded to 12,060 ASBs involving 152 Codebook TFs and 39 positive controls. Significant allelic imbalance was observed in reads supporting the reference and alternative alleles in ChIP–seq and GHT-SELEX datasets. Results are shown in Extended Data Fig. 5 and Supplementary Table 8. SNP calls and ASB results are available from Zenodo⁷³ at https://doi.org/10.5281/zenodo.18224872.

Variant Calling

For variant calling directly from ChIP–seq and GHT-SELEX data, raw ChIP–seq reads and pre-trimmed GHT-SELEX reads¹² were aligned to the hg38 human genome using bwa-mem (v.0.7.1) with default settings.

Reads with more than two mismatches or mapping quality below 10 were removed using filter_reads.py, originally obtained from stampipes (https://github.com/StamLab/stampipes/tree/encode-release; accessed September 2022).

SNV calling and read counting followed a previously described workflow²¹, available at GitHub (accessed December 2025):

1. samtools reheader (v.1.16.1) was used to assign the same sample SM field to all alignment files.
2. SNP calling was performed with bcftools mpileup (v.1.10.2)⁷⁴ using --redo-BAQ --adjust-MQ 50 --gap-frac 0.05 --max-depth 10000, followed by bcftools call with --keep-alts --multiallelic-caller.
3. Multiallelic records were split into biallelic records using bcftools norm with --check-ref x -m -. Variants were filtered with bcftools filter -i “QUAL>=10 & FORMAT/GQ>=20 & FORMAT/DP>=10” --SnpGap 3 --IndelGap 10 and bcftools view -m2 -M2 -v snps, retaining biallelic SNPs covered by at least 10 reads.
4. SNPs were annotated with bcftools annotate using --columns ID,CAF,TOPMED and dbSNP (v.151)⁷⁵.
5. Heterozygous variants on reference chromosomes were retained when genotype quality (GQ) was at least 20, read depth was at least 10 and each allele had at least five supporting reads. Filtering was performed with awk (v.5.0.1).
6. WASP (v.0.3.4)⁷⁶ was used with bwa-mem and filter_reads.py to correct for reference-mapping bias.
7. count_tags_pileup_new.py, an adapted version of count_tags_pileup.py from GitHub (https://github.com/vierstralab/nf-allelic-mapping/tree/main/bin), was used with pysam (v.0.20.0) to calculate allelic read counts. DNase-seq-specific settings were removed from the original script.
8. recode_vcf.py was used to convert the resulting BED files to VCF format.

Triallelic SNVs were split into two biallelic records.

ASB Calling and Annotation

ASB calling was conducted independently for GHT-SELEX and ChIP–seq data. To account for aneuploidy and copy-number variation, relative background allelic-dosage profiles were reconstructed with BABACHI (v.2.0.26) using default settings⁷⁷.

Allelic imbalance was estimated with MIXALIME (v.2.14.17)²¹, beginning with mixalime create. We then fitted a marginalized compound negative-binomial model (MCNB) with mixalime fit. The analysis used MCNB with --window-size 1,000 for GHT-SELEX and --window-size 10,000 for ChIP–seq, reflecting the lower coverage and smaller number of called SNPs in GHT-SELEX.

TF-specific ASB calls were obtained using mixalime test, followed by the TF-wise mixalime combine procedure. This analysis identified 12,060 ASBs at a 5% false-discovery rate (FDR), corrected for multiple tested SNPs.

SNP and ASB genotype qualities were substantially higher than the default thresholds, and most ASB SNPs were supported by at least two datasets, as shown in Extended Data Fig. 5d,e.

We identified 3,564 ASBs overlapping a PWM match with P < 0.001 for the associated TF. ASBs that did not overlap a PWM match may represent marker variants linked to causal SNVs.

For ASBs overlapping PWM matches, we calculated PWM scores for both alleles. Right-tailed P-values were estimated against a uniform background distribution using PERFECTOS-APE (v.3.0.6)⁷⁸. The log₂ fold change between uncorrected alternative-allele and reference-allele P-values, log₂[Alt/Ref], represented the PWM-predicted allelic preference. Positive and negative values indicated preference for the alternative and reference alleles, respectively.

ASBs with an absolute log₂[fold change] greater than 1 were classified as motif concordant or motif discordant. Classification was based on whether the PWM-predicted allelic preference agreed with the allele showing greater ChIP–seq or GHT-SELEX read coverage, either Ref > Alt or Alt > Ref, as shown in Extended Data Fig. 5c.

To provide global support for the Codebook ASB calls, we generated a joint set of 22,064 ASBs using MIXALIME multiple_combine (v.2.28.0), which aggregates allelic-imbalance P-values across all processed datasets (Supplementary Table 9).

For the resulting joint ASB–SNP set, we assessed overlap with GTEx (v.8)²⁴, ADASTRA (v.6.1)²² and the EBI GWAS Catalog (v.1, e115_r2025-12-03_full)²³ using two-sided Fisher’s exact tests.

Reporting Summary

Additional information about the research design is available in the Nature Portfolio Reporting Summary linked to this article.

Source: www.nature.com