Base-pair Resolution Research Articles

Small-molecule compounds boost genome-editing efficiency of cytosine base editor.

Cytosine base editor (CBE) enables targeted C-to-T conversions at single base-pair resolution and thus has potential therapeutic applications in humans. However, the low efficiency of the system limits practical use of this approach. We reported a high-throughput human cells-based reporter system that can be harnessed for quickly measuring editing activity of CBE. Screening of 1813 small-molecule compounds resulted in the identification of Ricolinostat (an HDAC6 inhibitor) that can enhance the efficiency of BE3 in human cells (2.45- to 9.21-fold improvement). Nexturastat A, another HDAC6 inhibitor, could also increase BE3-mediated gene editing by 2.18- to 9.95-fold. Ricolinostat and Nexturastat A also boost base editing activity of the other CBE variants (BE4max, YE1-BE4max, evoAPOBEC1-BE4max and SpRY-CBE4max, up to 8.32-fold). Meanwhile, combined application of BE3 and Ricolinostat led to >3-fold higher efficiency of correcting a pathogenic mutation in ABCA4 gene related to Stargardt disease in human cells. Moreover, we demonstrated that our strategy could be applied for efficient generation of mouse models through direct zygote injection and base editing in primary human T cells. Our study provides a new strategy to improve the activity and specificity of CBE in human cells. Ricolinostat and Nexturastat A augment the effectiveness and applicability of CBE.

Capillary electrophoresis of DNA with high resolution based on copoly(pentaerythritoltetra succinimidylcarboxypentyl/aminopropyl polyoxyethylene) hydrogel

Capillary gel electrophoresis is widely applied for determination of sequence and size of DNA, in which the sieving gel plays an unignorable role. Herein, a pore-size controllable hydrogel was synthesized in the capillary with two symmetrical tetrahedron-like macromonomers consisting of pentaerythritoltetra (succinimidylcarboxypentyl) polyoxyethylene (PS) and pentaerythritoltetra (aminopropyl) polyoxyethylene) (PA). By capillary electrophoresis of the DNA fragments with this hydrogel, it is found that a homogenous structure of hydrogel which is more suitable for the DNA separation can be achieved when the molecular weight of PA is approximate to that of PS. DNA fragments smaller than 1500 bp can be well resolved in this hydrogel within 13 min. More than 100 consecutive runs can be carried out in such a dynamically coated capillary before performance begins to degrade. Notably, such hydrogel can realize separation of dsDNA up to single base pair resolution and same length of dsDNA with 1 bp difference.

High-Resolution Architecture and Partner Genes of MYC Rearrangements in Lymphoma with DLBCL Morphology

High-Resolution Architecture and Partner Genes of MYC Rearrangements in Lymphoma with DLBCL Morphology

Chemical map-based prediction of nucleosome positioning using the Bioconductor package nuCpos

BackgroundAssessing the nucleosome-forming potential of specific DNA sequences is important for understanding complex chromatin organization. Methods for predicting nucleosome positioning include bioinformatics and biophysical approaches. An advantage of bioinformatics methods, which are based on in vivo nucleosome maps, is the use of natural sequences that may contain previously unknown elements involved in nucleosome positioning in vivo. The accuracy of such prediction attempts reflects the genomic coordinate resolution of the nucleosome maps applied. Nucleosome maps are constructed using micrococcal nuclease digestion followed by high-throughput sequencing (MNase-seq). However, as MNase has a strong preference for A/T-rich sequences, MNase-seq may not be appropriate for this purpose. In addition to MNase-seq-based maps, base pair-resolution chemical maps of in vivo nucleosomes from three different species (budding and fission yeasts, and mice) are currently available. However, these chemical maps have yet to be integrated into publicly available computational methods.ResultsWe developed a Bioconductor package (named nuCpos) to demonstrate the superiority of chemical maps in predicting nucleosome positioning. The accuracy of chemical map-based prediction in rotational settings was higher than that of the previously developed MNase-seq-based approach. With our method, predicted nucleosome occupancy reasonably matched in vivo observations and was not affected by A/T nucleotide frequency. Effects of genetic alterations on nucleosome positioning that had been observed in living yeast cells could also be predicted. nuCpos calculates individual histone binding affinity (HBA) scores for given 147-bp sequences to examine their suitability for nucleosome formation. We also established local HBA as a new parameter to predict nucleosome formation, which was calculated for 13 overlapping nucleosomal DNA subsequences. HBA and local HBA scores for various sequences agreed well with previous in vitro and in vivo studies. Furthermore, our results suggest that nucleosomal subsegments that are disfavored in different rotational settings contribute to the defined positioning of nucleosomes.ConclusionsOur results demonstrate that chemical map-based statistical models are beneficial for studying nucleosomal DNA features. Studies employing nuCpos software can enhance understanding of chromatin regulation and the interpretation of genetic alterations and facilitate the design of artificial sequences.

Defining genome architecture at base-pair resolution.

In higher eukaryotes, many genes are regulated by enhancers that are 104-106base pairs (bp) away from the promoter. Enhancers contain transcription-factor-binding sites (which are typically around 7-22bp), and physical contact between the promoters and enhancers is thought to be required to modulate gene expression. Although chromatin architecture has been mapped extensively at resolutions of 1kilobase and above; it has not been possible to define physical contacts at the scale of the proteins that determine gene expression. Here we define these interactions in detail using a chromosome conformation capture method (Micro-Capture-C) that enables the physical contacts between different classes of regulatory elements to be determined at base-pair resolution. We find that highly punctate contacts occur between enhancers, promoters and CCCTC-binding factor (CTCF) sites and we show that transcription factors have an important role in the maintenance of the contacts between enhancers and promoters. Our data show that interactions between CTCF sites are increased when active promoters and enhancers are located within the intervening chromatin. This supports a model in which chromatin loop extrusion1 is dependent on cohesin loading at active promoters and enhancers, which explains the formation of tissue-specific chromatin domains without changes in CTCF binding.

Independent variations in genome-wide expression, alternative splicing, and DNA methylation in brain tissues among castes of the buff-tailed bumblebee, Bombus terrestris

Caste differentiation in social hymenopterans is an intriguing example of phenotypic plasticity. However, the co-ordination among gene regulatory factors to mediate caste differentiation remains inconclusive. In this study, we determined the role of gene regulation and related epigenetic processes in pre-imaginal caste differentiation in the primitively eusocial bumblebee Bombus terrestris. By combining RNA-Seq data from Illumina and PacBio and accurately quantifying methylation at whole-genomic base pair resolution, we found that queens, workers, and drones mainly differentiate in gene expression but not in alternative splicing and DNA methylation. Gynes are the most distinct with the lowest global level of whole-genomic methylation and with the largest number of caste-specific transcripts and alternative splicing events. By contrast, workers exhibit few uniquely expressed or alternatively spliced genes. Moreover, several genes involved in hormone and neurotransmitter metabolism are related to caste differentiation, whereas several neuropeptides are linked with sex differentiation. Despite little genome-wide association among differential gene expression, splicing, and differential DNA methylation, the overlapped gene ontology (GO) terms point to nutrition-related activity. Therefore, variations in gene regulation correlate with the behavioral differences among castes and highlight the specialization of toolkit genes in bumblebee gynes at the beginning of the adult stage.

Nonequilibrium thermodynamics of the RNA-RNA interaction underlying a genetic transposition program.

Thermodynamic descriptions are powerful tools to formally study complex gene expression programs evolved in living cells on the basis of macromolecular interactions. While transcriptional regulations are often modeled in the equilibrium, other interactions that occur in the cell follow a more complex pattern. Here, we adopt a nonequilibrium thermodynamic scheme to explain the RNA-RNA interaction underlying IS10 transposition. We determine the energy landscape associated with such an interaction at the base-pair resolution, and we present an original scaling law for expression prediction that depends on different free energies characterizing that landscape. Then, we show that massive experimental data of the IS10 RNA-controlled expression are better explained by this thermodynamic description in nonequilibrium. Overall, these results contribute to better comprehend the kinetics of post-transcriptional regulations and, more broadly, the functional consequences of processes out of the equilibrium in biology.

Accurate modeling of DNA conformational flexibility by a multivariate Ising model

The sequence-dependent structure and deformability of DNA play a major role for binding of proteins and regulation of gene expression. So far, most efforts to model DNA flexibility are based on unimodal harmonic stiffness models at base-pair resolution. However, multimodal behavior due to distinct conformational substates also contributes significantly to the conformational flexibility of DNA. Moreover, these local substates are correlated to their nearest-neighbor substates. A description for DNA elasticity which includes both multimodality and nearest-neighbor coupling has remained a challenge, which we solve by combining our multivariate harmonic approximation with an Ising model for the substates. In a series of applications to DNA fluctuations and protein-DNA complexes, we demonstrate substantial improvements over the unimodal stiffness model. Furthermore, our multivariate Ising model reveals a mechanical destabilization for adenine (A)-tracts to undergo nucleosome formation. Our approach offers a wide range of applications to determine sequence-dependent deformation energies of DNA and to investigate indirect readout contributions to protein-DNA recognition.

Genome-wide analysis of DNA replication and DNA double-strand breaks using TrAEL-seq.

Faithful replication of the entire genome requires replication forks to copy large contiguous tracts of DNA, and sites of persistent replication fork stalling present a major threat to genome stability. Understanding the distribution of sites at which replication forks stall, and the ensuing fork processing events, requires genome-wide methods that profile replication fork position and the formation of recombinogenic DNA ends. Here, we describe Transferase-Activated End Ligation sequencing (TrAEL-seq), a method that captures single-stranded DNA 3' ends genome-wide and with base pair resolution. TrAEL-seq labels both DNA breaks and replication forks, providing genome-wide maps of replication fork progression and fork stalling sites in yeast and mammalian cells. Replication maps are similar to those obtained by Okazaki fragment sequencing; however, TrAEL-seq is performed on asynchronous populations of wild-type cells without incorporation of labels, cell sorting, or biochemical purification of replication intermediates, rendering TrAEL-seq far simpler and more widely applicable than existing replication fork direction profiling methods. The specificity of TrAEL-seq for DNA 3' ends also allows accurate detection of double-strand break sites after the initiation of DNA end resection, which we demonstrate by genome-wide mapping of meiotic double-strand break hotspots in a dmc1Δ mutant that is competent for end resection but not strand invasion. Overall, TrAEL-seq provides a flexible and robust methodology with high sensitivity and resolution for studying DNA replication and repair, which will be of significant use in determining mechanisms of genome instability.

Deep learning-based enhancement of epigenomics data with AtacWorks

ATAC-seq is a widely-applied assay used to measure genome-wide chromatin accessibility; however, its ability to detect active regulatory regions can depend on the depth of sequencing coverage and the signal-to-noise ratio. Here we introduce AtacWorks, a deep learning toolkit to denoise sequencing coverage and identify regulatory peaks at base-pair resolution from low cell count, low-coverage, or low-quality ATAC-seq data. Models trained by AtacWorks can detect peaks from cell types not seen in the training data, and are generalizable across diverse sample preparations and experimental platforms. We demonstrate that AtacWorks enhances the sensitivity of single-cell experiments by producing results on par with those of conventional methods using ~10 times as many cells, and further show that this framework can be adapted to enable cross-modality inference of protein-DNA interactions. Finally, we establish that AtacWorks can enable new biological discoveries by identifying active regulatory regions associated with lineage priming in rare subpopulations of hematopoietic stem cells.

Abstract PR-08: Accurate quantification of tumor DNA in liquid biopsies using deep learning

Abstract Circulating tumor DNA (ctDNA) fraction has been shown to be a prognostic factor of treatment response in patients with metastatic cancers. However, algorithms that accurately determine ctDNA fraction across purity ranges are still lacking. Commonly used liquid biopsy gene panels (~70 genes) often approximate the ctDNA fraction in a sample by using somatic variant allele frequencies, resulting in large confidence intervals due to the scarcity of somatic mutations detected per sample and uncertainty of their clonality and copy number. Genome wide signatures such as nucleosome positioning, measured by windowed protection scores (WPS), have been shown to reflect genomic and epigenomic abnormalities specific to cancer and proposed as biomarkers for cancer detection. However, genome wide assays sacrifice somatic mutation detection due to low sequencing coverage and in clinical practice a second assay is required to find actionable mutations. We leverage data from ~6,000 cell-free DNA samples to develop methods for accurately estimating the ctDNA fraction in clinical samples sequenced with a 554-gene panel covering 3MB of genomic content. First, we use Gradient Boosting and 17 features extracted from standard variant calling, copy number and QC pipelines. This approach has been shown to be highly accurate in samples of high ctDNA fraction and therefore provide reliable training labels for those samples. We suggest that more sophisticated methods, using genome wide signatures, can increase prediction accuracies in samples of lower ctDNA fraction. Second, we develop a deep learning model for predicting ctDNA fraction using base pair resolution WPS over the genomic territory covered by the gene panel. Using Morlet wavelet, we transformed each autosomal WPS unidimensional representation into a two-dimensional image. A 22 channels image per sample was generated and used as input to a Convolutional Neural Network (22-CNN) model. The 22-CNN was trained to predict the ctDNA fraction and to classify if a sample was from a ctDNA shedder (ctDNA fraction >0) or not. The dataset was equally divided into training and test sets of ~2800 samples each. We used five-fold cross-validation (CV) during training and evaluated each model in the unseen test dataset. The Mean Square Error of the 22-CNN ctDNA prediction model on the testing set was 0.0072 +- 0.0013 [0.006-0.009] and the accuracy of the classification model on the testing set was 0.79 +- 0.006 [0.78-0.79]. In this work, we present both machine learning and deep learning models to reliably estimate the ctDNA fraction in plasma samples from cancer patients. These methods address a key challenge in the field and are applicable to other liquid biopsy assays. Citation Format: Aamna M. Al-Shehhi, Chintan Parmar, Ben Gustafson, Angad Singh, Rebecca Leary, Markus Riester, O. Alejandro Balbin. Accurate quantification of tumor DNA in liquid biopsies using deep learning [abstract]. In: Proceedings of the AACR Virtual Special Conference on Artificial Intelligence, Diagnosis, and Imaging; 2021 Jan 13-14. Philadelphia (PA): AACR; Clin Cancer Res 2021;27(5_Suppl):Abstract nr PR-08.

Base-resolution models of transcription-factor binding reveal soft motif syntax.

The arrangement of transcription factor (TF) binding motifs (syntax) is an important part of the cis-regulatory code, yet remains elusive. We introduce a deep learning model, BPNet, that uses DNA sequence to predict base-resolution ChIP-nexus binding profiles of pluripotency TFs. We develop interpretation tools to learn predictive motif representations and identify soft syntax rules for cooperative TF binding interactions. Strikingly, Nanog preferentially binds with helical periodicity, and TFs often cooperate in a directional manner, which we validate using CRISPR-induced point mutations. Our model represents a powerful general approach to uncover the motifs and syntax of cis-regulatory sequences in genomics data.

Stability and folding pathways of tetra-nucleosome from six-dimensional free energy surface

The three-dimensional organization of chromatin is expected to play critical roles in regulating genome functions. High-resolution characterization of its structure and dynamics could improve our understanding of gene regulation mechanisms but has remained challenging. Using a near-atomistic model that preserves the chemical specificity of protein-DNA interactions at residue and base-pair resolution, we studied the stability and folding pathways of a tetra-nucleosome. Dynamical simulations performed with an advanced sampling technique uncovered multiple pathways that connect open chromatin configurations with the zigzag crystal structure. Intermediate states along the simulated folding pathways resemble chromatin configurations reported from in situ experiments. We further determined a six-dimensional free energy surface as a function of the inter-nucleosome distances via a deep learning approach. The zigzag structure can indeed be seen as the global minimum of the surface. However, it is not favored by a significant amount relative to the partially unfolded, in situ configurations. Chemical perturbations such as histone H4 tail acetylation and thermal fluctuations can further tilt the energetic balance to stabilize intermediate states. Our study provides insight into the connection between various reported chromatin configurations and has implications on the in situ relevance of the 30 nm fiber.

Base-pair resolution analysis of the effect of supercoiling on DNA flexibility and major groove recognition by triplex-forming oligonucleotides

In the cell, DNA is arranged into highly-organised and topologically-constrained (supercoiled) structures. It remains unclear how this supercoiling affects the detailed double-helical structure of DNA, largely because of limitations in spatial resolution of the available biophysical tools. Here, we overcome these limitations, by a combination of atomic force microscopy (AFM) and atomistic molecular dynamics (MD) simulations, to resolve structures of negatively-supercoiled DNA minicircles at base-pair resolution. We observe that negative superhelical stress induces local variation in the canonical B-form DNA structure by introducing kinks and defects that affect global minicircle structure and flexibility. We probe how these local and global conformational changes affect DNA interactions through the binding of triplex-forming oligonucleotides to DNA minicircles. We show that the energetics of triplex formation is governed by a delicate balance between electrostatics and bonding interactions. Our results provide mechanistic insight into how DNA supercoiling can affect molecular recognition, that may have broader implications for DNA interactions with other molecular species.

Whole-genome methylation analysis reveals novel epigenetic perturbations of congenital scoliosis

Congenital scoliosis (CS) is a congenital disease caused by malformations of vertebrae. Recent studies demonstrated that DNA modification could contribute to the pathogenesis of disease. This study aims to identify epigenetic perturbations that may contribute to the pathogenesis of CS. Four CS patients with hemivertebra were enrolled and underwent spine correction operations. DNA was extracted from the hemivertebrae and spinal process collected from the specimen during the hemivertebra resection. Genome-wide DNA methylation profiling was examined at base-pair resolution using whole-genome bisulfite sequencing (WGBS). We identified 343 genes with hyper-differentially methylated regions (DMRs) and 222 genes with hypo-DMRs, respectively. These genes were enriched in the mitogen-activated protein kinase (MAPK) signaling pathway, calcium signaling pathway, and axon guidance in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and were enriched in positive regulation of cell morphogenesis involved in differentiation, regulation of cell morphogenesis involved in differentiation, and regulation of neuron projection development in Biological Process of Gene Ontology (GO-BP) terms. Hyper-DMR-related genes, including IGHG1, IGHM, IGHG3, RNF213, and GSE1, and hypo DMR-related genes, including SORCS2, COL5A1, GRID1, RGS3, and ROBO2, may contribute to the pathogenesis of hemivertebra. The aberrant DNA methylation may be associated with the formation of hemivertebra and congenital scoliosis.

DNA methylome profiling of circulating tumor cells in lung cancer at single base-pair resolution

DNA methylation plays a pivotal role in regulating cellular processes, and altered DNA methylation pattern is a general hallmark of cancer. However, DNA methylome in circulating tumor cells (CTCs) is still a mystery due to the lack of proper analytical techniques. We introduced an efficient workflow, LCM–µWGBS, which can efficiently profile the DNA methylation of microdissected CTC samples. LCM–µWGBS combines the laser capture microdissection (LCM)-based CTC capture method and whole-genome bisulfite sequencing in very small CTC population (µWGBS) to gain insight into the DNA methylation landscape of CTCs. We herein profiled the DNA methylome of CTCs from lung cancer patients. Deriving from a comprehensive analysis of CTC methylome, a unique “CTC DNA methylation signature” that is distinct from primary lung cancer tissues was identified. Further analysis showed that promoter hypermethylation of epithelial genes is a hallmark of stable epithelial–mesenchymal transition process. Moreover, it has been suggested that CTCs are endowed with a stemness-related feature during dissemination and metastasis. This work constitutes a unique DNA methylation analysis of CTCs at single base-pair resolution, which might facilitate to propose noninvasive CTC DNA methylation biomarkers contributing to clinical diagnosis.

Structural and Genetic Determinants of Convergence in the Drosophila tRNA Structure\u2013Function Map

The evolution of tRNA multigene families remains poorly understood, exhibiting unusual phenomena such as functional conversions of tRNA genes through anticodon shift substitutions. We improved FlyBase tRNA gene annotations from twelve Drosophila species, incorporating previously identified ortholog sets to compare substitution rates across tRNA bodies at single-site and base-pair resolution. All rapidly evolving sites fell within the same metal ion-binding pocket that lies at the interface of the two major stacked helical domains. We applied our tRNA Structure–Function Mapper (tSFM) method independently to each Drosophila species and one outgroup species Musca domestica and found that, although predicted tRNA structure–function maps are generally highly conserved in flies, one tRNA Class-Informative Feature (CIF) within the rapidly evolving ion-binding pocket—Cytosine 17 (C17), ancestrally informative for lysylation identity—independently gained asparaginylation identity and substituted in parallel across tRNAAsn paralogs at least once, possibly multiple times, during evolution of the genus. In D. melanogaster, most tRNALys and tRNAAsn genes are co-arrayed in one large heterologous gene cluster, suggesting that heterologous gene conversion as well as structural similarities of tRNA-binding interfaces in the closely related asparaginyl-tRNA synthetase (AsnRS) and lysyl-tRNA synthetase (LysRS) proteins may have played a role in these changes. A previously identified Asn-to-Lys anticodon shift substitution in D. ananassae may have arisen to compensate for the convergent and parallel gains of C17 in tRNAAsn paralogs in that lineage. Our results underscore the functional and evolutionary relevance of our tRNA structure–function map predictions and illuminate multiple genomic and structural factors contributing to rapid, parallel and compensatory evolution of tRNA multigene families.

Detection of genetic variation and base modifications at base-pair resolution on both DNA and RNA

Accurate decoding of nucleic acid variation is critical to understand the complexity and regulation of genome function. Here we use a single-molecule magnetic tweezer (MT) platform to identify sequence variation and map a range of important epigenetic base modifications with high sensitivity, specificity, and precision in the same single molecules of DNA or RNA. We have also developed a highly specific amplification-free CRISPR-Cas enrichment strategy to isolate genomic regions from native DNA. We demonstrate enrichment of DNA from both E. coli and the FMR1 5’UTR coming from cells derived from a Fragile X carrier. From these kilobase-length enriched molecules we could characterize the differential levels of adenine and cytosine base modifications on E. coli, and the repeat expansion length and methylation status of FMR1. Together these results demonstrate that our platform can detect a variety of genetic, epigenetic, and base modification changes concomitantly within the same single molecules.

Anopheles gambiae Genome Conservation as a Resource for Rational Gene Drive Target Site Selection.

Simple SummaryMalaria is a huge public health burden that affects predominantly sub-Saharan Africa and is transmitted by Anopheles mosquitoes. As a measure for population control, a method called gene drive has been recently developed, which relies on genetic engineering to introduce specific genetic traits into mosquito populations. Gene drives are designed to insert at specific target sites in the mosquito genome. The efficacy of gene drives greatly depends on the selection of appropriate target sites that are functionally or structurally constrained and less likely to tolerate mutations that can hinder the spread of the desired trait in the population. The aim of this study was to perform a genome-wide analysis of highly conserved genomic regions in Anopheles gambiae and introduce a measure of conservation that could indicate sites of functional or structural constraint. The results of this analysis are gathered in a publicly available dataset that can support gene drive target selection and can offer further insights in the nature of conserved genomic regions.The increase in molecular tools for the genetic engineering of insect pests and disease vectors, such as Anopheles mosquitoes that transmit malaria, has led to an unprecedented investigation of the genomic landscape of these organisms. The understanding of genome variability in wild mosquito populations is of primary importance for vector control strategies. This is particularly the case for gene drive systems, which look to introduce genetic traits into a population by targeting specific genomic regions. Gene drive targets with functional or structural constraints are highly desirable as they are less likely to tolerate mutations that prevent targeting by the gene drive and consequent failure of the technology. In this study we describe a bioinformatic pipeline that allows the analysis of whole genome data for the identification of highly conserved regions that can point at potential functional or structural constraints. The analysis was conducted across the genomes of 22 insect species separated by more than hundred million years of evolution and includes the observed genomic variation within field caught samples of Anopheles gambiae and Anopheles coluzzii, the two most dominant malaria vectors. This study offers insight into the level of conservation at a genome-wide scale as well as at per base-pair resolution. The results of this analysis are gathered in a data storage system that allows for flexible extraction and bioinformatic manipulation. Furthermore, it represents a valuable resource that could provide insight into population structure and dynamics of the species in the complex and benefit the development and implementation of genetic strategies to tackle malaria.

Deep neural networks identify sequence context features predictive of transcription factor binding.

Transcription factors (TFs) bind DNA by recognizing specific sequence motifs, typically of length 6–12bp. A motif can occur many thousands of times in the human genome, but only a subset of those sites are actually bound. Here we present a machine learning framework leveraging existing convolutional neural network architectures and model interpretation techniques to identify and interpret sequence context features most important for predicting whether a particular motif instance will be bound. We apply our framework to predict binding at motifs for 38 TFs in a lymphoblastoid cell line, score the importance of context sequences at base-pair resolution, and characterize context features most predictive of binding. We find that the choice of training data heavily influences classification accuracy and the relative importance of features such as open chromatin. Overall, our framework enables novel insights into features predictive of TF binding and is likely to inform future deep learning applications to interpret non-coding genetic variants.