Transcription factor binding site clusters identify target genes with similar tissue-wide expression and buffer against mutations

Similarity between GTEx tissue-wide expression profiles of genes

We used data from the Genotype-Tissue Expression (GTEx, version 6p) project which measured expression levels for each gene in 53 tissues, in different numbers of individuals (5 to 564). For each tissue population, the median expression value is given in RPKM (Reads Per Kilobase of transcript per Million mapped reads) for each gene²⁰. Data are available on Zenodo²¹. To capture the tissue-wide overall expression pattern of a gene instead of within a single tissue, the tissue-wide expression profile of a gene was defined as its median RPKM across GTEx tissues, which is described by a 53-element vector (Equation 1). Note that different isoforms whose expression patterns may significantly differ from each other cannot be distinguished by this approach.

where EP^A is the tissue-wide expression profile of Gene A, $ME{V}_{\text{1}}^A$ is the median expression value of Gene A in Tissue 1, $ME{V}_{\text{2}}^A$ is the median expression value of Gene A in Tissue 2, etc.

To discover other genes whose tissue-wide expression profiles are similar to a given gene, we computed the Bray-Curtis Similarity (Equation 2) between the tissue-wide expression profiles of all gene pairs. Relative to other similarity measures (Table 1, Additional file 1²²), this function exhibits desirable properties, including: 1) being bounded between 0 and 1, 2) achieving maximal similarity of 1, if and only if two vectors are identical, and 3) larger values having a larger impact on the resultant similarity value.

Prediction of genes with similar tissue-wide expression profiles

The framework for predicting whether the tissue-wide expression profile of a gene resembles a particular gene is indicated in Figure 1 (panels A and B). The genomic locations of all DHSs in 95 cell types from the ENCODE project[23; hg38 assembly] were filtered for known promoters²⁴; these sequences were then scanned for TFBSs using 94 iPWMs corresponding to the primary binding motifs for 82 TFs³. Data are available on Zenodo²¹. Heterotypic TFBS clusters were detected with the IDBC algorithm by specifying a minimum threshold of 0.1*R_sequence for the R_i values of individual TFBS elements in potential clusters; this eliminated weak binding sites detected with iPWMs corresponding to false positive, non-functional TFBSs³.

Figure 1. The general framework for predicting genes with similar tissue-wide expression profiles and TF targets.

Red and blue contents are respectively specific to prediction of genes with similar tissue-wide expression profiles and prediction of TF targets. (A) An overview of the ML framework. The steps enclosed in the dashed rectangle vary across prediction of genes with similar tissue-wide expression profiles and TF targets. The step with a dash-dotted border that intersects promoters with DHSs is a variant of the primary approach. In the IDBC algorithm (Additional file 1²²), the parameter I is the minimum threshold on the total information contents of TFBS clusters. In prediction of genes with similar tissue-wide expression profiles, the minimum value was 939, which was the sum of mean information contents (R_sequence values) of all 94 iPWMs; in prediction of direct targets, this value was the R_sequence value of the single iPWM used to detect TFBSs. The parameter d is the radius of initial clusters in base pairs, whose value, 25, was determined empirically. The seven ML features derived from TFBS clusters are described in the Methods section. The performance of seven different classifiers was evaluated with ROC curves and 10-fold cross validation (Additional file 1²²). (B) Obtaining the positives and negatives for identifying genes with similar tissue-wide expression profiles to a given gene (Additional file 2²²). (C) Obtaining the positives and negatives for predicting target genes of seven TFs using the CRISPR-generated perturbation data in K562 cells (Additional file 3²²). (D) Obtaining the positives and negatives for predicting target genes of 11 TFs using the siRNA-generated knockdown data in GM19238 cells (Additional file 4²²).

The seven information density-related ML features (Additional file 1²²) derived from each TFBS cluster included: 1) The distance between this cluster and the transcription start site (TSS), 2) The length of this cluster, 3) The information content of this cluster (i.e. the sum of R_i values of all TFBSs in this cluster), 4) The number of binding sites of each TF within this cluster, 5) The number of strong binding sites (R_i > R_sequence) of each TF within this cluster, 6) The sum of R_i values of binding sites of each TF within this cluster, and 7) The sum of R_i values of strong binding sites (R_i > R_sequence) of each TF within this cluster.

Each of the Features 1–3 was defined in a gene as a vector, whose size equalled the number of clusters in the gene promoter; each cluster was mapped to a single value in the vector. In Features 4–7, each cluster itself was mapped to a vector corresponding to binding sites for 82 TFs (Additional file 1²²). If two genes contained different numbers of clusters, the maximum number of clusters among all instances was determined, and null clusters were added to the 5’ end of promoters with the smaller numbers of clusters; this enabled all genes to exhibit the same cluster counts. This allowed all genes to be used as training data for ML classifiers. Default parameter values for these classifiers were used to generate ROC curves with a built-in MATLAB function (Additional file 1²²).

Prediction of TF targets

Perturbed target gene expression after CRISPR-based mutation of TF genes. Dixit et al. performed CRISPR-based perturbation experiments using multiple guide RNAs for each of ten TFs in K562 cells, resulting in a matrix of coefficients indicating the effect of each guide RNA on each of 22,046 genes¹⁶. Data are available on Zenodo²¹. The coefficient of a guide RNA is defined as the log₁₀ fold change in gene expression level¹⁶. We previously derived iPWMs for the primary binding motifs of 7 of the ten TFs (EGR1, ELF1, ELK1, ETS1, GABPA, IRF1, YY1)³. The framework for predicting TF targets (Figure 1A and 1C) was applied to these TFs. We defined a positive TF target gene in a cell line as :

If the coefficients for all guide RNAs of the TF for a gene were zero, the gene was defined as a negative (i.e. a non-target gene). As the threshold ε increases, the number of positives strictly decreases. As ε decreases, we have lower confidence that the positives were DE as a result of the TF perturbation. We balanced the number of positives obtained against our confidence that they were DE by evaluating different values of ε (i.e. 1.01, 1.05 and 1.1; Additional file 1²²). For each TF, all ENCODE ChIP-seq peak datasets from the K562 cell line were merged to determine positives. Data are available on Zenodo²¹. To avoid imbalanced datasets that significantly compromised the classifier performance²⁷, the Bray-Curtis function was applied to compute the similarity values in the tissue-wide expression profile between all negatives and the positive gene with the largest average coefficient, and then negatives with the smallest values were selected, resulting in the same number of positive and negative genes (Figure 1C).

Because TFs may exhibit tissue-specific sequence preferences due to different sets of target genes³, the K562 cell line-derived iPWMs of EGR1, ELK1, ELF1, GABPA, IRF1, YY1 were used to detect binding sites in DHS-intersected intervals; for ETS1, we used the only available iPWM from GM12878³. Six features (Features 1-5 and 7) were derived from each homotypic cluster (i.e. Feature 6 became identical to Feature 3, because only binding sites from a single TF were used) (Figure 1A). The results of 10 rounds of 10-fold cross validation were averaged to evaluate the accuracy of the classifier.

Target gene expression changes after siRNA-based knockdown of TFs. Significant changes in expression of target genes upon knockdown of 59 TFs in the GM19238 line were identified from the probability of differences in expression level (P-value) relative to the null hypothesis of no change¹⁴. Data are available on Zenodo²¹. DE genes with larger changes exhibited smaller P-values. The distribution of ChIP-seq peaks were considered to be evidence of TF binding to the promoters of these genes¹⁴. Among these 59 TFs, we have previously derived iPWMs exhibiting primary binding motifs for 11 (BATF, JUND, NFE2L1, PAX5, POU2F2, RELA, RXRA, SP1, TCF12, USF1, YY1)³. For this reason, transcription targets of these TFs were predicted in the GM19238 cell line (Figure 1A, D).

For ML training, we defined a positive (i.e. a target gene) for a TF, if the P-value of this gene was ≤ 0.01, and the promoter interval up to 10 kb upstream of a TSS overlapped one or more ChIP-seq peaks of the TF in GM12878. Other genes with P-values > 0.01 exhibited insufficient support for being TF targets and were labeled as negatives.

The iPWMs from GM19238 were used to detect binding sites for all TFs except for RELA, RXRA and NFE2L1. The iPWM from the GM19099 line was used for RELA, and for RXRA and NFE2L1, the only available iPWMs were derived from HepG2 and K562 cells, respectively. The DHSs in GM19238 were first remapped to the hg38 assembly using liftOver, prior to being intersected with known promoters²⁸. Data are available on Zenodo²¹. Although the knockdowns were performed in GM19238, GM12878 and GM19099 are also lymphoblastic cell lines, with GM19099 and GM19238 both being derived from the Yoruban population. For this analysis, the iPWMs derived in GM12878 and GM19099 were more appropriate sources of accessible TFBSs than those from HepG2 and K562, since GM12878 and GM19099 are more likely comparable to GM19238 than HepG2 and K562. ML results were evaluated by averaging cross validation, as described above.

Mutation analyses on promoters of TF targets

To understand the significance of individual binding sites on the regulatory state of TF targets, we evaluated the effects of sequence changes in TFBSs that altered the R_i values of these sites, the definition of TF clusters, and whether consequential IDBC-related changes impacted the prediction of TF target genes. For example, mutations were sequentially introduced into the strongest binding sites in TFBS clusters of the EGR1 target gene, MCM7, to determine the threshold for cluster formation and whether disappearing clusters disabled MCM7 expression. For one target gene of each TF from the CRISPR-generated perturbation data, effects of naturally occurring TFBS variants present in dbSNP²⁹ were also evaluated to explore aspects of TFBS organization that enabled both clusters and promoter activity to be resilient to binding site mutations. This was done by analyzing whether the occurrence of individual or multiple single nucleotide polymorphisms (SNPs) lead to the loss of binding sites and the corresponding clusters, and resulted in changes in the predictions for these targets.

The Bray-Curtis Function can accurately quantify the similarity between tissue-wide gene expression profiles

We computed the similarity values (Equation 2) between the tissue-wide expression profiles of the glucocorticoid receptor (GR or NR3C1) gene and all other 18,812 PC genes to find genes with related profiles. NR3C1 is an extensively characterized TF with many known direct target genes³⁰. As a constitutively expressed TF activated by glucocorticoid ligands, the protein mediates up-regulation of anti-inflammatory genes by binding as homodimers to glucocorticoid response elements. Transcription of proinflammatory genes is down-regulated by complexing with other activating TFs (e.g. NFKB and AP1), thereby eliminating their ability to bind and activate targets³⁰. NR3C1 can bind its own promoter forming an auto-regulatory loop, which also contains functional binding sites of 11 other TFs (e.g. SP1, YY1, IRF1, NFKB) whose iPWMs have been developed and/or mutual interactions have been described previously^3,30. The tissue-wide expression profile of NR3C1 comprises all different splicing and translational isoforms (GRα-A, GRα-B, GRα-C, GRα-D, GRβ, GRγ, GRδ). However, the profile averages out tissue-specific preferences of some isoforms, for example, GRα-C isoforms are significantly higher in the pancreas and colon, whereas levels of GRα-D are highest in spleen and lungs³⁰. We found that SLC25A32 and TANK have the greatest similarity in expression to NR3C1 (0.880 and 0.877 respectively), based on their overall similar expression patterns across the 53 tissues (Figure 2).

Figure 2. GTEx tissue-wide expression profiles of NR3C1, SLC25A32 and TANK.

Visualization of the expression values (in RPKM) of these genes across 53 tissues from GTEx. For each gene, the colored rectangle belonging to each tissue indicates the valid RPKM of all samples in the tissue, the black horizontal bar in the rectangle indicates the median RPKM, the hollow circles indicate the RPKM of the samples considered as outliers, and the grey vertical bar indicates the sampling error. A comparison of the panels shows that the overall expression patterns of the three genes across the 53 tissues resemble each other (e.g. all three genes exhibit the highest expression levels in lymphocytes and the lowest levels in brain tissues).

The Decision Tree classifier performed best in prediction of genes with similar tissue-wide expression profiles

Several ML classifiers (Naïve Bayes, Decision Tree (DT), Random Forest and Support Vector Machines (SVM) with four different kernels) were evaluated to determine how well TFBS-related features could predict genes with tissue-wide expression profiles similar to NR3C1. Their performance were compared using ROC curves, for complete promoters or for accessible promoter sequences that were first intersected with DHSs (Figure 3). DT exhibited the largest AUC (area under the curve) under both scenarios, and was one of two most stable classifiers (i.e. ΔAUC < 0.01), with the other being the SVM with RBF kernel. Consistent with previous findings^10,31,32, inclusion of DHS information significantly improved AUC values of the other classifiers with the exception of Naïve Bayes. In many instances, all TFBSs in a contiguous DHS interval formed a single binding site cluster.

Figure 3. Comparison between the performance of different classifiers in prediction of genes with similar tissue-wide expression profiles to NR3C1.

(A) ROC curves and AUC of seven classifiers without intersecting promoters with DHSs. (B) ROC curves and AUC of seven classifiers after intersecting promoters with DHSs. The Decision Tree classifier exhibited the largest AUC under both scenarios, and inclusion of DHS information significantly improved other classifiers’ AUC except for Naïve Bayes.

The Decision Tree classifier was predictive of TF target genes

Based on its performance in distinguishing genes with NR3C1-like tissue-wide expression profiles, the DT classifier was used to predict TF targets respectively based on the CRISPR-¹⁶ and siRNA-generated¹⁴ perturbation data. Performance was assessed with 10 rounds of 10-fold cross validation (Tables 2 & Table 3). Gini importance values were also used to assess the relative contribution of the six ML features to the predictive power of the classifier (Additional file 5²²). For the CRISPR-perturbed TFs in the K562 cell line, clustering Features 1-3, particularly the TFBS cluster information densities (Feature 3), were most important. The TFBS-level Feature 5 comprising the distribution of strong binding sites accounts for the largest contribution to classifier performance for the 11 siRNA-perturbed TFs in the GM19238 cell line. To assess the value of all ML features in capturing the distribution and composition of CRMs in promoters, all but one (TFBS counts) were sequentially removed, and the impact on accuracy of the classifier was determined. In each instance, classifier performance decreased, except for CRISPR-perturbed GABPA, IRF1 and YY1 upon inclusion of DHS information (Additional file 5²²).

The DT classifier predicted TF targets with greater sensitivity and specificity, after eliminating TFBSs in inaccessible promoter intervals in the CRISPR-generated knockdown data (Table 2). Specifically, predictions for EGR1, ELK1, ELF1, ETS1, GABPA, and IRF1 were more accurate than for YY1, which itself represses or activates a wide range of promoters by binding to sites overlapping the TSS (Table 2). Accordingly, perturbation results showed that YY1 has ~4-22 fold more PC targets in the K562 cell line than the other TFs (ε = 1.05). Binding of YY1 more significantly impacts the expression levels of target genes. The ratio of the PC targets at ε = 1.1 vs ε = 1.01 was 0.334, which significantly exceeded those of other TFs in this study (0.017-0.082; Additional file 3²²). This is consistent with the extensive interactions of this factor with many other TF cofactors in K562 cells, and its central role in specifying erythroid-specific lineage development³.

We found that promoters of most TF targets contain accessible, likely functional binding sites that are significantly correlated with changes in gene expression levels. Despite a high accuracy of target recognition, sensitivity did not exceed specificity, except for IRF1 (Table 2), due to a relatively large number of false negative genes. Promoters of non-targets are either devoid of accessible binding sites, or these sites are non-functional. In these instances, the classifier was unable to distinguish between likely functional binding sites in targets and non-functional sites in non-targets. In vivo co-regulation mediated by interacting cofactors, which was excluded by the classifier, assisted in distinguishing these non-functional sites that do not significantly affect gene expression¹⁴.

As the minimum differential expression threshold ε increased, the accuracy of the classifier for all the TFs monotonically increased, as expected (Figure 4). In general, more significantly DE genes have been associated with a higher number of TFBSs in their promoters¹⁴. Thus, at greater ε, there are larger differences in the values of ML features derived from TFBS clusters between direct targets and non-targets.

Figure 4. Accuracy of the Decision Tree classifier when using three different values for.

Each accuracy value was averaged from 10 rounds of 10-fold cross validation. The minimum threshold ε of the average fold change in gene expression levels (for all guide RNAs) of the TF was determined for fold changes: 1.01, 1.05 and 1.1. The accuracy value of each individual round is indicated in Additional file 5²². As ε increased, accuracy for all seven TFs monotonically increased.

Some TF target genes also display similar tissue-wide expression profiles to the TFs, themselves

To determine how many TF targets have similar tissue-wide expression profiles, we intersected the set of targets with the set of 500 PC genes with the most similar tissue-wide expression profiles for each TF (Table 4, Additional file 6²²). For example, the B lymphocyte expression profiles of TFs PAX5 and POU2F2 are similar to their respective targets, IL21R and CD86. The intersected targets for YY1 include 21 and 7 nuclear mitochondrial genes (e.g. MRPL9 and MRPS10, which are subunits of mitochondrial ribosomes), respectively, in the K562 and GM19238 cell lines³³. YY1 is known to upregulates a large number of mitochondrial genes by complexing with PGC-1α in C2C12 cells³⁴, and genes involved in the mitochondrial respiratory chain in K562 cells¹⁶. Our results are consistent with the idea that YY1 may broadly regulate mitochondrial function within all 53 tissues, in addition to the erythrocyte, lymphocyte and skeletal muscle cell lines.

Between 0.4%–25% of genes with similar tissue-wide expression profiles to the TFs are actually their targets (Table 4). The majority are non-targets whose promoters contain non-functional binding sites that lack co-regulation by corresponding cofactors. For YY1 and EGR1, we contrasted the flanking cofactor binding site distributions and strengths in the promoters of the most similarly expressed target genes (YY1: MRPL9, BAZ1B; EGR1: CANX, NPM1) with non-target genes (YY1: ADNP, RNF25; EGR1: GUCY2F, AWAT1). In these target gene promoters, strong and intermediate TFBSs recognized by SP1, KLF1, CEBPB formed heterotypic clusters with adjacent YY1 sites. Additionally, TFBSs of SP1, KLF1, and NFY tended to be present adjacent to EGR1 binding sites (Additional file 7²²). These patterns contrasted with enrichment of CTCF and ETV6 binding sites in gene promoters of YY1 and EGR1 non-targets (Additional file 7²²). Previous studies have reported that KLF1 is essential for terminal erythroid differentiation and maturation³⁵. Direct physical interactions between YY1 and the constitutive activator SP1 synergistically induce transcription³⁶, through activating CEBPB which promotes differentiation and suppresses proliferation of K562 cells by binding the promoter of the G-CSFR gene encoding a hematopoietin receptor³⁷. EGR1 and SP1 synergistically cooperate at adjacent non-overlapping sites on EGR1 promoter but compete for binding at overlapping sites³⁸, whereas occupied CTCF binding sites often function as an insulator blocking the effects of cis-acting elements and preventing gene activation by mediating long-range DNA loops to alter topological chromatin structure^39,40. ETV6, a member of the ETS family, is a transcriptional repressor required for bone marrow hematopoiesis and associated with leukemia development⁴¹.

Transcription factor binding site clusters buffer against expression changes from mutations in single sites

These results and previous studies indicate that the promoters of direct target genes contain multiple binding site clusters. We used our ML models of TFBS organization to investigate the effects of mutations in individual binding sites on the predicted expression of TF targets. We hypothesized that alternative TF clusters in the same promoter might stabilize and compensate for the loss of a mutated TF cluster, enabling mutated promoters to retain some capacity to regulate gene transcription upon TF binding. First, we introduced artificial variants into binding sites in silico in the promoter of the target gene MCM7 of EGR1 and examined the effect on the output of the ML classifier (Figure 5). In the K562 line, MCM7 is upregulated by EGR1. Knockdown of MCM7 has an anti-proliferative and pro-apoptotic effect on K562 cells⁴² and the loss of EGR1 increases leukemia initiating cells⁴³, which suggests that EGR1 may act as a tumor suppressor in K562 cells through the MCM7 pathway.

Figure 5. Mutation analyses on the target MCM7 of EGR1.

This figure depicts the effect of a mutation in each EGR1 binding site cluster of the MCM7 promoter on the expression level of MCM7, which is a target of the TF EGR1. The strongest binding site in each cluster were abolished by a single nucleotide variant. Upon loss of all three clusters, only weak binding sites remained and EGR1 was predicted to no longer be able to effectively regulate MCM7 expression. Multiple clusters in the promoters of TF targets confer robustness against mutations within individual binding sites that define these clusters.

The strongest binding site (at position chr7:100103347[hg38], - strand, R_i = 12.0 bits) in the promoter was eliminated by a G->A mutation, resulting in the loss of Cluster 1 (Figure 5), which consists of two sites (the other site at chr7:100103339, -, 4.3 bits). The other two clusters comprising weaker binding sites of intermediate strength (chr7:100102252, +, 7.0 bits; chr7:100102244, +, 1.3 bits; chr7:100101980, +, 8.9 bits; chr7:100101977, +, 2.2 bits; chr7:100101984, +, 1.9 bits) were still predicted to compensate for this mutation, so that the promoter maintains the capacity to induce MCM7 expression (Figure 5). Adjacent sites within the same TFBS cluster, which may individually not have sufficient affinity to strongly bind TFs and activate transcription, are capable of stabilizing binding to adjacent sites². Weaker sites can direct TF molecules to strong sites and extend the duration of physical association, termed the funnel effect². Binding stabilization between adjacent sites and the funnel effect enable CRMs comprised of information-dense clusters to be robust to mutations in individual binding sites^2,44. In this example, Clusters 2 and 3 were also respectively removed by G->T and C->G mutations abolishing the strongest site in either cluster, altering the prediction of the classifier, that is, EGR1 is expected to fail to induce MCM7 transcription (Figure 5). The remaining four sparse weak sites do not form a cluster and cannot completely supplant the disrupted strong sites.

We then examined the predicted impacts of natural SNPs on binding site strengths, clusters and predicted the regulatory state of the promoter for a direct target of each of the seven TFs from the CRISPR-generated perturbation data (Table 5). We found that a single SNP (e.g. rs996639427 of EGR1) could affect the strengths of multiple binding sites within a cluster (Table 5). Apart from SNPs that are predicted to abolish binding (Figure 5), leaky variants that are expected to merely weaken TF binding are also common (Table 5). Neither mutations that abolish TFBSs nor leaky SNPs in flanking weak binding sites would be expected to inactivate TFBS clusters (e.g. rs1030185383 and rs5874306 of IRF1), whereas SNPs with large reductions in R_i values of strong sites are more likely to abolish clusters (e.g. rs865922947, rs946037930, rs917218063 and rs928017336 of YY1) (Table 5). Multiple TF clusters enable promoters to be resilient to the effects of these mutations; only the complete inactivation of all clusters by concurrent effects of multiple SNPs within TFBSs would be capable of making a promoter to be unresponsive to TF binding (e.g. rs997328042, rs1020720126 and rs185306857 of GABPA) (Table 5). Conversely, a small number of SNPs capable of strengthening TF binding and reinforcing the regulatory effect of the TF were also observed (e.g. rs887888062 of EGR1 and rs751263172 of ELF1) (Table 5). In addition to deleterious mutations, potentially benign variants may also be found in promoters, consistent with the expectations of neutral theory⁴⁵.

Underlying data

The median RPKM, TSS coordinate, DNase I hypersensitivity and ChIP-seq data were respectively obtained from the GTEx Analysis V6p release (www.gtexportal.org), Ensembl Biomart (www.ensembl.org) and ENCODE (www.encodeproject.org). The CRISPR- and siRNA-generated knockdown data were obtained from the supplementary information files of Dixit et al.¹⁶ and Cusanovich et al.¹⁴. The source datasets generated and/or analysed by this framework, along with sample results and compiled software are available from Zenodo. DOI: https://doi.org/10.5281/zenodo.1707423²¹.

Data are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).

Extended data

Additional files are available from Zenodo. DOI: https://doi.org/10.5281/zenodo.2611953²².

Additional file 1. The mathematical definitions of the four other similarity metrics, the workflow of the IDBC algorithm, an example feature vector, the mathematical definitions of five statistical variables to measure classifier performance, the default parameter values of classifiers in MATLAB, and histograms visualizing the first two criteria for selecting positives from the CRISPR-generated knockdown data.

Additional file 2: The lists of positives and negatives in the ML classifiers to predict genes with similar tissue-wide expression profiles to NR3C1.

Additional file 3: The lists of positives and negatives in the DT classifier to predict TF targets based on the CRISPR-generated knockdown data.

Additional file 4: The lists of positives and negatives in the DT classifier to predict DE direct targets based on the siRNA-generated knockdown data.

Additional file 5: The performance of the DT classifier using only TFBS counts, accuracy of each round of 10-fold cross validation, and Gini importance values of the ML features.

Additional file 6: The list of the most similar 500 PC genes to each TF in terms of tissue-wide expression profiles, and the intersection of these 500 genes and target genes of the TF.

Additional file 7: Cofactor binding sites adjacent to YY1 and EGR1 sites in the promoters of their targets and non-targets.

Additional file 8: The percentages of positives and negatives whose promoters do not overlap DHSs for the CRISPR-perturbed TFs.

Data are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).