Articles

[Link]

Similarity regression predicts evolution of

transcription factor sequence specificity
Samuel A. Lambert1, Ally W. H. Yang2, Alexander Sasse1, Gwendolyn Cowley3, Mihai Albu2,
Mark X. Caddick3, Quaid D. Morris 1,2,4,5,6, Matthew T. Weirauch 7,8 and Timothy R. Hughes *
1,2,9

Transcription factor (TF) binding specificities (motifs) are essential for the analysis of gene regulation. Accurate prediction of
TF motifs is critical, because it is infeasible to assay all TFs in all sequenced eukaryotic genomes. There is ongoing controversy
regarding the degree of motif diversification among related species that is, in part, because of uncertainty in motif predic-
tion methods. Here we describe similarity regression, a significantly improved method for predicting motifs, which we use to
update and expand the Cis-BP database. Similarity regression inherently quantifies TF motif evolution, and shows that previous
claims of near-complete conservation of motifs between human and Drosophila are inflated, with nearly half of the motifs in
each species absent from the other, largely due to extensive divergence in C2H2 zinc finger proteins. We conclude that diver-
sification in DNA-binding motifs is pervasive, and present a new tool and updated resource to study TF diversity and gene
regulation across eukaryotes.

T
o understand the function of noncoding DNA—for exam- Improving motif predictions would help to answer questions
ple, in gene regulation—it is essential to know the potential regarding the degree of diversity and evolution of eukaryotic TF
TFs that can bind to any sequence. Libraries of experimen- motifs. TF-binding specificities are thought to be highly conserved
tally derived TF motifs, most typically position weight matrices between Drosophila and mammals8; however, the specificity resi-
(PWMs)1, are widely used and encompass at most a few thousand dues for C2H2 ZF proteins are often very different even among
motifs and are oriented mainly towards well-studied TFs in human Drosophila species (and very different from those in mammals),
and model systems (for example, JASPAR)2. However, hundreds of suggesting a much faster pace of change for some TF families9,10.
eukaryotic genomes have now been sequenced, and the analysis of TF motif diversification has also been observed in other lineages,
gene expression and corresponding sequences in regulatory regions including plants and fungi, indicating that TF evolution occurs
can be performed in almost all of them. To enable such analyses, we broadly, in parallel to the more established cis-regulatory turn-
previously described Cis-BP, a database of measured and predicted over11. Measuring motif diversification cannot be addressed by sim-
TF motifs for 59,998 TFs from 340 sequenced eukaryotes3. The ply examining the orthology patterns of whole genes and proteins,
predictions in Cis-BP were made by simple amino acid sequence because there are examples in which one-to-one orthologs (which
identity between DNA-binding domains (DBDs), with cut-offs for would be expected to have conserved motifs) in fact have differ-
each DBD type established on the basis of replicate experiments and ent motifs10,12,13 and, conversely, there are entire families (in which
pairwise comparisons of motifs from different proteins with homol- members might be expected to diverge in their binding sites) that
ogous DBD types. The cut-off prediction method yielded an 89% have identical motifs (for example, the core sequence for E2F and
precision (with undetermined recall) on these data. regulatory factor X (RFX) TFs, as well as the rigid specificity of plant
It may be possible to improve both precision and recall by WRKY TFs)3,14. Indeed, how motif diversification is dictated by pro-
adapting predictions to specific DBDs. One approach is to use tein structure and mechanisms of DNA binding is largely unknown,
known ‘specificity residues’ or to prioritize DNA-contacting except in a few cases13,15,16.
residues when measuring amino acid sequence similarity. Other We reasoned that developing a system for determining both the
approaches including affinity regression4, which predicts affin- similarity and dissimilarity of TF motifs would provide uniform
ity to DNA/RNA k-mers on the basis of amino acid k-mer com- and unbiased estimates of the conservation of TFs and their DNA-
position of proteins. Affinity regression was applied to only two binding functions. Here, we describe such a system, its incorpora-
families, however: homeodomain TFs and RNA recognition motif tion into the Cis-BP database, validation experiments in several
(RRM)-containing RNA-binding proteins. DBD-specific ‘recogni- eukaryotes and use of the system to broadly describe TF motif evo-
tion codes’ have also been described, which predict binding motifs lution across eukaryotes.
on the basis of DNA-contacting residues, for Cys2His2 (C2H2) zinc
finger (ZF) and homeodomain proteins5–7. It is unclear whether Results
and how these methods will extend to the approximately 100 other Similarity regression predicts TF motif similarity. We devel-
types of DBDs. oped a homology-based system that, when calculating sequence

Department of Molecular Genetics, University of Toronto, Toronto, Ontario, Canada. 2Donnelly Centre for Cellular and Biomolecular Research,
1

University of Toronto, Toronto, Ontario, Canada. 3Institute of Integrative Biology, University of Liverpool, Liverpool, UK. 4Department of Computer Science,
University of Toronto, Toronto, Ontario, Canada. 5Canadian Institutes For Advanced Research (CIFAR) Artificial Intelligence Chair, Vector Institute, Toronto,
Ontario, Canada. 6Ontario Institute of Cancer Research, Toronto, Ontario, Canada. 7Divisions of Biomedical Informatics and Developmental Biology,
Center for Autoimmune Genomics and Etiology (CAGE), Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA. 8Department of Pediatrics,
University of Cincinnati College of Medicine, Cincinnati, OH, USA. 9CIFAR, Toronto, Ontario, Canada. *e-mail: [Link]@[Link]

NatUre Genetics | VOL 51 | JUNE 2019 | 981–989 | [Link]/naturegenetics 981

Articles NATure GeneTiCS

a
Ptx1
(D. melanogaster)

Hoxa7
(M. musculus)

1 10 20 30 40 50 57
Amino acid position

b
AA identity

BLOSUM62

1 10 20 30 40 50 57
Amino acid position
c
4
content (bits)
Information

3
2
1

~100,000 Homeodomain pairs

DNA-binding similarity

E-score
overlap
0 0.5 1.0 ~
AA %ID

1 10 20 30 40 50 57
Amino acid position

d
6
*
SR weight

2 *
**
** *
0 *
1 10 20 30 40 50 57
Amino acid position

Fig. 1 | Overview of the similarity regression method. Similarity regression uses TF protein similarity to predict the similarity in TF sequence specificities.
The procedure and results are outlined for homeodomain TFs as an example. a, The DBD sequence of each TF is aligned to the Pfam HMM as a common
reference to generate a global alignment. Amino acids (AA) shown are colored according to standard clustal colors for two homeodomains. b, For
each pair of TFs, the amino acid similarity is measured at each position of the alignment, recording whether the two residues are identical or similar
(BLOSUM62 substitution score). This procedure is repeated for every pair of homeodomain TFs with PBM data. c, Regression is performed for a matrix
in which each row is a pair of homeodomain TFs, with the similarity of their DNA sequence specificities (E-score overlap) as the dependent Y variables
(left) and positional protein similarity as independent X values (right, identity shown). Sequence diversity among the TFs is represented here for reference,
plotted as a sequence logo (in bits) above the protein similarity matrix. %ID, percent identity. d, The regression outputs a weight vector that indicates how
much amino acid similarity in each position of the DBD contributes to DNA-binding similarity. Known specificity residues37 are indicated by an asterisk.

similarity in DBDs, assigns greater importance to positions of the all pairs for each DBD class (Fig. 1d). We tested four variations of
DBDs (for example, base-contacting or specificity residues) that this scheme, including two different regression approaches (linear
have more impact on sequence specificity. To do this, we assigned and logistic) and two different representations of sequence similar-
a weight to each residue when calculating similarity between two ity (identity and BLOSUM62 substitution scores). Here, we trained
DBDs of the same class (for example, C2H2 ZF, ETS or Forkhead). the regression models to learn highly overlapping 8-base oligo-
Figure 1 displays the overall procedure. We use regression to assign mer E-score preferences obtained from universal protein-binding
the weights: for each pair of proteins, the independent variables are microarrays (PBMs)17, as cataloged in Cis-BP3. These scores are
the binary vector of amino acid similarity at each position of an comparable among different studies, thus circumventing the poten-
alignment to the Pfam hidden Markov model (HMM) (Fig. 1a–c), tially confounding impact of motif derivation18; however, we note
whereas the dependent variable is the similarity in DNA sequence that the model could be trained on any metric of motif identity
preference (Fig. 1c). The weights are the coefficients learned over or similarity. Each variation of the model generates a different set

982 NatUre Genetics | VOL 51 | JUNE 2019 | 981–989 | [Link]/naturegenetics

NATure GeneTiCS Articles
a b
Sox
1.00 1.00 HSF
ETS E2F
Nuclear receptor SBP DM
Target precison RFX
0.75 0.75

Similarity regression
T-box
Features
No. of proteins
SR (identity) Homeodomain

Precision
SR (BLOSUM62) Myb/SANT 100
0.50 Alignment %ID 0.50 200
NAC/NAM 300
AP2 GATA 400
bHLH Forkhead Pipsqueak
0.25 0.25 Features
Zinc cluster
Random C2H2 SR (identity)
ZF bZIP SR (BLOSUM62)
0 0

0 0.25 0.50 0.75 1.00 0 0.25 0.50 0.75 1.00

Recall Alignment %ID

c d
100
Predicted TF similarity

Quartiles 75
1 Observed TF

Percentage
2 similarity
3 Highly
50 similar
Dissimilar 4
Ambiguous

Dissimilar
Ambiguous 25

Highly
similar
0
0.00 0.25 0.50 0.75 1.00 Highly similar Ambiguous Dissimilar
E-score overlap Predicted TF similarity

Fig. 2 | Similarity regression classification of TFs for highly similar or dissimilar sequence specificities. a, Precision/recall curves for homeodomains
are shown for three prediction methods: alignment percent identity, similarity regression using amino acid identity (SR (identity)), and similarity
regression using BLOSUM similarity (SR (BLOSUM62)), on held-out data across all cross-validation folds. ‘Positives’ are pairs of TFs with highly similar
specificities (E-score overlap >twenty-fifth percentile of replicate experiments). ‘Negatives’ are all other pairs. b, Scatter plot comparing recall values
(predicting highly similar specificities at 75% precision threshold) for similarity regression versus percent identity, for each TF family. The best of the
four similarity regression models is shown. Points are sized according to the number of PBM experiments used for training and colored according to the
amino acid features used in each model. Domain abbreviations are taken from Pfam47, where full names can be found. c, Smoothed density estimates
for homeodomain E-score overlaps in each predicted TF similarity class. Densities are filled according to the quartiles of the data. Vertical dashed lines
indicate the E-score overlap thresholds used to define dissimilar (blue line) and highly similar (black line) TF specificities in the initial data. d, Percentage
of actual TF similarities within each predicted TF similarity class, for new PBM data. White dotted lines show expected percentages for the highly similar
and dissimilar classes (that is, thresholds were chosen to achieve these levels on training data).

of weights, which are selected by leave-one-out cross-validation among TFs. For Sox proteins, however, the weights are much higher
(Supplementary Fig. 1). Using the same cross-validation, we select in residues that contact the minor groove, consistent with struc-
the best model for each DBD class among these four models (linear tural data19, whereas GAL4/zinc cluster proteins, the dimerization
or logistic regression, using amino acid identity or similarity) or, of which is organized along the DNA backbone20,21, receive high
as a fifth possibility, the simple alignment percent identity that is weights in backbone-contacting residues (Supplementary Fig. 3c).
currently used by Cis-BP3. We refer to this procedure as similar- Similarity regression shows notable improvement in recall at
ity regression. Application of similarity regression to TF families identical precision values over percent identity alone. This improve-
for which DBDs are present in arrays (for example, C2H2 ZFs) is ment is notable for homeodomains (Fig. 2a) and other families with
explained in Supplementary Fig. 2. a large amount of PBM data (summarized in Fig. 2b). In these preci-
Similarity regression has several advantages over previous sion/recall curves, positives are those pairs of proteins with E-score
approaches. It identifies residues that are informative regarding overlap that exceeds the twenty-fifth percentile of experimental rep-
DNA sequence specificity. The weights obtained are highly biased licates (the same threshold was used previously3) and negatives are
towards DNA-contacting regions and specificity residues, if known. all other pairs. Thus, our precision/recall curves are likely underes-
Figure 1d illustrates the weights for the well-studied homeodomain timates, because our stringent negative definition includes highly
class, which has established specificity residues in DNA-contacting similar experiments that are just below the threshold.
positions5. Weights for all eukaryotic DBD families with similar- Importantly, similarity regression can also be used to predict
ity regression models are given in Supplementary Data 1 (models whether two proteins are highly unlikely to share DNA sequence
for homeodomains and C2H2 ZFs are shown in Supplementary specificity: using the same learned weights described above, a
Fig. 3a,b, respectively). These weights correspond to known mecha- threshold can be identified below which proteins will almost always
nisms of DNA recognition: there is a strong relationship between bind very different sequences. In this analysis, we defined differ-
similarity regression model weight and DNA contact frequency ent sequence preferences to be an overlap of 20% or less among
(Supplementary Fig. 3c). For example, similarity regression pin- the highly preferred 8-mers oligomers. We allowed some overlap
points known binding modes: for most TFs, weights are higher in because many families bind a characteristic sequence ‘core’. For
the residues that contact the major groove, which is predominant example, many homeodomains bind TAAT-like sequences, even

NatUre Genetics | VOL 51 | JUNE 2019 | 981–989 | [Link]/naturegenetics 983

Articles NATure GeneTiCS
a b
DBD AA identity (%) DBD AA identity (%)
20 40 60 80 100 C. sativa A. thaliana 20 40 60 80 100 N. crassa S. cerevisiae
PK27213.1 cre-1

PK27693.1 msn-1
PK19621.1
NCU03043
PK21166.1
NCU05061
PK28187.1
NCU09576
PK07741.1
NCU05909
PK12011.1
NCU06487
PK05732.1

PK26523.1 NCU01629

PK07924.1 NCU03421

PK21687.1 c A. nidulans S. cerevisiae

HSF ANIA_06998

Predicted motif similarity

Highly similar Dissimilar
Homeodomain ANIA_00885
Ambiguous DBD not shared
Rap1 ANIA_01906

Myb/SANT ANIA_04524

Pipsqueak ANIA_07516

Fig. 3 | PBM data from the plant C. sativa and model fungi A. nidulans and N. crassa for TFs with conserved and dissimilar motifs. NNs for each new TF
with PBM data were identified by finding the most similar TF (by similarity regression score) with a motif from either Arabidopsis thaliana (for C. sativa)
or Saccharomyces cerevisiae (for A. nidulans and N. crassa). a–c, Motifs for Myb/SANT TFs from C. sativa (a), C2H2 ZF TFs from N. crassa (b) and TFs
from five other TF families in A. nidulans (c) are shown, with a neighbor-joining tree scaled by DBD alignment percent identity in a and b. The colored bar
represents predicted motif similarity. See Supplementary Fig. 5 for a comparison between similarity regression-predicted similarity and NN TF similarity
for all new PBM data. Silhouettes of each species are displayed, adapted from Phylopic ([Link] under a Creative Commons license (https://
[Link]/licenses/by-sa/3.0/).

though their most highly preferred 8-base oligomers differ among that are dissimilar to TFs with known motifs, or both. We used these
family members. For each DBD type, we set a similarity regression data as a validation set to test how well similarity regression models
score threshold using a negative predictive value (NPV), at which measure the similarity of TF sequence specificity on unseen data
95% of pairs of proteins at that similarity score indeed have different (Fig. 2d). The TF similarity classifications for the newly analyzed
sequence preferences. As shown in Supplementary Fig. 4a, similar- proteins are correct for 81.2% of the highly similar and 95.2% of the
ity regression outperformed percent identity at discriminating these dissimilar pairs. These results held over a range of percent identity
dissimilar pairs by achieving a higher recall at the same NPV. to other proteins in Cis-BP, confirming that the models are accu-
In all subsequent analyses, we use similarity regression to classify rate with independent data. Indeed, there is an overall correlation
all pairs of proteins that share the same DBD type as ‘highly similar’ between similarity regression score and E-score overlap between the
if their similarity regression score is over the positive pair threshold, held-out data and the most similar training construct (by similar-
‘dissimilar’ if it is below the negative pair threshold or ‘ambiguous’ if ity regression score) for the similarity regression model of each TF
it is between the two thresholds. The ambiguous score range can, in family (Supplementary Fig. 5, median R2 = 0.63). Figure 3 provides
some cases, be quite large; Fig. 2c shows that it is predictive of inter- examples of conservation and divergence of motifs in the new data.
mediate 8-base oligomer overlap for homeodomain TFs. Similar
phenomena are observed in other TF families (data not shown). Comparison to alternative motif prediction methods. We also
Multi-class accuracy of the similarity regression models, and their investigated whether similarity regression could accurately predict
improvement over percent identity, is summarized by the Matthews motifs by comparing motif predictions with generalist methods
correlation coefficient in Supplementary Fig. 4b. Similarity regres- (for example, affinity regression4 applied to all TF families) and
sion outperforms percent identity in all but four TF families. domain-specific recognition codes. Similarity regression predicts
similarity in DNA sequence specificity of two proteins, whereas
New PBM data validate motif similarity classifications. To con- affinity regression directly predicts preferences of TFs or RNA-
firm that the models correctly classify previously unseen proteins, we binding proteins to individual DNA or RNA sequences on the basis
generated new PBM data for 340 TFs, including 15 human TFs and of their protein sequences. Nonetheless, the two can be compared
other TFs that represent multiple eukaryotic kingdoms, with a par- using similarity regression to predict 8-base oligomer preferences
ticular focus on Cannabis sativa (a medicinal plant), Caenorhabditis from proteins that should have highly similar sequence preferences,
briggsae (a nematode), Aspergillus nidulans and Neurospora crassa and similarity regression outperforms affinity regression in many
(model fungi) (Supplementary Table 1). These TFs were selected to cases. Using an identical training set (that is, the same experiments
increase the number of experimentally determined motifs for TFs in on the same proteins), similarity regression slightly outperformed
these species of interest, to obtain novel motifs by analyzing proteins affinity regression when predicting 8-base oligomer Z-score profiles

984 NatUre Genetics | VOL 51 | JUNE 2019 | 981–989 | [Link]/naturegenetics

NATure GeneTiCS Articles
for our 315 held-out constructs across 19 TF families (see Methods 100
for details) using either the single most similar protein (nearest
neighbor (NN), P < 0.05; Supplementary Fig. 6a) or by predict- 75
ing the Z-scores as a composite of up to five most similar proteins

Percentage TFs
(top 5, P < 0.01; Supplementary Fig. 6a). Similarity regression has

Conserved
the added benefit that it does not make predictions for dissimilar 50
proteins when the predictions are poor, whereas affinity regression
makes a prediction for every protein, without an associated quality 25
metric. In comparison to affinity regression, similarity regression
predictions that use only highly similar TF pairs (n = 104) have a
0
higher correlation to the measured Z-score profiles than affinity
regression using the NN (P < 0.01) or top 5 predictions (P < 0.0001)
100
(Supplementary Fig. 6). These outcomes hold for most (although
not all) individual TF families analyzed in isolation. For example,
whereas similarity regression performs equivalently to affinity

Potentially conserved
75

Percentage TFs
regression for zinc cluster TFs, it scores higher for the homeodo-
mains and C2H2 ZF families (Supplementary Fig. 6b–d). 50
We also compared similarity regression predictions to those
produced by state-of-the-art recognition code algorithms for
C2H2 ZF and homeodomain TFs, which directly predict PWMs 25
from DBD sequences. One of the C2H2 ZF predictors uses sup-
port vector machines22 whereas the other uses random forests 0
(ZFModels)23; the homeodomain predictor (PreMoTF)5 also
uses random forests. To compare to these specialist prediction 100
tools, we converted the similarity regression-predicted Z-score
profiles to PWMs using a simple alignment of the top-scoring
75
8-base oligomers (PWMalign)18. We compared the motif similar-
Percentage TFs

ity between predicted motifs to those derived from the held-out

Diverged
PBM experiments described above, which included 34 C2H2 ZF 50
and 17 homeodomain proteins. Supplementary Fig. 7a shows
that, in most cases, similarity regression motif predictions for 25
C2H2 ZFs are more similar to the experimental motifs than the
motifs predicted by either of two recognition codes, regardless of
whether similarity regression predictions were filtered to be highly 0

similar. For homeodomains, there was no significant difference 100 250 500 750 1,000 1,250
between predictions made by the recognition code and those Divergence time (Ma)
made by similarity regression based on multiple (that is top 5 or Kingdom Fungi Plants Metazoa
highly similar) neighbors (P > 0.05); however, PreMoTF outper-
forms similarity regression NN-based motif predictions (P < 0.05) Fig. 4 | Conservation of TF motifs within major eukaryotic kingdoms. The
(Supplementary Fig. 7b). average percentage of TFs for which the closest TF in the other species is
conserved (similarity regression classifies as highly similar), potentially
New TF similarity predictions improve Cis-BP. To capitalize on conserved (ambiguous) or diverged (dissimilar and unshared DBDs) was
the increased recall of similarity regression relative to percent iden- calculated for each pair of species from the same kingdom (48 metazoan,
tity, we implemented the method in Cis-BP, which compiles known 15 plant and 15 fungal species; Supplementary Fig. 7b). Each point represents
TF motifs and tracks homology relationships among similar TFs. the average percentage of TFs within each category, for each pair of species
Since Cis-BP was described in 2014, both the number of sequenced (that is, average of species X versus species Y and Y versus X), plotted
eukaryotes and the number of known motifs has roughly doubled. against divergence time in millions of years. Divergence time is plotted on
We therefore updated Cis-BP, which now includes 741 genomes a square root scale to visualize differences between closely related species.
(updated from 340) and 11,491 experimentally determined motifs Lines and shading show the LOESS regression fit and 95% confidence
that correspond to 4,559 distinct TFs (updated from 6,559 motifs interval, respectively. Ma, million years ago.
for 3,202 distinct TFs) and implemented similarity regression across
all 392,333 known and putative eukaryotic TFs. We also updated
many other properties of the database (for example, genome builds Cis-BP database can be found at [Link]
and DBD models) (Methods). where TF annotations, motifs and PBM data compiled from our
The incorporation of similarity regression in Cis-BP increases the laboratory and other public databases can be accessed and down-
number of TFs with predicted motifs by more than 6,000 compared loaded. In addition to increased coverage, the new build—which
to our previous method, at the same expected precision (a 4.2% over- contains many more genomes—also identifies many new families of
all increase, on identical genomes, DBDs and motifs). Coverage of TFs with still-unknown sequence specificity.
numerous TF families is increased markedly (Supplementary Fig. 8a).
For instance, ten TF families more than doubled their motif cover- Evolution of TF sequence specificity across Eukarya. Finally,
age, including zinc cluster (123% increase) and Sox (162% increase) we used the motif predictions and the Cis-BP update to gain an
TFs, the second and seventh most abundant families in Cis-BP, overview of TF motif conservation and divergence over eukary-
respectively. The average species now has 4% more TFs with motifs otic evolution. We focused on 84 species with well-annotated
(experimental and predicted), yielding an average motif coverage of genomes (present in Ensembl and/or Uniprot; species are listed in
37% (with 73% for human) (Supplementary Fig. 8b) and a total cov- Supplementary Fig. 8b). We included all TF DBDs in this analysis,
erage of 165,030 out of 392,333 eukaryotic TFs (42%). This updated regardless of conservation level or patterns of orthology, to gain

NatUre Genetics | VOL 51 | JUNE 2019 | 981–989 | [Link]/naturegenetics 985

Articles NATure GeneTiCS

or
a

he ept
ai
bH dom

T
F o r re
Homo sapiens

AN
ZF

a
eo

/S

A
2

LH

le

x
om

IP

bo

AT
2H

rk

S
yb
uc
x

ET
So
bZ

T-
M

G
C

N
Pan troglodytes

Macaca mulatta

Mus musculus

Canis familiaris TF comparisons

H. sapiens or
A. thaliana
Gallus gallus

Tetraodon nigroviridis Comparison

species

Drosophila melanogaster

Caenorhabditis elegans

Nematostella vectensis
800 600 400 200 0
TF similarity Highly similar Ambiguous Dissimilar DBD not shared
Divergence time (Ma)

n
ai
b

m
T

x
AM

do

bo
AN

Arabidopsis thaliana

ZF
eo
/N

KY

S
AS
/S

LH

1
AD
AC

IP

om

2H

B
yb

R
R

R
AP

LO
bH

B3

bZ

FA
W

M
M

G
N

C
Arabidopsis lyrata

Cannabis sativa

Vitis vinifera

Solanum lycopersicum

Zea mays

Oryza sativa

Amborella trichopoda

Physcomitrella patens

Selaginella moellendorffii
400 200 0
Divergence time (Ma)

Fig. 5 | Motif divergence of TF families in metazoans and plants. a, Nested pie charts showing the percentage of human TFs for which the closest TF in
other metazoans is highly similar, ambiguous, dissimilar or not shared, for the 11 most abundant metazoan DBDs. The outer ring of each pie chart shows
the proportion of human TFs in each similarity regression-predicted similarity class relative to the other species; the inner ring shows the proportion of TFs
for the other species, relative to humans. b, Motif similarity between A. thaliana and other plants, for the 13 most abundant plant DBDs. Silhouettes of each
species are displayed, adapted from Phylopic ([Link] and/or are available under a Creative Commons license ([Link]
publicdomain/zero/1.0/).

a complete picture of motif conservation. For each protein, we Figure 4 shows that eukaryotic kingdoms display qualitatively
identified the protein with the highest similarity regression model similar trends in the proportion of TFs within each of the categories,
score as described above in each other species and recorded the with respect to divergence time. Around 100 million years ago (for
classification (that is, highly similar, ambiguous or dissimilar). If example, origin of placental mammals and eudicot plants), approxi-
there is no protein with the same DBD type in the other species, mately 75% of motifs are conserved (highly similar) and an addi-
the TF is labeled ‘DBD not shared’ with the other species. Thus, tional 5–25% are potentially conserved (ambiguous category) for
there are four possible labels for each TF–species comparison that metazoans and plants, respectively. However, around 900 million
are mutually exclusive. years ago (the origin of metazoans), only around 60% are conserved

986 NatUre Genetics | VOL 51 | JUNE 2019 | 981–989 | [Link]/naturegenetics

NATure GeneTiCS Articles
a b BLAST NN
1−1 Other
D. melanogaster + H. sapiens
D. melanogaster

similar
SR dissimilar

Highly
CG7368 (B1H) ZNF777 (SMiLE−seq)
H. sapiens
l(3)neo38 (B1H) ZNF777 (SMiLE−seq)

Ambiguous
pad (B1H) ZNF468 (RCADE)
Species
vfl (B1H) ZNF682 (SMiLE−seq)

Aef1 (HT-SELEX) ZNF485 (RCADE)

Dissimilar
odd (HT-SELEX) OSR2 (HT-SELEX)

TF family CG4424 (HT-SELEX) ZNF436 (RCADE)

C2H2 ZF
Homeodomain gl (HT-SELEX) ZSCAN22 (RCADE)
Forkhead
MADF

DBD not
CG3065 (B1H) SP4 (HT-SELEX)

shared
Sox
FLYWCH
Myb/SANT Kr (HT-SELEX) BCL6 (PBM)
Pipsqueak
Other CG17802 (HT-SELEX) ZNF35 (RCADE)
0 10 20 30 40 50 0 10 20 30 40 50
sug (HT-SELEX) GLIS2 (HT-SELEX)
TFs (%) TFs (%)
c BLAST NN
Clamp (PBM) ZNF429 (RCADE)

CG3407 (B1H) ZNF558 (RCADE)

D. melanogaster + H. sapiens
mamo (B1H) ZNF121 (RCADE)
SR highly similar
bowl (HT-SELEX) OSR2 (HT-SELEX) CG10654 (HT-SELEX) ZNF324B (ChIP−seq)

sob (HT-SELEX) OSR2 (HT-SELEX) Cf2 (B1H) ZNF713 (HT-SELEX)

CG12605 (HT-SELEX) SCRT1 (HT-SELEX) CG8319 (B1H) ZNF84 (RCADE)

scrt (HT-SELEX) SCRT1 (HT-SELEX) CG14667 (HT-SELEX) ZNF584 (RCADE)

esg (B1H) SNAI2 (PBM) her (HT-SELEX) ZNF513 (RCADE)

Spps (B1H) SP4 (HT-SELEX) hkb (HT-SELEX) SP5 (Transfac)

Sp1 (HT-SELEX) SP8 (HT-SELEX) Sry-beta (Transfac) ZNF764 (ChIP−seq)

cbt (HT-SELEX) KLF11 (PBM) shn (B1H) HIVEP2 (Transfac)

dar1 (B1H) KLF5 (HT-SELEX) jim (B1H) ZNF250 (HT-SELEX)

luna (HT-SELEX) KLF7 (ChIP−seq) CG4360 (HT-SELEX) ZNF879 (RCADE)

CG42741 (B1H) KLF8 (Transfac) su(Hw) (ChIP−chip) ZNF25 (ChIP−seq)

sr (B1H) EGR3 (HT-SELEX) CTCF (ChIP−chip) CTCFL (ChIP−seq)

phol (B1H) YY1 (SMiLE−seq) CG12769 (HT-SELEX) ZNF524 (HT-SELEX)

pho (B1H) YY1 (SMiLE−seq) peb (B1H) RREB1 (SELEX)

sens-2 (B1H) GFI1B (HT-SELEX) lola (B1H) ZBTB20 (HT-SELEX)

sens (HT-SELEX) GFI1B (PBM) ttk (SMiLE−seq) ZBTB17 (Misc)

ab (B1H) ZBTB14 (HT-SELEX)

br (B1H) ZBTB20 (HT-SELEX)

CG12236 (B1H) ZBTB14 (HT-SELEX)

Trl (HT-SELEX) ZBTB37 (HT-SELEX)

D19B (B1H) ZNF778 (RCADE)

fru (B1H) ZBTB49 (HT-SELEX)

crol (B1H) ZNF84 (RCADE)

Fig. 6 | TF motif conservation between human and Drosophila melanogaster. a, Percentage of TFs in human or Drosophila (as indicated) that fall into each
similarity regression motif similarity class. TF conservation is partitioned by whether the TF has a one-to-one ortholog (reciprocal best BLAST hit) or all
other orthology relationships. Colors of the stacked bar plots indicate TF family. b,c, Experimentally determined motifs for individual Drosophila and human
C2H2 ZF TFs, shown in pairs that correspond to the BLASTP best hit (Drosophila query to human database). Reciprocal best BLASTP matches (putative
one-to-one orthologs) are indicated with bidirectional arrows. b, Pairs predicted to be dissimilar by similarity regression are shown. c, Pairs predicted to
be highly similar by similarity regression are [Link]–seq, chromatin immunoprecipitation with high-throughput sequencing; ChIP–chip, ChIP with
DNA microarray; SMiLE–seq, selective microfluidics-based ligand enrichment followed by sequencing; RCADE, recognition code-assisted discovery of
regulatory elements. Silhouettes of each species are displayed, adapted from Phylopic ([Link] under a Creative Commons license (https://
[Link]/licenses/by-sa/3.0/).

or potentially conserved; a similar proportion is obtained for the those that are comparable (that is, DBD families that are present in
origin of fungi (around 1,055 million years ago). Within the plant both), the majority have dissimilar or ambiguous motifs.
kingdom (approximately 1,160 million years ago), only slightly Much of the divergence in motifs occurs in a small number of TF
more motifs are conserved or potentially conserved (around 65%). families (Fig. 5, Supplementary Fig. 9); however, these families have
Across kingdoms (for example, between fungi and metazoan), most a large number of members and are in general already known for
DBDs are not shared24 and are thus not comparable. Even among their lineage-specific expansions: C2H2 ZF in metazoans9, nuclear

NatUre Genetics | VOL 51 | JUNE 2019 | 981–989 | [Link]/naturegenetics 987

Articles NATure GeneTiCS

hormone receptors in nematodes25,26 and Myb proteins in plants27. could, in turn, contribute to the long-term objective of developing
The similarity regression analysis thus underscores DNA sequence accurate recognition codes that do not rely on as much experimen-
specificity as a mode of diversification following duplication of tal data as similarity regression does.
these proteins. However, many other families appear rigid in their Similarity regression can also predict when proteins are
DNA-binding motifs and have presumably diversified in function unlikely to share sequence preferences, thus enabling system-
by other mechanisms (for example, basic-leucine zipper (bZIP) atic examination of the overall degree of trans-regulatory change
and basic helix–loop–helix (bHLH) proteins are able to diversify among eukaryotes. Our analyses lend strong support to the notion
through changes in heterodimerization partners)28,29. that cis-regulatory turnover is accompanied by alterations to
One notable example of C2H2 diversification is counter to a pre- trans-regulators between species, even over relatively short times-
vious claim in the literature, but is supported by extensive experi- cales (<100 million years). TFs with divergent motifs are concen-
mental data. A previous study8 claimed that only a few TF-binding trated in TF families with established patterns of lineage-specific
motifs have diverged in sequence specificity between human and expansions, although changes in one-to-one orthologs do occur
Drosophila. The discrepancy appears to be due to the fact that the as have previously been observed10,12,13. This study provides an
HT-SELEX data in this previous study were highly biased towards extensive analysis of TF sequence specificity for both Cannabis
TF families that have not diversified. In particular, it included only and Aspergillus, and both the outputs of similarity regression and
a small minority of C2H2 ZFs, which represent the largest class of the newly generated data highlight the diversity of DNA-binding
TFs in both species. Similarity regression predicts that the majority motifs in both the plant and fungal lineages. Despite lower diver-
of C2H2 ZF proteins do not have conserved motifs (Fig. 6a), even sity in the specificity residues of individual C2H2 ZF domains in
between TFs that have putative one-to-one orthology as defined fungi, relative to metazoa16, proteins that contain these domains
by BLAST. Experimental data confirm our similarity regression contribute substantially to diversification of motifs in fungi, pre-
predictions (motifs for C2H2 TFs are shown in Fig. 6b,c; motif sumably due to the fact that multiple C2H2 domains can be com-
similarities for all TF comparisons are shown in Supplementary bined in different ways. Myb domains also contribute substantially
Fig. 10). Examples of orthologous C2H2 TFs, even one-to-one to motif divergence in multiple lineages (both plants and fungi).
orthologs, differing substantially in their DNA-binding specificity The GAL4/zinc cluster domain proteins, which have expanded in
are shown in Fig. 6b, illustrating that simple orthology alone can fungi, have largely conserved monomeric binding specificity in
be a poor predictor of shared motifs. As a control, BLAST NNs their DBDs, and are thus more likely contribute to TF diversifica-
predicted by similarity regression to have highly similar motifs tion by alterations in spacing and orientation of dimeric sites as
between human and Drosophila do display more similar motifs homo- or heterodimers33.
in the experimental data than TFs with predicted ambiguous or Similarity regression also confirms the extreme diversity of
dissimilar specificities (Supplementary Fig. 10), even when they motifs in the C2H2 ZF family in metazoa. C2H2 ZFs are the fast-
were obtained using different techniques (primarily high-through- est evolving TF family in the recent human lineage34 and similarity
put SELEX (HT-SELEX)8,30 compared to bacterial one-hybrid regression indicates—and experimental data confirm—that their
(B1H) assays31,32). sequence specificities are largely distinct from those in Drosophila,
even among clear orthologs. Our findings differ from a previ-
Discussion ous conclusion that TF-binding specificities are highly conserved
We anticipate that similarity regression will contribute to under- between Drosophila and mammals8, mainly because the vast major-
standing of TF function in several ways. First, it presents several ity of the C2H2 ZF proteins were absent from the HT-SELEX data
advantages in the task of predicting motifs. Like simple homology in the previous study. This absence could be due to low success
(that is, alignment percent identity), the score it produces serves as rates in HT-SELEX (and other in vitro assays), due to long binding
a confidence measure that can be used to avoid incorrect predic- sites or other factors. Importantly, most C2H2 ZFs are bona fide
tions. At the same time, the increased recall (that is, coverage) of TFs, which bind specific DNA sequences in vivo and/or in vitro7,34.
similarity regression, relative to percent identity, provides a sub- Notably, among Drosophila species, even one-to-one orthologs of
stantial increase in the number of predicted motifs, which are now C2H2 TFs frequently differ in specificity residues and these differ-
included in our update of the Cis-BP database. Similarity regression ences are predicted to influence DNA sequence preferences10. The
can also be adapted to predict new motifs, by combining binding bulk of motif differences occur in TFs with more complex orthol-
data from related proteins. These motifs score favorably relative to ogy patterns, however. In mammals, there is strong evidence that
both a related general-purpose prediction method (affinity regres- the need to recognize new retroelements for silencing by KRAB-
sion), as well as state-of-the-art recognition codes for specific DBD containing C2H2 ZFs has played a part in their evolution35, but the
families (C2H2 ZFs and homeodomains). We do note that the KRAB domain is restricted to tetrapods and it is unclear what forces
motifs for the homeodomain-specific prediction tool PreMoTF5 are driving C2H2 ZF motif diversification in other lineages. Even in
scored more highly than similarity regression among the held-out humans, most C2H2 ZF proteins do not appear to bind to retroele-
homeodomains for which there was no highly similar protein (that ments36 and presumably have other functions.
is, a protein with a high similarity regression score) in the training Knowing the sequence specificities of TFs is an important first
set. Thus, although the use of the similarity regression score as a step in their characterization. Overall, we anticipate that similarity
confidence metric can be seen as an advantage, this outcome also regression and the results it produces will represent a major advance
highlights a disadvantage of making predictions solely on the basis in our understanding of the function and evolution of both TFs and
of similarity among proteins: a well-formulated recognition code gene regulatory mechanisms.
has the potential to make accurate predictions for completely novel
proteins. Unfortunately, such recognition codes do not exist for the Online content
vast majority of DBD types. Any methods, additional references, Nature Research reporting
Second, the weights (that is, coefficients) produced by similar- summaries, source data, statements of code and data availability and
ity regression are often highest for known specificity residues and associated accession codes are available at [Link]
DNA-contacting positions. Thus, unstudied positions with high s41588-019-0411-1.
weights represent candidates for new determinants of TF sequence
specificity. Together with structural data, these weights may also Received: 13 November 2018; Accepted: 4 April 2019;
shed new light on biophysical aspects of DNA binding—which Published online: 27 May 2019

988 NatUre Genetics | VOL 51 | JUNE 2019 | 981–989 | [Link]/naturegenetics

NATure GeneTiCS Articles
References 28. Grove, C. A. et al. A multiparameter network reveals extensive divergence
1. Stormo, G. D. DNA binding sites: representation and discovery. between C. elegans bHLH transcription factors. Cell 138, 314–327 (2009).
Bioinformatics 16, 16–23 (2000). 29. Reinke, A. W., Baek, J., Ashenberg, O. & Keating, A. E. Networks of bZIP
2. Mathelier, A. et al. JASPAR 2016: a major expansion and update of the protein–protein interactions diversified over a billion years of evolution.
open-access database of transcription factor binding profiles. Nucleic Acids Science 340, 730–734 (2013).
Res. 44, D110–D115 (2016). 30. Jolma, A. et al. Multiplexed massively parallel SELEX for characterization
3. Weirauch, M. T. et al. Determination and inference of eukaryotic of human transcription factor binding specificities. Genome Res. 20,
transcription factor sequence specificity. Cell 158, 1431–1443 (2014). 861–873 (2010).
4. Pelossof, R. et al. Affinity regression predicts the recognition code of nucleic 31. Noyes, M. B. et al. A systematic characterization of factors that regulate
acid-binding proteins. Nat. Biotechnol. 33, 1242–1249 (2015). Drosophila segmentation via a bacterial one-hybrid system. Nucleic Acids Res.
5. Christensen, R. G. et al. Recognition models to predict DNA-binding 36, 2547–2560 (2008).
specificities of homeodomain proteins. Bioinformatics 28, i84–i89 (2012). 32. Zhu, L. J. et al. FlyFactorSurvey: a database of Drosophila transcription factor
6. Persikov, A. V. et al. A systematic survey of the Cys2His2 zinc finger binding specificities determined using the bacterial one-hybrid system.
DNA-binding landscape. Nucleic Acids Res. 43, 1965–1984 (2015). Nucleic Acids Res. 39, D111–D117 (2011).
7. Najafabadi, H. S. et al. C2H2 zinc finger proteins greatly expand the human 33. MacPherson, S., Larochelle, M. & Turcotte, B. A fungal family of
regulatory lexicon. Nat. Biotechnol. 33, 555–562 (2015). transcriptional regulators: the zinc cluster proteins. Microbiol. Mol. Biol. Rev.
8. Nitta, K. R. et al. Conservation of transcription factor binding specificities 70, 583–604 (2006).
across 600 million years of bilateria evolution. eLife 4, e04837 (2015). 34. Lambert, S. A. et al. The human transcription factors. Cell 175,
9. Liu, H., Chang, L. H., Sun, Y., Lu, X. & Stubbs, L. Deep vertebrate roots for 598–599 (2018).
mammalian zinc finger transcription factor subfamilies. Genome Biol. Evol. 6, 35. Ecco, G., Imbeault, M. & Trono, D. KRAB zinc finger proteins. Development
510–525 (2014). 144, 2719–2729 (2017).
10. Nadimpalli, S., Persikov, A. V. & Singh, M. Pervasive variation of 36. Schmitges, F. W. et al. Multiparameter functional diversity of human C2H2
transcription factor orthologs contributes to regulatory network evolution. zinc finger proteins. Genome Res. 26, 1742–1752 (2016).
PLoS Genet. 11, e1005011 (2015). 37. Noyes, M. B. et al. Analysis of homeodomain specificities allows the
11. Lynch, V. J. & Wagner, G. P. Resurrecting the role of transcription factor family-wide prediction of preferred recognition sites. Cell 133,
change in developmental evolution. Evolution 62, 2131–2154 (2008). 1277–1289 (2008).
12. Baker, C. R., Tuch, B. B. & Johnson, A. D. Extensive DNA-binding specificity 47. Finn, R. D. et al. The Pfam protein families database. Nucleic Acids Res. 38,
divergence of a conserved transcription regulator. Proc. Natl Acad. Sci. USA D211–D222 (2010).
108, 7493–7498 (2011).
13. Sayou, C. et al. A promiscuous intermediate underlies the evolution of Acknowledgements
LEAFY DNA binding specificity. Science 343, 645–648 (2014). We thank Xiaoting Chen and Mario Pujato for computational support. S.A.L. was funded
14. Morgunova, E. et al. Structural insights into the DNA-binding specificity of by a Natural Sciences and Engineering Research Council of Canada Doctoral Fellowship.
E2F family transcription factors. Nat. Commun. 6, 10050 (2015). T.R.H. holds the Billes Chair of Medical Research at the University of Toronto. This
15. McKeown, A. N. et al. Evolution of DNA specificity in a transcription factor work was supported by a Canadian Institutes of Health Research grant (FDN-148403)
family produced a new gene regulatory module. Cell 159, 58–68 (2014). and a Natural Sciences and Engineering Research Council of Canada grant (RPGIN-
16. Najafabadi, H. S. et al. Non-base-contacting residues enable kaleidoscopic 2016-05643) to T.R.H., National Institutes of Health (NIH) grants R01 AR073228, R01
evolution of metazoan C2H2 zinc finger DNA binding. Genome Biol. 18, NS099068 and R01 GM055479, Lupus Research Alliance ‘Novel Approaches’, CCRF
167 (2017). Endowed Scholar and CCHMC CpG Award 53553 to M.T.W. and a Canadian Institutes
17. Berger, M. F. et al. Compact, universal DNA microarrays to comprehensively of Health Research Operating grant (MOP-125894) to Q.D.M. and T.R.H.
determine transcription-factor binding site specificities. Nat. Biotechnol. 24,
1429–1435 (2006).
18. Weirauch, M. T. et al. Evaluation of methods for modeling transcription Author contributions
factor sequence specificity. Nat. Biotechnol. 31, 126–134 (2013). S.A.L., M.T.W. and T.R.H. conceived the study and oversaw it to completion. S.A.L.
19. Love, J. J. et al. Structural basis for DNA bending by the architectural analyzed the data, made the figures and performed all computational analyses except for
transcription factor LEF-1. Nature 376, 791–795 (1995). experiments for which A.S. reimplemented the affinity regression pipeline and applied
20. Marmorstein, R., Carey, M., Ptashne, M. & Harrison, S. C. DNA recognition it to new data. Q.D.M. guided the computational and statistical analyses. M.A., S.A.L.
by GAL4: structure of a protein–DNA complex. Nature 356, 408–414 (1992). and M.T.W. maintained and updated the Cis-BP database. G.C. and M.X.C. produced
21. King, D. A., Zhang, L., Guarente, L. & Marmorstein, R. Structure of a the clones for Aspergillus PBM experiments. A.W.H.Y. produced the remainder of the
HAP1–DNA complex reveals dramatically asymmetric DNA binding by a clones and performed all PBM experiments. S.A.L. and T.R.H. wrote the manuscript with
homodimeric protein. Nat. Struct. Biol. 6, 64–71 (1999). feedback and approval from all authors.
22. Persikov, A. V. & Singh, M. De novo prediction of DNA-binding specificities
for Cys2His2 zinc finger proteins. Nucleic Acids Res. 42, 97–108 (2014). Competing interests
23. Gupta, A. et al. An improved predictive recognition model for Cys2-His2 The authors declare no competing interests.
zinc finger proteins. Nucleic Acids Res. 42, 4800–4812 (2014).
24. de Mendoza, A. et al. Transcription factor evolution in eukaryotes and the
assembly of the regulatory toolkit in multicellular lineages. Proc. Natl Acad. Additional information
Sci. USA 110, E4858–E4866 (2013). Supplementary information is available for this paper at [Link]
25. Narasimhan, K. et al. Mapping and analysis of Caenorhabditis elegans s41588-019-0411-1.
transcription factor sequence specificities. eLife 4, e06967 (2015). Reprints and permissions information is available at [Link]/reprints.
26. Robinson-Rechavi, M., Maina, C. V., Gissendanner, C. R., Laudet, V. & Correspondence and requests for materials should be addressed to T.R.H.
Sluder, A. Explosive lineage-specific expansion of the orphan nuclear receptor
HNF4 in nematodes. J. Mol. Evol. 60, 577–586 (2005). Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in
27. Stracke, R., Werber, M. & Weisshaar, B. The R2R3-MYB gene family in published maps and institutional affiliations.
Arabidopsis thaliana. Curr. Opin. Plant Biol. 4, 447–456 (2001). © The Author(s), under exclusive licence to Springer Nature America, Inc. 2019

NatUre Genetics | VOL 51 | JUNE 2019 | 981–989 | [Link]/naturegenetics 989

Articles NATure GeneTiCS

Methods Comparing similarity regression weights with known DNA-contacting residues.

Similarity regression. Similarity regression is formulated as a regression We used the DNAproDB database44 to compare the similarity regression weights
task in which the dependent variable (Y) is a metric of similarity in DNA with known protein–DNA contacts. DNAproDB catalogues DNA–protein
sequence specificity between pairs of proteins (see below) and the independent complexes present in the Protein Data Bank45, annotating the amino acid residues
variables (the feature vector X) are the identity and similarity of amino acid that contact the DNA backbone and bases in the major and minor grooves. We
residues at each individual position of the aligned DBDs, for the same pairs transferred these annotations to our models by first extracting all of the protein
of proteins. To make the alignment, each instance of a DBD is aligned to its sequences in DNAproDB and identified DBDs using hmmscan46 and the same
corresponding Pfam HMM using the semi-global method implemented in aphid38, Pfam HMM models47 and thresholds as the Cis-BP database. We then parsed
recording match positions (that is, positions that are present in the HMM). An the nucleotide–residue interactions for each structure into backbone, major and
example alignment of two homeodomain sequences is presented in Fig. 1a. At each minor groove interactions (using DNAproDB-recommended buried solvent
position of the aligned sequences, either identity (as binary values) or similarity accessible surface area, hydrogen bond and van der Waals interaction thresholds)
(BLOSUM62 substitution score39) is recorded (Fig. 1b), yielding the feature vector and associated them with the position of the residue in the DBD alignment. We
for each TF pair. For TF families that have DBDs present in arrays (mainly C2H2 represented the interactions as a contact frequency for each type of DNA contact,
ZFs and Myb) the best ungapped and overlapping pairwise alignment of DBD by normalizing the number of nucleotide–residue interactions that occurred in
arrays (Supplementary Fig. 2a) is found by selecting the alignment offset with the each position of the DBD by the number of protein–DNA structures that contained
maximum amino acid identity. For a multi-DBD alignment, the feature vector is that DBD. Correspondence between similarity regression weights and the three
generated by the average score (identity or similarity) in each position of the DBD classes of DNA contacts were evaluated using partial Pearson correlations, which
alignment from all DBD arrays, normalizing by the DBD length of the longest assessed the correlation between each contact type after removing the effects of the
protein (Supplementary Fig. 2b). other two contacts on the similarity regression weights.
In the analyses described, the metric of similarity in DNA sequence
specificity between pairs of proteins (Y) is calculated from the 8-base oligomer Comparison of similarity regression with affinity regression. Affinity
PBM data as the fraction of high-scoring 8-base oligomers (E > 0.45) that are regression predicts Z-scores of DNA 8-base oligomers from short peptides in the
shared between two TFs (that is intersection/union for two experiments, referred protein sequence. Here, we implemented a soft-coded Python version of affinity
to as ‘E-score overlap’). For each TF family, E-score overlaps that exceed the regression, ensuring similar performance on the previously published, original
twenty-fifth percentile of experimental replicates (the same threshold was used data4 and using identical constructions of the protein and DNA features. A single
previously3) are taken as having ‘highly similar’ specificities. The highly similar affinity regression model for each TF family was trained using the same data as
labels are not only used as positives for training logistic similarity regression the corresponding similarity regression model and the number of informative
models (see next paragraph), but also for evaluating the performance of components selected after dimensionality reduction was set to capture 90% of the
similarity regression (for example, precision/recall analysis). E-score overlaps weight of each singular value. To predict the Z-scores of uncharacterized/tested
that are less than 0.2 are taken as having ‘dissimilar’ specificities, allowing some transcription factors, affinity regression determines their protein k-mer vectors to
overlap because many families bind a characteristic sequence while their highest- predict the similarities of the held-out TF to all characterized protein profiles in the
affinity 8-base oligomers differ. The dissimilar labels are used as negatives to training set. Affinity regression uses these similarities to reconstruct the Z-score
define the score threshold below which TFs are unlikely to share specificities in profiles by a similarity-weighted sum of Z-score profiles using either the NN or
a NPV analysis. top 5 NNs and a geometrical reconstruction from the span of the training vectors,
For each TF family with sufficient PBM data, we trained four similarity proposed and applied in the previously published study that describes affinity
regression models that varied in the representation of protein similarity (identity or regression4. Affinity regression was applied to the new TFs that were present in the
BLOSUM substitution score) and in the representation of the data (either linear or PBM data from this study.
logistic regression models). Each regression model is trained in R40 using glmnet41 We used three means to predict the 8-base oligomer Z-score profile for each
constrained to fit positive regression coefficients, selecting the optimal ridge (L2) held-out TF using similarity regression by first copying the Z-score profile from
regularization strength using cross-validation. Because the data consist of pairs, the protein with the highest similarity regression score (that is, the single NN),
normal k-fold cross-validation is invalid, as random training and test splits would after which the Z-score profiles of the five proteins with the highest similarity
not be independent. To solve this problem, we train the models using leave-one- regression scores (top 5) were combined, weighting the Z-scores for each of
TF-out cross-validation (testing on data points made from the comparisons of a the five by the corresponding similarity regression or percent identity score
single TF), which we implemented using the caret package42. This performance and by then combining the Z-scores from all TFs in the training set that are
measure can be interpreted as how well a similarity regression model generalizes to predicted by Similarity regression to have highly similar specificities (similarity
unseen TFs and is used to select the optimal regularization parameters and score regression highly similar), weighting the Z-scores for each of these values by the
thresholds for each regression model. corresponding similarity regression score. We evaluated the accuracy of similarity
An outline of the similarity regression model generation and selection for regression, percent identity and affinity regression predictions using the Pearson
homeodomain TFs is presented in Supplementary Fig. 1. First, the optimal correlation coefficient between the predicted Z-score profile and the experimental
regularization strength is selected using the cross-validation procedure Z-scores. We used two-sided paired Wilcoxon signed-rank tests to identify
implemented in caret, yielding a selected model for each feature–output significant differences in mean Pearson correlation coefficient ranks between
combination. For each regression model, and the percent identity method, two Z-score reconstruction methods.
thresholds are derived to predict TFs with highly similar or dissimilar specificities.
To select these thresholds, the predictions on held-out data from each cross- Comparison of similarity regression with recognition codes. We used ‘PWM_
validation fold are combined and compared with their known TF similarity labels. align’18 to convert the similarity regression and percent identity-predicted Z-score
To identify TFs with highly similar sequence specificities (E-score overlap > the TF profiles (described above) into PWMs. We used published webservers with default
family replicate threshold (the twenty-fifth percentile of experimental replicates as settings to obtain predicted motifs (homeodomains: PreMoTF5, [Link]
previously described3)), a precision/recall curve is generated on the held-out data [Link]/PreMoTF/); C2H2 predictions, linear-expanded support vector machine
and a score threshold is selected from the curve such that it yields 75% precision method22 ([Link] and ZFModels5 ([Link]
(a heuristic identical to the one used in the previous study3). A threshold for ZFModels/). We measured motif similarity between the predicted motifs and
dissimilar specificities is derived by finding a NPV cut-off that classifies 95% of TFs experimentally determined motifs using MoSBAT energy scores48 (parameters:
below that score threshold as having truly dissimilar specificities (E-score overlap N = 100,000, L = 25 nt). We used two-sided paired Wilcoxon signed-rank tests to
of <0.2). For each threshold, the recall of positive and negative predictions was identify significant differences in mean MoSBAT energy score ranks between motif
recorded to evaluate the improvement of similarity regression models over percent prediction methods.
identity. If thresholds could not be derived for a TF family a global threshold of
70% identity for predicting highly similar TFs (identical to our previous study3) Updates to the Cis-BP database. We performed extensive updates to the Cis-BP
and a threshold of 25% identity for predicting dissimilar TFs (selected to yield database, encompassing changes to both the data and the methodologies. Build
95% NPV) were derived. The highly similar and dissimilar thresholds are then 2.0 of Cis-BP now contains data for 741 species (increased from 340) ([Link]
applied to the predictions to classify each TF pair in the held-out data as having [Link]/). In addition to adding new species, updated genome builds
highly similar, ambiguous or dissimilar specificities for each similarity regression were incorporated for all existing species, where available. Each of these updates
model (and for the percent identity method). The best similarity regression model includes the latest available protein sequences, protein and gene identifiers, gene
is then selected by comparing the three-class predictions to ground-truth labels names and gene aliases. Furthermore, the set of human TFs contained in Cis-BP
and selecting the model with the best Matthews correlation coefficient, a metric of now matches the set of 1,639 curated TFs provided in a recently published review34.
multi-class classification accuracy that is sensitive to class imbalance43. This process DBD scans were performed using updated Pfam HMM models47, including
yields a single final similarity regression model for each TF family, composed of a models for EBF1 (COE1_DBD), FLYWCH and ICP4 (Herpes_ICP4_N). We
weight vector (that is, coefficients for X values, which are the selected measure of also removed models for DP and SART-1, which are now known to not bind to
protein similarity), as well as two thresholds for the dependent variable (Y) that DNA with specificity. A total of 1,358 new motifs were obtained from 38 different
are used to predict whether two TFs have highly similar, ambiguous or dissimilar sources, including 541 HT-SELEX motifs obtained for human TFs from methylated
sequence specificities. and unmethylated DNA49, 534 DNA affinity purification sequencing (DAP-seq)

NatUre Genetics | [Link]/naturegenetics

NATure GeneTiCS Articles
motifs for Arabidopsis thaliana50, 248 HT-SELEX D. melanogaster motifs8 and 221 Reporting Summary. Further information on research design is available in the
ChIP-exo (ChIP–seq with increased resolution from exonuclease treatment) and Nature Research Reporting Summary linked to this article.
ChIP–seq-derived C2H2 ZF motifs51. Existing motif sources such as UNIPROBE52,
Transfac53, JASPAR54 and HOCOMOCO55 were also updated to include data from Data availability
the latest database builds. In addition to these improvements in the database New PBM data and motifs are deposited in GEO (accession number GSE121420)
contents, this update of Cis-BP incorporates several methodological advances. and the Cis-BP database (v.2.0; [Link]
First, when two predicted DBDs overlap in a given protein, only the DBD with
the most significant HMMER P value is retained. Second, matches to the Pfam
Myb/SANT domain are now further subclassified into Myb (which binds to DNA Code availability
specifically and also contains Myb-like sequences that are also likely to bind DNA) The Similarity Regression code, and examples, are available on GitHub (https://
or SANT (which does not bind to DNA specifically). In brief, we scored each Myb/ [Link]/smlmbrt/SimilarityRegression).
SANT domain with the Myb (PS51294), Myb-like (PS50090) and SANT (PS51293)
specific PROSITE56 models and annotated domains by the profile with the highest References
score. This procedure is now applied to remove SANT-only containing proteins 38. Wilkinson, S. P. aphid: an R package for analysis with profile hidden
(which are not TFs) and remove SANT domains from proteins that contain both Markov models. Bioinformatics [Link]
Myb and SANT domains. Third, we removed one-to-one orthologs (reciprocal best (2019).
BLAST hits) of metazoan proteins with false-positive human TFs derived from 39. Henikoff, S. & Henikoff, J. G. Amino acid substitution matrices from protein
a recent curation effort34. Finally, motif inferences in Cis-BP are now performed blocks. Proc. Natl Acad. Sci. USA 89, 10915–10919 (1992).
using the similarity regression approach described in this manuscript, as opposed 40. R Core Team. R: A Language and Environment for Statistical Computing (R
to the original method, which was based on percent identity. Foundation for Statistical Computing, 2013);[Link]
41. Friedman, J., Hastie, T. & Tibshirani, R. Regularization paths for generalized
Predicting TF motif conservation across species. To evaluate motif conservation linear models via coordinate descent. J. Stat. Softw. 33, 1–22 (2010).
between species, we used the TF annotations and DBD sequences from Cis-BP 42. Kuhn, M. Building predictive models in R using the caret package. J. Stat.
(version 2.0). For each pair of species analyzed, we used similarity regression to Softw. 28, 1–26 (2008).
predict TF similarity for all pairs of TFs from the same TF family. To calculate the 43. Gorodkin, J. Comparing two K-category assignments by a K-category
conservation of each TF in each species relative to a second species, we report the correlation coefficient. Comput. Biol. Chem. 28, 367–374 (2004).
maximum similarity regression score among all TFs in the second species, and the 44. Sagendorf, J. M., Berman, H. M. & Rohs, R. DNAproDB: an interactive tool
resulting similarity classification. If the TF was from a family that is not shared for structural analysis of DNA–protein complexes. Nucleic Acids Res. 45,
between species (for example, DBD families that are clade-specific), we assume W89–W97 (2017).
that the motif is not conserved, and report the TF as uncomparable with the label 45. Berman, H. M. et al. The protein data bank. Nucleic Acids Res. 28,
‘DBD not shared’. We obtained the time to the last common ancestor (divergence 235–242 (2000).
time) from the TimeTree database57. 46. HMMER: biosequence analysis using profile hidden Markov models (Howard
To identify the most similar proteins between human and Drosophila, we used Hughes Medical Institute, 2015); [Link]
BLASTP58 with default settings, using full-length TF sequences present in Cis-BP. 48. Lambert, S. A., Albu, M., Hughes, T. R. & Najafabadi, H. S. Motif comparison
The closest TF in each species (BLAST NN) was identified using the minimum based on similarity of binding affinity profiles. Bioinformatics 32,
E value and reciprocal best BLAST NNs (putative one-to-one orthologs) were 3504–3506 (2016).
recorded. We measured motif similarity between BLAST NNs with experimentally 49. Yin, Y. et al. Impact of cytosine methylation on DNA binding specificities of
determined motifs using MoSBAT energy scores48 (parameters: N = 200,000, human transcription factors. Science 356, eaaj2239 (2017).
L = 50 nt). 50. O’Malley, R. C. et al. Cistrome and epicistrome features shape the regulatory
DNA landscape. Cell 165, 1280–1292 (2016).
DBD cloning. In total, 350 novel A. nidulans TF DBDs were selected for analysis 51. Barazandeh, M., Lambert, S. A., Albu, M. & Hughes, T. R. Comparison of
and 180 were successfully cloned into the expression vector (pTH6838) and ChIP-seq data and a reference motif set for human KRAB C2H2 zinc finger
validated by sequencing. These were cloned using RNA that was extracted from the proteins. G3 (Bethesda) 8, 219–229 (2018).
wild-type A. nidulans strain (FGSC A4). cDNA was generated by RT–PCR using 52. Hume, M. A., Barrera, L. A., Gisselbrecht, S. S. & Bulyk, M. L. UniPROBE,
random hexamer primers. Proof-reading KOD Hot Start DNA polymerase was update 2015: new tools and content for the online database of protein-
used to amplify the DBD-coding region and flanking regions up to 50 amino acids binding microarray data on protein–DNA interactions. Nucleic Acids Res. 43,
long and products were extracted from a 1% agarose gel using a Silica Bead DNA D117–D122 (2015).
Gel Extraction Kit (Thermo Fisher Scientific). Double digests were performed 53. Matys, V. et al. TRANSFAC and its module TRANSCompel:
using the restriction endonucleases AscI (10 U μl−1) (Thermo Fisher Scientific) and transcriptional gene regulation in eukaryotes. Nucleic Acids Res. 34,
SbfI-HF (20 U μl−1) (New England Biolabs). The fragments were ligated into the D108–D110 (2006).
expression vector using T4 DNA ligase (New England Biolabs). Constructs were 54. Khan, A. et al. JASPAR 2018: update of the open-access database of
verified by Sanger sequencing (GATC Biotech). Other DBDs were cloned using transcription factor binding profiles and its web framework. Nucleic Acids
previously reported procedures3. Res. 46, D1284 (2018).
55. Kulakovskiy, I. V. et al. HOCOMOCO: towards a complete collection of
PBMs. PBM laboratory methods were performed as described previously18,59. Each transcription factor binding models for human and mouse via large-scale
DBD-encoding plasmid was analyzed in duplicate on two different arrays with ChIP-seq analysis. Nucleic Acids Res. 46, D252–D259 (2018).
differing probe sequences. The 8-base oligomer Z- and E-scores were calculated 56. Sigrist, C. J. et al. PROSITE: a documented database using patterns and
as previously described17. We deemed experiments successful if at least one 8-base profiles as motif descriptors. Brief. Bioinform. 3, 265–274 (2002).
oligomer had E > 0.45 on both arrays, the complementary arrays produced highly 57. Kumar, S., Stecher, G., Suleski, M. & Hedges, S. B. Timetree: a resource for
correlated E- and Z-scores and yielded similar PWMs based on the PWM_align timelines, timetrees, and divergence times. Mol. Biol. Evol. 34,
algorithm18. Motifs shown (and deposited in Cis-BP) for each TF are chosen by 1812–1819 (2017).
cross-replicate evaluation of three motif derivation methods (PWM_align, PWM_ 58. Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local
align_Z and BEEML-PBM)3,60. alignment search tool. J. Mol. Biol. 215, 403–410 (1990).
59. Lam, K. N., van Bakel, H., Cote, A. G., van der Ven, A. & Hughes, T. R.
Statistics and experimental design. Two-sided paired Wilcoxon signed-rank tests Sequence specificity is obtained from the majority of modular C2H2
were used to identify significant differences between motif prediction methods. zinc-finger arrays. Nucleic Acids Res. 39, 4680–4690 (2011).
Distributions were summarized with box plots where appropriate and described in 60. Zhao, Y. & Stormo, G. D. Quantitative analysis demonstrates most
the relevant figure legends. Additional details of the experimental design and data transcription factors require only simple models of specificity. Nat. Biotechnol.
are included in the Nature Research Reporting Summary. 29, 480–483 (2011).

NatUre Genetics | [Link]/naturegenetics

 nature research | reporting summary
Corresponding author(s): Hughes, TR

Reporting Summary
Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency
in reporting. For further information on Nature Research policies, see Authors & Referees and the Editorial Policy Checklist.

Statistical parameters
When statistical analyses are reported, confirm that the following items are present in the relevant location (e.g. figure legend, table legend, main
text, or Methods section).
n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement
An indication of whether measurements were taken from distinct samples or whether the same sample was measured repeatedly
The statistical test(s) used AND whether they are one- or two-sided
Only common tests should be described solely by name; describe more complex techniques in the Methods section.

A description of all covariates tested

A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons
A full description of the statistics including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND
variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals)

For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted
Give P values as exact values whenever suitable.

For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings
For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes
Estimates of effect sizes (e.g. Cohen's d, Pearson's r), indicating how they were calculated

Clearly defined error bars

State explicitly what error bars represent (e.g. SD, SE, CI)

Our web collection on statistics for biologists may be useful.

Software and code

Policy information about availability of computer code
Data collection Universal Protein Binding Microarray (PBM) Analysis Suite ([Link]
python Dependancies: pandas, biopython
R Dependancies (packages at: [Link] caret (v6), glmnet (v 2.0-13), PRROC (v1.3), aphid (v 1.0.1), seqinr (v3.4.5)

Data analysis The SR code, and examples, are made available on GitHub ([Link]
hmmscan (v3.2.1, part of the HMMER toolkit, [Link]
Affinity Regression ([Link]
PreMoTF ([Link]
ZFModels ([Link]
C2H2 SVM Method ([Link]
April 2018

For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors/reviewers
upon request. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information.

1
Data

nature research | reporting summary

Policy information about availability of data
All manuscripts must include a data availability statement. This statement should provide the following information, where applicable:
- Accession codes, unique identifiers, or web links for publicly available datasets
- A list of figures that have associated raw data
- A description of any restrictions on data availability
New PBM data and motifs are deposited in GEO (accession number: GSE121420), and Cis-BP ([Link]

Field-specific reporting
Please select the best fit for your research. If you are not sure, read the appropriate sections before making your selection.
Life sciences Behavioural & social sciences Ecological, evolutionary & environmental sciences
For a reference copy of the document with all sections, see [Link]/authors/policies/[Link]

Life sciences study design

All studies must disclose on these points even when the disclosure is negative.
Sample size We successfully generated new PBM data for 340 TFs representing multiple eukaryotic kingdoms, with a particular focus on Cannabis sativa (a
medicinal plant), Caenorhabditis briggsae (a nematode), Aspergillus nidulans and Neurospora crassa (model fungi), and also 15 human TFs.
These TFs were selected on the basis of at least one of two different criteria: first, to increase the number of experimentally determined
motifs for TFs in these species of interest, and second, to obtain novel motifs by analyzing proteins that are dissimilar to TFs with known
motifs. For A. nidulans TF 350 novel DBDs were selected for analysis and 180 were successfully cloned into the expression vector (pTH6838)
and validated by sequencing. Sample size was not predetermined, but is large enough to test the SR and other prediction algorithms on
unseen data from relevant/highly-abundant TF families.

Data exclusions No data passing our PBM success criteria was excluded from the analysis.

Replication Successful replication of SR model accuracy was evaluated using the new PBM data, and motif similarity calculated from other motifs included
in the human vs. fly analysis.

Randomization Randomization was not necessary in this study as samples were allocated to different groups (TF similarity) based on uniform rules learned
during cross-validation.

Blinding Blinding was not necessary as SR models are developed automatically using labeled training data based on established thresholds.
Performance metrics were calculated without human intervention.

Reporting for specific materials, systems and methods

Materials & experimental systems Methods

n/a Involved in the study n/a Involved in the study
Unique biological materials ChIP-seq
Antibodies Flow cytometry
Eukaryotic cell lines MRI-based neuroimaging
Palaeontology
Animals and other organisms
Human research participants
April 2018

Unique biological materials

Policy information about availability of materials
Obtaining unique materials TF clones used for PBMs are available upon direct request to the corresponding author.