Multiscale footprints reveal the organization of cis-regulatory elements

Experimental methods

Cell culture

HepG2 cells (ATCC HB-8065; authenticated by short tandem repeat profiling and tested for mycoplasma by ATCC) were cultured in DMEM (catalogue no. 11965092, ThermoFisher), with the addition of 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin. Cells were incubated at 37 °C in 5% CO₂ and maintained in exponential phase, followed by digestion with TrypLE express (catalogue no. 12604013, ThermoFisher) for preparation of single-cell suspensions.

BMMC sample processing

Frozen human bone marrow mononuclear cells (BMMCs, AllCells; commercial donor cell experiment conducted under approval from Harvard Institutional Review Board) were thawed in a 37 °C water bath for 1 min and transferred to a 15-ml centrifuge tube. Next, 10 ml of prewarmed DMEM with 10% FBS was added dropwise to cells, which were then spun at 400g for 3 min at room temperature. Following removal of supernatant, cells were washed twice in 0.5 ml of PBS with 0.04% bovine serum albumin (BSA). To deplete neutrophils, cells were resuspended in 100 μl of chilled DPBS with 0.2% BSA and 10 μl of human TrueStain FcX (BioLegend, catalogue no. 422302) and incubated on ice for 10 min to reduce non-specific labelling. Cells were then incubated on ice for a further 30 min following the addition of 0.5 μl of biotin antihuman CD15 antibody (BioLegend, catalogue no. 301913). Following immunostaining, 25 μl of MyOne T1 beads was added to the sample for 5 min at room temperature, to capture neutrophils. We then added 900 μl of DPBS with 0.2% BSA to dilute the sample. The sample was placed on a magnet for 3 min, and 1 ml was transferred to a new tube while the sample was on the magnet. Cells were then ready for fixation and the SHARE–seq experiment.

Fixation

Cells were centrifuged at 300g for 5 min and resuspended at 1 million cells ml⁻¹ in PBSI. Cells were fixed by the addition of formaldehyde (catalogue no. 28906, ThermoFisher) to a final concentration of 1%, and incubated at room temperature for 5 min. Fixation was stopped by the addition of 56.1 μl of 2.5 M glycine, 50 μl of 1 M Tris-HCl pH 8.0 and 13.3 μl of 7.5% BSA on ice. The sample was incubated at room temperature for 5 min and then centrifuged at 500g for 5 min to remove supernatant. All centrifugation was performed on a swing-bucket centrifuge. The cell pellet was washed twice with 1 ml of PBSI and centrifuged at 500g for 5 min between washings. Cells were resuspended in PBS with 0.1 U μl⁻¹ Enzymatics RNase Inhibitor and aliquoted for transposition.

SHARE–seq

Following fixation, SHARE–seq was performed as previously described²⁷, with the following modifications. To improve transposition, this was performed using preassembled Tn5 (seqWell, Tagify(TM) SHARE-seq Reagent). To improve RNA capture, we added poly-A to transcripts before reverse transcription. To do this, transposed cells (60 μl) were mixed with 240 μl of poly(A) mix (final concentration of 1× Maxima RT buffer, 0.25 U μl⁻¹ Enzymatics RNase Inhibitor, 0.25 U μl⁻¹ SUPERase RI, 0.018 U μl⁻¹ Escherichia coli poly-A enzyme (New England Biolabs catalogue no. M0276L) and 1 mM ribonucleic ATP). The sample was aliquoted to 50 μl per PCR tube and incubated at 37 °C for 15 min.

Quantification and sequencing

Both scATAC–seq and scRNA-seq libraries were quantified with the KAPA Library Quantification Kit and pooled for sequencing. Single-cell libraries were sequenced on the Nova-seq platform (Illumina) using a 200 cycle kit (read 1, 50 cycles; index 1, 99 cycles; index 2, 8 cycles; read 2, 50 cycles). Bulk libraries were sequenced on the Nova-seq platform (Illumina) using a 100 cycle kit (read 1, 50 cycles; index 1, 8 cycles; index 2, 8 cycles; read 2, 50 cycles).

SHARE–seq data preprocessing

SHARE–seq data were processed using the SHARE-seqV2 alignment pipeline (https://github.com/masai1116/SHARE-seq-alignmentV2/) and aligned to hg38. Open chromatin region peaks were called on individual samples using MACS2 peak caller (v.2.2.9.1)⁵⁶ with the following parameters: –nomodel –nolambda –keep-dup -call-summits. Peaks from all samples were merged, and those overlapping with ENCODE blacklisted regions (http://mitra.stanford.edu/kundaje/akundaje/release/blacklists/hg38-human/hg38.blacklist.bed.gz) were filtered out. Peak summits were extended by 150 bp on each side and defined as accessible regions (for footprinting analyses, these peaks were later resized to 1,000 bp in width). Fragment counts in peaks and TF scores were calculated using chromVAR (v.1.24.0)¹⁰. Cell barcodes with under 30% reads in peaks (fraction of reads in peaks, FRiP) or under 250 unique fragments were removed. Aligned reads were then intersected with peak window regions, producing a matrix of chromatin accessibility counts in peaks (rows) by cells (columns). To examine cell identity, cisTopic (50 topics)⁵⁷ was used for dimension reduction, followed by Louvain clustering. Progenitor populations were subclustered to obtain finer cell identity. The data were projected into two-dimensional space by UMAP⁵⁸. Seurat v.3 (v.5.0.3)⁵⁹ was used to scale the digital gene expression matrix by total UMI count, multiplied by the mean number of transcripts, and values were log transformed.

Generation of BAC naked DNA data

We selected 25 chromatin regions based on overlap with a manually selected set of key transcription factors and differentiation related genes. Bacterial artificial chromosome (BAC) clones (BACPAC Resources) were cultured in lysogeny broth for 14 h. BAC DNA was extracted using a ZR BAC DNA Miniprep Kit (Zymo, catalogue no. D4048), following the manufacturer’s instructions. Purified DNA was quantified using Qubit (ThermoFisher). BAC DNA was tagmented similarly to the SHARE–seq ATAC–seq experiment. Briefly, 50 ng of BAC DNA from multiple clones was pooled for tagmentation following the SHARE–seq transposition condition. Tagmented DNA was purified using a Qiagen Minelute PCR clean-up kit and then PCR amplified for seven cycles. To minimize batch effect, we generated five biological replicates and pooled all materials for sequencing. The library was sequenced on a Nova-seq platform (Illumina) using a 100 cycle kit (read 1, 50 cycles; index 1, 8 cycles; index 2, 8 cycles; read 2, 50 cycles). Sequencing data were processed as for SHARE–seq ATAC–seq data.

Generation of human genomic DNA data

Human (female) genomic DNA was obtained from Promega (catalogue no. G1521). Genomic DNA was digested with the restriction enzyme SbfI-HF (New England Biolabs, catalogue no. R3642L). For each of two replicates, 25 μg of DNA was digested with 200 units of SbfI-HF in a 500 μl reaction at 37 °C overnight. The digested DNA was run on a 1% agarose gel and fragments corresponding to 2–2.5, 2.5–3, 3–4 and 4–5 kb were excised. All fragments from replicate 1 and the 3–4 kb fragment from replicate 2 were purified using the QIAquick gel extraction kit (Qiagen). Purified products were concentrated using the Zymo DNA clean and concentrator kit. Finally, tagmentation and library preparation were performed following the above protocol for BAC DNA tagmentation, using the same ratio of DNA mass to Tn5 when a lower amount of DNA was recovered following size selection.

In vitro footprinting

BAC DNA (50 ng) was incubated with either recombinant c-MYC/MAX (Active Motif, catalogue no. 81087) or CEBPA (OriGene, catalogue no. TP760685), tagmentation buffer (20 mM Tris, 10 mM MgCl₂ and 20% dimethylformamide) and water in a 40-μl volume at room temperature for 1 h. Next, a master mix consisting of 0.15 μl of preassembled Tn5 (seqWell, Tagify(TM)), 4.85 μl of dilution buffer (50 mM Tris, 100 mM NaCl, 0.1 mM EDTA, 1 mM DTT, 0.1% NP-40 and 50% glycerol) and 5 μl of tagmentation buffer was added to the samples for tagmentation in a 50-μl final volume (final TF concentrations of 0, 50 and 100 nM) for 30 min at 37 °C. Tagmented DNA was purified using a Qiagen Minelute PCR clean-up kit and subsequently PCR amplified for five cycles. Samples were then pooled and sequenced on a Nova-seq platform (Illumina). Sequencing data were processed in the same way as bulk ATAC–seq data.

Ageing multiome experiment

Mouse experiments were approved and performed in compliance with Harvard University’s Institutional Animal Care and Use Committee. C57BL6 mice were obtained from either Jackson Laboratory or the National Institute on Aging Aged Rodent Colony (Charles River Laboratory), housed at a density of between two and five per cage in standard ventilated racks and provided with food and water ad libitum in a specific-pathogen-free facility accredited by the Association and Accreditation of Laboratory Animal Committee. Mouse cages contained Anderson’s Bed-o’Cob bedding (The Anderson Inc.), two nestlets (Ancare, 2 × 2-inch² compressed cotton squares) and a red mouse hut (Bioserv). For HSC isolation and flow cytometry, cells from the bone marrow of long bones (two femurs and two tibias per mouse) from young (n = 10, 11 weeks old) and aged (n = 5, 24 months old) male C57BL/6 mice were flushed with a 21-gauge needle into staining medium (HBSS/2% FBS), pelleted and resuspended in ACK lysis buffer for 5 min on ice. Mouse numbers were determined by the anticipated cell yield and input needs for single-cell assay; cells from mice were pooled so that no blinding or randomization was performed. Cells were then washed with staining medium, filtered through a 40-micrometre cell strainer, pelleted and incubated with the following cocktail of rat anti-mouse, biotin-conjugated lineage antibodies on ice for 30 min: CD3 clone C145-2c11 (BioLegend, catalogue no. 100304; 1:100), CD4 clone GK15 (BioLegend, catalogue no. 100404; 1:400), CD5 clone 53-7.3 (eBioscience, catalogue no. 13-0051-85; 1:400), CD8 clone 53-6.7 (BioLegend, catalogue no. 100704; 1:400), CD19 clone 6D5 (BioLegend, catalogue no. 115504; 1:400), B220 clone RA3-6B2 (BioLegend, catalogue no. 103204; 1:200), GR1 (Ly6-G/Ly6-C) clone RB6-8C5 (eBioscience, catalogue no. 13-5931-82; 1:400), Mac1/CD11b clone M1/70 (BioLegend, catalogue no. 101204; 1:800) and Terr119 clone TERR-119 (BioLegend, catalogue no. 116204; 1:100). Cells were then washed in staining medium, with a small aliquot reserved for each sample to serve as a non-depleted control, and lineage depleted using sheep anti-rat Dynabeads (Invitrogen, catalogue no. 1135) on a magnet. Cells were washed, pelleted and incubated with the following cocktail of anti-mouse antibodies on ice for 45 min to identify HSCs: Pacific Orange Streptavidin (Invitrogen, catalogue no. S32365; 1:500), PE/Cy7 Sca1(Ly-6a/E) clone D7 (eBioscience, catalogue no. 25-5981-82; 1:200), APC cKit clone 2B8 (BD Pharmingen, catalogue no. 553356; 1:200), FITC CD48 clone HM48-1 (BioLegend, catalogue no. 103403; 1:200) and PE CD150 clone Tc15-12F12.2 (BioLegend, catalogue no. 115904; 1:200). Following incubation, cells were washed and resuspended in staining medium, and 7-AAD (BD Pharmingen, catalogue no. 559925; 1:50) added immediately before flow cytometry. Cell sorting of HSCs (Live Lin⁻ Sca1⁺ cKit⁺ CD48⁻ CD150⁺) was performed on a BD FACS Aria II. Cells within the same age group were sorted into the same tube for subsequent sequencing. Data analysis was performed using BD FACS Diva (v.8.0.2) and FlowJo (v.10.8.2) software. Data processing was performed using Cell Ranger ARC 2.0.0.

Following sorting, nuclei were isolated following the 10x Genomics demonstrated protocol Low Cell Input Nuclei Isolation, which is described in the CG000365 User Guide. Nuclei were then processed using the Chromium Single Cell Multiome ATAC + Gene Expression kit (10x Genomics), following the manufacturer’s instructions, to obtain between 2,000 and 10,000 cells per sample. Libraries were sequenced on an Illumina Next-seq system using the following sequencing formats: read 1, 28; i7 index, 10; i5 index, 10; read 2, 44 (scRNA-seq); read 1, 30; i7 index, 8; i5 index, 24; read 2, 30 (scATAC–seq). Data processing was performed using CellRanger ARC 2.0.0 software from 10x Genomics.

Tn5 sequence bias modelling

Getting Tn5 insertion counts

The ends of fragment files were shifted by +4/−4 to obtain the centre of the 9 bp staggered end created by Tn5 transposition. The number of insertions at each single base-pair position within each cCRE from each sample was then quantified and stored in a sample-by-cCRE-by-position three-dimensional tensor for rapid data retrieval.

Data preprocessing

The model takes local DNA sequence context as input and predicts single base-pair resolution Tn5 bias. To this end, the ±50-bp DNA sequence surrounding each position of interest is encoded by one-hot encoding into a 101 × 4 matrix and used as model input. For the prediction target, we use local relative Tn5 bias as the target value. More specifically, the raw Tn5 insertion count at each position is divided by the average Tn5 insertion count within a ±50-bp window. Positions with low local coverage (fewer than 20 insertions per base pair) were removed to guarantee quality of training data. To facilitate model training, the resulting observed Tn5 bias values were log₁₀ transformed and rescaled. For dataset partition, we randomly split all BACs into 80%, 10% and 10% for the training, validation and test sets, respectively; in other words, all data originating from the same BAC belonged to the same partition. This is to prevent overlapping local sequence contexts being placed in both training and testing datasets, which might lead to overestimation of performance. To guarantee equal coverage of examples with different bias levels, we binned all training examples into five bins based on their Tn5 bias value, and upsampled each bin so that all bins ended up with the same number of examples. In addition, given the symmetric nature of Tn5 insertion, we generated reverse-complement sequences of training examples as data augmentation. The original and reverse-complement data were combined, shuffled and then used for model training.

Model architecture

The convolutional network consists of three convolution and max-pooling layers and two fully connected layers. Each convolution and max-pooling layer performs convolution, rectified linear unit (ReLU) nonlinear activation⁶⁰ and max pooling sequentially. We used 32 filters of width 5 for each layer, along with the ‘same’ padding mode and stride size of 1. The two following fully connected layers have output dimension of 32 and 1, respectively; ReLU activation is used by the first fully connected layer and linear activation by the second (that is, the final output layer).

Model training and evaluation

The model was trained on the training set, and hyperparameters were optimized based on performance on the validation set. Final performance of the frozen model was evaluated on the test set. The model was implemented using Keras⁶¹, trained with mean-square error as loss function and optimized using the Adam optimizer⁶² with default parameters. Training was performed with a batch size of 64 and early stopping based on model loss on the validation set.

Benchmarking with other Tn5 bias models

Methods including k-mer models (k = 3, 5, 7) and PWM methods (single nucleotide and dinucleotide) were included in benchmarking. For k-mer methods, the foreground and background frequencies for all possible k-mer sequences were quantified. Foreground/background frequency ratio was used as estimated Tn5 bias for the corresponding k-mer. For single-nucleotide PWM, we calculated foreground and background base frequencies within a ±10-bp window (total length, 21 bp) and computed the PWM of Tn5 insertion. Dinucleotide PWM scores were calculated using TOBIAS¹⁸ with default settings. It is suggested to train custom ChromBPNet bias models on inaccessible chromatin regions for each dataset, to represent Tn5 sequence bias and achieve the highest-quality models⁶³. Accordingly, we trained ChromBPNet on HepG2 and K562 ATAC–seq data to evaluate performance at its recommended setting.

Genome-wide Tn5 bias reference tracks

Sequences of reference genomes for Homo sapiens (hg38), Mus musculus (mm10), Drosophila melanogaster (dm6), Saccharomyces cerevisiae (sacCer3), Caenorhabditis elegans (ce11), Danio rerio (danRer11) and Pan troglodytes (panTro6) were downloaded from the UCSC genome browser website⁶⁴ (https://hgdownload.soe.ucsc.edu/goldenPath/). The aforementioned Tn5 bias neural network model was applied to each position in the reference genomes to generate genome-wide Tn5 bias tracks.

Computation of footprint scores

For detection of DNA–protein interactions at different scales within cCREs, we implemented a framework for computing footprint scores for each base-pair position in the cCRE. In short, for each single base-pair position, we define a centre footprint window and flanking windows on both sides (Fig. 1a). We then calculated the observed ratio of centre/(centre + flanking) Tn5 insertion counts. The foreground observed ratio was compared with background distribution to calculate statistical significance, which was then converted to a footprint score.

Estimation of background dispersion

Given a specific combination of centre bias, flanking bias and local coverage, we expect a certain distribution of centre/(centre + flanking) insertion ratio when no protein is bound; this is defined as background distribution. Such background distribution can be estimated using BAC naked DNA Tn5 insertion data. To this end, we first randomly sampled 100,000 positions from the BAC dataset and retrieved their local coverage (defined as the total insertion number in centre and flanking areas), centre bias and flanking bias. Then, for each sampled position A, we identified the 500 nearest-neighbour positions (NN₁–NN₅₀₀) in the three-dimensional space of centre bias, flanking bias and local coverage. To ensure that each dimension was weighed equally, the values of each dimension were first normalized to zero mean and unit variance. The 500 nearest-neighbour observations can be considered as background observations with nearly identical bias and coverage, with the centre/(centre + flanking) ratio of NN₁–NN₅₀₀ forming the background distribution of position A. Therefore, for each of the 100,000 sampled positions, we can calculate the mean and standard deviation of its background ratio distribution. This allows us to train a background dispersion model that takes the tuple (centre bias, flanking bias and local coverage) as input and efficintly predicts the mean and standard deviation of background distribution. To ensure that the model is exposed to training examples with a wide range of local coverage, we downsampled the BAC dataset to 50, 20, 10, 5 and 1% of the original sequencing depth. Finally, we trained a neural network with a single hidden layer (32 nodes, ReLU activation⁶⁰) and linear output layer activation. The dataset was randomly split into 80% training, 10% validation and 10% test. The model was implemented using Keras⁶¹, and trained on the training dataset with mean-squared error loss using the Adam optimizer^61,62. Early stopping was determined using loss on the validation set, and performance of the final model was evaluated on the test set. In addition, we trained separate models for each footprint radius due to the marked differences in total centre or flank bias when footprint radius varied (Supplementary Notes).

Calculation of footprint scores

For each position in the cCRE, we define a centre footprint window and flanking windows on both sides. We first calculate the foreground observed centre/(centre + flanking) ratio of Tn5 insertion counts, then we apply the pretrained background dispersion model to calculate the mean and standard deviation of its background distribution. We next use a lower-tailed z-test to calculate the P value for footprinting. If the observed ratio is significantly lower than background distribution, this position is likely to be bound by a protein. More specifically, to avoid calling footprints at positions at which only one flanking side shows higher Tn5 insertion than the centre window but not the other, we perform centre-versus-left and centre-versus-right tests separately and retain the higher P value (Supplementary Notes). The −log₁₀(P values) are smoothed by running-max and running-mean smoothing and then used as the final footprint scores.

Aggregate footprinting

For calculation of aggregate footprints, Tn5 insertions surrounding TF or nucleosome binding sites across the genome are first aggregated and then used to calculate footprint scores. For TFs, we selected sites with a matched TF motif using motifmatchr version 1.24.0 (ref. ⁶⁵) (P cutoff 1 × 10⁻⁵) and those overlapping with a ChIP–seq peak of the corresponding TF. For motif matches on the reverse strand, the Tn5 insertion profile surrounding the motif is inverted so that insertions for different sites are aligned in the same direction.

Footprint-to-TF prediction

Note that seq2PRINT was chosen as our primary TF binding predictor, but we still provide this lightweight footprint-to-TF prediction model for its speed advantages. For comparison between footprint-to-TF prediction and seq2PRINT-based TF binding prediction, see ‘Multiple methods to predict TF binding from multiscale footprinting’ in Supplementary Notes.

To predict the landscape of TF binding, we trained a binary classifier that predicts whether any TF motif site is bound by the corresponding TF. Motif sites are identified by the matchMotifs function in the motifmatchr package⁶⁵. All sites with a matching P value below 5 × 10⁻⁵ are retained. For any TF motif site, we use multiscale (20, 40, 60, 100, 160 or 200-bp diameter) footprints within a ±100-bp local area centred around the motif, as well as a motif match score as input to the model. The motif match score returned by the matchMotifs function is quantile transformed to uniform distribution. As a result, by combining of the 201-dimensional footprint vectors from six different scales with a single motif match score, we create a 1,207-dimensional vector as the final model input. The first 1,206 dimensions of footprint scores are standardized individually to zero mean and unit variance. For the prediction target, we assign a label of 1 to all sites overlapping with a ChIP peak of the same TF, and a label of 0 to those not overlapping with ChIP. Some TFs were found to have a very low percentage of motif sites overlapping with ChIP (10% or under), potentially due to low quality of the motif or ChIP dataset, and these were removed from model training and testing. We also added additional negative examples as well as reverse-complement examples for data augmentation.

For data partition, HepG2 data were used as training data and GM12878 data (previously published in the original SHARE–seq paper²⁷) as validation. After fixing model hyperparameters, HepG2 and GM12878 data were combined into a larger training dataset to train a final footprint-to-TF model. The final model was tested on K562 scATAC data⁶⁶, as well as on three cell types (naive B cells, CD14 monocytes and late-erythroid cells) in the human BMMC SHARE–seq dataset.

Model architecture and training

The TF binding prediction model is a neural network model with two hidden layers (128 + 32 nodes). ReLU activation⁶⁰ is used by both hidden layers, and sigmoid activation by the final output layer. The model was implemented using Keras⁶¹, and trained on the training dataset with a batch size of 128 using the Adam optimizer⁶⁷. Binary cross-entropy is used as the loss function. Early stopping was used based on model loss on the validation set.

ChIP validation and benchmarking with previous methods

For evaluation of model performance, we used ChIP–seq as ground truth and validated predicted binding events. HepG2 and GM12878 data were downloaded from UniBind⁶⁸ (https://unibind.uio.no/search) for model training. ChIP–seq for BMMC cell types was downloaded from cistromeDB⁶⁹. For benchmarking with previous methods, to ensure that we included only high-quality TF binding sites, we downloaded K562 ChIP-based TF binding data from UniBind. For cistromeDB datasets, we applied quality control filters as specified on the cistromeDB website (http://cistrome.org/db/#/about). More specifically, we included the following filters: FRiP 0.01 or over, FastQC 0.25 or over, uniquely mapped ratio 0.6 or over, peaks with fold change above ten, 500 or more, peaks union DNase I hypersensitive site ratio 0.7 or above and PCR bottleneck coefficient 0.8 or above. Datasets with the following cell type labels were included: Monocyte, B Lymphocyte, Erythroid cell, Erythroid Progenitor Cell and Erythroid progenitor. If there was more than one dataset for the same TF, the intersection of all datasets for the same TF was retained as the final list of high-confidence binding sites for model training.

The K562 datasets from UniBind were used for benchmarking with previous methods, including DNase I footprinting, TOBIAS and sequence attribution scores obtained from ChromBPNet. In short, the same ATAC–seq data were used as input to all ATAC footprinting methods and ChromBPNet, whereas DNase I footprinting results for K562 were obtained from ENCODE datasets ENCLB253REF, ENCLB843GMH and ENCLB096YUZ from ref. ¹¹. To guarantee fair comparison, we used the same set of motif match positions (previously published in ref. ⁷⁰) as candidate binding sites, and mapped the predicted scores of each method to these sites for comparison. For TFs with multiple UniBind datasets, their intersection was used for benchmarking. Then, for each method, we ranked candidate sites by predicted binding score and evaluated the precision of prediction using the top 10% of sites. Visualization of predicted and ground truth binding sites was done with the Gviz package (v.1.46.1)⁷¹. Furthermore, to evaluate the false-positive rate of each model, we also tested all three ATAC-based models on our BAC naked DNA data. The same data were used as input to each model, and the number of predicted binding events was used to represent false-positive predictions.

Footprint-to-nucleosome prediction

Input data

For prediction of nucleosome occupancy, we trained a regression neural network model. For any genomic position, we used multiscale (20, 40, 60, 100 and 160 bp in diameter) footprints within a ±100-bp local area as input to the model. A scale of 200 bp was not included, to prevent the model from learning the co-occupancy of adjacent nucleosomes. To train this model, we used previously published chemically mapped nucleosome occupancy data on S. cerevisiae⁷² as the training label, and computed multiscale footprints using previously published S. cerevisiae ATAC–seq data²⁰ as training input. We retained observations in regions with local ATAC–seq coverage above 10, and rescaled the 5 and 95% percentiles of nucleosome occupancy values to 0 and 1, respectively. For data partition, we randomly split all data by chromosomes into training (chrVII, chrXI, chrIX, chrI, chrV, chrX, chrVIII and chrXII), validation (chrIV and chrII) and test (chrVI, chrXVI, chrXIII, chrIII, chrXIV and chrXV) sets.

Model architecture and training

The nucleosome prediction model is a neural network model with two hidden layers (64 + 16 nodes). ReLU activation is used by both hidden layers, and linear activation by the final output layer. The model was implemented using Keras. The model was trained on the training dataset with a batch size of 128 using the Adam optimizer. Mean-squared error was used as the loss function. Early stopping was used based on model loss on the validation set.

Performance evaluation

Model performance was evaluated using data on the test yeast chromosomes mentioned above. In total, 859 regions of length 5 kb on the test yeast chromosomes (chrVI, chrXVI, chrXIII, chrIII, chrXIV and chrXV) were used for testing. We detected summits of predicted nucleosome signal and ground truth nucleosome occupancy as predicted and observed nucleosomes, respectively. Precision was calculated as the percentage of predicted nucleosomes having a ground truth nucleosome within a certain cutoff distance (50 or 75 bp in this study). Recall was calculated as the percentage of ground truth nucleosomes having a predicted nucleosome within the same cutoff distance.

Seq2PRINT

Model architecture and training

The seq2PRINT model is a convolutional neural network that takes one-hot encoded DNA sequences (a DNA sequence of length L encoded into an L × 4 matrix, where each row has one element set to 1 representing the specific nucleotide) as input and predicts the corresponding multiscale footprints at base-pair resolution. The architecture is depicted in Extended Data Fig. 4a. To facilitate articulation of the architecture, we divide the seq2PRINT model into three parts. Although the architecture of seq2PRINT can be flexible depending on the available computational resources and the depth and scale of the training data, for all results in this work we build the model and choose the hyperparameters as follows.

The first convolutional layer consists of 1,024 filters of width 21 bp with Gaussian error linear unit (GELU) activation, aiming to capture informative sequence patterns from the input DNA sequences (that is, sequence motifs). The output is subsequently passed to eight layers of convolutional blocks with residual connection. Each convolutional block consists of one grouped dilated convolutional layer (n filters = 1,024, width = 3, groups = 8, dilation = 2ⁱ, i = 1,…, 8), followed by a position-wise feed-forward layer (implemented as a convolutional layer with n filters = 1,024, width = 1). Batch normalization layers with GELU activations are inserted between these convolutional blocks. The increased dilation rate results in an expanding receptive field for the neural network, capturing the relationship of the sequence patterns and their context. With eight layers of the convolutional block, seq2PRINT has a receptive field of ±920 bp for each cCRE. The use of grouped convolutional layers enables the construction of wider models that capture richer information with limited parameters, providing a regularization effect and reduced computational complexity. Finally, the output of the stacks of convolutional blocks is passed to the output layers.

In this paper we designed two output layers, a multiscale footprint layer (a convolutional layer of filter width 1) that outputs the multiscale footprints and an accessibility layer (a global average pooling layer followed by a fully connected layer) to predict the number of Tn5 insertions in a specific cCRE.

To facilitate the training of the seq2PRINT model, we implemented batch-efficient multiscale footprint calculation on graphics processing units, which follows the same mathematical models as the described footprint calculation, the only difference being that it outputs the z-statistics rather than the P value calculated from the z-test.

During training, model weights are updated to minimize the following loss function:

where MSE represents mean-squared error \({\rm{MSE}}(x,y)=\sum {(x-y)}^{2}\), \({{\rm{footprint}}}_{{\rm{PRINT}}}\) and \({{\rm{footprint}}}_{{\rm{pred}}}\) represent multiscale footprints calculated by the PRINT framework and seq2PRINT model, respectively, and n_obs and n_pred represent the observed and predicted Tn5 insertions, respectively, in this region.

Notably, the gradient back-propagation for the accessibility layer is broken before the convolutional blocks: in other words, the sequence patterns and relationships among them learned by the preceding layers of the seq2PRINT model are driven solely by the multiscale footprint objective. The accessibility output layer and corresponding loss function reweight these learned sequence features for interpretation purposes only (see the following section for a detailed explanation). This design alleviates the need to choose weights between footprint loss and accessibility loss, and also makes seq2PRINT a footprint-driven sequence model, differentiating it from previous accessibility-driven models (for example, Basset⁷³).

The seq2PRINT model is optimized with the Adam optimizer at a learning rate of 3 × 10⁻⁴, and uses exponential moving averages to stabilize and improve model convergence. In this study we used chromosome-based, fivefold cross-validation, with outputs across folds averaged for use as the final predictions.

Deriving sequence attribution scores

We use the DeepLIFT⁷⁴ method to calculate sequence attribution scores, which represents how each base pair in a given input DNA sequence contributes to a specific scalar output from the trained seq2PRINT model. The output of the accessibility layer is a scalar for each region, making it naturally suitable as the target for DeepLIFT to calculate attribution scores. However, the multiscale footprint layer is not a scalar but a matrix of size L × number of scales. Therefore, we designed two strategies to summarize output footprints into a scalar value.

Both strategies involve conversion of the predicted z-statistics to log(P value) footprint scores as the PRINT framework. The first strategy involves manual inspection, as demonstrated in Fig. 2b, where we locate the region × scale of interest from the observed and predicted multiscale footprints, sum footprint scores within the region and calculate the corresponding sequence attribution score. The second strategy sums footprint scores within the whole peak region, and is more suitable for genome-wide calculations. Without further specification, we refer to sequence attribution scores calculated from the accessibility layer as count sequence attribution scores, and those calculated from the footprint layer as footprint sequence attribution scores.

For all results in this manuscript, we use 20 dinucleotide-shuffled input sequences as reference sequences for the DeepLIFT algorithm. We adapted implementation of the DeepLIFT algorithm from the DeepSHAP package to accommodate the custom nonlinear functions used in this framework; custom implementation is available as part of https://github.com/buenrostrolab/scPrinter.

De novo motif identification based on sequence attribution scores

TF-MoDISco (tfmodisco-lite v.2.2.1)⁷⁵ was utilized to infer de novo motifs based on sequence attribution scores. Briefly, TF-MoDISco identifies local regions in input sequences with high sequence attribution scores (seqlets), then aligns and clusters similar seqlets into groups of de novo motifs. We set the number of seqlets as 1,000,000, and the remaining parameters as default. The de novo discovered motifs are assigned to known motifs using TomTom (meme suite v.5.5.7)⁷⁶. To infer the matching of de novo motifs at cCREs, we use the software finemo (commit no. 830d7f3)⁷⁷, which accepts both de novo motifs and sequence attribution scores.

TF binding prediction using attribution scores

The calculated sequence attribution scores highlight TF binding sites affecting footprints and accessibility. We thus trained a binary classifier that is similar to the footprint-to-TF model but, rather than multiscale footprints, it accepts both count and footprint sequence attribution scores.

The training and validation TF binding sites remain the same as in the previous footprint-to-TF model. For each motif site, the features include the count and footprint sequence attribution scores within a ±100-bp area centred around the motif, the motif-matching score and also the Pearson correlation between the motif of interest and sequence attribution score at the motif-matching site, in a 405-dimensional vector.

The TF binding model consists of three hidden layers (256, 128 and 64 neurons), with GELU activations and 0.25 dropout rates in between. The model is trained using the Adam optimizer, with binary cross-entropy as the objective function.

This TF binding model is the final TF binding prediction model that we used for Figs. 2–4, due to its superior performance. For a comparison between footprint-to-TF prediction and seq2PRINT-based TF binding prediction, see ‘Multiple methods to predict TF binding from multiscale footprinting’ in Supplementary Notes.

LoRA enables efficient sequence modelling for scATAC–seq data

To make the sequence model much more scalable on scATAC–seq with diverse cell types or cell states, we used the LoRA technique for parameter-efficient fine-tuning of neural network models. Specifically, for given scATAC–seq data, we first train a seq2PRINT model (referred to as the pretrained model) by aggregating Tn5 insertions over all cells in the dataset. Next, for each pseudo-bulk aggregating cells over specific cell states, we use the LoRA fine-tuning technique to derive a pseudo-bulk-specific seq2PRINT model. We describe this fine-tuning process as follows.

For any neural network layer parameterized by initial weights W₀, where \({W}_{0}\in {R}^{{d}_{{\rm{p}}}}\), the LoRA model would learn an updating parameter \(\varDelta W\in {R}^{{d}_{{\rm{p}}}}\), where d_p is the number of learnable parameters of the layer, and the sum of these two parameters is used as the weight parameter W for the resulting fine-tuned neural network \(W={W}_{0}+\varDelta W\). For the residual grouped convolutional layer we used in seq2PRINT,

If we fine-tune the model individually for n pseudo-bulks of interest, this results in a total of n × d_p parameters to be learned.

Motivated by the LoRA model, we instead learn a low-rank decomposition of this updating parameter, guided by the intrinsic low ranking of cell states captured by cell embeddings. Specifically, with d_e as the dimension of single-cell embeddings and r as the rank hyperparameter, we learn two sets of weights, A and B, where \(A\in {R}^{{d}_{{\rm{e}}}\times r+r\times \frac{n\_{\rm{filter}}}{{\rm{group}}}}\)\(B\in {R}^{{d}_{{\rm{e}}}\times r+r\times n\_{\rm{filter}}\times {\rm{kernel}}\_{\rm{size}}}\). This design reduces the number of parameters to be learned from n × d_p to \(r\times ({d}_{{\rm{p}}}^{* }+{d}_{{\rm{e}}})\), where \({d}_{{\rm{p}}}^{* }=\frac{n\_{\rm{filter}}}{{\rm{group}}}\,+\)\(n\_{\rm{filter}}\times {\rm{kernel}}\_{\rm{size}} < \,\frac{n\_{\rm{filter}}}{{\rm{group}}}\times n\_{\rm{filter}}\times {\rm{kernel}}\_{\rm{size}}\). In this study we chose r = 32, which is much smaller than the number of pseudo-bulks that we studied.

Predicted impact of TF motifs on footprints

To probe the relationships between multiscale footprints and DNA sequences learned by the seq2PRINT model, we generated the marginalized prediction from seq2PRINT given a known or de novo discovered motif. For a de novo discovered motif of interest, we first identified its consensus sequence by taking the nucleotide with highest probability at each position. We then randomly selected 25,000 cCREs from the dataset, planted the consensus sequence at the centre and averaged model predictions around the sampled cCREs. For a known motif, we gathered its motif-matching positions in cCREs, scrambled these and averaged model predictions. In both approaches, we calculated marginalized prediction as the difference between the presence and absence of a given motif .

Tracking TF binding dynamics across human haematopoiesis

Generation of pseudo-bulks

Single cells in the human BMMC dataset were first embedded into lower dimensional space using cisTopic⁵⁷, and then grouped into 1,000 pseudo-bulks based on their spatial proximity in the cisTopic space. More specifically, we first sampled 1,000 cells as pseudo-bulk centres and then identified the k-nearest neighbours (KNN, k = 5,000) of each centre cell in the cisTopic space as other members of the same pseudo-bulk. We reasoned that sampling of centre cells with low local connectivity could help increase coverage of the state space by prevention of oversampling of densely connected local neighbourhoods. Therefore, we first randomly sampled 10,000 scaffold cells and used these to construct a KNN graph (k = 10), then selected the 1,000 cells with the lowest in-degree in the KNN graph as pseudo-bulk centres.

Computing pseudo-time

Pseudo-time along human haematopoietic lineages was computed using the Palantir package (v.1.0.0)⁷⁸. To reduce computing time, we randomly sampled 100,000 cells from the human BMMC dataset as scaffold cells. The cisTopic embedding of both scaffold and pseudo-bulk centre cells was used as input to Palantir.

Footprinting repressor TFs

Because BHLHE40 has been studied extensively for its role in T cells, the BHLHE40 footprinting results in Extended Data Fig. 2 were obtained using CD8⁺ T cells from the human bone marrow dataset.

TF binding in erythroid and lymphoid cell cCREs

For erythroid and lymphoid development trajectories, we first identified relevant genes by selecting those with a correlation between RNA level and pseudo-time above 0.5, respectively. We then filtered cCREs within ±50 kb of the transcriptional start sites of related genes and correlation between cCRE accessibility and pseudo-time above 0.5. For all cCREs, we located candidate TF binding sites activated during development by retaining those at which the correlation between seq2PRINT TF binding score and pseudo-time was above 0. The remaining TF binding sites were used for dynamic quantifications. To quantify the inferred timing of binding, we rescaled the TF binding score of each binding site to [0,1] and calculated the area under the curve, with higher values representing earlier rise in binding signals and vice versa.

PCA measures TF binding pattern complexity at each cCRE

To reveal the complexity of TF binding patterns within each cCRE across diverse cell populations in the human BMMC dataset, we used a principal component analysis (PCA)-based method. To reduce computational complexity, we first collapsed the 1,000 LoRA fine-tuned seq2PRINT models (corresponding to 1,000 pseudo-bulks) into 20 models corresponding to 20 cell types. Model weights of pseudo-bulks corresponding to the same cell type were averaged during this process. We then generated cCRE-wide TF binding scores for these 20 cell types. We tiled each cCRE with 10-bp windows, and averaged TF binding scores within each window. PCA was then used as a dimension reduction method on this 20 × 10-bp window matrix for each cCRE, with the minimum number of PCs explaining 98% of the variance being used to quantify the complexity of TF binding patterns.

Characterizing age-related cCRE reorganizations

Data preprocessing

Cells with FRIP lower than 0.3 and depth less than 300 were first removed. In addition, we used ArchR (https://www.archrproject.com/)⁷⁹ to calculate doublet scores for each single cell and removed those with the top 5% of doublet scores. The remaining cells were then processed with the Seurat V3 package (v.5.0.3)⁸⁰. Cells were embedded into lower dimensional space using latent semantic indexing⁸¹ and then clustered. Seurat clusters corresponding to HSCs were selected for pseudo-bulking and downstream differential testing. Cells with the LinNeg FACS sort label were excluded from HSC-specific analyses. To identify representative cell states, we used SEACells (https://github.com/dpeerlab/SEACells/tree/main)⁸² to identify 100 representative cell states across HSCs. Representative cells were used as centres to form pseudo-bulks. Each pseudo-bulk was generated by serial inclusion of nearest-neighbour cells from the centre cell in order of increasing distance, until we reached a total of 5 million fragments.

Defining HSC subpopulations

The scRNA expression data obtained from 10x Multiome was first filtered by removal of the top 5% of highly expressed genes, as well as ribosomal and mitochondrial genes; cells with fewer than 100 RNA counts were removed. We then ran SCTransform (v.0.4.1) with 5,000 variable features to normalize the data. Normalized values were used as input to Spectra⁸³ (https://github.com/dpeerlab/spectra). We next ran Spectra initialized with default HSC and global pathways, and including additional gene sets obtained from the published literature; gene sets with fewer than five genes covered were removed. We scored the expression of each Spectra program in each pseudo-bulk. To this end, we used DESeq2 to normalize the pseudo-bulk-by-gene RNA count matrix, and then rescaled per-gene values. For any specific Spectra program, we generated 100 background programs consisting of genes with matched overall expression levels. The average gene expression in the Spectra program of interest was then compared wih that in sampled background programs to derive a z-score. Finally, the pseudo-bulk-by-program z-score matrix was used to cluster pseudo-bulks into HSC subpopulations.

Differential testing

Differential RNA and cCRE accessibility testing was performed using DESeq2 (v.1.42.1)⁸⁴. For each pseudo-bulk, we quantified total RNA read counts for each gene and Tn5 insertion counts in each peak (resized to 1 kb), and used DESeq2 to identify significant differential genes and cCREs, with age as the covariate.

Predicting Ets/Runx dimer structures using Alphafold3

For each de novo identified motif representing the Ets/Runx composite motif, we generated the consensus DNA sequence. We then input the same segments of the amino acid sequences of Ets and Runx as those in PDB-4L0Z (Runx, 54–212; Ets, 332–432), along with the consensus DNA sequence and its reverse complements, into Alphafold3. For the structure corresponding to the validated dimer structure, we aligned the AF3-predicted structure to the PDB structure using PyMol (v.2.6) and calculated root mean-squared deviation. All structures were visualized using ChimeraX v.1.8. Solvent-inaccessible bases were identified using the ChimeraX interface function with default parameters.

Reporting summary

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