{'fraction_non_alphanumeric': 0.031279790276453764, 'fraction_numerical': 0.017218779790276455, 'mean_word_length': 4.811980609418282, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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}}
Identifying Cell Type-Specific Transcription Factors by Integrating ChIP-seq and eQTL Data–Application to Monocyte Gene Regulation

We describe a novel computational approach to identify transcription factors (TFs) that are candidate regulators in a human cell type of interest. Our approach involves integrating cell type-specific expression quantitative trait locus (eQTL) data and TF data from chromatin immunoprecipitation-to-tag-sequencing (ChIP-seq) experiments in cell lines. To test the method, we used eQTL data from human monocytes in order to screen for TFs. Using a list of known monocyte-regulating TFs, we tested the hypothesis that the binding sites of cell type-specific TF regulators would be concentrated in the vicinity of monocyte eQTLs. For each of 397 ChIP-seq data sets, we obtained an enrichment ratio for the number of ChIP-seq peaks that are located within monocyte eQTLs. We ranked ChIP-seq data sets according to their statistical significances for eQTL overlap, and from this ranking, we observed that monocyte-regulating TFs are more highly ranked than would be expected by chance. We identified 27 TFs that had significant monocyte enrichment scores and mapped them into a protein interaction network. Our analysis uncovered two novel candidate monocyte-regulating TFs, BCLAF1 and SIN3A. Our approach is an efficient method to identify candidate TFs that can be used for any cell/tissue type for which eQTL data are available.


Introduction
Mapping human cell type-specific gene regulatory networks is a significant and longstanding challenge due to the large number (.1,200) of transcription factors (TFs) and the expansiveness of noncoding regions of the genome. 1 Expression quantitative trait locus (eQTL) studies, in which expression single nucleotide polymorphisms (eSNPs) that are statistically associated with population variation in transcript abundance are systematically mapped, are an increasingly used tool for revealing functional regulatory regions in cell types of interest. 2 The genomic regions that are in linkage disequilibrium (LD) with eSNPs include cell type-specific enhancers (and other types of cis-regulatory elements) that contain TF binding sites (TFBS). However, in order to map cell type-specific regulatory regions to specific TFs, additional information is needed. One approach is to bioinformatically scan-within regulatory regions-genomic sequence for matches to known TFBS sequence patterns 3,4 ; however, this approach has a significant false-positive rate 5 and TFBS sequence preferences may vary depending on the cellular context. 6 Following the recent availability of published binding site locations (based on chromatin immunoprecipitation-to-tag-sequencing [ChIP-seq]) for 133 TFs in various human cell lines through the Encyclopedia of DNA Elements (ENCODE) project, 7 we devised a novel method to directly analyze eQTL regulatory regions using ChIP-seq-based cellular TFBS measurements. Conceptually, our approach builds on previous efforts in which eQTL data were leveraged in motif-based identification of TFBS 8 and in identifying genomic correlates of regulatory SNPs. 9 We applied the method to search for TF regulators using eQTLs from two population studies of human monocytes. 10,11 Monocytes are innate immune cells that are vital in host defense, have important roles in many infectious and chronic inflammatory diseases, and are of research interest as models of cellular gene regulatory responses to microbial challenge. [12][13][14][15] We hypothesized that binding sites of monocyteregulating TFs are more enriched within monocyte eQTLs than the TFBS that are not monocyte regulators. To investigate this hypothesis, we ranked ENCODE-profiled TFBS ChIP-seq experiments based on the degrees of overlap of their genome-wide binding site profiles (in human cell lines) with monocyte eQTLs. We then tested the enrichment of myeloid cell-derived data sets in the ranked list and the enrichment of known monocyte-regulating TFs in the resulting list of TFs. We found that TFs that are known to have regulatory functions in monocytes are concentrated toward significant P values and are well connected to other candidate monocyte TFs in a protein interaction network. Our study uncovered two candidate monocyte-regulating TFs, BCLAF1 and SIN3A, which are previously known to have functions related to immunity.

results and discussion
To investigate whether anchoring the analysis of TFs using monocyte eQTL data would enable prioritizing TFBS data sets by their relevance to monocytes, we ranked the 397 ChIP-seq-derived TF or cofactor binding site (hereafter, "TFBS") location data sets based on their overlap (as measured by a z-score of the number of TFBS within the eQTLs; z was computed using a background model based on genomewide-randomized eQTL placement; "Methods" section) with eQTL regions defined by LD with 80,000 monocyte eSNPs. Then, for two cell types, a cell type related to monocytes (the human myelogenous leukemia cell line K562 16 ) and a negative control cell type (the human embryonic stem cell line H1-hESC), we measured the degree to which the data sets of the cell types are enriched at low ranks (corresponding to high overlap with monocyte eQTLs), using the GSEA statistical test. 17 The ChIP-seq data sets from K562 cells were biased toward low ranks (corresponding to high overlap with monocyte eQTLs), with a GSEA score E = 0.359 ( Fig. 1A; E ≠ 0 at P , 0.001); by comparison, the ChIP-seq data sets from H1-hESC cells were biased toward higher ranks (corresponding to low overlap with monocyte eQTLs), with E = −0.316 (Fig. 1B). From this analysis, we concluded that our method enables the prioritization of cell type-specific TFBS data sets based on eQTL data. Next, to investigate the power of our method to uncover specific TFs that regulate gene expression in a specific cell type, for each TFBS data set, we scored the significance of the number of TFBS that were within eQTLs using a one-tailed binomial test in which the background probability is based on the genome-widerandomized background model ("Methods" section). We adjusted the 397 resulting P values for multiple hypothesis testing and for each of the 133 unique TF types that we analyzed (Supplementary Table 1), we selected the TFBS ChIPseq experiment for that TF that had the largest effect size (enrichment ratio) based on the binomial test. a This analysis produced a list of 133 TFs that we ranked by adjusted P value (P adj ) (Supplementary Table 2). We hypothesized that TFs that regulate gene expression in monocytes would be biased toward low ranks. To test this hypothesis, we reviewed the lite rature in order to compile a list of 71 TFs that are known to regulate gene expression in monocytes or monocyte-derived immune cells (see Supplementary Table 3), of which 26 TFs a It is acknowledged that this selection criterion will somewhat bias the P adj values; however, we note that the highest P adj among the top 20% of the TFs is less than 10, and thus, there is little risk of a false rejection of the null hypothesis assuming the binomial test assumptions hold.  Table 1 for a complete list of the 397 ChIP-seq experiments whose data were analyzed, including cell type and TF target for each experiment). Horizontal coordinates of the data points show rank (z) for ChIP-seq data sets for the indicated cell line, with lowest rank corresponding to highest z. Here, consistent with the definition in the original GSEA study, 17 P hit is the cumulative density function (CDF) of an indicator function that is 1 for a K562 data set and 0 for a non-K562 data set, and P miss is the CDF for an indicator function that is 0 for a K562 data set and 1 for a non-K562 data set.
were among the 133 TFs profiled by ENCODE. We then tested for enrichment of the known monocyte-regulating TFs at low ranks in the list, using Gene Set Enrichment Analysis (GSEA). 17 We found that the known monocyte-regulating TFs are strongly biased toward low ranks ( Fig. 2; E = 0.329; P , 0.01).
Several of the NMTFs have functions related to immunity or interact with previously identified monocyte-regulating TFs, supporting their potential roles in monocyte gene regulation. For example, the NMTF with the smallest P adj , BCLAF1, is highly expressed in various immune cell types (based on data from ImmGen and BioGPS; "Methods" section), and it plays a crucial role in signal transduction between NF-κB and C/EBPβ (both of which are known to regulate gene expression in monocytes) in multiple cell types. 22 Another significant NMTF, SIN3A, plays a role (along with the known monocyte-regulating TF, Myc) in a signaling circuit implicated in the pathogenesis of acute myeloid leukemia. 23 The NMTFs TAF1 and TAF7 may function together in monocytes, as they are coexpressed 24 and are both components of the transcription factor II D (TFIID) complex. TAF7 is highly expressed in immune cells and in monocytes in particular (ImmGen and BioGPS). The NMTF SP2 (along with other known monocyte-regulating TFs, STAT1, STAT3, IRF1, and EGR1) is among the top 10 TFs associated with response to parasitic infection. 25 The NMTF YY1 regulates the activity of the promoter of murine IRF3, which plays a key role in inflammation and immune response. 26 THAP1, ETS1, TR4, and NRF1 have been previously identified to have specific immune regulatory functions. THAP1 is known to play a role in inducing T-cell apoptosis. 27 STAT3 is a known regulator of ETS1 in inflammation control in mouse macrophages. 28 In addition, a TR4-CD36 pathway is used by mouse miR-133a to regulate lipid accumulation in macrophages. 29 One general transcription factor, TBP whose transcript abundance is thought to be stable in human monocytes, 30 is also in the top 20%; one possible mechanism   could be TFIID recruitment to TATA-less promoters by cell type-specific regulators. 31

discussion and conclusions
We have demonstrated an efficient method to identify candidate TFs that regulate gene expression in a specific cell type of interest, by leveraging publicly available human genome-wide location data sets and eQTL data sets. A key advantage of this approach is that it can leverage the extensive array of publicly accessible ChIP-seq data from the ENCODE project (1,485 TF ChIP-seq data sets as of the current date). We note that one limitation of the approach is that human ChIP-seq data sets for many TFs are at present primarily available from cancerous cell lines, 7 whose gene regulatory networks may be altered compared to noncancerous cells of the same lineage type. Nevertheless, the strong result for enrichment of known monocyte TF regulators in the list of TFs ranked using our method suggests that this method is useful to screen for TF regulators for noncancerous cell types. Undoubtedly, confirmatory evidence regarding a specific function for a candidate TF in a given cell type of interest, (such as gene expression as we used in this work), would have an important role in avoiding false-positive identification of a TF due to a cancer cell line-specific artifact. From our analysis using monocyte eQTL data, we conclude that TFs with known functions in regulating gene expression in monocytes tend to have binding sites located in the vicinity of monocyte eQTLs. Of the NMTFs that were identified in our analysis, the two most significant TFs were BCLAF1 and SIN3A. As revealed by protein network reconstruction, the top 20% of TFs identified by our approach are highly interconnected. Further targeted experiments would be useful to elucidate their specific roles in monocyte gene regulation. The approach is generally applicable to other cell types or tissues in that it could be used to identify candidate TFs for any cell type or tissue type for which eQTLs have been mapped from human population studies.

Methods
Software: We used the R statistical computing software (v3.2.2) and Bioconductor (v3.2) packages on OS X 10.11.4. databases. We obtained 397 TFBS location data files (derived from ChIP-seq analysis of 133 TFs in 68 human cell lines) from the public repository of the ENCODE project 7 via the University of California Santa Cruz Genome Browser portal. Additionally, we obtained data files containing eSNPs and associated eQTL P values from online supplementary data for two large-cohort monocyte eQTL studies involving 1,773 unrelated individuals. 10,11 We carried out all genome analyses using coordinates from the GRCh37/hg19 human reference genome assembly. We identified mouse orthologs expressed in monocytes and other immune cell types of the top TFs using the Immunological Genome (ImmGen) project database. 32 We used the gene annotation portal, BioGPS, 33 to obtain human microarray-based gene expression measurements for the top-ranked TFs. eQtL intervals. We mapped genomic intervals containing SNPs that are in LD with each eSNP using the SNP Annotation and Proxy (SNAP) tool 34 using haplotype data from the 1,000 Genomes Project, 35 including only SNPs with R 2 $ 1.0 and distance #10 kbp.
tFbs/eQtL overlap count significance testing. We created a set of allowed regions for random eQTL LD block placement (for use as our background model for TFBS/eQTL overlap counts) by combining peak regions from all 397 TF ChIP-seq data sets with the monocyte eQTLs and merging together any regions whose edges were within 1.9 kbp of one  (Table 1); shading denotes a known role for the TF in monocyte gene regulation (Supplementary Table 3). Edges denote interactions or coexpression ("Methods" section) as indicated in the legend. A dashed edge means that the two TFs are predicted (by the GeneMANIA database) to physically interact ("predicted interactions"). (b) Hive plot 21 representation of the network. TFs are arranged into "source" (red), "intermediate" (blue), and "sink" (green) axes based on the graph topology. Edge thickness represents the interaction type (thick arc = physical interaction; medium = predicted interaction; thin = coexpression). For example, in the hive plot, NF-κB is identified as an upstream regulator in the monocyte gene regulatory network.
another (selected so that the resulting set of allowed regions covered ∼25% of the genome). As a background model for the counts of TFBS/eQTL overlaps, we placed all eQTL LD blocks into random nonoverlapping locations within the set of allowed regions. We carried out 1,000 iterations of the LD block randomization. For each TF ChIP-seq data set, the count of peak center coordinates that overlapped with an eQTL LD block was computed, for both the human eQTL LD blocks and (as a background model) 1,000 iterations of the randomly placed LD blocks. For each ChIP-seq data set, we computed z-scores for the overlap counts using the mean and standard deviation of the background model-based counts (where z .. 0 indicates that the count of regions of overlap between the eQTLs and the ChIP-seq peaks is higher than would be expected based on the overlap count distribution of the null model where the eQTL LD blocks are randomly placed). For scoring the enrichment of cell types among ChIP-seq data sets ranked by z (in descending order) as well as scoring the enrichment of monocyte-regulating TFs among TFs ranked by adjusted P value (P adj ) (in ascending order), we used the gene set enrichment analysis (GSEA) statistical test 17 with P = 0 (unweighted) to obtain an enrichment score (E score). For testing cell types among ChIP-seq data sets, a higher E score means that the ChIP-seq experiments corresponding to the cell type or TF have lower ranks for their z-scores (ie, higher z-scores) than would be expected by chance. For testing monocyte-regulating TFs among all TFs ranked by P adj , a higher E score means that the monocyte-regulating TFs have lower P adj value ranks than would be expected by chance. To test for E ≠ 0 significance, we used a permutation-based approach 36 with at least 1,000 iterations. For filtering TFs, we computed the average (over the 1,000 iterations) fraction of ChIP-seq peaks whose centers overlap one of the randomly placed LD blocks; the probability was used in a one-tailed binomial test. We adjusted the binomial P values for multiple hypothesis testing using the Benjamini-Hochberg false discovery rate method. 37 tF network construction and visualization. We used GeneMANIA 38 (4/16/2016 release) to construct the interaction network for the top 27 TFs (by P value). We used Cytoscape 39 v3.4.0 for network layout and visualization. We used the HiveR R package 40 for generating the hive plot.
supplementary information. The software source code for this method is available under a free-software, open-source (Apache 2.0) license at github.com/ramseylab/eqtlchiptest.