{'fraction_non_alphanumeric': 0.028196003284971256, 'fraction_numerical': 0.04641977239842009, 'mean_word_length': 5.0354024073637005, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Clustered, information-dense transcription factor binding sites identify genes with similar tissue-wide expression profiles

Background The distribution and composition of cis-regulatory modules (e.g. transcription factor binding site (TFBS) clusters) in promoters substantially determine gene expression patterns and TF targets, whose expression levels are significantly regulated by TF binding. TF knockdown experiments have revealed correlations between TF binding profiles and gene expression levels. We present a general framework capable of predicting genes with similar tissue-wide expression patterns from activated or repressed TF targets using machine learning to combine TF binding and epigenetic features. Methods Genes with correlated expression patterns across 53 tissues were identified according to their Bray-Curtis similarity. DNase I HyperSensitive region (DHS) -accessible promoter intervals of direct TF target genes were scanned with previously derived information theory-based position weight matrices (iPWMs) of 82 TFs. Features from information density-based TFBS clusters were used to predict target genes with machine learning classifiers. The accuracy, specificity and sensitivity of the classifiers were determined for different feature sets. Mutations in TFBSs were also introduced to examine their impact on cluster densities and the regulatory states of predicted target genes. Results We initially chose the glucocorticoid receptor gene (NR3C1), whose regulation has been extensively studied, to test this approach. SLC25A32 and TANK were found to exhibit the most similar expression patterns to this gene across 53 tissues. Prediction of other genes with similar expression profiles was significantly improved by eliminating inaccessible promoter intervals based on DHSs. A Random Forest classifier exhibited the best performance in detecting such coordinately regulated genes (accuracy was 0.972 for training, 0.976 for testing). Target gene prediction was confirmed using CRISPR knockdown data of TFs, which was more accurate than siRNA inactivation. Mutation analyses of TFBSs also revealed that one or more information-dense TFBS clusters in promoters are required for accurate target gene prediction. Conclusions Machine learning based on TFBS information density, organization, and chromatin accessibility accurately identifies gene targets with comparable tissue-wide expression patterns. Multiple, information-dense TFBS clusters in promoters appear to protect promoters from the effects of deleterious binding site mutations in a single TFBS that would effectively alter the expression state of these genes.


ABSTRACT 21
Background: The distribution and composition of cis-regulatory modules (e.g. transcription 22 factor binding site (TFBS) clusters) in promoters substantially determine gene expression 23 patterns and TF targets, whose expression levels are significantly regulated by TF binding. TF 24 knockdown experiments have revealed correlations between TF binding profiles and gene 25 expression levels. We present a general framework capable of predicting genes with similar 26 tissue-wide expression patterns from activated or repressed TF targets using machine learning 27 to combine TF binding and epigenetic features. promoter intervals of direct TF target genes were scanned with previously derived information 31 theory-based position weight matrices (iPWMs) of 82 TFs. Features from information density-32 based TFBS clusters were used to predict target genes with machine learning classifiers. The 33 accuracy, specificity and sensitivity of the classifiers were determined for different feature sets. 34 Mutations in TFBSs were also introduced to examine their impact on cluster densities and the 35 regulatory states of predicted target genes. 36 Results: We initially chose the glucocorticoid receptor gene (NR3C1), whose regulation has 37 been extensively studied, to test this approach. SLC25A32 and TANK were found to exhibit the 38 most similar expression patterns to this gene across 53 tissues. Prediction of other genes with 39 similar expression profiles was significantly improved by eliminating inaccessible promoter 40 intervals based on DHSs. A Random Forest classifier exhibited the best performance in 41 detecting such coordinately regulated genes (accuracy was 0.972 for training, 0.976 for testing). 42 Target gene prediction was confirmed using CRISPR knockdown data of TFs, which was more 43 accurate than siRNA inactivation. Mutation analyses of TFBSs also revealed that one or more 44 information-dense TFBS clusters in promoters are required for accurate target gene prediction. 45 TF binding profiles, either derived from in vivo ChIP-seq peaks [6][7][8] or computationally 70 detected binding sites and CRMs [9], have been shown to be predictive of absolute gene 71 expression levels using a variety of tissue-specific machine learning classifiers and regression 72 models. Because signal strengths of ChIP-seq peaks are not strictly proportional to TFBS 73 strengths [3], representing TF binding strengths by ChIP-seq signals may not be appropriate; 74 nevertheless, both achieved similar accuracy [10]. CRMs have been formed by combining two 75 or three adjacent TFBSs [9], which is inflexible, as it arbitrarily limits the number of binding sites 76 contained in a module, and does not consider differences between information densities of 77 different CRMs. Chromatin structure (e.g. histone modification (HM) and DNase I 78 hypersensitivity) were also found to be highly redundant with TF binding in explaining tissue-79 specific mRNA transcript abundance at a genome-wide level [7,8,11,12], which was attributed to 80 the heterogeneous distribution of HMs across chromatin domains [8]. Combining these two 81 types of data explained the largest fraction of variance in gene expression levels in multiple cell 82 lines [7,8], suggesting that either contributes unique information to gene expression that cannot 83 be compensated for by the other. 84 The number of genes directly bound by a TF significantly exceeds the number of genes 85 whose expression levels significantly change upon knockdown of the TF. Only a small subset of 86 genes whose promoters overlap ChIP-seq peaks were differentially expressed (DE) after 87 individually knocking 59 TFs down using small interfering RNAs (siRNAs) in the GM19238 cell 88 line [13]. Correlation between TFBS counts and gene expression levels across 10 different cell 89 lines among 8,872 genes from these knockdown data were more predictive of DE targets than 90 setting a minimum threshold on TFBS counts [14]. Their TFBS counts were defined as the 91 number of ChIP-seq peaks overlapping the promoter, though it was unknown how many binding 92 sites were present in these peaks; true positives might not be direct targets in the TF regulatory 93 cascade, as the promoters of these targets were not intersected with ChIP-seq peaks. verified predicted direct TF targets whose promoters overlap tissue-specific ChIP-seq peaks, in 116 contrast with correlation-based approaches [14]. 117

METHODS 118
To identify genes with similar tissue-wide expression patterns, we formally define gene 119 expression profiles and pairwise similarity measures between profiles of different genes. A 120 general machine learning framework relates features extracted from the organization of TFBSs 121 in these genes to their tissue-wide expression patterns. True positives (TPs) and negatives 122 (TNs) for predicting direct DE TF targets were validated using CRISPR-and siRNA-generated 123 knockdown data (see below). 124

Similarity between gene expression profiles 125
The median RPKM (Reads Per Kilobase of transcript per Million mapped reads) of 56,238 126 genes across 53 tissues were obtained from the Genotype-Tissue Expression (GTEx) project 127 [18]. To capture the tissue-wide overall expression pattern of a gene instead of within a single 128 tissue, the expression profile of a gene was defined as its median RPKM across the 53 tissues, 129 which forms a vector of size 53 and does not distinguish between different isoforms whose 130 expression patterns may significantly differ from each other. To obtain ground-truth genes that 131 have similar expression profiles to a given gene, the Bray-Curtis Similarity (Equation 1) was 132 used to compute the similarity value between the expression profiles of two genes, because it 133 takes both the directions and lengths of the vectors into account while maintaining strict bounds 134 of 0 and 1. 135

Prediction of genes with similar expression profiles 140
The framework for identifying genes that have similar expression profiles to a specific gene is 141 shown in Figure 1A  The information density-related features derived from each TFBS cluster include: 1) The 148 distance between this cluster and the transcription start site (TSS); 2) The length of this cluster; 149 3) The information content of this cluster (i.e. the sum of Ri values of all TFBSs in this cluster); 4) 150 The number of binding sites of each TF within this cluster; 5) The number of strong binding sites 151 TFs, we have previously derived iPWMs exhibiting primary binding motifs for seven (EGR1, 167 ELF1, ELK1, ETS1, GABPA, IRF1, YY1) [3]. Therefore, the framework for predicting direct TF 168 targets in the K562 cell line ( Figure 1A and 1C) was applied to these TFs. The criteria for 169 defining a TP (i.e. a DE direct target), of a TF was: 170 1) The fold change in the expression level of this gene for each guide RNA of the TF was > 171 (or <) 1, consistent with the possibility that the gene was regulated by the TF, and 172 2) The average fold change in the expression level of this gene for all guide RNAs of the TF 173 was > threshold (or < 1⁄ ), and 174 3) The promoter interval (10 kb) upstream of a TSS of this gene overlaps a ChIP-seq peak of 175 the TF in the K562 cell line. 176 If the coefficients of all guide RNAs of the TF for a gene are zero, the gene was defined as a 177 TN. As the threshold ε increases, the number of TPs strictly decreases; as ε decreases, we 178 have increasingly lower confidence in the fact that the TPs were indeed differentially expressed 179 because of the TF perturbation. To achieve a balance between sensitivity and specificity, we 180 test sets ( Figure 1D). 212 The DHSs in the GM19238 cell line mapped from the hg19 genome assembly were first 213 remapped to the hg38 assembly using liftOver (available at genome.ucsc.edu) [23]. Aside from 214 RELA and NFE2L1, the iPWMs of TFs from the GM12878 cell line were used to detect binding 215 sites. For RELA, the iPWM from the GM19099 cell line was used; for NFE2L1, the only 216 available iPWM was derived from K562 cells and was applied. Although the knockdown was 217 performed in GM19238, GM12878 and GM19099 are also lymphoblastic cell lines, with 218 GM19099 and GM19238 both being derived from Yorubans. For this analysis, the iPWMs 219 derived in GM12878 and GM19099 were more appropriate than the iPWM from K562, since 220 GM12878 and GM19099 are of the same tissue type and are thus more likely comparable to 221 GM19238 than to K562. 222

Mutation analyses on promoters of differentially expressed direct targets 223
To better understand the significance of individual binding sites for information-dense 224 clusters and the regulatory state of direct targets, we evaluated the effects of sequence changes 225 that altered the Ri values of these sites on cluster formation and whether a gene was predicted 226 to be a TF target. Mutations were sequentially introduced into the strongest binding sites in 227 TFBS clusters of the EGR1 target gene, MCM7, to determine the threshold for cluster formation 228 after disappearing clusters disabled induction of MCM7 expression. For one target gene of each 229 TF from the CRISPR-generated perturbation data, effects of naturally occurring TFBS variants 230 present in dbSNP (https://www.ncbi.nlm.nih.gov/projects/SNP/) [24] were also evaluated to 231 explore aspects of TFBS organization that enabled both clusters and promoter activity to be 232 resilient to binding site mutations. This was done by analyzing whether the occurrence of 233 individual or multiple single nucleotide polymorphisms (SNPs) lead to the loss of binding sites 234 and the clusters that contain them, and result in changes in the predictions of these targets. 235

Similarity between gene expression profiles 237
To confirm that the Bray-Curtis Similarity can indeed effectively measure how akin the 238 expression profiles of two genes are to each other, it was applied to compute the similarity 239

Prediction of genes with similar expression profiles 255
The framework for predicting genes with similar expression profiles was based on promoter 256 scans with each TFBS, followed by the derivation of the spatial density-and information density-257 related features from clusters in each promoter for genes with an NR3C1-like expression pattern 258 (as shown in Figure 1A and 1B). We investigated two versions of this framework, depending on 259 whether promoter sequences were first intersected with DHSs. Under both scenarios, all 260 classifiers (Naïve Bayes, two types of Decision trees and three types of Support vector 261 machines (SVM)) were applied to both the training and test sets, successfully distinguishing 262 similar from dissimilar genes in terms of expression profiles (i.e. accuracy, sensitivity and 263 specificity all > 0.5) ( Table 1, Additional file 5). We found, however, that generally all TFBSs in a 264 DHS formed a binding site cluster, and the performance of all classifiers were significantly 265 improved by inclusion of DHS information (i.e. accuracy, sensitivity and specificity were all 266 increased) ( Table 1, Additional file 5). The SVM classifier with the RBF kernel and the Random 267 Forest classifier were the only two classifiers with accuracies exceeding 0.97, and each 268 performed equally well on both the training and test sets (Table 1). 269

Prediction of differentially expressed direct TF targets 270
Between the two classifiers with the best performance in distinguishing genes with similar 271 expression profiles to NR3C1 from others (i.e. SVM with RBF kernel and Random Forest), we 272 After eliminating TFBSs in inaccessible promoter intervals, i.e. those excluded from tissue-276 specific DHSs, the RF classifier predicted direct targets with greater accuracy and specificity 277 (Table 2 and 3, Additional file 5). Specifically, predictions based on CRISPR-generated 278 knockdown data for TFs: EGR1, ELK1, ELF1, ETS1, GABPA, and IRF1 were more accurate 279 than for YY1, which itself represses or activates a wide range of promoters by binding to sites 280 overlapping the TSS (Table 2, Additional file 5). Accordingly, the perturbation data indicated that 281 YY1 has ~3-23 times more targets in the K562 cell line than the other TFs (ε = 1.05), and its 282 binding has a more significant impact on the expression levels of target genes (for YY1, the ratio 283 of the target counts at ε = 1.1 vs ε = 1.01 was 0.328, which significantly exceeded those of the 284 other TFs (0.019-0.081); Additional file 3). This is concordant with our previous finding that YY1 285 extensively interacts with 11 cofactors (e.g. DNA-binding IRF9 and TEAD2; non-DNA-binding 286 DDX20 and PYGO2) in K562 cells, consistent with a central role in specifying erythroid-specific 287 Despite a high accuracy of target recognition, sensitivity was consistently higher than 289 specificity (Table 2, Additional file 5), implying that the classifier more effectively identified direct 290 targets compared to non-targets. This is attributable to the fact that the promoters of false 291 positive target genes also contain accessible, but non-functional TFBSs. In vivo co-regulation 292 mediated by interacting cofactors, which was excluded by the classifier, assisted in 293 distinguishing these non-functional sites that do not significantly affect gene expression [3,13]. between direct targets and non-targets (Additional File 1). Note that this inference holds valid 303 only when taking all direct targets and non-targets of a TF into account; it may not be true for a 304 specific pair of genes (i.e. the promoter of a gene that is not a DE target may contain a greater 305 number of accessible, but non-functional TFBSs) (Additional File 1). We noted this trend on the 306 test sets of only ELF1 and IRF1 ( Figure 3); for the other five TFs (EGR1, ELK1, ETS1, GABPA, 307 YY1), differences in the clustered TFBS counts between targets and non-targets did not 308 necessary increase with larger values of , since the test set consists of both targets and non-309 targets in equal proportions (Additional File 1). However, the classifier performed well in each 310 instance, because the count differences were still sufficiently large to discriminate between 311 targets and non-targets (Figure 3). 312 With the siRNA-generated knockdown data, the performance of the RF classifier was 313 compared to an approach inferring DE targets by correlating TF binding with gene expression 314 levels across ten cell types [14]. In this correlation-based approach, three measures (i.e. the 315 absolute Pearson correlation coefficient (PC), the absolute Spearman correlation coefficient 316 (SC), and the absolute combined angle ratio statistic (CARS)), whose performance was 317 evaluated with precision-recall curves, were alternatively used to compute a correlation score 318 between the number of ChIP-seq peaks overlapping the promoter and gene expression values. 319 Genes predicted to be DE targets had scores above the threshold resulting in a 1.5-fold 320 increase compared to the background precision. For example, in the case of the TF YY1, which 321 was analyzed by both approaches, the performance of the RF classifier on the training set was 322 0.66 (precision) and 0.456 (recall), and the test set was 0.672 and 0.396 (Table 3) intersected ChIP-seq peaks over shorter 5kb promoter intervals upstream of the TSS. 327

Intersection of genes with similar expression profiles and direct targets 328
To determine how many direct targets have similar tissue-wide expression profiles, we 329

Mutation analyses on promoters of direct targets 361
Because the promoters of most direct targets contain multiple binding site clusters, we 362 anticipate that this enables these genes' expression to be naturally robust against binding site C->G mutations abolishing the strongest site in either cluster, which altered the prediction, that 385 is, EGR1 lost the capability to induce MCM7 transcription (Figure 4). The remaining four sparse 386 weak sites do not form a cluster and cannot completely supplant the disrupted strong sites. 387 Further, we examined the impacts of known natural SNPs on binding site strengths, clusters 388 and the regulatory state of the promoter for a direct target of each of the seven TFs from the 389 CRISPR-generated perturbation data (Table 5). Often a single SNP (e.g. rs996639427 of EGR1) 390 can affect the strengths of multiple binding sites (Table 5). Apart from SNPs that are predicted 391 to abolish binding (Figure 4), leaky variants that merely weaken TF binding are common ( Table  392 5). Binding stabilization between adjacent sites and the funnel effect enable the CRMs 393 comprised of information-dense clusters to be robust to mutations in individual binding sites. In 394 this way, neither mutations that abolish TFBSs nor leaky SNPs in flanking weak sites can 395 destroy functional clusters (e.g. rs1030185383 and rs5874306 of IRF1), whereas SNPs with 396 large reductions in Ri values of central strong sites are more likely to abolish clusters (e.g. 397 rs865922947, rs946037930, rs917218063 and rs928017336 of YY1) (Table 5). More generally, 398 the presence of multiple clusters enables promoters to be effectively resilient to the effects 399 binding site mutations; only the complete abolishment of all clusters resulting from the 400 simultaneous occurrence of multiple SNPs can transform the promoter to be unresponsive to TF 401 binding to residual weak sites (e.g. rs997328042, rs1020720126 and rs185306857 of GABPA) 402 (Table 5). Furthermore, a relatively small number of SNPs that strengthen TF binding and 403 eventually amplify the regulatory effect of the TF on the gene expression level are also present 404 (e.g. rs887888062 of EGR1 and rs751263172 of ELF1) ( Table 5), suggesting that, in addition to 405 deleterious mutations, benign variants may also be found in promoters, consistent with the 406 expectations of neutral theory [35]. 407

DISCUSSION 408
In this study, the Bray-Curtis Similarity function was initially shown (for the NR3C1 gene) to 409 measure the relatedness of overall expression patterns between genes across a diverse set of The relatively poor performance of the classifier on YY1 (Table 2) is attributable to its smaller 433 percentage of TNs with inaccessible promoters (Additional file 8). Additionally, the Random 434 Forest classifier was more predictive of functional TF binding on the CRISPR-generated 435 knockdown data than the siRNA-generated ones (Table 2 and 3). This larger discrepancy in 436 feature vectors between TPs and TNs from CRISPR-based perturbations is also attributable to 437 the greater differences in the percentages between TPs and TNs with inaccessible promoters 438 (Additional file 8). Among the 22,046 genes whose expression levels were measured in the 439 CRISPR-based perturbations, most of the TNs with inaccessible promoters merely have one 440 transcript and specific functions (e.g. VENTXP1 for the TF, EGR1), whereas many such TNs 441 were excluded from the 8,872 genes whose knockdown data were generated by siRNA 442

inactivation. 443
Our mutation analyses revealed that some deleterious TFBS mutations could be 444 compensated for by other information-dense clusters in a promoter; thus predicting the effects 445 of mutations in individual binding sites would not be sufficient to interpretation of downstream 446 effects. Though compensatory clusters may maintain gene expression, the promoter will provide 447 lower levels of activity than the wild-type promoter could, which is a recipe for achieving natural 448 Apart from these TFs, the Bray-Curtis Similarity can be directly applied to identify the ground-455 truth genes with overall similar tissue-wide expression patterns to any other gene whose 456 expression profile is known. Further studies could investigate the biological significance 457 underlying the phenomenon that all these genes share a common expression pattern, including 458 the similarity between other regulatory regions besides proximal promoters in terms of TFBSs 459 and epigenetic markers. This machine learning framework can also be applied to predict direct 460 DE targets for other TFs and in other cell lines, depending on the availability of corresponding 461 knockdown data. 462 There are a number of limitations of our approach. The Bray-Curtis function seems unable to 463 accurately measure the similarity between gene expression profiles of a ubiquitously expressed 464 gene (e.g. NR3C1) and a tissue-specific gene (e.g. stomach-specific PGA3), which exhibit quite 465 different tissue-wide expression patterns (i.e.

Ethics approval and consent to participate 536
Not applicable 537

Consent for publication 538
Not applicable 539

Availability of data and materials 540
The median RPKM, TSS coordinate, DNase I hypersensitivity and ChIP-seq data are 541 respectively available from the GTEx Analysis V6p release (www.gtexportal.org), Ensembl 542 Biomart (www.ensembl.org) and ENCODE (www.encodeproject.org). The CRISPR-and siRNA-543 generated knockdown data are available from the supplementary information files of Dixit   In the IDBC algorithm (Additional file 1), the parameter I is the minimum threshold on the total 656 information contents of TFBS clusters. In prediction of genes with similar expression profiles, 657 the minimum value was 939, which was the sum of mean information contents (Rsequence values) 658 of all 94 iPWMs; in prediction of direct targets, this value was the Rsequence value of the single 659 iPWM used to detect TFBSs in each promoter. The parameter d is the radius of initial clusters in 660 base pairs, whose value, 25, was determined empirically. Eight types of three different 661 classifiers were evaluated with statistics (accuracy, sensitivity and specificity) to measure the 662 classifier performance (Additional file 1). B) Formation of the training and test sets for identifying 663 genes with similar expression profiles to a given gene (Additional file 2). C) Formation of the 664 training and test sets for predicting direct targets of seven TFs using the CRISPR-generated 665 perturbation data in K562 cells (Additional file 3). D) Formation of the training and test sets for 666 predicting direct targets of 11 TFs using the siRNA-generated knockdown data in GM19238 667 cells (Additional file 4). When genes with single transcripts were more than the TPs, those with 668 the largest P-values were selected as TNs (null hypothesis of differential expression cannot be 669 rejected); when genes with single transcripts were fewer than the TPs, those genes with two 670 transcripts and the largest P-values were also selected. This step was iterated until the number 671 of TNs equaled that of TPs. 672