Computational Identification of Diverse Mechanisms Underlying Transcription Factor-DNA Occupancy

Summary

ChIP-based genome-wide assays of transcription factor (TF) occupancy have emerged as a powerful, high-throughput method to understand transcriptional regulation, especially on a global scale. This has led to great interest in the underlying biochemical mechanisms that direct TF-DNA binding, with the ultimate goal of computationally predicting a TF's occupancy profile in any cellular condition. In this study, we examined the influence of various potential determinants of TF-DNA binding on a much larger scale than previously undertaken. We used a thermodynamics-based model of TF-DNA binding, called “STAP,” to analyze 45 TF-ChIP data sets from Drosophila embryonic development. We built a cross-validation framework that compares a baseline model, based on the ChIP'ed (“primary”) TF's motif, to more complex models where binding by secondary TFs is hypothesized to influence the primary TF's occupancy. Candidates interacting TFs were chosen based on RNA-SEQ expression data from the time point of the ChIP experiment. We found widespread evidence of both cooperative and antagonistic effects by secondary TFs, and explicitly quantified these effects. We were able to identify multiple classes of interactions, including (1) long-range interactions between primary and secondary motifs (separated by ≤150 bp), suggestive of indirect effects such as chromatin remodeling, (2) short-range interactions with specific inter-site spacing biases, suggestive of direct physical interactions, and (3) overlapping binding sites suggesting competitive binding. Furthermore, by factoring out the previously reported strong correlation between TF occupancy and DNA accessibility, we were able to categorize the effects into those that are likely to be mediated by the secondary TF's effect on local accessibility and those that utilize accessibility-independent mechanisms. Finally, we conducted in vitro pull-down assays to test model-based predictions of short-range cooperative interactions, and found that seven of the eight TF pairs tested physically interact and that some of these interactions mediate cooperative binding to DNA.

Introduction

A major challenge in the analysis of genomic sequences is the annotation of cis-regulatory elements. Significant progress has been made towards this goal through high throughput methods such as ChIP-chip and ChIP-SEQ that describe the locations where specific transcription factors (TFs) bind to the genome in vivo [1]–[3]. ChIP-based characterization of TF binding profiles can help elucidate specific regulatory interactions between TFs and genes [4]. A number of genome-wide ChIP data sets, corresponding to diverse TFs and cellular conditions, have been generated through the efforts of various laboratories and consortia [1], [3]. Such data sets also offer the opportunity to apply computational and statistical methods to understand the determinants of TF-DNA binding at a quantitative level [5]–[7]. Given the central role of the TF-DNA binding process in the regulatory activity of a TF, such an understanding can provide a holistic view of transcriptional regulation and also set the stage for future computational methods for predicting cell type-specific TF-binding profiles.

The most extensively studied determinant of TF occupancy is the DNA binding specificity of the TF. Various experimental approaches [8]–[11] have been successful in obtaining motifs representing the diversity and relative affinities of DNA sequences bound by an individual TF. An initial expectation is that a TF's motif will allow prediction of its binding levels genome-wide. On the other hand, it is clear that interactions with other TFs can significantly influence binding to regulatory sequences. For example, interaction of Hox proteins with a cofactor results in greater DNA binding specificity [12] and the tramtrack (TTK) protein can regulate transcription independent of its own DNA binding domain via its interaction with the Trithorax-like (TRL, also known as GAGA binding factor) [13]. Furthermore, TF occupancy at a genomic location in a given cell type also depends on the concentration of that TF, as well as the motifs and concentrations of other TFs that might facilitate or inhibit DNA-binding at the location [14].

A number of recent studies have used genome-wide datasets to characterize parameters that correlate with TF occupancy. In several studies, genome-wide measurements of in vivo DNA accessibility were tested for the ability to help describe TF ChIP data. These studies clearly demonstrate that TF occupancy has a close relationship with in vivo DNA accessibility [6], [7], with both factors likely influencing each other [6], [15]–[19]. While these studies reveal that experimental analysis of accessibility can improve modeling of ChIP data, they do not reveal the underlying genomic sequence features that contribute to accessibility. In another study [5], sequence motifs experimentally and computationally identified in Drosophila were shown to contribute to context-specific TF occupancy. Application of discriminative motif analysis to a TF assayed across multiple conditions can successfully identify predictive motifs associated with context-specific binding. However, whether TFs bound to these discriminative motifs contribute to occupancy by direct interaction with the primary TF, accessibility or other mechanisms is not assessed.

In this work, we test the influence of various potential sequence determinants of in vivo TF-DNA binding – the TF's binding motif, as well as the positive or negative influence of other TFs binding in the vicinity – on each of 45 TF-ChIP data sets in Drosophila. For this analysis, we took advantage of over 300 distinct DNA binding specificity motifs determined for individual TFs [20], which encompasses approximately 40% of all predicted Drosophila TFs, and relied upon stage-specific whole-genome RNA-SEQ data [21] to determine which secondary TFs are expressed at the time of the ChIP experiment. We follow the general framework proposed by Kaplan et al. [6], which involves: (1) building computational models that predict TF binding at a location, and (2) assessing how well a baseline model that only uses the “primary motif” (i.e., binding motif of the “ChIP'ed” TF) fits the ChIP data, as compared with more complex models that incorporate additional determinants such as motifs for additional secondary TFs (i.e., TFs other than the ChIP'ed TF). We use the biophysical model STAP [22] to perform these tests. Improvements in the goodness-of-fit measure are evaluated statistically, and a cross-validation framework is adopted to account for differing model complexity in the comparisons. We evaluate each potential determinant separately in order to limit the number of free parameters in the models. For each identified secondary TF, we performed statistical tests to categorize the mechanistic basis of its contribution. In particular, we asked if a secondary motif's influence is likely to be (a) through long-range (≤150 bp) versus short-range (≤30 bp) interactions with the primary motif, (b) through synergistic or antagonistic interactions, and (c) through modulation of local DNA accessibility or direct interactions between TFs.

We find widespread evidence of the effect of secondary TFs on the primary (ChIP'ed) TF's binding levels, including both enhanced occupancy (“cooperativity”) and reduced occupancy (“antagonism”). Cooperative and antagonistic influences of secondary motifs can act through: 1) long-range interactions between primary and secondary motifs, suggesting indirect effects such as chromatin remodeling, 2) short-range interactions with specific inter-site spacing biases, suggesting a direct association, or 3) through overlapping binding sites, suggesting competition for site occupancy.

Two types of experimental evidence support our computational assignments of secondary TFs that influence occupancy via local chromatin architecture or cooperative DNA binding. Extending previous observations [6], [7], we find that DNA accessibility is the primary genomic feature correlated with TF occupancy across the majority of the 45 data sets examined here. We then use accessibility data to re-examine secondary TF motifs that improved prediction of ChIP data in our accessibility-agnostic analysis. We identify several secondary motifs whose contribution is reduced or lost when accessibility information is part of the model, suggesting that the secondary TF influences binding mainly by modulating accessibility patterns. The TFs vielfaltig (VFL, also known as Zelda) and TRL (also known as GAGA factor) appear to synergistically influence the binding of several primary TFs in early and mid-stage embryonic development respectively. Interestingly, the influence of VFL is sometimes imposed through accessibility, while in other cases it is independent of accessibility. In contrast, the influence of TRL is imposed exclusively through accessibility. The TF motifs for extradenticle (EXD), retained (RETN), jing interacting gene regulatory 1 (JIGR1) and homothorax (HTH) commonly antagonize TF occupancy through accessibility-mediated and accessibility-independent mechanisms. We find many cases where the secondary motif's influence remains significant upon accounting for accessibility, thus suggesting alternative mechanisms such as cooperative or antagonistic DNA-binding by the primary and secondary TFs. We identify eight examples where the arrangement of primary and secondary motifs implies cooperative binding via physical interaction, and demonstrate that for all but one of these cases the TFs do, in fact, directly interact in vitro and that several bind cooperatively in vitro to sequences that are occupied in vivo. Overall, our analysis demonstrates that a biophysical model for the combinatorial action of primary and secondary TFs used with an extensive collection of binding motifs for known TFs can describe the mechanistic basis for in vivo patterns of TF occupancy.

Results

Model evaluation and baseline results using the ChIP'ed transcription factor's motif

We began our analysis with fifty-five TF-ChIP data sets obtained from diverse sources (see Methods). Each TF-ChIP data set was represented by 1000 peaks and 1000 random non-coding sequence windows, all of length 500 bp. This representation was selected with the goal of identifying TF motifs that improve the ability to properly rank the occupancy within the peak group and/or improve the ability to discriminate between peaks and random sequences. The average ChIP score of each window was treated as the TF occupancy level in that window (Methods), and is henceforth called the “ChIP score”. For each data set a Position Weight Matrix (PWM or motif) representing the DNA binding specificity of the ChIP'ed TF was identified (Methods) and designated as the “primary motif”. We used the STAP program [22] (Figure 1A) to predict a binding level, henceforth called “STAP score”, for each sequence window in a data set, using the primary motif from that data set. We then computed the Pearson Correlation Coefficient (CC) between ChIP scores and STAP scores across the 2000 windows in each data set, and call this the “baseline CC” for the data set. This value captures the ability of the primary TF's binding motif to determine that TF's relative occupancy levels both within the most highly bound regions and in peak versus non-peak regions. Since STAP has one free parameter for which it requires training data (sequences and their binding levels), we performed 4-fold cross-validation to obtain STAP scores for all 2000 windows, with 500 test windows in each fold. Out of the 55 data sets, seven did not show a sufficiently high correlation (CC ≥0.15 and p-value<1E-11) here or in any other test that we report in the following sections, and three data sets presented technical problems in the training phase, e.g., inconsistent parameter values learned over different folds of cross-validation. These 10 data sets (Supplementary Table S1) are excluded from the rest of our report. In all of these examples, the TF motif is broadly confirmed by similarity to motifs for the same TF obtained by other methods or to motifs for homologous TFs. Thus, the low correlation may reflect a high degree of recruitment to DNA by other proteins, technical problems with this group of ChIP datasets, or with the model as applied to these datasets.

The results of this first exercise are shown in Table 1 and Figure 2A. We noted the baseline CC in this test to be ≥0.15 (p-value<1E-11) for 39 data sets, with the highest CC reported for the data set “TRL_Cchip_s5_14”, i.e., ChIP-chip data for the TF TRL in stage 5–14 embryo, obtained from the (C)avalli laboratory (see Table 1 legend for data set nomenclature scheme). We repeated this exercise, for all 55 data sets, using a second program, TRAP, also based on a biophysical model of DNA binding [23] with default parameter settings, and noted that CC values from STAP were generally better (Figure 2B), although there were several data sets where the two methods gave almost identical CC values. We also observed from Figure 2B that the ten data sets that we exclude from most of the analyses in this work (red symbols) received poor CC values from both STAP and TRAP. The purpose of this exercise was not to identify a superior method for occupancy prediction; such an attempt would have been biased since we have more experience with STAP than TRAP, and our TRAP analysis was run without training of free parameters. Our goal was to provide evidence that STAP-based predictions provide a reasonable baseline for more advanced models that will be examined below.

Figure 2C provides a scatter-plot visualization of the STAP results on the data set “TRL_Cchip_s5_14”, which has the best baseline CC value (CC = 0.765). Figure 2D provides an alternative visualization of the same results, as an ROC curve showing how the false positive rate of calling a ChIP peak based on STAP scores varies as we vary the STAP score threshold. We see that 89.7% of the 1000 ChIP peaks can be detected using STAP scores while making 10.3% false positive predictions; the Area Under Curve (AUC) is 0.96. Next to TRL, the TF with the highest CC is biniou (BIN), with the data set “BIN_Fchip_s14” exhibiting a CC of 0.654 and an AUC of 0.895 (Supplementary Figure S1). We note that this data set had been previously observed, in [24], to have a high enrichment of the TF motif in ChIP peaks. The ROC for a data set with a more typical value of CC is shown in Supplementary Figure S2 (CC = 0.305, AUC = 0.679). Figure 2E provides a different visualization of the accuracy of STAP predictions, as genome browser tracks of ChIP and STAP scores for the TF BCD on a single gene locus.

The CC values reported above can arise from differences in STAP scores of peaks and non-peaks in a data set, and from correctly modeling the ChIP scores within peaks and/or non-peaks. To examine the contribution of these two types of agreement between data and model, we separately calculated CC values among the peaks and non-peaks (Supplementary Table S10, Supplementary Figures S5A,B). We found several data sets where a significant overall CC was accompanied by a lower but significant CC within peaks, e.g., BCD_Bseq_s5, where the overall CC is 0.560, and the CC within peaks is 0.466 (Supplementary Figure S5A,C). These are examples where the goodness-of-fit arises from discrimination of peaks and non-peaks as well as from quantitative modeling of binding levels. In a few data sets, the signal appears to arise mainly from separation of peaks and non-peaks, e.g., BIN_Fchip_s14, where the overall CC is 0.654 but the CC within peaks and non-peaks is 0.233 and 0.185 respectively. By and large, the CC values within peaks were higher than those within non-peaks, as expected (Supplementary Figure S5B). Interestingly, for a few data sets the correlation within non-peaks was much greater than within peaks. These include UBX_Mchip_s5_14, EN_Mchip_s5_14, DISCO_Mseq_s5_11, EVE_Mseq_s14, and MAD_Bchip_s5, with peaks of the latter two exhibiting significant negative correlation between STAP predictions and ChIP scores (Supplementary Figure S5D). For these TFs, the motif score may be uniformly high in the top peaks, but help discriminate between very low and intermediate occupancy levels in the randomly selected regions, leading to stronger correlation within non-peaks.

We also repeated the evaluation of the above “baseline” model with a modified definition of data sets: now, the 1000 non-peaks of each data set was replaced with 1000 non-peaks randomly chosen from the ChIP peaks of other TFs. CC values analogous to those of Table 1 (Column “CC(M1)”) were computed and compared to those from Table 1 (Supplementary Table S15 and Figure S9). We observed that for a few data sets the new CC value is lower, suggesting that the primary TF motif in these cases may represent common features of TF bound regions. (We visit these cases in a later section.) We also noted several cases where the CC values were significantly higher when using other TFs' peaks as the non-peaks of a data set (e.g., Supplementary Figure S10). We believe such examples better reveal the role of the primary TF motif in determining the TF-DNA binding strength within accessible regions, since all segments considered in the newly defined data sets are ChIP peaks of some TF. Overall, our analysis of primary TF motif scores in different sets of genomic regions supports the idea that the 2000 regions selected for further study can provide insight into diverse types of mechanisms contributing to in vivo TF binding.

ChIP data supports concentration-dependence of TF-DNA binding

STAP uses a simple thermodynamic model to define TF-DNA occupancy for a genomic region based on binding site affinities, the equilibrium constant of the TF for its optimal site, and the TF's concentration level [22]. While the binding site affinity relative to that of the optimal site can be quantified using the PWM [25], the latter two quantities are formally unknown. The formula used by STAP (see Methods) features these two quantities as a mutual product, which is treated as a free parameter in the model. This TF-specific free parameter, henceforth denoted by γ, may lead to less of a difference in the contributions of high and moderate affinity sites. That is, at higher γ values, as would result from a high effective TF concentration, both high and moderate affinity sites may be fully occupied (saturated occupancy) whereas a stronger bias for high affinity sites will be observed at lower values (Figure 1A). As noted above, we use a cross-validation setting where the parameter is trained on three-quarters of the data and used to predict STAP scores in the left-out quarter, and the process is repeated four times. We examined the role of this parameter in the accuracy of the STAP model by varying it in the broad range 10⁻¹ to 10⁵ and recording the CC at each value of γ. As shown in Figures 3A–C, the optimal parameter value varies across data sets, between 10⁰ to 10⁴, with a roughly equal split into low, medium and high regions of the allowed range. All experiments reported in the rest of this paper were constrained to use γ in the range 10⁰–10⁴. We note that a value of γ = 10⁰ indicates that the optimal site has a fractional occupancy of 0.5 at cellular levels of TF concentration, while a value of γ = 10⁴ indicates a fractional occupancy of ∼1.

Figures 3A–C also reveal that for any given data set there is a substantial variation in the accuracy of STAP scores as we vary the TF-specific γ parameter. For instance, the CC value in the data set “TWI_Fchip_s9” (TWI at stage 9, source: Furlong lab ChIP-chip data) is about 0.25–0.30 at the two extreme values of γ (10⁻¹ and 10⁵ respectively), but reaches a much higher value of 0.42 at γ = 10². This dependence on the γ parameter, along with the variability of optimal γ across data sets underscores the importance of this parameter in the model. The parameter is analogous to the motif transition probability parameter in HMM-based models used in motif scanning, and our observation highlights the need for data set-specific training of this parameter in order to achieve the most accurate predictions. More generally, we conclude that simply adding the strengths of motif matches in a window is not necessarily the best way to predict TF occupancy in that window.

For a given TF, the γ parameter is proportional to the TF's concentration level in the experimental conditions. Therefore, if we have ChIP data on the same TF from two different stages, the optimal γ values ought to reflect the relative concentration levels in those stages. The examined collection of data sets included eight such pairs of data sets comprising ChIP data for the same TF from two different developmental stages. We therefore plotted the ratio of the trained γ values in the two stages versus the ratio of the TF's expression levels in those stages. We noted (Figure 3D) that the ratios of γ values were roughly consistent with ratios of expression levels, in that if one ratio is >1, the other ratio is also greater than or close to 1, and not <1. Expression levels were obtained from RNA-SEQ data from whole-embryos and may therefore be only a crude approximation of cell type-specific protein concentrations. This, and the fact that all ChIP experiments were performed on whole-embryo extracts, are expected to affect the sensitivity of this analysis, and may be the reason why we did not see a more quantitative agreement between the two ratios (i.e., points always close to the diagonal).

VFL and TRL binding sites frequently influence occupancy through long-range interactions

Our tests so far examined how different aspects of the primary TF, such as its binding specificity and concentration, affect its DNA-binding profile. In the next set of tests, we sought to evaluate the role of TFs other than the primary TF in determining the latter's occupancy. To this end, we used STAP with two motifs – the primary motif and one secondary motif at a time – and allowed cooperative interaction between TF molecules bound at sites within a certain distance, called the “distance threshold”, of each other (Figure 1B). There are now three free parameters: the two γ parameters corresponding to the primary and secondary motif, and a parameter representing the interaction energy between bound molecules of the primary and secondary TF. Evaluations performed under a cross-validation scheme ensured that CC values here are comparable to those in the baseline results from Table 1.

In the first set of tests of cooperative effects, we set the distance threshold to be 150 bp, therefore allowing long-range interaction that is similar to the length of DNA in one nucleosome. (We use “long-range” here to contrast with “short-range” interactions inferred from site pairs with ≤30 bp spacing in the next subsection, but note that “long-range” has different connotations in other contexts, e.g., to refer to interactions beyond enhancer boundaries [26], [27].) For each data set, we tested a secondary motif for every TF among the most highly expressed genes in the appropriate developmental stage, based on RNA-SEQ data [21]. We compared the CC of a (primary motif, secondary motif) pair to that from the primary motif (Table 1), and examined all cases where the improvement in CC was ≥0.04 (see Methods). The improvement, henceforth called ΔCC, was subjected to two different assessments of statistical significance. First, we recomputed the ΔCC with one hundred random variants of the secondary motif (see Methods), and asked what fraction of these random ΔCC values were better than the original ΔCC, thus obtaining a “ΔCC p-value”. Second, we utilized the ΔCC values from every tested secondary motif to compute a Z-score (see Methods). This mimics standard outlier detection procedures and designates a ΔCC value as significant if it appears to be an outlier compared to other observed ΔCC values for this data set. This is analogous to a multiple hypothesis correction, since we test over 50 candidate secondary motifs per data set. Additionally, we required that the cooperative interaction model has a greater CC than a model where the secondary motif alone is used by STAP. Thereby, we identified data sets where the combination of the primary and secondary motifs, through cooperative interactions, can describe the primary TF's occupancy better than either motif in isolation.

This analysis revealed 25 cases of significant improvements (ΔCC ≥0.04, p-value≤0.05 and Z-score ≥3), spread over 18 data sets (Supplementary Table S2). Table 2 tabulates the secondary motif with the most pronounced effect for each of these data sets. We noted that these effects arise mainly from an improved ability to discriminate peaks from non-peaks, and in only 4 (respectively 2) of these 18 cases the cooperativity model improves the CC even within peaks (respectively non-peaks) (Supplementary Table S11, Supplementary Figure S6). Remarkably, for 15 of these 18 data sets, the most influential secondary motif was either VFL (8 cases) or TRL (7 cases). Figure 4A shows an example where the use of VFL as a secondary motif significantly improves the ability to discriminate ChIP peaks from non-peaks. Overall, the VFL motif significantly improves primary TF occupancy predictions for 10 data sets (Supplementary Table S3, Figure S4A), of which 9 were from an early developmental stage (stage 5), and the tenth was from a broader span of developmental stages including stage 5. We noted that VFL is highly expressed in later stages as well and its motif was tested as a secondary motif in the corresponding data sets, but significant influences were not detected in those data sets. VFL has been proposed to play a “pioneer factor” role [28] in early development [29], [30], and its motif has been found to be highly over-represented in so-called “HOT” regions that represent the most accessible regions of the genome [31], [32]. Yanez-Cuna et al. [5] recently showed the VFL motif to be required for DNA-binding by the TF TWI, as well as for regulatory activity of TWI-bound enhancers, and to be enriched in early binding sites of other TFs such as MEF2. Our findings support these strong lines of evidence for an important facilitative role of VFL in determining TF-binding, and explicitly quantify this role for 10 different TF-ChIP data sets.

We found the TRL motif to influence the binding levels of primary TFs in eight data sets overall, of which six are from the later developmental stages 9–14 (Supplementary Table S4). As discussed in [33]–[35], TRL plays an important role in regulating the chromatin structure and packaging large segments of the chromosome into active (euchromatic) or inactive (heterochromatic) domains. The TRL motif was also prominent among sequence signatures of context-specific TF-DNA binding reported in [5], although this previous study did not explicitly quantify its facilitative effect on various primary TFs. It is interesting to note that TRL is the TF with the highest baseline CC (Table 1), reflecting the possibility that TRL-DNA binding is largely dependent on the TF itself and does not require facilitative effects of secondary TFs. This is consistent with speculation that TRL is a “pioneer factor” [28], [36].

In these initial tests, STAP was configured to allow interaction between primary and secondary motif as long as their bound sites were within 150 bps. We next asked if the promiscuous effects of VFL and TRL could be observed when reducing this distance threshold to 30 bps, which would suggest that short-range mechanisms of interaction might be involved. We found that in most cases the effects of these two motifs were not significant at the shorter distance range (Supplementary Tables S17, S18), and in the four cases where significant effects were detected at this range, the magnitude of the effect was lower than that at 150 bp range. A possible interpretation of this finding, especially in light of available knowledge about these two proteins, is that they act as chromatin remodelers over relatively long scales (150 bp or greater) and facilitate TF binding by making binding sites of the primary TF more accessible. (We revisit this point in a later section, by directly examining accessibility data.) Notably, the data sets for VFL and TRL themselves did not reveal any secondary motifs with significant effects, once again supporting a possible pioneer factor role for these two TFs.

Short-range interactions with secondary TF sites influence primary TF occupancy and can predict physical interactions between TFs

While VFL and TRL clearly show the most frequent effects on TF binding, a number of other influential secondary TF motifs were also revealed by our analysis; these are shown in Table 3. For each of these cases we report the ΔCC values at both distance thresholds (30 bp and 150 bp). Of particular interest were the (primary motif, secondary motif) pairs where the ΔCC was significant only at the 30 bp threshold, since this may reflect direct interactions. (These significant short-range interactions were reflected in a better ability to discriminate peaks from non-peaks rather than an improved ranking of the peaks; see Supplementary Table S12 and Supplementary Figure S7.) A case in point is the data set HB_Bchip_s9, for the TF hunchback (HB), where the secondary motif Adh transcription factor 1 (ADF1) improves the baseline CC of 0.204 to 0.303 when modeling heterotypic cooperativity at distance threshold 30 bp. The ΔCC of 0.099 is highly significant (empirical P-value = 0, i.e., no shuffled motif yielded better ΔCC), while that at the 150 bp threshold does not meet our significance criteria. A similar effect was observed for the ADF1 motif on HB ChIP data in stage 5 embryos. We hypothesized that this is evidence for direct physical interaction between HB and ADF1 resulting in modulation of HB binding levels. We searched for sequence signatures of such a hypothesized interaction in the relative spacing of HB and ADF1 binding sites. Examination of the 250 highest ChIP peaks in the data set showed a statistically significant bias (P-value 3E-4, see Methods) for spacing in the range 18–23 bps (Figure 4B). A similar test on 250 non-peaks from the data set showed no bias for this range or any other. This analysis suggests that proximally located binding sites of HB and ADF1 result in increased HB occupancy in ChIP peaks. We examined other data sets where the ΔCC was significant, and found similar evidence of biased inter-site spacing in ChIP peaks (Figure 4B), supporting the hypothesis that direct cooperative interactions may be a key factor in determining TF binding profiles in these cases. In some cases, e.g., the pair (D, TTK) and (GT, TTK), we noticed more than one preferred spacing range, separated by 11 bp, as might be expected due to proper phasing requirements between physically interacting TFs [37]. We also tested for biased inter-site spacing between the TFs Distal-less (DLL) and Zif Zinc-finger protein (ZIF, also called CG10267) (Figure 4B), because the ΔCC was found to be significant for this pair (empirical p-value 0), although the z-score of this ΔCC was 2.435, slightly below our chosen threshold of 3.0.

For each of the predicted heterotypic interactions shown in Figure 4B, we assayed for direct physical interactions between the TFs using a modification of the LUMIER method [38], [39]. In these experiments, one partner is expressed as a fusion to Maltose Binding Protein (MBP) and the other partner as a fusion to luciferase (luc). To avoid possible bridging interaction by other eukaryotic proteins, the proteins were expressed using a purified prokaryotic in vitro expression system and then combined for analysis. MBP-tagged proteins were isolated using amylose beads and the luciferase activity retained on the beads (via primary TF-secondary TF interaction), relative to a negative control with unfused luc, was used to calculate a luminescence intensity ratio (LIR, see Methods). A value of seven or greater was selected as a cutoff for positive interactions. This threshold is based on a set of positive and negative control interactions among bHLH protein dimers examined using this assay (HNP and MHB, unpublished) as well as additional negative controls using luc fused to the TF CLK or MBP without a fusion partner (Supplementary Figure S12). This threshold is more than twice as stringent than those used in previous studies examining protein interactions in cell culture [38], [39] and consequently may exclude some weaker interactions, including some that may only be significant in the context of cooperative binding to DNA. Each predicted interaction pair was examined in both configurations (e.g., the primary TF was fused to MBP in one experiment and to luc in the second). In addition, since Mothers against dpp (MAD) and Medea (MED) are known to bind DNA as a heteromeric complex [40]–[43], it is possible that any interaction computationally identified for one of these proteins is the result of an interaction with the other one. In our in vitro experiments, only direct physical interactions between two proteins are tested. Therefore, each of the predicted interactions for either MAD or MED was also tested with the other.

For five of the eight tested pairs (i.e., those not involving MAD or MED), a clear in vitro interaction was observed in both configurations (Figure 4C, Supplementary Table S9). For the two predicted interactions involving MAD, one of the two configurations gave a signal while the other was just below our selected cutoff. In one additional case, no physical interaction was observed between ribbon (RIB) and MED, but RIB was observed to interact with the MED binding partner, MAD. None of our tested negative controls was near the threshold and the interaction signal for most of the tested pairs was similar to two, well-established positive control interactions for this set of proteins, a MAD-MED heterodimer and a homodimer of giant (GT), which is a member of the bZIP family of TFs that bind DNA as homodimers [44]. Thus, all of the tested predictions are supported by a moderate to strong in vitro interaction, demonstrating that at least some of the short range cooperative interactions identified by our computational model reflect actual physical interactions that were previously unrecognized in large scale protein-protein interaction screens.

The physical interaction of TFs suggests that they may use cooperativity to increase binding of the primary TF to DNA sites with properly spaced binding sites for both TFs [45]. We tested this prediction for three of the above TF pairs using a variation of a previously described microwell assay [46] (Figure 5). The primary TF fused to luciferase and a secondary TF fused to MBP are used in an in vitro pull down assay with biotinylated dsDNA oligos containing a sequence from the ChIP peaks that contains binding sites for both TFs. The TFs are mixed with the biotinylated target site and an excess of unlabeled wild type or mutant competitor DNA. The competitor sequences used to examine cooperative DNA binding of ZIF and DLL are shown in Figure 5A and all sequences are shown in Supplemental Table S16. Streptavidin-mediated recovery of luc-TF/biotin-DNA complexes in the presence of excess wild type competitor (wt) indicates the background signal. In experiments with both TFs present (Figure 5B), the recovery of the luciferase-tagged primary TF in the presence of a competitor with mutations in both TF binding sites (e.g. ΔZIFΔDLL) increased 8–18 fold over the background in the presence of wt competitor (Figure 5B, upper panels). In contrast, little increase was observed when this experiment was repeated without the secondary TF (Figure 5B, lower panels), indicating that the secondary TF facilitates binding of the primary TF to these sites. The specificity of this interaction was confirmed by testing mutant competitor DNAs that disrupt the individual TF binding sites (e.g., ΔZIF or ΔDLL) or that increase the intersite distance by five base pairs (e.g., “+5”). Each of these alterations in the DNA sequence results in reduced competition by the mutant DNA competitor relative to wild type and increased recovery of the primary TF (Figure 5B). Furthermore, reduced competition is observed even when adding two competitors with mutations in one or the other individual TF binding site and each present at the same concentration as the wild type control; thus, high affinity binding requires the two TF binding sites to be present on the same DNA molecule with the proper spacing. These results indicate that the physical interactions detected for each of these pairs mediate cooperative DNA binding to an endogenous sequence from one of the top ChIP peaks.

In light of the possibility that the influence of short range cooperative interactions may be more pronounced when the interacting TFs are at relatively modest concentration levels, we extended the tests reported in Table 3 to include all candidate secondary TFs with expression in the top 50%. The results, shown in Supplementary Table S14, reveal that for several data sets stronger influences are detectable when allowing lower expression levels of the secondary TF. On the one hand, this means that the list of interactions identified in Table 3 is likely incomplete. On the other hand, the list shown in Table S14 must be interpreted with caution since testing more candidate secondary TFs may lead to spurious interactions being reported due to similarity of motifs between two candidates.

Competitive binding by secondary TFs frequently influences occupancy

Cooperative interactions are not the only manner in which one TF's binding may influence another's. Two TFs competing for overlapping binding sites can modulate each other's binding levels at the location [47]. Our next set of tests searched for evidence of this phenomenon in ChIP data sets. We used a two motif STAP model with no interaction terms, and compared the cross-validation CC from this model to the baseline CC of Table 1. The only way in which a secondary TF site can influence the binding prediction for the primary TF in the two-motif model is if their sites overlap (Figure 1B). The results, shown in Table 4, comprise 17 cases of significant ΔCC over the baseline model (ΔCC ≥0.04, P-value≤0.05, Z-score ≥3). In at least 10 of these cases, the secondary motif's presence is strongly anti-correlated with the primary TF's ChIP score, i.e., the competing motif is more frequent in non-peaks or in lower ranking peaks than in strong peaks. This may imply that the strong peaks exhibit selection against sites of the secondary TF competing with the primary motif. Figure 6A shows three examples of the pattern of overlap between sites of a primary TF and a secondary TF, observed in sequences with high STAP scores and low ChIP scores. We noted that in all of these cases, the overlapping sites tended to be suboptimal matches to either motif.

Two different data sets involving HB, one from stage 5 and the other from stage 9, were influenced by overlapping sites of the RETN motif (Table 4). RETN is a well-known repressor that acts through competitive binding when inhibiting activation by the TF engrailed (EN) [48]. Two other secondary motifs that seem to influence multiple data sets are EXD and HTH. Both of these homeodomain proteins play prominent roles during development as cofactors in repressor complexes with both Hox and other homeodomain proteins. Interestingly, in all three cases where EXD influences binding, there is no correlation between EXD sites and the primary TF occupancy, while in all three cases where HTH exerts an influence, there is a strong negative correlation (∼−0.18) between HTH motif presence and primary TF binding (see Discussion).

Antagonistic influence of a secondary motif can manifest even with non-overlapping sites

The next set of tests was directed at detecting evidence of antagonistic binding at non-overlapping sites. A possible mechanism for such a phenomenon is that of the secondary TF upon binding rendering the local DNA inaccessible, e.g., through recruitment of HDACs [49], as is speculated to be the case with some short-range repressors in Drosophila [50]. We used a two-motif STAP model with a TF-TF interaction term that is fit on training data, and compared the resulting CC to that from the primary motif alone (Table 1). This interaction term was constrained to be <1, corresponding to an unfavorable energy of interaction in the underlying thermodynamics model (Figure 1B). Note that this model incorporates both competitive binding and antagonistic influence from non-overlapping sites. Comparing the CC achieved by this model at either the 30 bp or the 150 bp distance threshold to the baseline (Table 5, Supplementary Figure S4B), we found 35 cases of significant improvements (ΔCC ≥0.04, P-value<0.05, Z-score ≥3). These included 6 data sets influenced by the EXD motif, 4 data sets by the HTH and RETN motifs, and 3 data sets by the JIGR1 motif. We noted that these four motifs were also observed to influence binding through competitive binding to overlapping sites (Table 4) above. However, in such cases where a secondary motif had significant effect on binding levels both in the competitive binding mode as well as the antagonistic binding mode, the magnitude of the effect was always stronger in the latter mode. The strongest case of antagonistic influence at the 30 bp distance threshold was estimated for the data set CAD_Bseq_s5, for the TF caudal (CAD), where the RETN motif improves the CC from 0.178 to 0.401. On the other hand, the strongest influence at the 150 bp threshold was by the EXD motif, also on the CAD_Bseq_s5 data set, where the baseline CC of 0.178 improved to 0.412, and this effect was exclusive to the 150 bp range. In fact, a large majority of the antagonistic binding influences were significant exclusively at either the short (30 bp) or the long (150 bp) range (Table 5). This may suggest that the underlying mechanisms of short and long-range antagonistic influences are different, although we did not observe any motif-specific preferences for one range versus the other.

We searched for inter-site spacing biases that might provide additional insights into the significant antagonistic influences identified above. It was commonly the case that ChIP peaks had a significant bias towards specific spacing values while non-peaks tended to avoid that range (Figure 6B, e.g., D-EXD). Interestingly, though less commonly, such spacing biases were also observed in non-peaks (Figure 6B, Supplementary Tables S5, S19). Even when examining antagonistic influences of the same secondary TF, e.g., HTH, we found some data sets where the spacing bias was exclusive to ChIP peaks and others where the bias was present in non-peaks.

Separate examination of peaks and non-peaks for effects of antagonistic influence revealed that such effects are manifested in a better discrimination of peaks versus non-peaks as well as a better modeling of ChIP scores with peaks alone or, more commonly, within non-peaks (Supplementary Table S13 and Supplementary Figure S8).

DNA accessibility data provides clues about mechanisms of secondary TF action

Recent work [6], [7] has shown that DNA accessibility data, which reflects nucleosome positioning and other chromatin-related effects, has a very strong correlation with TF occupancy, and when used in conjunction with the primary TF's motif can lead to highly accurate predictions of occupancy. This has been demonstrated in the context of five TFs in Drosophila (data from whole embryo) and six TFs in human (data from two cell lines). These prior results motivated us to examine the same hypothesis for the much larger collection of TF-ChIP data sets studied here. In all of our tests in this section we used DNaseI hypersensitivity data from [19]. In the first tests, we used a high threshold (90^th percentile) on developmental stage-specific accessibility to designate “accessible regions”, predicted zero occupancy in inaccessible regions, and used STAP and the primary motif to predict occupancy in accessible regions. Accessibility-filtered STAP scores computed in this manner correlated very highly with ChIP data (Supplementary Table S6), and led to substantial improvements upon the baseline results of Table 1, for 38 of the 45 data sets. This confirms that the observations made by Kaplan et al. and Pique-Regi et al. are manifest over a larger dataset.

The test above showed that motif and accessibility information together provide highly accurate predictions of ChIP scores. A natural question that arises then is: how strong is the influence of the primary TF's motif in determining its occupancy, beyond the influence of accessibility? To answer this question we computed the “semi-partial correlation coefficient” (SPCC) between ChIP and STAP scores, which subtracts or “partials out” the contribution of accessibility information. Technically, this amounts to first predicting ChIP scores using accessibility alone, and then correlating the residual ChIP scores with STAP scores (see Methods). We found that for the majority of data sets the SPCC values (Table 1, column SPCC(m1)) were comparable to the baseline CC values, demonstrating that, as expected, the primary motif plays a major role in shaping TF binding profiles. For ten data sets, SPCC was better than baseline CC, most notably for the data set TIN_Fchip_s9 where the primary motif's correlation improves from 0.428 to 0.507 upon partialing out accessibility. In these cases, factoring out the accessibility effects better reveals the expected relationship between primary motif presence in the sequence and occupancy. In contrast, five data sets showed a dramatically lower SPCC than CC (Table 1); these were related to the TFs VFL, TRL and MED. This is consistent with hypothesis emerging in this work (also see next paragraph) and in recent literature that VFL and TRL have direct influence on accessibility patterns, and partialing out the correlation with accessibility results in much reduced correlation between primary motif and TF occupancy. The third of the trio of TFs identified here, MED, is also believed to direct the co-factor CBP to the genome [51] and thus influence accessibility profiles. The SPCC was lower than CC also for TWI, D and SLP1, though not as dramatically. Sandmann et al. [52] have previously found TWI to bind to a large number of mesodermal enhancers and speculated that its role may be to facilitate chromatin remodeling. D is a SOX domain protein and there has been suggestion that this family of TFs may function as chromatin remodelers [53]. Interestingly, independent evidence in support of the accessibility-mediated effect of VFL, TRL, TWI, SLP1 and MED emerged when we repeated the evaluation of the single motif STAP model (baseline, Table 1) on data sets composed of the top 1000 ChIP peaks and 1000 random non-peaks selected from ChIP peaks of other TFs (Supplementary Table S15 and Supplementary Figure S9). We found the CC on these data sets to be conspicuously below that on the default data sets where the non-peaks were random genomic segments. This implies that the primary motif in these cases is better able to discriminate peaks of the primary TF from random non-peaks than from other accessible regions (peaks of different TFs). This in turn suggests that the motifs of VFL, TRL, TWI, SLP1 and MED may be common features of many ChIP peaks that discriminate them from random non-coding sequences irrespective of the bound TFs.

Our next tests examined the effect of cooperative binding with secondary TFs in the light of accessibility information. Recall that the VFL and TRL motifs had emerged as the most promiscuous influences in our tests above (Table 2), and that their influence was noted as being predominantly long-range (Supplementary Tables S3, S4), leading us to speculate that they may be mediated through modulation of local accessibility. We therefore asked if the improvements in CC due to either of these motifs are observed after removing the effects of accessibility information. We computed SPCC values of the cooperative interaction model after partialing out accessibility (Figure 7A), similar to that described in the previous paragraph. We found that the effects of TRL disappear in all 8 data sets where it had been significant before considering accessibility, adding evidence in favor of our hypothesis that TRL's influence is mediated by accessibility. In contrast, VFL was found to exhibit a more diverse behavior: in 7 of 10 data sets its effects vanished after considering accessibility, while in 2 data sets (CAD_Bseq_s5 and HB_Bseq_s5), a pronounced influence (ΔSPCC ≥0.04) remained even after partialing out accessibility (Supplementary Table S7). These two data sets also showed evidence of an inter-site spacing bias between VFL and the primary motif (Supplementary Figure S3). These findings suggest that VFL's influence on TF binding may involve distinct mechanisms, including not only a general effect on local accessibility, but also more TF-specific mechanisms potentially involving direct interactions with the primary TF.

We repeated the above analysis on data sets where secondary motifs other than VFL and TRL had led to significant improvements in CC through a cooperative binding model (Table 3). The results, shown in Table 6 and Figure 7B, reveal that in most cases the influence of the secondary motif is pronounced even after partialing out accessibility information. This suggests that most of these secondary TFs operate through primary TF-specific interactions rather than by only influencing accessibility. Similar results were obtained when examining the cases of antagonistic influence by secondary motifs (Figure 7C).

Discussion

We studied mechanistic determinants of TF-DNA binding by computationally modeling genomic occupancy from over 40 ChIP data sets obtained from four different stages of embryonic development, in conjunction with over 300 TF motifs and stage-specific DNA accessibility and RNA-SEQ data. Our ultimate goal is to use the insights revealed here, both general and data set-specific, to develop improved computational tools that can quantify functional TF-DNA interactions genome-wide. Such tools can potentially inform models of TF regulatory networks in the same way that ChIP data is beginning to be used today [1], [4]. We note that characterizing hundreds of TFs by the whole-genome ChIP-SEQ in the vast number of different cellular conditions is not currently feasible. Computational tools therefore offer an attractive alternative, especially if they can be shown to predict cell type-specific occupancy. TF motifs are already being characterized through high throughput technologies such as Bacterial 1-Hybrid [9], SELEX [11], [54], and Protein-Binding Microarrays [55]. Cell type-specific DNA accessibility profiles and TF expression levels only need to be characterized once for a given cell state, and can then be used to predict binding profiles for all TFs. Our work provides initial evidence for the feasibility of this vision. At the same time, we note that the CC values reported here should not be interpreted as correlation coefficients between genome-wide predictions and observed levels of TF binding. The manner in which we chose to evaluate various models, i.e., by examining agreement with ChIP scores on 1000 bound regions and 1000 randomly selected non-peaks, was dictated primarily by the goal of detecting significant influences on primary TF occupancy. We also note that the CC values varied substantially across data sets, from 0.765 for TRL to 0.062 for Dorsal (DL) (Table 1). This variation in model performance may reflect weaknesses of certain data sets or PWMs, or a variable reliance of ChIP scores on the primary TF's binding.

Despite a general appreciation of the potential role of various determinants of TF binding, there have been very few systematic studies of the extent of their influence across a large number of TFs. We review three such studies that set the stage for our own work and explain the main goals and contributions of our work in the backdrop of these important prior studies.

Kaplan et al. [6] studied ChIP-SEQ data on five TFs in early Drosophila development, and concluded that the TF motif and DNA accessibility are the most informative correlates of TF-DNA binding, as determined by the agreement between measured and predicted occupancy profiles. They also used TF sequence signatures to examine the role of competitive and cooperative interactions with other TFs with similar developmental roles and concluded that these interactions do not play a significant role overall. Their negative finding regarding secondary motifs may be limited to the small number of data sets examined, or be a limitation of the specific methodology adopted in the study (including the use of a more limited set of motifs that were available then). Here, we perform much more extensive tests of the role of the above-mentioned binding determinants of TF binding, by analyzing 45 TF-ChIP data sets spanning multiple stages of embryonic development in D. melanogaster. We primarily consider the influence of a large number of secondary TFs that are highly expressed in that developmental stage. In contrast to the earlier findings, we find many cases where the primary TF's binding levels are significantly influenced by the presence or absence of binding sites for other TFs.

In a related study, Pique-Regi et al. [7] considered the problem of classifying primary motif matches within ChIP peaks versus those outside of ChIP peaks, in the context of six ChIP-SEQ data sets from two human cell lines. They found accessibility and specific histone modifications to be the most useful features in this classification task, but did not consider the influence of secondary TFs. However, there are fundamental differences in the goals of our study from that of Pique-Regi et al. Their objective was to build a computational tool for annotating TF-bound sites genome-wide, and therefore their algorithm integrates several variables that correlate with binding, including evolutionary conservation, transcription start site proximity, DNA accessibility and histone marks. On the other hand, our focus is on the influence of variables that are expected to be mechanistic determinants of binding, and whose influence can be reasonably understood within an intuitive biophysical framework. We therefore focus specifically on testing whether and how binding sites of secondary TFs shape the primary TF's binding profile. In this pursuit, we rely upon motif, sequence and TF expression data, treating these as the “predictor variables” with which to model ChIP data. We do not include other variables such as evolutionary conservation (which is not a mechanistic determinant) or start site proximity (whose influence cannot be easily modeled biophysically) as predictors in this statistical exercise. DNA accessibility data is used in our analysis, not to improve occupancy prediction per se, but to answer a specific mechanistic question about how secondary TFs influence binding. Also, there is a fundamental technical difference between the data types modeled in the two studies: the variable we propose to model is not tied to TF-DNA interaction at an individual binding site as in [7], but to the aggregate effect of all binding events within a 500 bp window. For the simplicity, we ask whether a model can predict the actual ChIP score at a genomic position, rather than ask whether a model can predict whether a putative motif match falls within a significant ChIP peak or not.

A recent study by Yanez-Cuna et al. [5] searched for motif signatures of context specific binding of TFs. In particular, they analyzed ChIP data sets for the same TF from two different cellular conditions and asked if peaks exclusive to either condition could be discriminated on the basis of motif presence. They showed that such motif signatures do exist for the seven TFs examined and that general-purpose machine learning methods such as support vector machines can accurately classify context-specific binding sites using tens of motifs. In the same vein, they showed that bound and non-bound regions of a TF can be discriminated using a combination of tens of motifs, for most of the 21 TF-ChIP data sets examined. Additionally, they performed a closer examination of the binding determinants of one particular TF, twist (TWI), and demonstrated that binding sites for the secondary TFs VFL and TTK significantly affect the correct prediction of many context-specific TWI binding sites. While Yanez-Cuna et al. mostly focused on demonstrating that accessory motif signatures can distinguish TF-DNA binding regions in different developmental stages, our primary goal was to precisely identify the most influential secondary motifs for each of 45 different TF-ChIP data sets. To this end, we focused largely on quantifying the influence of secondary motifs and assessing their statistical significance rigorously. By performing our analysis over many data sets, we were able to gain more general insights about the widespread or TF-specific roles of particular secondary TFs. In particular, our statistical tests are geared towards explaining the mechanistic basis of such roles: short- versus long-range effects, synergistic versus antagonistic effects, chromatin mediated versus direct interactions, etc.

The review by Biggin [56] uses findings from recent studies to argue that accessibility is more important than the role of secondary TFs in determining primary TF binding levels. However, we do not attempt here to characterize the effect of accessibility as being stronger or weaker than the effect of interacting TFs. Integrating perspectives from Biggin and others [15]–[17], [57], [58], DNA accessibility in vivo can be considered the result of multiple factors playing out simultaneously, possibly including innate sequence preferences of nucleosome location, a conglomerate of chromatin remodeling activities and displacement of nucleosomes by competition with TF binding. Under this view, there are practical limitations in the approach of directly comparing the improvement in occupancy prediction due to accessibility information to that due to secondary motif information alone. Moreover, while it may be possible to make broad statements regarding the influence of accessibility or other chromatin-related information on TF binding, secondary TFs , due to the combinatorial nature of gene regulation, will be factor-specific in their effects and thus will only be detectable on a few data sets. Accordingly, our goal is to characterize as many of these determinants of TF occupancy, from each ChIP data set, rather than assign any one number to the overall influence of, say, interactions between the primary and secondary TFs, which will be factor dependent by definition.

A related study that examined the effects of secondary TFs on ChIP data is that of Gordan et al. [59] who reported on TF-ChIP data sets in yeast where a secondary motif alone was a better correlate of peak location than the primary motif. In some cases, this may be due to a problem with the primary motif (H.N.P. and M.H.B. unpublished results). In other cases, such a situation may reflect indirect binding of the primary TF to the peak, via physical interaction with the bound secondary TF. It suggests an alternative model of ChIP data, where binding is predicted to be a sum or linear combination of the occupancy values of the primary TF (direct binding) and a secondary TF (indirect binding). We have not explored this model here, and believe that it is an important goal for future studies.

Our approach to including accessibility data in the analysis was to use partial correlations to examine secondary TF influences before and after factoring out the effect of accessibility on ChIP scores. Alternative approaches may directly include accessibility data in the occupancy models, as was done by Kaplan et al [6], who changed prior probabilities of binding in their probabilistic model based on accessibility, and Pique-Regi et al. [7], who included DHS and histone modification data as features in their classifier. Future modifications of our approach will attempt to include accessibility within the biophysical framework of STAP, and may potentially reveal the role of accessibility even more accurately.

An intriguing observation from our analyses was the influence of competitive binding by the secondary TF EXD despite there being no correlation between EXD sites and the ChIP scores of the primary TF. It is puzzling because it suggests that the frequency of EXD sites does not differ between peaks and non-peaks, yet these sites somehow make a significant difference to binding predictions. However, it is possible that the frequency of EXD sites overlapping with primary TF sites is different between peaks and non-peaks, and the advanced model uses the competition for overlapping sites to predict lower occupancy in certain sequences than that predicted by the baseline model, leading to improved agreement with ChIP scores (Supplementary Figure S11).

Our work opens up several important directions of future research into TF-DNA interaction on a genomic scale. While the models we explored used at most one secondary motif in one interaction mode, a more realistic model will require integration of more than one underlying mechanisms influencing primary TF occupancy. Accessibility information will play a crucial role in the predictive ability of such models. In the longer term, an important goal will be to develop integrative models where sequence, TF gene expression and developmental history is sufficient to predict, at least to a good approximation, both accessibility patterns and TF-DNA binding profiles. With the future availability of large collections of TF motifs, such computational surrogates for cell type-specific ChIP data will enable global studies of gene regulatory networks and provide specific regulatory assignments that can be experimentally confirmed.

Materials and Methods

Data

We used 55 TF-ChIP data sets on 37 TFs active in early stages of Drosophila embryonic development. These include five ChIP-seq data sets and 20 ChIP-chip data sets from BDTNP [60], seven ChIP-chip data sets from the Furlong lab [24], and 21 normalized ChIP-chip and ChIP-seq data sets from the ModEncode project [1], [61]. ChIP data of VFL and TRL were obtained respectively from [29] and [62]. Stage-specific genome-wide DNaseI hypersensitivity (chromatin accessibility) data, which is mapped to genome release 4 coordinates, was downloaded from the first replicate in the BDTNP web site and converted to release 5 coordinates using the liftOver tool and chain files from the UCSC web site (http://hgdownload.cse.ucsc.edu/downloads.html). We used 614 Drosophila transcription factor motifs, corresponding to 322 distinct TFs, from the FlyFactorSurvey database [20]. The motifs were ranked based on expression of the associated TF gene, using RNA-SEQ data [21] for the appropriate developmental stage. In cases where a TF-ChIP data set corresponded to a range of stages, expression values were stage-normalized and averaged before ranking. Motifs corresponding to heterodimeric complexes (such as HLH TFs in complex with DNA) were not considered. Motifs in the top 10% of the expression-based ranked list for the appropriate developmental stage were tested as candidate secondary motifs. The one exception to this are results in Table 3 where the top 25% of the ranked list was considered.

We smoothed each TF-ChIP data set and each DNase I data set by assigning scores to each 500 bps window over the genome, with a 50 bps shift. First, raw “read scores” in a data set were mapped to the nearest genomic position that is a multiple of 50. The score of a 500 bp window was then computed by averaging over all read scores mapped to positions in that window; we refer to this as the “ChIP score” of the window. After this transformation, we selected 1000 non-overlapping, highest scoring windows as “peaks” and randomly extracted 1000 non-exonic, non-overlapping windows without replacement from the remaining genome as “non-peaks”. This set of 2000 windows and their ChIP scores constitutes a TF- and stage-specific data set in our analyses. A “primary” motif was designated for the data set, based on the availability of motifs for the ChIP'ed TF. In cases where there were multiple motifs available for the ChIP'ed TF, the motif with the highest correlation between STAP scores and ChIP scores over all 2000 windows (see below) was selected. “Secondary” motifs tested for potential effects on the primary TF's binding were selected based on expression data, as mentioned above.

Evaluations of STAP-based TF-DNA occupancy prediction by 4-fold cross-validation

We used the STAP program [22] to predict the ChIP score of a window, using the primary motif and optionally a secondary motif for that TF. STAP has one or more free parameters that require training data – a set of sequences and their ChIP scores. Hence, we used cross-validation to train and test various models of TF-DNA occupancy that are encoded by STAP. We randomly divided the 1000 peaks into 4 equal partitions and also the 1000 non-peaks into 4 equal partitions. In each fold of cross-validation, three partitions from the peaks and non-peaks were used as the training set and one partition (i.e., 250 peaks and 250 non-peaks) was the test set. Predicted ChIP-scores on each of the test sets of windows were collected together, and the resulting set of 2000 real and predicted ChIP score pairs were subjected to evaluations. Evaluations on a data set were considered a failure if the STAP parameter values learned in the four folds were widely different; this happened for one data set.

The STAP model

This was described in [22]. STAP considers each molecular configuration σ that specifies which sites in the given sequence are bound by their respective TFs. Following standard statistical physics, the “Boltzmann weight” of σ, denoted by W(σ), represents the relative probability of the system being in configuration σ, and is calculated based on TF concentration and the estimated binding affinity of every bound site in σ. The Boltzmann weight is a product of terms contributed by each TF-bound site in the configuration. This corresponds to the assumption that each bound TF interacts independently with the DNA, with energy contributions that add up [63]. See Figure 1A for an example where the sequence has two sites (‘A’ and ‘B’) for TF ‘A’, or Figure 1B where there is one site for each of two TFs ‘A’ and ‘C’. A site's contribution, q(S), depends on the TF concentration and the strength of site S, and is given by:

where [TF] is the concentration of the TF (in arbitrary units), LLR(⋅) is the log likelihood ratio score of a site, computed based on the known position weight matrix (PWM) of the TF [25], S_max is the strongest binding site of the TF, and K(S_max) is the equilibrium constant of the TF binding to this site. The product K(S_max)[TF] is a TF-specific free parameter denoted by γ_TF. Let N_k(σ) denote the number of bound sites of TF k in configuration σ. The STAP model predicts the occupancy of TF k as:

Note that while N_k(σ) counts the number of bound sites for TF k only, the Boltzmann weight W(σ) depends on bound sites for all TFs.

Assessment of statistical significance by P-value and Z-score

The accuracy of STAP predictions was assessed by computing the Pearson correlation coefficient (CC) between real and predicted ChIP scores of 2000 windows in a data set. To assess the impact of a secondary motif M₂ in modeling a data set whose primary motif is M₁, we tested STAP in a single motif mode (“STAP(M₁)”) and in two-motif mode (“STAP(M₁,M₂)”) and compared the difference in their accuracies: ΔCC = CC(STAP(M₁,M₂)) – CC(STAP(M₁)).

A secondary motif M₂ was deemed as a significant influence on the data set if the following conditions were met:

• ΔCC ≥0.04. This is used as a basic criterion before performing further assessments. It is a heuristic for reducing the number of tests we had to perform, thus speeding up the analysis.

• ΔCC′ = CC(STAP(M₁,M₂)) – CC(STAP(M₂)) ≥0.04. This is to ensure that the model when using both the primary and the secondary motif fits the data better than a model using the secondary motif alone (without the primary motif).

• P-value(ΔCC) ≤0.05. To compute this, we permuted the secondary motif by randomly shuffling rows and columns of the PWM, which keeps the motif information content intact. ΔCC was computed for each permuted version, and by repeating this 100 times, we estimated a p-value for the original ΔCC.

• ≥3. Here, is a “Z-score” that compares the given ΔCC to corresponding values obtained from every candidate secondary motif that was tested, and is analogous to a multiple hypothesis testing correction. It is a measure borrowed from outlier detection theory [64], and reflects if the observed ΔCC is an outlier compared to a given set of ΔCC values. It is defined as

We also evaluated the best secondary motif effect for each data set by computing an “Area Under ROC” (AUC) value for the interaction model (Supplementary Table S8).

Tests of inter-site spacing bias

For each significant case of cooperative or antagonistic influence by a secondary motif, we searched for biases in the inter-site spacing between the primary and secondary TFs. Let us assume a pair of motifs (M1, M2) represents the binding specificities of the primary and secondary TFs. To test for a specific spacing bias, say ‘d’ base pairs, between (M1, M2) in a given set of segments, we grouped all pairs of adjacent heterotypic binding sites (located by FIMO program with threshold of e⁻⁷ [65]) into those having or not having inter-site distance of d. We counted the number of site pairs in each group and compared these counts to the corresponding counts in a “background” data set using one-tailed Fisher's exact test. The “background” data set was constructed by shuffling the locations of predicted sites in each segment, thus preserving the number of binding sites in each segment, and pooling together 10 such randomized data sets (Kazemian et al., manuscript in review). Tests of spacing bias were conducted on a set of top 250 scoring ChIP peaks and separately on a set of bottom 250 non-peaks.

Semi-partial correlation

Semi-partial correlation is a statistical technique generally employed to assess the association of one random variable X with the other random variable Y after eliminating the effect of a third random variable Z on Y [66]. In our tests, X represents predicted TF-DNA binding, Y the experimental TF-DNA occupancy from ChIP, and Z the accessibility. The semi-partial correlation score between X and Y, after “partialing out” Z from Y, is computed as , where is the correlation coefficient between A and B.

In vitro analysis of TF-TF interactions

Protein interactions were measured in a modification of the previously described LUMIER or LuMPIS methods [38], [39] except that each protein was expressed in vitro rather than in cell culture. Open reading frame (ORF) clones for transcription factors were part of the Berkeley Drosophila Genome Project the collection of universal donor clones [67]. ORFs were transferred into two vectors, pHPT7-FlRluc-BD and pHPT7-MBP-BD (HNP and MHB, unpublished), using Cre Recombinase (New England Biolabs, M0298L). For one TF, Mad, the ORF was PCR amplified ligated into AscI and PmeI restriction sites in each vector. These vectors contain a T7 promoter for in vitro transcription, a loxP site for cloning and either maltose binding protein (MBP) or Renilla luciferase (luc) coding regions. Clone names and primer sequences are provided in the supplementary information (Table S9).

Proteins were made by coupled in vitro transcription/translation using the PURExpress In Vitro Protein Synthesis Kit (NEB, E6800S). All samples were analyzed by Western Blot to confirm that some full-length product was obtained. Luciferase input was measured using the Renilla Luciferase Assay System (Promega, E2820). The proteins were diluted with IP Buffer (150 mM NaCl, 50 mM Tris pH 7.4) such that roughly 10⁶ luciferase counts were added to each sample and an equivalent amount of MBP protein were mixed. Proteins were incubated with gentle rocking for 4°C for 2 hours. Amylose Resin (NEB, E8021S) blocked with 5 percent BSA was added to proteins and incubated with rocking at 4°C for 2 hours. The samples were washed twice with IP buffer and transferred to 96-well plates (Corning, 07-200-589) for luciferase measurements. The luminescence intensity ratio was measured using as follows:

Each experiment was performed in duplicate, the experiments were averaged, and the standard deviation was calculated. Source cDNAs, amplification primers and luciferase data are compiled in Supplementary Table S9.

In vitro analysis of TF binding

In vitro synthesis of tagged TFs and luciferase assays were performed as described above. Target sequences were identified from the top ChIP peak regions that contained strong matches to the primary and secondary motifs with a spacing and orientation that was most frequently observed. Other criteria used in selecting target sequences included whether the ChIP peak lies within a known enhancer, and whether its predicted occupancy under STAP's cooperativity model is higher (in rank) than that under the baseline model without cooperativity. Double stranded DNA oligonucleotides were synthesized that contained wild type or altered sequences. One oligonucleotide containing the wild type sequence is biotinylated on the first base. The genomic coordinates for the wild type sequences and all mutant sequences are shown in Supplementary Table S16. Protein-DNA interactions were measured in a modification of a previously described microwell-based assay [46]. Tagged TFs were expressed in vitro rather than in cell culture and diluted with low-stringency binding buffer (140 mM KCl, 5 mM NaCl, 1 mM K₂HPO₄, 2 mM MgSO₄, 20 mM HEPES (pH 7.05), 100 µM EDTA, 1 µM ZnSO4) +1% BSA. Oligonucleotides were annealed and diluted using annealing buffer (50 mM Tris-HCl, 0.2 mM MgSO₄, pH 7.0). Annealed oligo mixes were prepared with 5 ul of 1.2 uM biotinylated oligos, 5 ul of 24 µM competitor oligo, 2 ul of 500 ng/ul Poly(dI-dC)*Poly(dI-dC), and 8 ul of annealing buffer (final volume 20 ul) and incubated for 1 hour. 10⁶ luciferase counts of the luc-tagged primary TF and (if appropriate) an equivalent amount of MBP-tagged secondary TF were mixed (30 ul volume). The diluted proteins were added to the DNAs and incubated with gentle rocking at 4°C for 2 hours. Streptavidin coated 96 well plates (ThermoScientific # 15502) were blocked with 5% BSA and low stringency binding buffer. The protein/oligo mixture was added to the plates and incubated for 2 hours at 4°C. The samples were washed twice with low stringency binding buffer. Recovered luciferase activity was measured directly in the plates. All values were normalized by dividing by the luciferase counts recovered in the sample containing an excess of wild type competitor DNA.

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10

Attachment 11

Attachment 12

Attachment 13

Attachment 14

Attachment 15

Attachment 16

Attachment 17

Attachment 18

Attachment 19

Attachment 20

Attachment 21

Attachment 22

Attachment 23

Attachment 24

Attachment 25

Attachment 26

Attachment 27

Attachment 28

Attachment 29

Attachment 30

Attachment 31