Beyond accessibility: ATAC-seq footprinting unravels kinetics of transcription factor binding during zygotic genome activation

While footprinting analysis of ATAC-seq data can theoretically enable investigation of transcription factor (TF) binding, the lack of a computational tool able to conduct different levels of footprinting analysis has so-far hindered the widespread application of this method. Here we present TOBIAS, a comprehensive, accurate, and fast footprinting framework enabling genome-wide investigation of TF binding dynamics for hundreds of TFs simultaneously. As a proof-of-concept, we illustrate how TOBIAS can unveil complex TF dynamics during zygotic genome activation (ZGA) in both humans and mice, and explore how zygotic Dux activates cascades of TFs, binds to repeat elements and induces expression of novel genetic elements. TOBIAS is freely available at: https://github.com/loosolab/TOBIAS.

true binding sites for each TF, which we used for validating each method after application to 106 the matched ATAC-seq data (see Methods; Validation). In terms of Tn5 bias correction, as 107 well as visualization of aggregate footprints, we found that TOBIAS clearly outperforms other 108 bias-correction tools in uncovering footprints and thereby distinguishing between 109 bound/unbound sites (Supp. Figure 2a, Supp. File 1). For the task of protein binding prediction, 110 we found that TOBIAS significantly outperformed the other de novo tools HINT-ATAC, PIQ 111 and Wellington (Supp. Figure 2b) and performed equally well as msCentipede overall (Supp. 112 Figure 2c). Notably, TOBIAS also showed robust performance across individual cell types 113 (Suppl. Figure 2d). Looking at individual TFs, TOBIAS outperforms msCentipede for factors 114 with a notable gain of footprints after Tn5 bias correction (Supp. Figure 2e), once again 115 highlighting the advantage of taking Tn5 bias into account. Although msCentipede implements 116 a motif centric learning approach, which can take TF specific binding patterns into account, it 117 did not yield overall higher accuracy in comparison to TOBIAS. Additionally, the approach of 118 building individual TF models took 300 times longer to compute than performing footprinting 119 using TOBIAS (Supp. Figure 2f and Table 1). Such learning approaches are therefore greatly 120 limited in the number of TFs and conditions that can realistically be included in an analysis. In 121 conclusion, we find that the TOBIAS framework shows unprecedented accuracy and 122 performance in the field of ATAC-seq footprinting. 123 In order to confirm the improvement of footprint detection after Tn5 bias correction, we made 124 use of another exemplary dataset derived from hESC 20 . Importantly, besides cases where the 125 footprint was hidden by Tn5 bias (Supp. Figure 3a; ZSCAN4), we also identified TFs for which 126 the motif itself disfavors Tn5 integration, thereby creating a false-positive footprint in 127 uncorrected signals, which disappears after Tn5 bias correction (Supp. Figure 3a; HLX). 128 Utilizing a footprint depth metric as described by 16 (Supp. Figure 3b) across uncorrected, 129 expected and corrected Tn5 signals, we found a high correlation between uncorrected and 130 expected footprinting depths (Supp. Figure 3c). In contrast, this effect vanished after TOBIAS 131 correction (Supp. Figure 3d), effectively uncovering TF footprints which were superimposed 132 by Tn5 bias. From a global perspective, taking 590 TFs into account, TOBIAS generated a 133 measurable footprint for 64% of the TFs (Supp. Figure 3e). This is in contrast to previous 134 reports wherein it has been suggested that only 20% of all TFs leave measurable footprints 135 16 . To summarize, we found that TOBIAS exceeded other tools in terms of uncovering 136 footprints and correctly identifying bound TF binding sites. 137

Footprinting uncovers transcription factor binding dynamics in mammalian ZGA 138
To demonstrate the full potential of TOBIAS, in particular in the investigation of processes 139 involving only few cells, we analyzed a series of ATAC-seq datasets derived from both human 140 and murine preimplantation embryos at different developmental stages ranging from 2C, 4C, 141 8C to ICM in addition to embryonic stem cells of their respective species 20, 21 . Altogether, 142 TOBIAS was used to calculate footprint scores for a list of 590 and 464 individual TFs across 143 the entire process of PD of human and mouse embryos, respectively. After clustering TFs into 144 co-active groups within one or multiple developmental timepoints, we first asked whether the 145 predicted timing of TF activation reflects known processes in human PD. Intriguingly, we found 146 10 defined clusters of specific binding patterns, the majority of which peaked between 4C and 147 8C, fully concordant with the transcriptional burst and termination of ZGA (Figure 2a and Supp. 148 Table 1). 149 Two clusters of TFs (Cluster 1+2; n=83) displayed highest activity at the 2-4C stage and 150 strongly decreased thereafter, suggesting that factors within these clusters are likely involved 151 in ZGA initiation. We set out to classify these TFs, and observed a high overlap with known 152 maternally transferred transcripts 22 (LHX8, BACH1, EBF1, LHX2, EMX1, MIXL1, HIC2, 153 FIGLA, SALL4, ZNF449), explaining their activity before ZGA onset. Importantly, DUX4 and 154 DUXA, which are amongst the earliest expressed genes during ZGA 5, 6 , were also contained 155 in these clusters. Additional TFs included HOXD1, which is known to be expressed in human 156 unfertilized oocytes and preimplantation embryos 23 and ZBTB17, a TF mandatory to generate 157 viable embryos 24 . Cluster 6 (n=67) displayed a particularly prominent 8C specific signature, 158 that harbored well known TFs involved in lineage specification such as PITX1, PITX3, SOX8, 159 MEF2A, MEF2D, OTX2, PAX5 and NKX3.2. Furthermore, overlapping TFs within Cluster 6 160 with RNA expression datasets ranging from the germinal vesicle to cleavage stage 5 , 12  161   additional TFs (FOXJ3, HNF1A, ARID5A, RARB, HOXD8, TBP, ZFP28, ARID3B, ZNF136,  162 IRF6, ARGFX, MYC, ZSCAN4) were confirmed to be exclusively expressed within this time 163 frame. Taken together, these data show that TOBIAS reliably uncovers massively parallel TF 164 binding dynamics at specific time points during early embryonic development. 165

Transcription factor scores correlate with footprints and gene expression 166
To confirm that TOBIAS-based footprinting scores are indeed associated with leaving bona 167 fide footprints we utilized the ability to visualize aggregated footprint plots as implemented 168 within the framework. Indeed, bias corrected footprint scores were highly congruent with 169 explicitly defined footprints (Figure 2b) of prime ZGA regulators at developmental stages in 170 which these have been shown to be active 7 . For example, footprints associated with DUX4, 171 a master inducer of ZGA, were clearly visible from 2C-4C, decreased from 8C onwards and 172 were completely lost in later stages, consistent with known expression levels 20 and ZGA onset 173 in humans. Footprints for ZSCAN4, a primary DUX4 target 5 , were exclusively visible at the 174 8C stage. Interestingly, GATA2 footprints were exhibited from 8C to ICM stages which is in 175 line with its known function in regulating trophoblast differentiation 25 . As expected, CTCF 176 creates footprints across all timepoints. Strikingly, we observed that these defined footprints 177 were not detectable without TOBIAS mediated Tn5 bias correction (Supp. Figure 3f). These 178 data show that footprint scores can be reliably confirmed by footprint visualizations, which 179 further allow to infer TF binding dynamics. 180 To test if the global footprinting scores of individual TFs correlate with the incidence and level 181 of their RNA expression, we matched them to RNA expression datasets derived from 182 individual timepoints throughout zygotic development, taking TF motif similarity into account. 183 Indeed, we found that TOBIAS scores for the majority of TFs either correlated well with the 184 timing of their expression profiles or displayed a slightly delayed activity after expression 185 peaked (Supp. Figure 4a). This is important because it shows that in conjunction with 186 expression data, TOBIAS can unravel the kinetics between TF expression (mRNA) and the 187 actual binding activity of their translated proteins. The value of this added information becomes 188 particularly apparent when analyzing activities of TFs that did not correlate with the timing of 189 their RNA expression (Supp. Figure 4a;  The timing of ZGA varies between mice (2C) and humans (4C to 8C) (reviewed in 27 28 ). By 203 integrating the TOBIAS scores from human and mouse (Supp. Figure 4b and Supp. Table 3), 204 and instrumentalizing the capability of TOBIAS to generate differential TF binding plots for all 205 time points automatically, we investigated similarities and differences of PD between these 206 species. Firstly, reflecting the shift of ZGA onset, we identified 30 TFs which appeared to be 207 ZGA specific in both human and mouse (Figure 2c) including several homeobox factors which 208 already have described functions within ZGA 29 as well as ARID3A which has been shown to 209 play a role in cell fate decisions in creating trophectoderm 30 . 210 Next, we used the differential TF binding plots to display differences in ZGA at the transition 211 between 2C and 4C in mouse (Figure 3a), and human 8C and ICM (Figure 3b) (Supp. File 2 212 + 3 for all pairwise comparisons). In mice, we observed a shift of Obox-factor activity in 2C to 213 an activation of Tead (Tead1-4) and AP-2 (Tfap2a/c/e) motifs in 4C. Notably, AP-2/Tfap2c is 214 required for normal embryogenesis in mice 31 and was also recently shown to act as a 215 chromatin modifier that opens enhancers proximal to pluripotency factors in human 32 . We 216 observed a similar shift of TF activity for homeobox factors such as PITX1-3, RHOXF1, CRX 217 and DMBX1 at the human 8C stage towards higher scores in ICM for known pluripotency 218 factors such as POU5F1 (OCT4) and other POU-factors. Taken together, these results 219 highlight the ability of TOBIAS to capture differentially bound TFs, not only across the whole 220 timeline, but also between individual conditions and species. 221 Throughout the pairwise comparisons, we observed that TFs from the same families often 222 display similar binding kinetics within species, which is not surprising since they often possess 223 highly similar binding motifs (Figure 3a Table 1+3). Interestingly, 247 despite the fact that Dux has already been proved to play a prominent role in ZGA 5-7, 33, 34 , 248 there is still a poor understanding of how Dux regulates its primary downstream targets, and 249 consequently its secondary targets, during this process. We therefore applied TOBIAS to 250 identify Dux binding sites utilizing an ATAC-seq dataset of Dux overexpression (DuxOE) in 251 mESC 5 . 252 As expected, the differential TF activity predicted by TOBIAS showed an increase in activity   Many of the TOBIAS-predicted Dux targets encode TFs themselves. Therefore, we applied 279 the TOBIAS network module to subset and match all activated binding sites to TF target genes 280 with the aim of inferring how these TF activities might connect. Thereby, we could model an 281 intriguing pseudo timed TF activation network. This directed network uncovered a TF 282 activation cascade initiated by Dux, resulting in the activation of 7 primary TFs which appear 283 to subsequently activate 32 further TFs (first three layers depicted in Figure 4e). As Dux is a 284 regulator of ZGA, we asked how the in vitro activated Dux network compared to gene 285 expression throughout PD in vivo. Strikingly, the in vivo RNA-seq data of the developmental 286 stages 21 confirmed an early 2C specific expression for Dux, followed by a slightly shifted 287 activation pattern for all direct Dux targets except for Rxrg ( Figure 4f). However, it is of note 288 that Rxrg is significantly upregulated in the in vitro DuxOE from which the network is inferred 289 (Supp . Table 5), pointing to both the similarities and differences between the in vivo 2C and 290 in vitro 2C-like stages induced by Dux. In conclusion, these data show that beyond identifying 291 specific target genes of individual TFs, TOBIAS can infer biological insight by predicting entire 292 TF activation networks. 293 Notably, many of the predicted Dux binding sites (40%) are not annotated to genes ( Figure  294 4g), raising the question what role these sites play in ZGA. Dux is known to induce expression 295 of repeat regions such as LTRs 5 and consistently, we found that more than half of the DUX-296 bound sites without annotation to genes are indeed located within known LTR sequences 297 ( Figure 4g) which were transcribed both in vitro and in vivo ( Figure 4h). Interestingly, we 298 additionally found that 28% of all non-annotated Dux binding sites overlap with genomic loci 299 encoding LINE1 elements. Although LINE1 expression does not appear to be altered in mESC 300 cells, there is a striking pattern of increasing LINE1 transcription from 4C-8C ( Figure 4h) in 301 vivo, pointing to a possible role of LINE1 regulation throughout PD. Finally, we found a portion 302 of the Dux binding sites which do not overlap with any annotated gene nor with putative 303 regulatory repeat sequences, even though transcription clearly occurs at these sites ( Figure  304 4h; bottom). One example is a predicted Dux binding site on chromosome 13, which coincides 305 with a spliced region of increased expression between control mESC/DuxOE and comparable 306 high expression in 2C, 4C and 8C (Supp. Figure 6). These data clearly indicate the existence 307 of novel transcribed genetic elements, the function of which remains unknown, but which are 308 likely controlled by Dux and could play a role during PD. 309 In conclusion, TOBIAS predicted the exact locations of Dux binding in promoters of target 310 genes, and could unveil how Dux initiates TF-activation networks and induces expression of 311 repeat regions. Importantly, these data further show that TOBIAS can identify any TFBS with 312 increased binding, not only those limited to annotated genes, which aids in uncovering novel 313 regulatory genetic elements. 314

Footprint scores reveal true characteristics of protein binding 316
To the best of our knowledge, this is the first application of a DGF approach to visualize gain 317 and loss of individual TF footprints in the context of time series, TF overexpression, and TF-318 DNA binding for a wide-range of TFs in parallel. Importantly, we found that these advances 319 could in large part be attributed to the framework approach we took in developing TOBIAS, 320 which enabled us to simultaneously compare global TF binding across samples and quantify 321 changes in TF binding at specific loci. The modularity of the framework also allowed us to 322 apply a multitude of downstream analysis tools to easily visualize footprints and gain even 323 more information about TF binding dynamics as exemplified by the discovery of the Dux TF-324 activation network. 325 The power of this framework to handle time-series data becomes especially apparent when 326 correlating the TOBIAS-based prediction of TF binding to RNA-seq data from the same time 327 points. For instance, TOBIAS could infer when the maternally transferred TF SALL4 is truly 328 active while its gene expression pattern alone does not allow to make such conclusions. Along 329 this line, TOBIAS is also powerful in circumstances where gene expression of a particular TF 330 appears to be anticorrelated with its binding activity. It is tempting to speculate that TFs for 331 which footprinting scores are low, even though their RNA expression is high, might act as 332 transcriptional repressors, because footprinting relies on the premise that TFs will increase 333 chromatin accessibility around the binding site. In support of this hypothesis, recent 334 investigations have suggested that repressors display a decreased footprinting effect in 335 comparison to activators 37 . Therefore, the integration of ATAC-seq footprinting and RNA-seq 336 is an important step in revealing additional information such as classification TFs into 337 repressors and activators, as well as the kinetics between expression and binding. 338

Species-specific TFs use common ZGA motifs in mice and human 339
By integration of human and murine TF activities using both differential footprinting and 340 species-specific TFBS overlaps, our analyses revealed that the majority of TF motifs are active 341 at corresponding timepoints of human and mouse ZGA. This is not necessarily surprising since 342 homologous TFs that exert the same functions usually use similar motifs (e.g Pou2f1/POU2F1, 343 Otx1/OTX1 and/or Foxa3/FOXA3). Interestingly though, we found that this is not the case for 344 all TF motifs. We found that the human RHOXF1 motif (Figure 2b) is likely not utilized by Rhox 345 proteins in mice even though more than 30 Rhox genes exist. Evidently, throughout multiple 346 duplications, Rhox genes seem to have obtained other functionalities in mouse 38 in 347 comparison to the two human RHOX genes that are expressed in reproductive tissues 39 . 348 Therefore, although we found the human RHOXF1 motif to be highly active in mice, this motif 349 is most likely utilized by other proteins such as the mouse specific Obox proteins. In support 350 of this conclusion, expression patterns of Obox proteins appear to be tightly regulated during 351 PD 40 ( 21 ). High expression of Obox 1/2/5/7 is observed from the zygote to 4C stage, while 352 Obox3/6/8 are expressed and peak at later stages (Supp . Table 4). Notably, there is a 353 significant sequence similarity of the homeobox domains but not in the other parts of the 354 RHOXF1 and Obox protein sequences, which supports the similarity in binding specificity. 355 Although the potential functional overlap of RHOXF1 and Obox factors remains unresolved, 356 our inter-species analysis suggests an unappreciated function of these factors and their 357 targets during PD, warranting an in depth investigation. 358 In the context of TF target prediction, the power of TOBIAS was particularly highlighted by the 359 fact that the analysis could identify almost all known Dux targets. In addition to coding genes, 360 our analysis disclosed novel Dux binding sites and significant footprint scores at LINE1 361 encoding genomic loci, which appear to be activated at the 4C/8C stage. This finding is 362 especially interesting because a recent study has shown that LINE1 RNA can interact with 363 Nucleolin and Kap1 to repress Dux expression 41 . Therefore, our findings give rise to a kinetics 364 driven model in which Dux not only initiates ZGA but also regulates its own termination by a 365 temporally delayed negative feedback loop. Exactly how this feedback loop is controlled 366 remains to be determined. 367

Limitations and outlook of footprinting analysis 368
Despite the striking capability of DGF analysis, some limitations and dependencies of this 369 method still remain. Amongst these is the need of high-quality TF motifs for matching footprint 370 scores to individual TFs with high confidence. In other words, while the binding of a TF might 371 create an effect that can be interpreted as a footprint, without a known motif, this effect cannot 372 be matched to the corresponding TF. This becomes evident in the context of DPPA2/4, a TF 373 described by several groups to act in PD and even upstream of Dux 34 . DPPA2/4 targets GC 374 rich sequences 34 , but its canonical binding motif remains unknown. It also needs to be noted 375 that footprinting analysis cannot take effects into account that arise from heterogeneous 376 mixtures of cells wherein TFs are bound in some cells and in others not. Therefore, if not 377 separated, the classification of differential binding will be an observation averaged across 378 many cells, possibly masking subpopulation effects. Recent advances have enabled the 379 application of ATAC-seq in single cells 42 , but this generates sparse matrices, rendering 380 footprinting approaches on single cells elusive. However, we speculate that by creating 381 aggregated pseudo-bulk signals from large clustered SC ATAC datasets, DGF analysis might 382 also become possible in single cells. 383

384
Here, we have illustrated the TOBIAS framework as a versatile tool for ATAC-seq footprinting 385 analysis which helps to unravel transcription factor binding dynamics in complex experimental 386 settings that are otherwise difficult to investigate. We showed that entire networks of TF 387 binding, which have previously been explored using a combination of omics methods, can be 388 recapitulated to a great extent by DGF analysis, which requires only ATAC-seq and TF motifs. 389 From a global perspective, we provided new insights into PD by quantifying the stage-specific 390 activity of specific TFs. Furthermore, we highlighted the usage of TOBIAS to study specific 391 transcription factors as exemplified by our investigations on Dux. Finally, we used the specific 392 TF target predictions to gain insights into the local binding dynamics of Dux in the context of 393 TF-activation networks, repeat regions and novel genetic elements. 394 In conclusion, we present TOBIAS as the first comprehensive software that performs all steps 395 of DGF analysis, natively supports multiple experimental conditions and performs visualization 396 within one single framework. Although we utilized the process of PD as a proof of principle, 397 the modularity and universal nature of the TOBIAS framework enables investigations of 398 various biological conditions beyond PD. We believe that continued work in the field of DGF, 399 including advances in both software and wet-lab methods, will validate this method as a 400 resourceful tool to extend our understanding of a variety of epigenetic processes involving TF 401   to individual repeat elements (as seen in figure 4G) was performed using "Bedtools intersect". 491 The sum of overlaps were counted by repeat class (LINE1/LTR). 492 Visualization 493 All TF-score heatmaps were generated by R Version 3.5.3 and complex heatmap package 494 version 3.6 58 . Individual gene views were generated by loading TOBIAS output tracks into 495 IGV version 2.6.2 59 or using the TOBIAS PlotTracks module, which is a wrapper for the 496 svist4get visualization tool 60 . TF networks were drawn with Cytoscape version 3.7.1 61 . 497 Heatmaps of genomic signal density were generated using Deeptools version 3.3.0 62 . All 498 other figures, such as footprint plots, volcano plots and motif clustering dendrograms were 499 generated by the TOBIAS visualization modules as described below. 500 The TOBIAS framework The term will be negative and will therefore raise the score if there is a high depletion of 552 cuts in the footprint (middle). If there is no depletion, the score will simplify to the mean of cuts 553 in the flanking regions, representing accessibility. It is therefore not necessary to see a 554 canonical footprint shape for the footprint score to be high. The footprint score can be 555 interpreted as higher scores being more evidence that a protein was bound at a given position. 556

(TOBIAS BINDetect module) 558
To match the calculated footprint scores to potential binding sites, TOBIAS BINDetect 559 integrates genomic sequence, footprint scores from several conditions and motifs to identify In the second step of the algorithm, every TF binding site found (for each motif given as input) 578 is split into bound and unbound sites based on a score threshold per condition. The threshold 579 is set at the level of significance of a normal-distribution fit to the background distribution of 580 scores (user-defined p-value). As well as the per-condition split, each site is assigned a 581 log2FC (fold change) per comparison, which represents whether the binding site has 582 larger/smaller footprint scores in comparison. The global distribution of log2FC's per TF is 583 compared to the background distributions to calculate a differential binding score, which is 584 calculated as: 585 586 where and are the means and standard deviations of the observed and 587 background log2FC distributions respectively. A p-value is also calculated by subsampling 588 100 log2FCs from the background and calculating the significance of the observed change 589 (Python's scipy.stats.ttest_1samp). By comparing the observed log2FC distribution to the 590 background log2FC, the effects of any global differences due to sequencing depth, noise etc. 591 are controlled. 592 The differential binding scores and p-values are visualized as a volcano plot per condition-593 comparison. All TFs with -log10(p-value) above the 95% quantile or differential binding scores 594 smaller/larger than the 5% and 95% quantiles (top 5% in each direction) are colored and 595 shown with labels. Below the plot, hierarchical clustering of the TFBS-distance matrix is shown 596 and all TFs with distances less than 0.5 (overlap of 50% of bp) are colored as separate 597

clusters. 598
The result of BINDetect is a folder-structure containing an overview of all potential binding 599 sites (as .bed as well as excel-files), the predicted split into bound and unbound sites, and a 600 global overview of differentially bound TFs per condition-comparison. 601

module) 603
Footprints are visualized using the subtool "TOBIAS PlotAggregate". Aggregate footprints are 604 created by aligning genomic signals centered on all binding sites (taking into account 605 strandedness), to create a matrix of (n sites) x (n bp). The aggregate signal is calculated as 606 the mean of each column (each bp). The default of +/-60bp from the motif center was used 607 throughout this manuscript. 608 The aggregate footprinting depth (FPD), which is applied in Supp. Figure 2c- Figure  2b). 615 Similarly to the investigations in previous literature 16 , we applied a mixture model from the 616 Mixtools R package 67 to estimate the fractions of TFs with/without measurable footprints 617 (Supp. Figure 2e). 618 619 Transcription factor binding network (TOBIAS CreateNetwork module) 620 The TF-TF network for Dux was built by subsetting all binding sites on the following 621 characteristics: Bound in the promoter of a target gene, labeled "Unbound" in Control, labeled 622 "Bound" in DuxOE, and log2FC footprint score increasing for DuxOE vs. Control. All targets 623 were further reduced to only include genes encoding TFs with available motifs. Motifs were 624 matched to genes as explained in the methods section "Processing of transcription factor 625 motifs". The network was then created using "TOBIAS CreateNetwork". The result is a network 626 of source and target nodes with directed edges, which in words can be described as: Source 627 TF binds in the promoter of Target TF. 628

TOBIAS framework output structure 629
The output generated by the TOBIAS framework is organized in a hierarchical folder structure, 630 which increases clarity of all steps of the analysis. The folder structure specifically organizes 631 input data, pre-processing output like peak-calling and annotation, genomic tracks such as 632 bias correction and footprints, as well as the local and global TF predictions. Particularly, the 633 output for every individual TF investigated is arranged into separate folders containing TF 634 specific plots, annotations and binding predictions. This structure makes it simple to use the 635 output for further downstream analysis, as was showcased in this work. An

Comparison of TOBIAS to existing methods 641
Although footprinting tools for DNase-seq exist 68-70 65, 71-73 74 , not all can be applied to paired-642 end ATAC-seq data. We have focused our comparison on tools which are easily obtainable 643 and installable, do not require ChIP-seq training-data, and are explicitly supporting ATAC-seq. 644 We have additionally added two metrics for "Accessibility" and "PWM score" to compare 645 TOBIAS to other footprinting-free metrics. The validation datasets and usage of existing tools 646 are described in the following sections. 647

Datasets 648
The TOBIAS framework was benchmarked using ATAC-seq data from four human cell types: 649 GM12878 (GEO: GSE47753), A549 (GEO: GSE114202), K562 (ENA: PRJNA288801) and 650 HEPG2 (ENA: PRJEB30461). ATAC-seq data was trimmed using cutadapt 44 and mapped 651 using Bowtie2 75 . All reads with a quality score <30 as well as non-proper paired reads were 652 removed. All replicates were merged to one joined .bam-file of reads. Peaks were called using 653 MACS2 47 with parameters "--nomodel --shift -100 --extsize 200 --broad --qvalue 0.01 --broad-654 cutoff 0.01". ChIP-seq peak regions (narrowPeak format) were downloaded from ENCODE 655 and associated to motifs from Jaspar CORE 2018 using "MEME Centrimo" 76 . Only ChIP-seq 656 experiments with motif enrichment > 1.0e-10 (Centrimo E-value) were kept. In case of more 657 than one ChIP-seq experiment for the same target in the same cell type, the one with the 658 highest motif enrichment was chosen. After filtering, there were 12 TFs for A549, 54 TFs for 659 GM12878, 64 TFs for HepG2, and 87 TFs for K562 for a total of 217 ChIP-seq experiments 660 matched to ATAC-seq. Bound binding sites per TF were defined as any TFBS within +/-50bp 661 from the paired ChIP-seq peak summit. In case of two or more binding sites per peak, the one 662 closest to the summit was set to bound, and others were excluded from the analysis. Unbound 663 binding sites were defined as any TFBS not overlapping any ChIP-seq peak, as well as not 664 overlapping bound sites from any other factors for this cell type. Bound and unbound sites 665 were further filtered to only include TFBS falling within ATAC-seq peaks for the cell type in 666 question. 667

Bias correction approaches 668
TOBIAS was compared to the existing bias correction methods as follows: 669 • seqOutBias ( 77 ) 670 The seqOutBias software was downloaded from GitHub 671 (https://github.com/guertinlab/seqOutBias). Following the vignette for ATAC-seq, 672 mappability files were created and ATAC-seq reads were corrected for plus/minus 673 strand reads separately. After correction, we further shifted the positive and negative 674 tracks +5 and -4bp respectively, as this was not performed by the tool itself. 675 • HINT-ATAC ( 14 ) 676 The HINT software was downloaded from PyPI as part of the RGT software suite. Bias-677 correction was performed from the ATAC-seq reads using the command "rgt-hint 678 tracks --bc --bigWig <bam>". 679 680 Aggregate footprints for each method across all (within peaks), bound and unbound binding 681 sites (see explanation above) were visualized using "TOBIAS PlotAggregate". 682 Footprinting 683 Existing footprinting tools were applied as follows: 684 • msCentipede ( 78 ) 685 The msCentipede software was downloaded from GitHub 686 (https://github.com/rajanil/msCentipede). For each TF, the binding model was built 687 using the 5000 TFBS with the highest PWM score genomewide. For model learning, 688 the "--mintol" parameter was set to 1e-3 as a tradeoff between accuracy and speed. 689 The resulting models were then used to infer the posterior binding-probability of TFBS 690 in peaks. 691 • Wellington ( 70 ) 692 The pyDNase software was downloaded from PyPI. Footprints in ATAC-seq peaks 693 were estimated using "wellington_footprints.py" with the "-A" option for ATAC-seq 694 mode. 695 • PIQ ( 79 ) 696 The PIQ software was downloaded from Bitbucket (https://bitbucket.org/thashim/piq-697 single/). The script bam2rdata.r was used to bring the input .bam-file into the correct 698 data format. Likewise, the script pwmmatch.exact.r was used to predict genomewide 699 TFBS. Finally, footprinting scores for each TF were obtained using the script pertf.r for 700 each motif/cell type pair. The purity score was taken as the probability for a certain 701 TFBS to be bound. 702 • HINT-ATAC ( 14 ) 703 The HINT software was downloaded from PyPI as part of the RGT software suite. 704 Footprints were identified using the command "rgt-hint footprinting --atac-seq --paired-705 end --organism=hg38 <bam> <peaks>". The output of HINT-ATAC footprinting is a 706 .bed-file of footprint ranges ranked by tag count. All TFBS overlapping a footprint with 707 more than 2/3 of the TFBS bases was assumed to be bound and scored using the tag 708 count of the footprint. The rest of the TFBS (within peaks) were set to score 0 (low 709 chance of protein binding). The auROC was calculated based on the ability of these 710 scores to predict true protein binding. It should be noted that this affects the shape of 711 the ROC curve, as all TFBS without overlaps are assumed to have the same probability 712 of being bound. However, this is a characteristic of the method, and HINT-ATAC was 713 therefore evaluated on the same premise as other tools. 714 • Accessibility 715 The "Accessibility" metric is defined as the sum of Tn5 insertions in a 300 basepair 716 window centered at the binding site. This score represents the accessibility of a certain 717 region not taking into account local footprint information. 718 • PWM score 719 The score of the motif-sequence match at the specific TFBS. As this is based on 720 sequence alone, the PWM-score is independent of chromatin accessibility. 721 Due to high computational times for some tools, the validation was limited to binding sites on 722 human chromosome 1. On the basis of the ChIP-seq labels, the area under the ROC curve 723