Article https://doi.org/10.

1038/s41467-024-45191-5

Predicting DNA structure using a deep

learning method
1 1 1,2,3,4
Received: 25 April 2023 Jinsen Li , Tsu-Pei Chiu & Remo Rohs

Accepted: 17 January 2024

Understanding the mechanisms of protein-DNA binding is critical in compre-

hending gene regulation. Three-dimensional DNA structure, also described as
Check for updates
DNA shape, plays a key role in these mechanisms. In this study, we present a
deep learning-based method, Deep DNAshape, that fundamentally changes
1234567890():,;
1234567890():,;

the current k-mer based high-throughput prediction of DNA shape features by

accurately accounting for the influence of extended flanking regions, without
the need for extensive molecular simulations or structural biology experi-
ments. By using the Deep DNAshape method, DNA structural features can be
predicted for any length and number of DNA sequences in a high-throughput
manner, providing an understanding of the effects of flanking regions on DNA
structure in a target region of a sequence. The Deep DNAshape method pro-
vides access to the influence of distant flanking regions on a region of interest.
Our findings reveal that DNA shape readout mechanisms of a core target are
quantitatively affected by flanking regions, including extended flanking
regions, providing valuable insights into the detailed structural readout
mechanisms of protein-DNA binding. Furthermore, when incorporated in
machine learning models, the features generated by Deep DNAshape improve
the model prediction accuracy. Collectively, Deep DNAshape can serve as
versatile and powerful tool for diverse DNA structure-related studies.

Binding interactions between DNA binding proteins such as tran- homeodomain family utilize shape readout in specific DNA regions
scription factors (TFs) and their DNA target sites are crucial for gene within the core binding site5–9. Narrower minor groove regions can
regulation; thus, it is vital to fully understand TF-DNA binding result from diverse DNA sequences, providing an additional signal for
mechanisms. These mechanisms can be viewed from two perspectives: binding specificity. Another example is the myocyte enhancer factor-2
base readout and shape readout1,2. Base readout occurs through direct (MEF2) family, which uses shape readout in the center of its degenerate
contacts between amino acids and DNA bases, and a relatively strict motif 10. In general, TFs use a combination of base and shape readout
pattern is followed. For instance, the zinc finger TF family heavily relies modes to achieve DNA binding specificity. The mismatch-induced DNA
on base readout to recognize DNA targets3,4. Shape readout happens structural alterations that influence TF-DNA binding affinity prove the
when proteins interact with the double helical conformation, rather importance of shape readout11. However, the influence of the base and
than interacting directly through hydrogen bonds with certain shape readout mechanisms can be difficult to parse.
atoms on nucleobases. As an example, proteins containing positively The influence of shape readout in DNA binding mechanisms
charged amino acids (e.g., arginine, protonated histidine, and lysine) extends to various aspects of DNA structural features, also known as
may favor a certain DNA structural feature, such as a narrower minor DNA shape features. To measure shape readout quantitatively, one
groove that is more negatively charged1. Some TFs from the must first define DNA shape. This article focuses on the B-form

1
Department of Quantitative and Computational Biology, University of Southern California, Los Angeles, CA 90089, USA. 2Department of Chemistry,
University of Southern California, Los Angeles, CA 90089, USA. 3Department of Physics and Astronomy, University of Southern California, Los Angeles, CA
90089, USA. 4Thomas Lord Department of Computer Science, University of Southern California, Los Angeles, CA 90089, USA. e-mail: [email protected]

Nature Communications | (2024)15:1243 1

Article https://doi.org/10.1038/s41467-024-45191-5

representation of DNA molecules - the canonical right-handed double highlights the limitations of the pentamer query table. Although it is
helix. Three-dimensional (3D) DNA shape features describe the global possible to generate a query table for longer k-mers, the number of
and local relationships between nucleobases, base pairs, and sugar- simulations that would be necessary to cover all longer k-mers is
phosphate backbones12. Although DNA oligonucleotides are flexible exponentially higher; meanwhile, there are not enough existing
molecules with varying DNA shape features, certain optimal con- experimental structures28. Therefore, there is a need for a method that
formations may be intrinsically favorable, as seen in their static or can accurately predict DNA shape features in a high-throughput
sometimes equivalent time-averaged representation of DNA shape manner using only limited data while considering longer-range effects.
features13. Proteins that favor certain intrinsic DNA shape features may To develop such a model, we began with assumptions about how
bind DNA targets with other DNA structures; however, they may the 3D DNA structure is affected by sequence. Firstly, we assumed that
require a higher energy to maintain the binding, thus exhibiting a DNA shape features at each bp are mainly influenced by their neigh-
lower binding affinity6. boring bp, and that this influence weakens with distance. Secondly, we
The energy required to alter an unfavored static DNA structure assumed that this influence can be statistically inferred. Therefore, we
may depend on the DNA flexibility. Therefore, it is important to con- designed a specialized deep learning architecture to deal with variable-
sider the conformational flexibility of the DNA and fluctuations in its length DNA sequences and computed the neighboring effects of
shape features. As revealed by, for instance, molecular dynamics (MD) flanking regions in a layer-by-layer manner (Fig. 1a–d, Supplementary
simulations, DNA shape features fluctuate in various ways14–16. For Fig. 1). We trained the model on DNA shape features that were pre-
example, CpG, CpA and TpA base-pair (bp) steps are generally more viously analyzed and compiled from MC simulations (Fig. 1b, Supple-
flexible17, whereas A-tracts consisting of consecutive ApA, ApT, or TpT mentary Figs. 2–3) and which had been experimentally validated27. The
bp steps are conformationally rigid17,18. The length of A-tracts has also model now considers longer range neighboring effects compared to
been observed to impact the flexibility of neighboring regions19. current data source limitations (Fig. 1e) and can be used for predicting
Methylated cytosine induces DNA flexibility changes that in turn DNA shape features from any given sequence (Fig. 1f, g). We then
influence nucleosome stability20,21. Another recent study using MD evaluated the predicted DNA shape features with a tetramer query
simulations showed that conformational flexibility contributed by the table derived from MD simulations (Supplementary Table 1)28. We
flanks plays an important role in homeodomain binding22,23. However, compared predictions from our resulting model, Deep DNAshape, to
no high-throughput methods have been proposed yet to predict DNA predictions from our original pentamer-based DNAshape (DNA-
shape fluctuations. shapeR) method39 (Figs. 2–4, Supplementary Table 2). To thoroughly
Methods such as MD or Monte-Carlo (MC) simulations and benchmark the model, we also trained it on alternative DNA shape
X-ray crystallography (XRC) can be used to acquire DNA structures sourced from experimentally solved structures40 and MD simulations41
for short DNA fragments. Methods such as Curves12,24,25 and 3DNA26 (see Supplementary Information), leading to the generation of Deep
have been developed to derive DNA shape features from computa- DNAshape (Expt) and Deep DNAshape (MD). These models were then
tional trajectories or experimentally solved structures. High- compared with the MC-trained Deep DNAshape method, along with
throughput methods such as DNAshape27–29 have been introduced comparative analysis across all model variants (Supplementary Fig. 4,
to circumvent the difficult and sometimes impossible task of simu- Supplementary Table 3).
lating or experimentally solving structures. Numerous studies have In addition, the design of the model unlocks the potential to
successfully employed the DNAshape method27, demonstrating the examine systematically how DNA shape fluctuations are affected by
effectiveness of using high-throughput methods to predict DNA extended flanking regions. We previously generated a query table
shape features30–34. containing standard deviation (SD) values for 13 DNA shape features
Nevertheless, although DNAshape27–29 is an efficient method, it and used them in a machine learning study28. Although these values
relies entirely on a pentamer query table containing all possible pen- were statistically computed, nevertheless, the model performance was
tamers compiled from extensive MC simulations35. The pentamer significantly improved when the values were included28. We assumed
length limits this method because only the nearest and next-nearest that these SD values were highly correlated with true fluctuation
neighbors are accounted for when considering the influence of the values. Therefore, although conformational flexibility of DNA is
sequence environment on the center of the pentamer. As for other important42, it is the fluctuation of shape features that is difficult to
data sources, such as MD simulation data and experimental structures access and frequently overlooked in research. Here, we used the same
in the Protein Data Bank (PDB)36, only tetramer query tables could be approach as in predicting the static DNA shape to predict DNA shape
generated due to data paucity28. Therefore, effects from longer-range fluctuation (FL) values: specifically, we directly calculated fluctuation
neighbors are totally neglected in this query table setup. from MC simulations using Curves12. We investigated whether these
Our approach, Deep DNAshape, overcomes the limitation of high-throughput-predicted fluctuation values aligned with previous
DNAshape, particularly its reliance on the query table search key. This findings on DNA flexibility42, and we considered the insights that these
advancement is pivotal, given that the limitation was only caused by fluctuation values might provide.
the available amount of data. Deep DNAshape enhances the capability Next, we compared the Deep DNAshape model with the original
to discern how the shape at the center of a pentamer region is influ- pentamer-based DNAshape method. We tested our model on data
enced by its extended flanking regions, providing a model that offers a from TF-DNA binding assays for quantifying the relative binding affi-
more accurate representation of DNA. nity of DNA sequences for any given TF. These binding assays con-
The pentamer query table contains all possible pentamers for sisted of multiple in vitro experimental methods, such as protein
DNA shape features of the central bp, considering up to a 2-bp flanking binding microarray (PBM)43, HT-SELEX44 and SELEX-seq7. We pre-
region, which consists of the nearest-neighbors and next-nearest- viously used these datasets in conjunction with an expanded set of 13
neighbors. Despite this definition, flanking regions exceeding two bp DNA shape features including groove features, inter-bp features, and
may still be influential, which has been shown for the central TpA intra-bp features (see Methods), demonstrating the effectiveness of
step37. For some TFs having a long core motif, a more complete view of DNA shape features in predicting TF-DNA binding specificities using
DNA shape considering longer-range flanking regions is necessary. In L2-regularized multiple linear regression28. We tested DNA shape
addition, some DNA shape parameters are bimodal15, and statistically predicted by Deep DNAshape against the same datasets to determine if
calculating an average value for the query table may not capture the improvements could be detected compared to the original pentamer
whole picture. The ability to approximate DNA shape features using DNAshape model (Fig. 5). Finally, we showed the potential of Deep
only sequence information such as mono- and dinucleotides38 also DNAshape by processing large genomic-level data.

Nature Communications | (2024)15:1243 2

Article https://doi.org/10.1038/s41467-024-45191-5

...
er
a b e f DNA sequence

Required k-mer occurance

lay 5
ape 4
Training structural data from simulation/experiment Sh

for query table

3 Deeper
Deep G
G
G
G
A
T
G
G
A
T
T
C
C
A
A
T
T
T
T
A
A
A
A
T
T C A T T T C C
C A T T T T C C
Deep
DNAshape
model
2 structural insights
DNAshape C G C G A C G C C C G C G Deep
MC DNAshape
A T C T G G A A T G A T T T A C G A C T T T C C A XRC
model MD
model
2Å Minor groove width (MGW) 7Å

=
c DNA sequence Layer
Self convolution
g DNA shape DNA shape fluctuation
Dropout and average Minor groove
Low High
EP
DNA shape MGW
Shape layer 1
Minor groove width
(single nucleotide)
Intra-base pair
y 
Shear Low High
Shape layer 2
Stretch
(up to 1-bp of flanking influence) x
Stagger
Buckle
ProT 
(up to 2-bp of flanking influence) y Low High
Opening
Shape layer k
x
Inter-base pair
(up to k-bp of flanking influence) Shift

d Slide
y

3’ 5’ Rise Low High
5’ 3’ Tilt
Roll x

HelT
Shape layer Collect and compute neighboring effect

Fig. 1 | Deep DNAshape schematic. Deep DNAshape model (a) and training data Deep DNAshape model and current data limitations. The shown labels ‘MC’, ‘XRC’
for the model (b). Sequences used in simulations are pooled together, and corre- and ‘MD’ are current limitations to generate pentamer or tetramer query tables
sponding DNA shape values are assigned to each position of each sequence. Some sourced from Monte Carlo simulations, X-ray crystallography, and molecular
values, such as MGW at the terminal bp, are not defined and therefore excluded. dynamics simulations, respectively. f Deep DNAshape model usage. Deep DNA-
c Schematic representation of the Deep DNAshape model. DNA sequences go shape can process a given DNA sequence as a string of characters (A, C, G, and T)
through several ‘shape layers’. Each time the model goes through a shape layer, and predict any specific DNA shape feature for each nucleotide position of a
features for each individual position are passed to its nearby two positions, sequence. g All DNA shape predicted by the Deep DNAshape model. In addition to
allowing one additional consideration for the flanking region. d Simple diagram of static DNA shape features, the model can predict DNA shape fluctuations (Cartoons
‘shape layer’. In each layer, features from neighboring nodes are collected and do not represent real values of low or high. Values can be negative or positive for
computed with information in each current node. Features are updated for each angular shape parameters.). Shown are graphical explanations of the four most
node on the sequence. See Methods for details. e Capacity comparison between the frequently-used DNA shape features (MGW, ProT, Roll, and HelT).

Results We validated our predictions against tetramers compiled from

Deep DNAshape predicts DNA shape and shape fluctuations MD simulations15 (Supplementary Table 1) and our previous pentamer
considering extended flanking influences without biases query table (Supplementary Table 2). Additionally, we validated our
MC simulations were used to generate 3D structural predictions for predictions by calculating and comparing the average inter-bp fea-
numerous DNA samples with analyzed DNA shape values for 2121 tures for 10 dinucleotides and intra-bp features for A-T and C-G bp,
different sequences27. For each individual DNA shape feature, a considering all possible flanking regions predicted by the Deep DNA-
training file containing all sequences and their sequence–position- shape model and its variants using different data sources (Supple-
wise DNA shape values was compiled from the simulation data. mentary Fig. 4, Supplementary Table 3; also refer to Supplementary
Despite the varying lengths of these sequences, the Deep DNAshape Text). This approach eliminated the use of a query table and enabled us
model was designed to accommodate such variation (Fig. 1a–d, to predict DNA shape features affected by longer flanking regions,
Supplementary Fig. 1). Following hyperparameter searches on compared to using a forcibly generated hexamer or heptamer query
training and validation data (Supplementary Figs. 2, 3), each model table with large numbers of missing values (Fig. 2a–d, Supplementary
was trained on the entire training data and, when used together, was Figs. 5–7, Supplementary Fig. 8 for Deep DNAshape (MD)). Using our
able to predict any DNA shape feature for any length of DNA Deep DNAshape method, the inferred DNA shape values based on
sequence. The predictions consider only the local bp information of effects of extended flanking regions (Fig. 2a–d, Supplementary
nearby bp (Supplementary Fig. 1). The models can predict DNA Figs. 5–7) are almost perfectly aligned (Supplementary Table 2) with
shape features by considering up to 7-bp flanking regions, and the statistically computed values. Compared to an interpolated query
number of bp of flanking regions considered can be selected for table, our method corrects the biases from different distributions of k-
optimal accuracy (Fig. 1e–g). Models were trained with minimal mers and artifacts from molecular simulations or experiments using
overfitting, while maintaining high accuracy in deeper layers, as DNA fragments (Supplementary Fig. 9).
evident from the training and validation split samples in the By performing high-throughput prediction of DNA shape features
hyperparameter searches (Supplementary Figs. 2, 3). considering extended flanking regions, Deep DNAshape can be used to

Nature Communications | (2024)15:1243 3

Article https://doi.org/10.1038/s41467-024-45191-5

a c e

n=16,384

MGW (Å)
b d f

Flanking Flanking Flanking Flanking

n=16,384
5’ 3’ 5’ 3’
N-N-N-A-N-N-N N-N-N-C-N-N-N
Fig. 2 | Minor groove width predicted by Deep DNAshape using extended G, and T in sequential order for flanking 5’ and in reverse order for flanking 3’.
flanking regions. a–d Heatmaps showing MGW values at central position for all e–f Boxplots showing MGW values predicted by Deep DNAshape for random
possible 7-mers. a, c MGW values as if we constructed a 7-mer query table from all sequences with fixed 5’ and 3’ caps. e Predictions are capped by ‘GCGC’, and
available MC simulations directly. b, d MGW predictions generated by the Deep f Predictions are capped by ‘AAAA’. Center line indicates the median. Box limits are
DNAshape method for all possible 7-mers. a, b All sequences with ‘A’ in the center. 75th and 25th percentiles. The whiskers extend 1.5 times the interquartile range (IQR)
c, d All sequences with ‘C’ in the center. Flanking regions from 5’ and 3’ are sorted from the top and bottom of the box. Outliers are removed in boxplots.
based on distance to the center and bp character in the heatmap. ‘_’ represents A, C,

investigate the effects of flanking regions on the core DNA structure prolongability of DNA double helical molecules. The effects of A-tracts
without requiring MD simulations or X-ray crystallography (XRC) on fluctuations can also be visualized (Supplementary Fig. 11). We
experiments. We initially examined four DNA shape features, minor observed that flanking A-tracts greatly elevate the MGW fluctuations of
groove width (MGW), propeller twist (ProT), Roll, and helix twist the cores (Supplementary Fig. 11), whereas GCGC ends slightly
(HelT), at the core of a DNA fragment for all sequence combinations, decrease the MGW fluctuations of the cores.
with the goal of examining the general neighboring effect on different
cores (Supplementary Fig. 10). The static shape of the core can show Deep DNAshape is superior to pentamer DNAshape in eluci-
bimodal or even trimodal distributions, affected by different flanking dating TF-DNA binding
sequences (Supplementary Fig. 10). High-throughput prediction of DNA shape can be applied to data from
One DNA fragment worth studying is the A-tract, which contains experimental binding assays on TF-DNA binding (Fig. 3a). The
at least three consecutive ApA, ApT, or TpT bp steps without a flexible pentamer-based DNAshape method27–29,39,47 lacked the ability to utilize
TpA bp step, and has been associated with a narrower MGW1,9. We longer-range flanking regions of sequences. Therefore, this method
permutated all k-mers with certain sequences as their cap, such as cannot diversify DNA shape changes in the core binding site affected
AAAA-NNNNNNN-AAAA, for 15-mers with A-tract caps. Next, we pre- by extended flanking regions (Fig. 3b) that may still contribute to TF
dicted the MGW using Deep DNAshape compared to counterparts binding specificity. For example, the DNA binding affinity of basic-
capped by GCGC ends. This approach permitted us to see effects from helix-loop-helix (bHLH) TFs binding to enhancer boxes (E-boxes)
different flanking cap combinations (Fig. 2e, f). The result indicated (Fig. 3c) is greatly affected by the flanking regions48. The pentamer-
that A-tracts increase the flanking MGW on their 5’ end and decrease based DNAshape method cannot access a flanking region-associated
the MGW on their 3’ end. In other words, the MGW for a short strand of change in the DNA shape of the central ‘CG’ dinucleotide in the most
DNA will be upregulated at the 3’ end and downregulated at the 5’ end, common E-box motif ‘CACGTG’. In the Deep DNAshape model, a
if accompanied by two A-tracts at the 5’ and 3’ ends. The A-tracts ‘shape layer’ capable of utilizing extended flanking regions (Fig. 3b)
showed a general trend of decreasing the MGW from 5’ to 3’, consistent may provide additional insights, such as into how bHLH proteins bind
with previous individual structural studies45,46. to their target DNA sequences.
The dynamics of DNA sequences can be difficult to describe. In We investigated genomic-context PBM data for three human
Deep DNAshape, DNA shape fluctuations can be predicted to repre- bHLH TFs, Max, c-Myc, and Mad2. We plotted MGW predictions for the
sent conformational flexibility for an individual DNA molecule. The top 25% of aligned binding data (Fig. 3d, Supplementary Figs. 12, 13)
training data are based on a compiled shape fluctuation dataset from and compared them to the predictions by the pentamer-based DNA-
MC simulations. Different bp or bp steps have different intrinsic flex- shape method. Subtler differences can be seen with Deep DNAshape in
ibilities. These fluctuation values are comparable to values that are the central motif cores than was possible with the pentamer-based
directly computed from MD simulations (Supplementary Table 1)15. DNAshape method, even though most of the cores had the same
The values contribute to the global bendability, twistability, and sequence identity ‘CACGTG’ (Fig. 3e–h). When we filtered out the TF-

Nature Communications | (2024)15:1243 4

Article https://doi.org/10.1038/s41467-024-45191-5

a High-throughput in vitro experiments b

Deep DNAshape
PBM
Neighboring effect:
Compute neighboring effect from flanking regions Shape layer
Enlarge core DNA shape
HT-SELEX
Shrink core DNA shape
No effects
SELEX-seq
3’ 5’
Analyze and align 5’
C 3’
TF-DNA binding data Pentamer-based DNAshape
Sequence Rel. binding affinity
AACCACGTGGTT 1.0 MGW
GACCACGTGGTC 0.95 AAAAA 3.38 GCACG
AAACACGTGTTT 0.95
CACGT Query
AACCACGTGGTC 0.80
AACCACGTGTTT 0.78
GCACG 5.74
CACCACGTGGTG 0.77
ACGTG
AAACACGTGGTC 0.76 TTTTG 4.76 CGTGA
... ...

TF: Max homodimer

Deep
DNAshape
DNA shape d Pentamer-based DNAshape Deep DNAshape
model considering longer flanks
Low-affinity binding site

c
Flanking Flanking
region region

MGW (Å)
High-affinity binding site

Core region
DNA shape

E-box motif:

e f g h
Pentamer-based DNAshape MGW (Å)

Deep DNAshape MGW (Å)

Fig. 3 | Comparison of Deep DNAshape and pentamer-based DNAshape meth- binding affinity. Color represents DNA shape values. Data are aligned based on core
ods on TF-DNA binding data. a General pipeline to apply Deep DNAshape on TF- binding site. Only top 25% of binding data are used. e–h Scatterplots of central four
DNA binding assay data. b Comparison of methodology to predict DNA shape core MGW values predicted by Deep DNAshape and pentamer-based DNAshape.
features from Deep DNAshape and pentamer-based DNAshape methods. c Co- Gaussian kernel density estimation plot is added to show the contour of scatter
crystal structure of DNA-bound Max protein homodimer (PDB ID: 1AN2). Flanking density. Histograms along axes are shown to visualize the distribution of DNA
regions indicate regions without protein contacts. Core is the region with protein shape values predicted by both methods. Values predicted by Deep DNAshape
contacts. d MGW predicted by pentamer-based DNAshape method vs. Deep reveal structural variations at the center of the E-box that the pentamer-based
DNAshape method for Max protein and DNA binding data, in order of relative DNAshape method was unable to reveal.

Nature Communications | (2024)15:1243 5

Article https://doi.org/10.1038/s41467-024-45191-5

DNA binding dataset to include only sequences with ‘CACGTG’ in the nearest- and next-nearest neighbors. Through Deep DNAshape, we can
core, we observed a consistent negative correlation between the propose a hypothesis for the potential binding mechanism used by the
binding affinities and the Roll and Roll-FL values predicted by Deep bHLH family of TFs to distinguish identical binding cores (Supple-
DNAshape (Supplementary Fig. 14), regardless of the in vitro experi- mentary Fig. 16), although this hypothesis will require further investi-
ment platform. These correlations are also observable using Deep gation to confirm.
DNAshape (Expt) (Supplementary Fig. 15a–c), but not when using Deep We also investigated some of the well-studied homeodomain
DNAshape (MD) (Supplementary Fig. 15d–f). Such visualizations had TFs from the homeobox (Hox) family (Fig. 4). Extradenticle (Exd)
been unachievable with the pentamer-based DNAshape method and Sex combs reduced (Scr) heterodimerize to bind to DNA
because this method only considers 2-bp flanking regions (which, in sequences. The resulting Exd-Scr heterodimer (Fig. 4a–c) prefers a
the present case, remain constant). These results highlight the effect of narrower MGW at the two ‘AY’ steps (Y: pyrimidine) in its motif
the flanks and the need for a method that accounts for more than the ‘NGAYNNAY’7. Another heterodimer, Exd with Ultrabithorax (Ubx),

Fig. 4 | Evaluations of Deep DNAshape on extended flanking regions using Hox- method was unable to predict. d Co-crystal structure of DNA-bound Exd-Ubx het-
TF binding data. a Co-crystal structure of DNA-bound Exd-Scr heterodimer (PDB erodimer (PDB ID: 4CYC). Arrows indicate locations of insertion from Exd (cyan)
ID: 2R5Z). Arrows indicate locations of insertion from Exd (cyan) and Scr (green) and Ubx (green) arginine residues into the minor groove. e, f MGW predicted by
arginine residues into the minor groove. b, c MGW predicted by the pentamer- the pentamer DNAshape and Deep DNAshape methods, for Exd-Ubx SELEX-seq
based DNAshape and Deep DNAshape methods, for Exd-Scr SELEX-seq data. Lines data. Lines are calculated based on cutoff values of relative binding affinities.
are calculated based on cutoff values of relative binding affinities. Darker colors in Darker colors in plot correspond to sequences with higher binding affinities. The
plot correspond to sequences with higher binding affinities. The comparison of the comparison of the two panels shows that Deep DNAshape predicts the MGW of the
two panels shows that Deep DNAshape predicts the MGW of the flanking regions flanking regions and the additional Exd minimum in MGW that the pentamer-based
and the additional Exd minimum in MGW that the pentamer-based DNAshape DNAshape method was unable to predict.

Nature Communications | (2024)15:1243 6

Article https://doi.org/10.1038/s41467-024-45191-5

does not have this preference at the second ‘AY’ step7 (Fig. 4d–f). Exd With the introduction of Deep DNAshape, which considers longer-
itself has another arginine in the N-terminal tail, located on the 5’ range effects, we can now replace the DNA shape features in these
side immediately upstream of the motif that inserts into the DNA machine learning models with the ones predicted by Deep DNAshape.
minor groove9. Compared to the pentamer-based DNAshape, Deep Models using DNA shape features predicted by Deep DNAshape out-
DNAshape can predict DNA shape features while considering longer perform models using features predicted with the pentamer-based
flanking regions and in the terminal regions of DNA target sites. DNAshape version (Fig. 5b, c). Furthermore, by using Deep DNAshape
Using Deep DNAshape, we revealed a narrower MGW at the correct derived features in combination with fluctuation values, we were able
locations according to co-crystal structures using the aligned 16-mer to surpass the performance of the 3-mer model with fewer degrees of
SELEX-seq data (Fig. 4b, c, e, f). Furthermore, for Scr, hypotheses can freedom (Fig. 5d). Additionally, the updated fluctuation values greatly
be made that the first ‘AY’ dip in MGW is required for binding, which improved the performance of the model compared to the previous SD
is seen across high- to low-binding-affinity sites. The binding affinity values28 (Supplementary Fig. 17). From a machine learning perspective,
is strengthened by a second ‘AY’ dip and by an additional dip for Exd the DNA shape and fluctuation values predicted by Deep DNAshape
binding for its N-terminus. For Ubx, unlike the pentamer-based contain more relevant information. We subsequently compared Deep
DNAshape method, Deep DNAshape predicted the first ‘AY’ dip in DNAshape with its variants trained by different underlying data sour-
the MGW at both A and Y bases (Fig. 4e, f). This observation is well- ces, finding that their performances were relatively similar (Supple-
aligned with other Hox protein research on the influence of DNA mentary Fig. 18). Deep DNAshape retains its performance in deeper
shape6. layers that consider longer flanking regions, while the variants peak at
layer 2 or 3 (See Supplementary Text).
Deep DNAshape exhibits improved prediction accuracy for TF-
DNA binding specificity Deep DNAshape reveals a more conserved relationship of DNA
The mechanisms of TF-DNA readout are complex, leading to sub- shape in transcription start sites between Drosophila species
stantial research on the use of machine learning to improve the pre- One key advantage of the DNAshape method is its ability to perform
diction accuracy for TF-DNA binding specificity28,49,50. In previous high-throughput predictions that can be easily applied to genomic-
studies, we have successfully improved the performance of multiple level data51. To assess the performance of Deep DNAshape on a large
linear regression tasks on TF-DNA binding affinity data by incorpor- dataset, we used a dataset of TSSs from four Drosophila species (D.
ating DNA shape features27,28 (Fig. 5a). This DNA shape encoding melanogaster, D. simulans, D. sechellia, and D. pseudoobscura)52. Pre-
approach reduced the degrees of freedom when using k-mer encod- vious research highlighted the importance of structural features in
ing, while maintaining the same level of performance. These models protein binding at TSSs53. Our analysis of DNA shape features at tran-
utilized both DNA sequence and shape features to achieve improved scription start sites (TSSs) for these four fly species revealed conserved
performance. relationships in DNA shape features among these genomic regions

Fig. 5 | Workflow of utilizing Deep DNAshape in machine learning modeling and 1mer+4shape models between the pentamer-based DNAshape and Deep DNAshape
performance comparison of binding specificity predictions. a Schematic methods. c R2 performance comparison of 1mer+13shape models between the
representation of predicting TF-DNA binding specificity with L2-regularized mul- pentamer-based DNAshape and Deep DNAshape methods. d R2 performance
tiple linear regression model, using encoded sequence features and DNA shape comparison between 1mer+13shape+FL models (representing the Deep DNAshape
features. b–d In vitro TF-DNA binding experimental data used in these models approach with all available features) and 3mer models without using DNA shape
include gcPBM, SELEX-seq, and HT-SELEX data. b R2 performance comparison of features.

Nature Communications | (2024)15:1243 7

Article https://doi.org/10.1038/s41467-024-45191-5

despite sequence differences51. Data for this analysis were derived Figs. 2 and 3). The Deep DNAshape method establishes a general
from gene expression experiments52, and Deep DNAshape was able to framework for predicting response variables at the single-
process these data in a matter of minutes on a machine equipped with nucleotide level, considering local environments for variable
an NVIDIA A100. Our results showed more evolutionarily conserved lengths of DNA, and provides valuable improvements upon the
relationships in the genomic structural features across these four fly widely used DNAshape method.
species, as evidenced by the MGW, ProT, Roll, and HelT values (Sup- DNA shape information is valuable in understanding TF-DNA
plementary Fig. 19, Supplementary Table 4). binding events1,2,28,48,54–56. One crucial application of the Deep DNA-
shape method is to investigate DNA shape preferences in TF-DNA
Discussion binding. With the ability to predict accurate DNA shape features con-
In this study, we introduce Deep DNAshape, a high-throughput deep- sidering long-range flanking regions, our results provide insight into
lerarning method that accurately predicts 3D DNA structural para- TF-DNA binding interactions that were previously challenging to
meters, or DNA shape, for any DNA sequence. Our pentamer-based investigate (Figs. 3 and 4). To quantify the extent to which Deep
high-throughput prediction method, DNAshape, successfully pre- DNAshape provides more information content than our pentamer-
dicted DNA shape features, including MGW, inter-bp features, and based DNAshape method, we applied it to the prediction of TF-DNA
intra-bp features28,39, and was useful for studying structural TF-DNA binding specificity and found a significant improvement compared to
binding mechanisms without the need of extensive molecular simu- pentamer-based DNAshape (Fig. 5). The Deep DNAshape method still
lations or structural biology experiments. However, it had several operates in a high-throughput manner, allowing for the investigation
major drawbacks and limitations. of data on a genomic scale. In a brief application, we demonstrated
First, DNAshape relied entirely on a pentamer query table to closer evolutionary relationships in DNA shape parameters between
predict DNA shape features at central positions of a pentamer. Thus, to four Drosophila species in TSS regions (Supplementary Fig. 19). Other
predict the DNA shape features for the pentamer core, only two bp of recent studies57–62 incorporating DNA shape features derived from the
adjacent nucleotides in each direction were considered. The query pentamer model will very likey benefit from using Deep DNAshape.
table could not cover the full two bp of possible flanking combinations Thus, the Deep DNAshape method unlocks a whole new level of pos-
for inter-bp features, and any sequences in the flanks beyond two bp sibilities for genomic studies of DNA shape.
were not included in predictions. This query table may be sufficient if Besides training the Deep DNAshape model using MC simulations
flanking regions over two bp away do not contribute to DNA shape at as the underlying data, the model is capable of being re-trained using
the core of the target sequence at first approximation. However, evi- alternative data sources. Our benchmark analysis revealed that the
dence has shown48 that such extended flanking regions may affect the Deep DNAshape variants Expt (trained using experimentally solved
binding affinity and DNA shape at the core of the DNA fragment of structural data) and MD (trained by MD simulation data) provided
interest. Such effects are evident in the training curve (Supplementary similar DNA shape predictions and performance metrics across mul-
Figs. 2, 3) of Deep DNAshape. tiple applications. However, the underlying data of these variants were
Second, sequence combinations of the flanking regions outside affected by noise, artifacts, or low coverage (see Supplementary Text),
the pentamers were not evenly distributed in the simulation data, making MC simulations the current optimal choice for studying effects
which introduced biases towards each pentamer calculation (Supple- of longer flanking regions. Future advancements in MD simulations
mentary Fig. 9). For example, if a pentamer is ‘ACGTA’ and the only could allow Deep DNAshape to be easily transitioned to using such
heptamer in the simulation data is ‘CACGTAG’, assuming the third bp data as the underlying source.
away from the center is still affecting the DNA shape in the core, the In addition to predicting DNA shape, Deep DNAshape is also a
prediction of DNA shape features for ‘ACGTA’ is always biased towards general framework for processing variable DNA lengths in a layer-by-
‘C’ and ‘G’ flanking the pentamer. layer manner, by using expanded neighboring-bp information in DNA
Our proposed Deep DNAshape method (Fig. 1) addresses the sequences and by inferring response variables (e.g., nucleosome
limitations of our previous method and significantly improves the positioning or 3D genome interaction) at the single-nucleotide level.
ability to predict DNA shape. Unlike the pentamer-based DNAshape Other research that utilizes DNA shape as a feature would likely benefit
method that relied on pentamer query tables and sliding-window from the use of the Deep DNAshape method. To further improve Deep
algorithms to assign DNA shape values, Deep DNAshape eliminates the DNAshape, one could focus on generating additional simulation data,
use of these tools. All k-mers needed to occur at least once in the optimizing the DNA shape inference equation or network architecture,
simulation data to complete the pentamer query table, but table and expanding the model to include chemically modified bp63 or
completion was not sufficient to prevent statistical bias or artifacts. mismatched bp11, if such data could be acquired on a large scale.
Even with pentamers that occurred more than 250 times28 in the
simulation data, the data still did not fully cover all hexamers (Sup- Methods
plementary Figs. 6 and 7), let alone heptamers (Fig. 2a, b, Supple- DNA structural simulations and DNA shape analyses
mentary Fig. 5), octamers, or longer k-mers. Different flanking regions DNA sequences with variable lengths are initialized and simulated by
have different effects on the DNA shape features at the core of DNA a Monte-Carlo (MC) method27. After filtering out artifacts, the total
target sites, and these influences may be relevant to sequence com- number of valid DNA simulations is 2,12127. DNA shape features are
ponents and are inferable. Deep DNAshape infers the missing values then calculated by Curves (version 5.3)12,25 and minor groove width
with the learnt influence of flanking regions on the core, while not (MGW) is symmetrized with respect to each bp27. Data are pooled
harming the prediction considering short flanking regions (Supple- into a training file that includes all DNA sequences and their corre-
mentary Table 2). sponding DNA shape values. DNA shape fluctuation values are also
Because the method learns the effects of flanking regions and calculated by Curves through analyzing the sampling process from
predicts DNA shape features in a layer-by-layer manner, it is self- the MC simulations. Fluctuation values represent the variance that
constrained. Each new layer infers the DNA shape features based on occurred during simulations. These values can be viewed as a cor-
the DNA shape features inferred by its previous layer, while con- related measure of DNA conformational flexibility of an individual
sidering one more bp of flanking regions. Outputs of all layers in the DNA molecule. In this study, the addition of ‘-FL’ to a DNA shape
training are used in calculating the loss function, enabling the model feature indicates the fluctuations (FL) of that specific DNA shape
to predict DNA shape features considering any length of flanking feature. Note that DNA shape can be derived from other data sour-
regions, with near-zero overfitting in deeper layers (Supplementary ces (see Supplementary Methods).

Nature Communications | (2024)15:1243 8

Article https://doi.org/10.1038/s41467-024-45191-5

Definitions of DNA shape features Dropout layers and average layers

DNA shape features are used to describe DNA structure numerically, in For each DNA shape layer, a feature vector is generated for every node
a bp-to-bp manner. The set of DNA shape features used includes six in the DNA sequence. To prevent overfitting, the feature vector passes
inter-bp features, six intra-bp features, and two minor groove features. through a dropout layer. During prediction, the dropout layer is not
Additional DNA shape features include helical axis features and back- used. Finally, the feature vectors are averaged into a single value and
bone torsion features, which can be analyzed and predicted if supplied post-processed to remove the effects of normalization before making
to the model. Inter-bp features are used to describe translational dis- a prediction.
tances (in Å) and rotational angles (in °) between adjacent bp. Speci-
fically, the six inter-bp features are ‘Shift’, ‘Slide’, ‘Rise’, ‘Tilt’, ‘Roll’, and DNA sequence encoding
‘Helix twist (HelT)’ (Fig. 1g). The intra-bp features describe the trans- To represent the DNA sequences, four characters, ‘A’, ‘C’, ‘G’ and ‘T’ are
lational distances and rotational angles between two bases in a single used, and are one-hot encoded into four arrays as [1,0,0,0], [0,1,0,0],
bp. The six intra-bp features are ‘Shear’, ‘Stretch’, ‘Stagger’, ‘Buckle’, [0,0,1,0] and [0,0,0,1]. For example, a sequence, ‘ACGTGCG’, is
‘Propeller twist (ProT)’, and ‘Opening’ (Fig. 1g). The minor groove represented as
features describe the groove geometry and electrostatic potential in 2 3 2 3
the center of the minor groove. The groove features used in this study A 1 0 0 0
are ‘Minor groove width (MGW)’ (Fig. 1g) and ‘Electrostatic Potential 6C7 60 1 0 07
6 7 6 7
(EP)’29. For detailed information, refer to28. 6 7 6 7
6G7 60 0 1 07
6 7 6 7
6 7 6 7
6T7=60 0 0 17 ð4Þ
Pre- and post-processing of DNA shape values 6 7 6 7
6G7 60 0 1 07
To ensure robustness in our DNA shape analysis, we utilized a nor- 6 7 6 7
6 7 6 7
malization method to account for different ranges of minimum and 4C5 40 1 0 05
maximum values of DNA shape features. We compensated for extreme G 0 0 1 0
values (e.g., from simulation artifacts) by using the following normal-
ization equation: To represent the DNA sequences as dinucleotides, one-hot
encoding for dinucleotides was used. Each dinucleotide will be
^ = ðS  e
S SÞ=ðS1st  S99th Þ ð1Þ assigned a vector with 16 binary values to represent 16 possible dinu-
cleotide combinations. The same example ‘ACGTGCG’ will be encoded
Here, S represents the DNA shape feature analyzed from the MC as
simulation. e
S denotes the median of the DNA shape feature values
2 3 2 3
within the dataset, while S1st and S99th mark the first percentile and last AC 0100,0000,0000,0000
percentile of the sorted DNA shape feature values, respectively. After 6 CG 7 6 0000,0010,0000,0000 7
6 7 6 7
normalization, S^ may still exceed the range of −1 to +1 for extreme 6 7 6 7
6 GT 7 6 0000,0000,0001,0000 7
6 7 6 7 ð5Þ
values. However, during the training process, the tanh activation 6 TG 7 = 6 0000,0000,0000,0010 7
6 7 6 7
function enforces a strict upper and lower bound of S1st and S99th , 6 7 6 7
4 GC 5 4 0000,0000,0100,0000 5
ensuring that the output remains within a normalized range of −1 to +1.
In postprocessing of the model output, given the normalized DNA CG 0000,0010,0000,0000
shape values predicted by our model, S can be directly computed from
^ using the reverse of the above equation given knowledge of e
S S, S1st , An unknown nucleotide ‘N’ can also be represented this way,
and S99th . where ‘N’ will be a vector by taking the average values upon all
possible values. In the implementation of Deep DNAshape, to bal-
Definition of the DNA shape layer ance the degree of freedom of the terminal bases, ‘N’ caps are
The DNA shape layer is used to treat linear DNA sequences as double- added to both terminals, but they are removed in the final predic-
linked nodes or a graph, where each node can be a single nucleotide or tion. These ‘N’ caps have a 3’ bond and a 5’ bond connected to
dinucleotide. Each node is connected by its 5’ node and 3’ node themselves.
through a forward and a backward ‘bond’ (edge). To compute the
output features Xei of the shape layer, features are gathered for each Deep DNAshape model design and learning objectives
node from the previous node (X i1 ), current node (X i ), and next node The model is designed for step-by-step, expandable learning of DNA
(X i + 1 ). Xei is computed through the following equation, shape features for any given DNA sequence, with input as one-hot-
encoded DNA sequences and output as predicted DNA shape features
 
f = f σðλ  X + αÞ, X
X ð2Þ for each position on the sequence. Inter-bp DNA shape features are
i i i
encoded as dinucleotides, and sequences are represented as linear
where, double-linked nodes. Mean absolute error (MAE) is used as the loss
function, calculated between predicted and MC-simulated DNA shape
λ = ω1  X i1 + ω2  X i + 1 + B1 values from output layers. A postprocessing step is used to recover
DNA shape values from normalization. Individual models are trained
for each DNA shape feature, with one ‘self’ convolutional layer for
α = θ1  X i1 + θ2  X i + 1 + B2 ð3Þ
input and seven ‘shape layers’ following it.

f is a trainable gated recurrent unit (GRU) cell with a sigmoid recurrent Hyperparameter search
activation function and tanh activation function. σ is ReLu activation Hyperparameters for the Deep DNAshape model include learning rate,
followed by batch normalization function. ω1 , ω2 , θ1 , θ2 , B1 and B2 are optimizer, filter size, and others. These hyperparameters are grid
trainable variables to be learnt from the dataset by the model. The self- searched for each DNA shape feature and evaluated on a separate
shape layer is a one-dimensional convolutional layer that transforms training and validation dataset. The best-performing hyperparameters
the feature number to match later DNA shape layers. The feature are selected and applied to all DNA shape features. After hyperpara-
dimension remains the same for all DNA shape layers. meter searches, the model parameters are set as follows: number of

Nature Communications | (2024)15:1243 9

Article https://doi.org/10.1038/s41467-024-45191-5

shape layers: 7, learning rate: 0.05, number of epochs: 1500, optimizer: 4. Siggers, T. W. & Honig, B. Structure-based prediction of C2H2 zinc-
stochastic gradient descent (SGD) with momentum 0.95, dropout finger binding specificity: sensitivity to docking geometry. Nucleic
ratio: 0.5 and filter size: 64. Acids Res. 35, 1085–1097 (2007).
5. Abe, N. et al. Deconvolving the recognition of DNA shape from
Data of TF-DNA binding assays sequence. Cell 161, 307–318 (2015).
We collected and used relative binding affinity data captured in mul- 6. Zeiske, T. et al. Intrinsic DNA shape accounts for affinity differences
tiple experiments, the same as were used in28. The genomic-context between Hox-cofactor binding sites. Cell Rep. 24, 2221–2230
PBM data include human TFs c-Myc, Max, and Mad264, with ‘Max’ (2018).
representing a Max-Max homodimer, ‘c-Myc’ a c-Myc-Max hetero- 7. Slattery, M. et al. Cofactor binding evokes latent differences in DNA
dimer, and ‘Mad2’ a Mad2-Max heterodimer. HT-SELEX data include binding specificity between Hox proteins. Cell 147, 1270–1282
many TFs in multiple TF families65. SELEX-seq data include TFs in the (2011).
Hox family from Drosophila5,7. 8. Kribelbauer, J. F. et al. Context-dependent gene regulation by
homeodomain transcription factor complexes revealed by shape-
L2-regularized multiple linear regression model for TF-DNA readout deficient proteins. Mol. Cell 78, 152–167.e11 (2020).
binding prediction 9. Rohs, R. et al. The role of DNA shape in protein-DNA recognition.
The L2-regularized multiple linear regression model is designed to Nature 461, 1248–1253 (2009).
predict relative TF-DNA binding specificity from the data of the TF- 10. Dantas Machado, A. C. et al. Landscape of DNA binding signatures
DNA binding assays28. The model encodes DNA sequences as k-mer of myocyte enhancer factor-2B reveals a unique interplay of base
(k = 1,2,3) sequence features and any number of DNA shape features. and shape readout. Nucleic Acids Res. 48, 8529–8544 (2020).
Unless otherwise specified, ‘4shape’ indicates four shape features 11. Afek, A. et al. DNA mismatches reveal conformational penalties in
including MGW, ProT, Roll, and HelT, and ‘13shape’ indicates MGW protein–DNA recognition. Nature 587, 291–296 (2020).
plus the 6 inter-bp features and 6 intra-bp features. The DNA shape 12. Lavery, R. & Sklenar, H. The definition of generalized helicoidal
features are normalized according to the minimum, maximum, and parameters and of axis curvature for irregular nucleic acids. J.
standard deviation seen in the dataset. Input data are separated into Biomol. Struct. Dyn. 6, 63–91 (1988).
10 folds. In each fold of the training and test data, another 10-fold 13. Pérez, A., Luque, F. J. & Orozco, M. Frontiers in molecular dynamics
cross validation is used in the training data to select lambda values simulations of DNA. Acc. Chem. Res. 45, 196–205 (2012).
for the L2 term. The model then uses the selected lambda to fit the 14. Pérez, A., Lankas, F., Luque, F. J. & Orozco, M. Towards a molecular
training data and predict the values for the test data. In the end, 10 dynamics consensus view of B-DNA flexibility. Nucleic Acids Res.
folds of predictions are combined to assess the model performance 36, 2379–2394 (2008).
as R2 . The number of sequence features, for sequences with length n, 15. Pasi, M. et al. μABC: a systematic microsecond molecular dynamics
is 4n for 1-mers, 16  ðn  1Þ for 2-mers, and 64  ðn  2Þ for 3-mers. study of tetranucleotide sequence effects in B-DNA. Nucleic Acids
The number of shape features, for sequences with length n, is n for Res. 42, 12272–12283 (2014).
each bp feature and groove feature, and n  1 for each bp step 16. Heddi, B., Oguey, C., Lavelle, C., Foloppe, N. & Hartmann, B.
feature. Intrinsic flexibility of B-DNA: the experimental TRX scale. Nucleic
Acids Res. 38, 1034–1047 (2010).
Reporting summary 17. Marin-Gonzalez, A., Vilhena, J. G., Perez, R. & Moreno-Herrero, F. A
Further information on research design is available in the Nature molecular view of DNA flexibility. Q. Rev. Biophys. 54, e8 (2021).
Portfolio Reporting Summary linked to this article. 18. Haran, T. E. & Mohanty, U. The unique structure of A-tracts and
intrinsic DNA bending. Q. Rev. Biophys. 42, 41–81 (2009).
Data availability 19. Nikolova, E. N., Bascom, G. D., Andricioaei, I. & Al-Hashimi, H. M.
TF-DNA binding datasets were derived from public resources (see Probing sequence-specific DNA flexibility in A-tracts and
Methods, Data of transcription factor (TF)-DNA binding assays). Raw pyrimidine-purine ssteps by nuclear magnetic resonance 13C
data for the underlying training datasets (MC, MD and Expt) were relaxation and molecular dynamics simulations. Biochemistry 51,
sourced from reference27 and public databases (see Supplementary 8654–8664 (2012).
Methods for details). PDB IDs used in the main text are 1AN2, 2R5Z 20. Ngo, T. T. M. et al. Effects of cytosine modifications on DNA flex-
and 4CYC. Process data of TF-DNA binding and underlying training ibility and nucleosome mechanical stability. Nat. Commun. 7,
data are deposited at Zenodo <https://doi.org/10.5281/zenodo. 10813 (2016).
10403307>. Source data are provided with this paper. 21. Li, S., Peng, Y., Landsman, D. & Panchenko, A. R. DNA methylation
cues in nucleosome geometry, stability and unwrapping. Nucleic
Code availability Acids Res. 50, 1864–1874 (2022).
All code was implemented in Python. All code related to training, 22. Ghoshdastidar, D. & Bansal, M. Flexibility of flanking DNA is a key
prediction, and the pre-trained models – as well as an executable determinant of transcription factor affinity for the core motif. Bio-
package for predicting DNA shape features from any DNA sequence – phys. J. 121, 3987–4000 (2022).
can be found at https://github.com/JinsenLi/deepDNAshape66. All 23. Chiu, T. P., Rao, S. & Rohs, R. Physicochemical models of
source code is provided under an BSD-3-Clause software license. protein–DNA binding with standard and modified base pairs. Proc.
Natl. Acad. Sci. USA 120, e2205796120 (2023).
References 24. Lavery, R., Moakher, M., Maddocks, J. H., Petkeviciute, D. & Zakr-
1. Rohs, R. et al. Origins of specificity in protein-DNA recognition. zewska, K. Conformational analysis of nucleic acids revisited:
Annu. Rev. Biochem. 79, 233–269 (2010). Curves+. Nucleic Acids Res 37, 5917–5929 (2009).
2. Inukai, S., Kock, K. H. & Bulyk, M. L. Transcription factor–DNA 25. Lavery, R. & Sklenar, H. Defining the structure of irregular nucleic
binding: Beyond binding site motifs. Curr. Opin. Genet. Dev. 43, acids: Conventions and principles. J. Biomol. Struct. Dyn. 6,
110–119 (2017). 655–667 (1989).
3. Paillard, G., Deremble, C. & Lavery, R. Looking into DNA recognition: 26. Lu, X. J. & Olson, W. K. 3DNA: A software package for the analysis,
zinc finger binding specificity. Nucleic Acids Res. 32, 6673–6682 rebuilding and visualization of three-dimensional nucleic acid
(2004). structures. Nucleic Acids Res 31, 5108–5121 (2003).

Nature Communications | (2024)15:1243 10

Article https://doi.org/10.1038/s41467-024-45191-5

27. Zhou, T. et al. DNAshape: A method for the high-throughput pre- 48. Gordân, R. et al. Genomic regions flanking E-box binding sites
diction of DNA structural features on a genomic scale. Nucleic Acids influence DNA binding specificity of bHLH transcription factors
Res. 41, W56–W62 (2013). through DNA shape. Cell Rep. 3, 1093–1104 (2013).
28. Li, J. et al. Expanding the repertoire of DNA shape features for 49. Zhou, T. et al. Quantitative modeling of transcription factor binding
genome-scale studies of transcription factor binding. Nucleic Acids specificities using DNA shape. Proc. Natl. Acad. Sci. USA 112,
Res. 45, 12877–12887 (2017). 4654–4659 (2015).
29. Chiu, T. P., Rao, S., Mann, R. S., Honig, B. & Rohs, R. Genome-wide 50. Alipanahi, B., Delong, A., Weirauch, M. T. & Frey, B. J. Predicting the
prediction of minor-groove electrostatic potential enables biophy- sequence specificities of DNA- and RNA-binding proteins by deep
sical modeling of protein–DNA binding. Nucleic Acids Res. 45, learning. Nat. Biotechnol. 33, 831–838 (2015).
12565–12576 (2017). 51. Chiu, T. P. et al. GBshape: a genome browser database for DNA
30. Barissi, S., Sala, A., Wieczór, M., Battistini, F. & Orozco, M. DNAf- shape annotations. Nucleic Acids Res. 43, D103–D109 (2015).
finity: A machine-learning approach to predict DNA binding affi- 52. Main, B. J., Smith, A. D., Jang, H. & Nuzhdin, S. V. Transcription start
nities of transcription factors. Nucleic Acids Res. 50, 9105–9114 site evolution in Drosophila. Mol. Biol. Evol. 30, 1966–1974 (2013).
(2022). 53. Bansal, M., Kumar, A. & Yella, V. R. Role of DNA sequence based
31. Wang, S. et al. Predicting transcription factor binding sites using structural features of promoters in transcription initiation and gene
DNA shape features based on shared hybrid deep learning archi- expression. Curr. Opin. Struct. Biol. 25, 77–85 (2014).
tecture. Mol. Ther. Nucleic Acids 24, 154–163 (2021). 54. Mathelier, A. et al. DNA shape features improve transcription factor
32. Demirci, S., Peters, S. A., Ridder, D. & Dijk, A. D. J. DNA sequence and binding site predictions in vivo. Cell Syst. 3, 278–286.e4 (2016).
shape are predictive for meiotic crossovers throughout the plant 55. Yang, J. & Ramsey, S. A. A DNA shape-based regulatory score
kingdom. Plant J. 95, 686–699 (2018). improves position-weight matrix-based recognition of transcription
33. Zhang, Q., Shen, Z. & Huang, D.-S. Predicting in-vitro Transcription factor binding sites. Bioinformatics 31, 3445–3450 (2015).
Factor Binding Sites Using DNA Sequence + Shape. IEEE ACM Trans. 56. Slattery, M. et al. Absence of a simple code: how transcription
Comput. Biol. Bioinf. 18, 667–676 (2021). factors read the genome. Trends Biochem. Sci. 39, 381–399 (2014).
34. Yang, J. et al. Prediction of regulatory motifs from human Chip- 57. Liu, Z. & Samee, M. A. H. Structural underpinnings of mutation rate
sequencing data using a deep learning framework. Nucleic Acids variations in the human genome. Nucleic Acids Res 51, 7184–7197
Res. 47, 7809–7824 (2019). (2023).
35. Rohs, R., Sklenar, H. & Shakked, Z. Structural and energetic origins 58. Zhang, Y. et al. A novel convolution attention model for predicting
of sequence-specific DNA bending: Monte Carlo simulations of transcription factor binding sites by combination of sequence and
papillomavirus E2-DNA binding sites. Structure 13, 1499–1509 shape. Brief. Bioinform. 23, bbab525 (2021).
(2005). 59. Ding, P. et al. DeepSTF: predicting transcription factor binding sites
36. Berman, H. M. et al. The nucleic acid database. A comprehensive by interpretable deep neural networks combining sequence and
relational database of three-dimensional structures of nucleic shape. Brief. Bioinform. 24, bbad231 (2023).
acids. Biophys. J. 63, 751–759 (1992). 60. Wang, Z., Xiong, S., Yu, Y., Zhou, J. & Zhang, Y. HAMPLE: deci-
37. Balaceanu, A. et al. Modulation of the helical properties of DNA: phering TF-DNA binding mechanism in different cellular environ-
next-to-nearest neighbour effects and beyond. Nucleic Acids Res. ments by characterizing higher-order nucleotide dependency.
47, 4418–4430 (2019). Bioinformatics 39, btad299 (2023).
38. Rube, H. T., Rastogi, C., Kribelbauer, J. F. & Bussemaker, H. J. A 61. Bhimsaria, D. et al. Hidden modes of DNA binding by human nuclear
unified approach for quantifying and interpreting DNA shape receptors. Nat. Commun. 14, 4179 (2023).
readout by transcription factors. Mol. Syst. Biol. 14, e7902 (2018). 62. Khan, S. R., Sakib, S., Rahman, M. S. & Samee, Md. A. H. DeepBend:
39. Chiu, T. P. et al. DNAshapeR: an R/Bioconductor package for DNA An interpretable model of DNA bendability. iScience 26, 105945
shape prediction and feature encoding. Bioinformatics 32, (2023).
1211–1213 (2015). 63. Jiang, Y., Chiu, T. P., Mitra, R., & Rohs, R. Probing the role of the
40. Young, R. T., Czapla, L., Wefers, Z. O., Cohen, B. M. & Olson, W. K. protonation state of a minor groove-linker histidine in Exd-Hox–DNA
Revisiting DNA sequence-dependent deformability in high- binding. Biophys. J. 123, 248–259 (2024).
resolution structures: Effects of flanking base pairs on dinucleotide 64. Mordelet, F., Horton, J., Hartemink, A. J., Engelhardt, B. E. & Gordân,
morphology and global chain configuration. Life 12, 759 (2022). R. Stability selection for regression-based models of transcription
41. Ivani, I. et al. Parmbsc1: A refined force field for DNA simulations. factor-DNA binding specificity. Bioinformatics 29, i117–i125 (2013).
Nat. Methods 13, 55–58 (2016). 65. Yang, L. et al. Transcription factor family-specific DNA shape
42. Chiu, T. P., Li, J., Jiang, Y. & Rohs, R. It is in the flanks: conformational readout revealed by quantitative specificity models. Mol. Syst. Biol.
flexibility of transcription factor binding sites. Biophys. J. 121, 13, 910 (2017).
3765–3767 (2022). 66. Li, J. DeepDNAshape: code release. https://doi.org/10.5281/zenodo.
43. Berger, M. F. & Bulyk, M. L. Universal protein-binding microarrays for 10403299. (Zenodo, 2023).
the comprehensive characterization of the DNA-binding specifi-
cities of transcription factors. Nat. Protoc. 4, 393–411 (2009). Acknowledgements
44. Jolma, A. et al. DNA-binding specificities of human transcription The authors acknowledge Tianyin Zhou for the development and vali-
factors. Cell 152, 327–339 (2013). dation of the pentamer-based DNAshape method. This work was sup-
45. MacDonald, D. et al. Solution structure of an A-tract DNA bend. J. ported by the National Institutes of Health [grant R35GM130376 to R.R.]
Mol. Biol. 306, 1081–1098 (2001). and the Human Frontier Science Program [grant RGP0021/2018 to R.R.].
46. Stefl, R., Wu, H., Ravindranathan, S., Sklenář, V. & Feigon, J. DNA
A-tract bending in three dimensions: Solving the dA4T4 vs. dT4A4 Author contributions
conundrum. Proc. Natl. Acad. Sci. USA 101, 1177–1182 (2004). J.L., T.C. and R.R. conceived the project. J.L. and T.C. designed the
47. Rao, S. et al. Systematic prediction of DNA shape changes due to model architecture. J.L. performed all training, evaluations, and ana-
CpG methylation explains epigenetic effects on protein–DNA lyses. J.L. wrote the manuscript with help from T.C. and R.R. The project
binding. Epigenet. Chromatin 11, 6 (2018). was supervised by R.R.

Nature Communications | (2024)15:1243 11

Article https://doi.org/10.1038/s41467-024-45191-5

Competing interests Publisher’s note Springer Nature remains neutral with regard to jur-
The authors declare no competing interests. isdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons

Additional information Attribution 4.0 International License, which permits use, sharing,
Supplementary information The online version contains adaptation, distribution and reproduction in any medium or format, as
supplementary material available at long as you give appropriate credit to the original author(s) and the
https://doi.org/10.1038/s41467-024-45191-5. source, provide a link to the Creative Commons licence, and indicate if
changes were made. The images or other third party material in this
Correspondence and requests for materials should be addressed to article are included in the article’s Creative Commons licence, unless
Remo Rohs. indicated otherwise in a credit line to the material. If material is not
included in the article’s Creative Commons licence and your intended
Peer review information : Nature Communications thanks Maria Popt- use is not permitted by statutory regulation or exceeds the permitted
sova and the other, anonymous, reviewer(s) for their contribution to the use, you will need to obtain permission directly from the copyright
peer review of this work. A peer review file is available. holder. To view a copy of this licence, visit http://creativecommons.org/
licenses/by/4.0/.
Reprints and permissions information is available at
http://www.nature.com/reprints © The Author(s) 2024

Nature Communications | (2024)15:1243 12