A large proportion of TSSs are found in regions of accessible chromatin

To understand how regions of open chromatin vary between cell types, we inspected the degree to which DHSs were shared in the 19 cell types. A DHS was classified as being specific to a cell line if it was only present in a single cell type or overlapped <50% of its length with a DHS from any of the other 18 cell types (Fig. 2A). Across all DHSs, ~14% were specific to a single cell line (Fig. 2B). Intergenic DHSs showed the highest percentage of being cell-type–specific (~17%). Conversely, TSS DHSs were largely not cell-type–specific with <1% being open specifically in a single cell type. Despite the broad panel of cell lines that vary in expression, the chromatin state at the TSS of these genes was open and largely invariant across multiple cell lines. This is in agreement with a recent study analyzing a subset of the cell types used here (Song et al. 2011).

Open in a separate window

Figure 2.

Cell-type specificity of hypersensitive regions. (A) Example (Chr1: 201,890,462–201,938,914) showing cell-type–specific DHSs across two cell lines (pink boxes). Note that we called a DHS cell-type–specific if it did not overlap another DHS by more than half in any of the 18 other cell lines. (B) Bar graph showing the proportions of cell-type–specific DHSs across different genomic locations averaged across all cell lines. (C) TSSs were divided by the number of cell lines that they overlapped in a region of open chromatin. For each set of TSSs, normalized CG content in the promoter regions (−900,100) of the TSSs are shown.

We determined the normalized CG content in the proximal promoter region of the gene, defining the proximal promoter as −900 to +100 bp around the TSS. If a gene had multiple TSSs, the average of the normalized CG content from each TSS was used. There was a steady positive trend in the number of cell lines in which a DHS overlapped a TSS and the CG content around the TSS (Fig. 2C). Previous studies have reported that gene expression can be predicted from the CG content in the proximal promoter region (Yamashita et al. 2005; Carninci et al. 2006; Zhu et al. 2008). Our result indicates that higher levels of CG dinucleotide content, and thus more frequent presence of CpG islands, are positively correlated with, and could be functioning to preserve, an open chromatin state surrounding the TSS. There were fewer genes with a TSS open in only one cell type (976 genes) and many with an open TSS across all 19 cell types (8393) (Supplemental Table 2).

Cell-type–specific expressed genes show differing patterns of accessible chromatin at their TSS

Gene expression data for the 19 cell lines were generated using Affymetrix exon arrays. Expression values for each gene were transformed to Z-scores across all of the cell lines. Genes with large positive or negative Z-score values thus showed a larger deviation from the mean expression across cell types. The Z-score transformed expression values were used to select subsets of genes with specific expression patterns (Fig. 3A–C; Supplemental Table 3). Up-regulated genes, exemplified by GCM1 (Fig. 3A), had a particularly high expression in one cell line, but expression close to the mean in the other cell lines. To identify genes exhibiting this type of expression pattern, we sorted the Z-score expression for the genes in each cell line. The top 200 genes in this sorted list were classified as being up-regulated in that cell type (UR genes). Down-regulated genes exhibit low expression levels in one cell type but are otherwise constitutively expressed in other cell lines (Fig. 3B; Thorrez et al. 2011). We classified the last 200 genes in the sorted Z-score expression list as being the cell-type–specific down-regulated genes (DR genes). Constitutively expressed genes (Fig. 3C) were identified by filtering all genes that were not in UR and DR gene sets in any cell line and had absolute expression Z-score values < 1.7 in all cell lines. Using this cutoff, 168 genes displayed a pattern of constant expression levels across all cell lines.

Open in a separate window

Figure 3.

Cell-type–specific gene expression and definition of gene classes. (A–C) Representative examples of different patterns of gene expression. Note that Z-score values are calculated from expression across all 19 cell lines. (A) A gene where the expression is specifically up-regulated in the first cell line (UR gene). (B) A gene that is specifically down-regulated in the first cell line (DR gene). (C) A gene that has low variability in expression (constitutively expressed gene). (D) Median expression Z-scores for the genes in each set in each cell line. (E) Normalized CG content from the promoter regions of genes. (F) The fraction of TSS in each gene set that were in a region of open chromatin. E and F share the same color map.

To address how up-regulated genes are expressed in one particular cell type, we grouped UR genes from all other cell types and denoted this group as UR-Other genes (Fig. 3A). We imposed the additional constraint that such genes would show an expression Z-score < 0 in the cell type of consideration, i.e., had expression below its mean expression. As an example, GCM1 (Fig. 3A) was highly expressed in the first cell type and in none of the others shown. It was therefore grouped into the UR class for the first cell type and into the UR-Other class in each of the other cell types. Similarly, genes denoted as DR-Other had to be classified as down-regulated in another cell line and had an expression Z-score > 0 in the cell type of consideration (Fig. 3B). In this way, we defined different classes of transcriptionally active (UR, DR-Other, and Constitutive) and transcriptionally inactive genes (UR-Other, DR) from the point of view of each cell line in comparison to other cell lines.

By definition, UR and DR genes displayed the highest and lowest Z-score gene expression, respectively (Fig. 3D). UR genes were consequently enriched in functions related to the tissue type of origin (Supplemental Table 4). DR-Other and Constitutive genes showed similar expression values and had higher expression values than genes in the UR-Other class. To understand whether these different classes had different properties in sequence composition and chromatin state, we first inspected normalized CG content in the proximal promoter region (Fig. 3E). UR and UR-Other genes had the lowest average CG content compared with the other classes of genes. Constitutively expressed genes displayed a particularly high CG content in their proximal promoter regions, as previously reported (Yamashita et al. 2005; Carninci et al. 2006; Zhu et al. 2008); however, this observation clearly extended to the DR and DR-Other gene classes, which had a CG content slightly lower than constitutively expressed genes.

Constitutively expressed and DR-Other genes had the highest proportion of their TSSs in regions of accessible chromatin (Fig. 3F). DR genes displayed slightly lower chromatin accessibility compared with DR-Other, indicating that repression of DR genes largely occurs while maintaining chromatin accessibility at the TSS. UR genes also showed a high proportion of genes containing an accessible TSS, at similar levels to DR and DR-Other. In stark contrast, UR-Other genes had the lowest fraction of TSSs that overlapped a DHS. These results indicate that even though UR-Other and DR genes are both transcriptionally inactive in a cell type of consideration compared with other cell types, they are likely to be regulated via different chromatin-remodeling mechanisms. Specifically, genes that are up-regulated in a small number of cell types likely maintain a closed chromatin conformation until cellular processes require up-regulation. In contrast, down-regulated genes may be viewed as constitutively expressed genes that are repressed in a single cell type. UR genes had intergenic and gene body DHSs associated with them (Supplemental Fig. 1A,B), in agreement with previous results indicating that cell-type–specific expression is mediated by distal cis-regulatory regions (Song et al. 2011). Overall, these results indicate that different classes of transcriptionally active and inactive genes have different CG content and chromatin accessibility at their TSS.

Classifying tissue-specific expression from sequence features in open chromatin

To predict gene expression patterns from sequence, approaches have frequently used features contained within fixed-size proximal promoter sequences. We used DHS data from a large number of cell types to determine whether using both proximal and distal regulatory regions with open chromatin would improve predictive models for cell-type–specific expression patterns. Position weight matrices (PWMs) for TFs in vertebrates were compiled from TRANSFAC, JASPAR, and UniProbe databases (Matys et al. 2006; Bryne et al. 2008; Newburger and Bulyk 2009). For each DHS, 789 PWMs were used to calculate TFBS scores that accounted for local dinucleotide composition. The maximum sliding window score for each PWM was used as the TFBS score for that DHS (Fig. 4A; Methods). To associate DHSs with specific genes that they are likely to regulate, we applied a simple approach of associating each DHS with the closest TSS (closest gene DHS). For each TF, we then chose the maximum TFBS score across all DHSs associated with a gene (Fig. 4B). As an alternative approach, we split DHSs into distal sites (a set including both Gene-Body and Intergenic DHSs) and TSS DHS sites and used the maximum TFBS in each set as individual features (split DHSs). This doubled the number of features and allowed us to identify different characteristics of TSS-overlapping versus distal DHSs. To compare our models with previous approaches, we also used TFBS features calculated in proximal promoters, defined here as −900 to +100 nt surrounding the TSS (Landolin et al. 2010).

Open in a separate window

Figure 4.

Transcription factor binding site features. (A) DHS and promoter sequences are scanned with PWMs. TFBS scores are log-likelihood ratios of PWM over the background model. A sliding window is used to identify the score for each DHS or promoter. (B) Example to show association of DHSs with genes. Numbers in the brackets are example TFBS scores for the DHS for a specific DHS. Two methods of association were used. In closest gene DHS, DHSs 1–4 from the GM12878 cell line are associated with the gene MAFB. For the TF in consideration, the maximum of all TFBS scores is 2.3. In Split DHS, we separated DHSs overlapping the TSS and other DHSs. This resulted in two features for each gene for each TF.

We used the TFBS scores as features for sparse logistic regression classifiers to discriminate between different gene classes. These classifiers balance the use of many available features against model complexity, effectively selecting a small subset of informative features that are used in the classification. We trained cell-type–specific classifiers on the task to discern whether a gene belonged to a specific expression pattern (e.g., UR vs. UR-Other, UR vs. DR, UR vs. Constitutive, etc.). The area under the receiver operating characteristic curve (AuROC) metric was used to evaluate the performance of a model, where a value of 0.5 indicates random assignments and 1.0 indicates perfect classification (see Methods). To not bias results due to different amounts of training data, the positive sets of up- and down-regulated genes were all of the same size.

The performance of the classifier using only proximal promoter information is close to that of a random classifier, across all tasks. All of the classifiers using DHS sequences display strong improvements in performance over this baseline in discriminating genes that are up-regulated in different cell types (UR vs. UR-Other) (Fig. 5A), with a greater improvement in performance coming from the Split DHS approach with separate features for the TSS and Distal DHSs (median AuROC ~0.73). Similar results were obtained when training classifiers to distinguish between specifically up- and down-regulated genes from the same cell types (UR vs. DR) (Fig. 5B), and to distinguish up-regulated from constitutively expressed genes (UR vs. Const.) (Fig. 5C). Discriminating down-regulated genes from different cell types (DR vs. DR-Other) and down-regulated from constitutively expressed genes (DR vs. Const.) resulted in lower accuracies but still showed the trend of better performance with DHSs compared with proximal promoter sequence (Supplemental Fig. 2A,B). All results clearly indicate that strong performance improvement is achieved by scanning for TFBS matches in open chromatin regions.

Open in a separate window

Figure 5.

Classifier performance for various classification tasks. (A–C) Performance of the classifier using all PWMs. Each figure compares the performance of two methods of associating DHSs to genes (Closest Gene DHS and Split DHS) with the proximal promoter. The solid black lines across the dots indicate the median. Across all figures, the promoter sequence classifier does not perform as well as the performance achieved by using Closest Gene DHS and Split DHS and is significant at the 0.05 level (paired t-test). (D–F) Impact of normalized CG dinucleotide content on classifier performance. Results using the Split DHS and promoter sequence are shown. Without CG, columns are the same as in A–C. All figures show average results from five iterations of fourfold cross-validation. The dotted line indicates an AuROC of 0.5, which is the performance of a random classifier.

Evaluating the influence of CG dinucleotide content

CG dinucleotide content in the proximal promoter sequence of genes is a common sequence feature that is directly or implicitly used to distinguish various classes of genes. Adding CG dinucleotide content as an additional feature led to a variable impact on classifier performance depending on the classification being considered (Fig. 5D–F). Specifically, when using open chromatin information, adding CG content did not substantially improve the performance of classifying UR genes from UR-Other genes (Fig. 5D). In the case of the Split DHS, only six of the cell lines had a significant coefficient for the CG dinucleotide coefficient (mean across all cell lines = −0.66, SD = 1.40; the coefficient was set to 0 when not significant). Due to the means being close to zero and the standard deviations being large, the effect of using CG content to discriminate between UR and UR-Other genes was largely negligible. Only for a few cell types, such as hepatocytes, did we observe a negative regression coefficient of significant magnitude (−5.11 for Split DHS and −3.31 for Promoters). This is in agreement with previous results showing that liver-specific genes have promoters with lower CpG content (Smith et al. 2005).

As anticipated by the trends observed in Figure 3E, UR vs. DR classification tasks benefited more by the addition of CpG content. Here, this feature was deemed to be significant in 17 of the cell lines for the Split DHSs (Fig. 5E). Furthermore, the regression coefficients were largely negative, which indicated the higher CG content among the DR genes (mean = −2.88, SD = 1.55). As has been shown before (Yamashita et al. 2005; Carninci et al. 2006; Zhu et al. 2008), we observed that high CpG content is predictive of constitutively expressed genes when compared with UR genes (Fig. 5F). The regression coefficient for the feature was significant in all cases (mean = −3.25, SD = 1.44). The CpG content feature had almost no impact in classifying DR from DR-Other genes (Supplemental Fig. 3). Finally, CpG content had a significant coefficient in classifying DR from constitutively expressed genes in only one cell line (mean = −0.09, SD = 0.39).

Adding CG content to the baseline proximal promoter models reconciled the apparent discrepancies between previous studies and the results reported in Figure 5A–C, because all classification tasks were improved upon for the proximal promoter. However, the DHS models with CG content outperformed baseline proximal promoter models with the inclusion of CG content (paired t-test < 0.05). In fact, in all cases except UR vs. Constitutive genes, DHS models even without CG content perform significantly better than both proximal promoter models (paired t-test < 0.05). Note that while adding CG content provided enormous performance gains for certain classification tasks (UR genes vs. Constitutive genes), this could be considered misleading. If TFBS scores are not explicitly normalized for local nucleotide composition, as we have done here, decent performance results can be achieved based solely on the different CG content observed for down-regulated and constitutively expressed genes compared with up-regulated genes. CG content is predictive in the case of classifying constitutive and DR genes from UR genes, but is not very useful in differentiating between genes that are up-regulated or down-regulated in different cell types. It is notable that the categories that are less aided by CG content are exactly those where our classifiers displayed the most predictive value.

Identifying candidate regulators

In addition to classifying genes belonging to different groups, we inspected the classifiers to identify motifs that were most informative in the classification task, i.e., those PWMs that had large regression coefficients (Supplemental Table 4). This identified several TFs with known impact on transcriptional output in the cell line of interest. For example, YY1, SPI1, and IRF8 are crucial in the specification of B-cells (GM12878 cell line) (Lu et al. 2003; Liu et al. 2007; Sokalski et al. 2011). We also identified the REST motif as a positive regulator of UR genes in the medulloblastoma cell line that is of neural origin (Supplemental Table 6). REST specifically down-regulates neuron-specific genes in many non-neuronal cell lines, and its expression is suppressed in neurons (Schoenherr and Anderson 1995). As a result, the model identified the cis-elements that are present in the DHSs associated with neuron-specific genes as the factor that separates these genes from the genes up-regulated elsewhere. This example illustrates that the inactivation of a repressor can also explain up-regulation of genes. Other well-characterized factors included ETS1 in HUVEC cells and HNF4A for HepG2 cells (Cereghini 1996; Oda et al. 1999; Yordy et al. 2005).

The feature set described thus far was comprehensive in that it used available PWM information from multiple sources, independent of the expression levels of transcription factors or the potential redundancy of features. To assess how much cell-type–specific regulation can be explained by the cell-type–specific expression of transcription factors themselves, we selected the top 10 TFs with highest absolute Z-scores from each cell line and had PWMs that were not similar to each other (Supplemental Table 7).

Using sparse logistic regression classifiers trained on these small sets of variables, we observed similar predictive trends, which indicated that a subset of cell-type–specific TFs were predictive of tissue-specific expression (Supplemental Fig. 4). Using only promoter CG content or the status of the chromatin at the TSS as features for a baseline comparison shows that motifs in DHS regions significantly contribute to the performance improvement across all comparisons. In addition, we used genomic regions identified as conserved in the 46-way placental mammal phastCons track from the UCSC Genome Browser. We note that using conserved sequences and particularly conserved noncoding regions improved performance compared with the promoter. However, the AuROC was still highest when DHS sequences were used, indicating that the presence of motifs in weakly conserved DHS regions contributes to the performance improvement. Finally, to assess the potential influence of insulators, we excluded DHSs that overlapped CTCF binding sites for classifiers trained specifically for the nine cell types for which genome-wide CTCF ChIP data were available (Supplemental Fig. 5). While this did not impact classification of UR genes, it reduced the accuracy of identifying DR genes, demonstrating that regions containing insulator sites are likely to contain regulatory information for the repression of genes.

Knowing both the regression coefficient in our model and the expression level of a potential regulator provided clues as to whether the TF in question is an activator or a repressor in the cell line, as highlighted for REST in medulloblastoma cells (Table 1; Supplemental Table 7). As another example, NR2F2 was identified as a positive predictor of up-regulated genes for embryonic stem cells. However, NR2F2 is a known negative regulator of POU5F1, a critical gene involved in pluripotency (Rosa and Brivanlou 2011). As expected, NR2F2 is down-regulated in ES cells (Supplemental Table 7). We also identified other known positive regulators, such as GATA1 in K562 cells (Huang et al. 2005) and MYF6 in myotubes (Fan et al. 2011). Note that genes that have both positive and negative coefficients have different effects when in TSS DHSs and Distal DHSs.

Table 1.

Candidate TFs identified by the classifier for each cell line using Split DHS from the top 10 highest absolute Z-score of expression and nonredundant TFs

Open in a separate window

For HNF4A in HepG2 and GATA1 in K562 cells, ChIP data are available from the ENCODE Project. To validate the predictions made by our model, we looked for overlap of these ChIP sites with DHS sites associated with different sets of genes. In HepG2 cells, 19% of all genes with an associated DHS overlapped an HNF4A binding site. Strikingly, 64.5% of the UR genes had a DHS overlapping an HNF4A ChIP peak (P-value < 1 × 10⁻¹², binomial test). Conversely, only 10.5% of DR genes had a DHS that overlapped an HNF4A site (P-value < 1 × 10⁻³). In K562 cells, we found that 6% of all genes had an associated DHS with a GATA1 ChIP peak. However, 31.5% (P-value < 1 × 10⁻¹²) of UR genes and only 3.5% (P-value < 0.1) of DR genes had a DHS with a GATA1 ChIP peak. The ChIP binding data provided strong and independent evidence that our models identify relevant factors that regulate the transcriptional program in these cells.

In addition, we investigated the accuracy of our predictions of TFBS locations in DHSs. In HepG2, 5215 of the 6597 HNF4A ChIP peaks overlapped the predicted TFBS in DHSs. Furthermore, using TFBS scores led to a high accuracy on discriminating between positive and negative sets defined by ChIP peaks (AuROC of 0.79). In K562 cells, only 315 of the 1704 GATA1 ChIP peaks overlapped the predicted TFBS in the DHS; yet, the AuROC still remained high at a value of 0.88. This indicated that high-scoring TFBSs accurately predicted binding of GATA1 to these sites. We note that the low percent overlap may arise from nonspecific or indirect binding of GATA1.

To assess the presence of additional sequence motifs not accounted for by the sets of known PWMs, we used the discriminative version of MEME to perform motif finding (Bailey et al. 2010), identifying motifs differentially enriched between UR and UR-Other, respectively, DR genes (Supplemental Table 8). While some of the identified motifs corroborated the importance of features from the set of top 10 TFs (FOXA2 [formerly HNF3B] in HepG2), others corresponded to TFs that were not in this list. These are candidate TFs that are not among the most differentially expressed, but still might be involved in the transcriptional program, potentially through other steps of activation. We note that we largely did not recover the motifs recently identified in a subset of seven of the 19 cell lines (Song et al. 2011). In contrast to this study, which used the sequences from cell-type–specific DHSs as foreground and the subsets of cell-type–specific DHSs in other cell types as background, we analyzed the sequences from all DHSs associated with a gene and defined the background according to the classification tasks.

Classifying DHSs based on genomic location

The RefSeq hg19 database was downloaded from the UCSC Genome Browser and used to classify DHSs based on their genomic location. If a DHS overlapped the TSS of any transcript variant of a gene, it was classified as being a TSS DHS for that gene. Other DHSs were similarly classified as Gene Body DHSs if they overlapped any region of the gene excluding the TSS. All other DHSs were classified as Intergenic DHSs.

Microarrays

We used Affymetrix Human Exon 1.0ST microarrays to measure gene expression following ENCODE protocols. We normalized 110 microarrays (measuring 40 cell lines) together, then extracted the subset (19 cell lines) used in the present study. Probesets flagged as cross-hybridizing were first removed from the analysis (Salomonis et al. 2010). Although these arrays provide exon-level probesets, we sought gene-level expression estimates, so we grouped probesets by gene for normalization (Bemmo et al. 2008). To normalize, we used Affymetrix Power Tools (APT) with the chipstream command “rma-bg, med-norm, pm-gcbg, med-polish.” This chipstream calls for an RMA normalization with gc-background correction using antigenomic background probes. After normalizing, we noticed an effect due to a switch in microarray reagents. Some of our arrays were processed differently, because our earlier array reagents become unavailable partway through the experiment. Using hierarchical clustering and multidimensional scaling, we found the arrays to group on the basis of reagent used, rather than by biological relatedness. To make the arrays comparable, we used an R script (ComBat) to correct for this batch effect (Johnson et al. 2007). After correction with ComBat, the arrays grouped according to expected biological similarity.

Cell-type–specific expression

We identified the genomic location of genes based on matching gene symbols to the RefSeq hg19 database. If a gene from the array did not have a matching gene symbol, it was dropped from the analysis. If a gene did not have expression above background (>4.8) in at least one cell type, it was dropped from the analysis. For the remaining genes, expression values across the 19 cell lines were Z-transformed. The Z-scores for expression in each cell line were sorted. The top 200 genes were classified as UR genes and the bottom 200 genes as DR genes in that cell type. For each cell type, the UR-Other genes were compiled as follows: We first made a set comprising all UR genes from all cell types. We then removed the current cell-type UR genes from this global UR gene set. We further removed genes from this set that had expression Z-score ≥ 0 in the current cell type to exclude genes that had higher than mean expression in that cell type. This ensured that this set of genes was purely composed of genes that were up-regulated in other cell types and had lower than mean expression in the current cell type. A similar procedure identified DR-Other genes. To identify constitutively expressed genes, we selected genes in neither UR or DR sets across all cell lines that were expressed above background in all cell lines and had a maximum |Z-score| < 1.7, which resulted in 168 genes. This size compared well with the other positive sets of UR and DR genes. This gave us a list of genes that did not have a significant variance in their expression.

GO analysis

DAVID was used to identify functional categories of genes up-regulated in each cell type (Huang et al. 2009a,b).We used the GO categories and also the UP_Tissue category to identify the tissue showing the closest gene expression profile to the cell type in question. A P-value < 0.05 was used as the significance threshold.

PWM scans of DHS and promoter sequences

We collected PWMs for vertebrate TFs from the TRANSFAC, JASPAR, and UniProbe databases. This gave us a collection of 789 PWMs, some of which refer to the same TF. Note that we allowed multiple PWMs for each TF in our data set since it is generally not known which one reflects binding affinities more accurately. Furthermore, multiple PWMs may also reflect different binding modalities for the same factor. This could be due to, for example, the presence of specific cooperative binding partners in one cell type but not in another.

We scanned the sequence from each DHS site and (proximal) promoter sequence (−900 to +100 relative to each TSS of a gene) using these PWMs. A score was assigned to each location in the sequence based on the log-likelihood ratio of the PWM score (probability of the PWM generating the specific sequence) versus the probability that the sequence was generated by a background model. The background model used was a first-order Markov Model trained over a 500-bp window centered at the base pair being scored. This effectively corrects for the underlying dinucleotide composition and allows us to separate signal from noise (Megraw et al. 2009). Scores with a log-likelihood ratio less than 0 were not included in further analysis.

After these scores were generated for each base pair, we summed scores across a sliding window of size 60 to account for local clusters of multiple, potentially overlapping binding sites. Clusters of binding sites have been shown to be more likely to be bound by TFs as opposed to single binding sites (Gotea et al. 2010). For each DHS or promoter region, we assigned the maximum sliding window score as the TFBS score for that TF for that region.

Associating DHS TFBS scores with genes

To associate each DHS with a gene that it was potentially regulating, we found the closest TSS to the midpoint of the DHS. The DHS was then associated with that gene. In general, there were several associated DHSs for each gene, and these were assumed to be the putative regulatory regions for that gene. To assign one score for a TFBS to each gene, we picked the maximum TFBS score across all of the associated DHSs (closest gene DHS). In Split DHS, we separated DHSs into two groups—overlapping TSSs and those in other parts of the genome. We selected the maximum for each set, therefore having two features per gene per TF.

Sparse logistic regression classifier

We used a sparse logistic regression classifier that minimizes an objective function that is a linear combination of the sum of squared residuals and the ₁-norm of the weights (Koh et al. 2007). We divided our data into four parts to perform a fourfold cross-validation. Three parts of the divided data were used as the training set. This training set was further divided into six parts, one of which was used as the validation set to learn the hyperparameter for the contribution of the ₁-norm in the objective function. The optimal hyperparameter was selected from 10 values (0.001, 0.011, …, 0.091). This value was then used to evaluate the performance of the model on the original test set, and the area under the receiver operating characteristics (AuROC) was computed as a measure of performance. We performed five iterations of each fourfold cross-validation, with the data shuffled before each iteration. Results for each classification task are therefore averages over 20 different partitions of the data, which makes our results more robust to chance arrangements of the data.

Cell-type–specific expressed nonredundant TFs

To identify TFs that were cell-type–specifically expressed, we used the absolute Z-score and extracted gene symbols for TFs with available PWMs. We again used the gene symbol from the Affymetrix arrays to match expression to TF names in our PWM list. By using absolute Z-scores, we picked out genes that were cell-type–specifically up- and down-regulated. To ensure that we were picking a nonredundant set of TFs, we used STAMP (Mahony and Benos 2007) with the default settings. Starting from the TF with the highest absolute Z-score, we only added a TF when it had a significantly different PWM (E-value > 0.25) from the TFs already chosen. We stopped adding to the set when we had 10 TFs.

Conservation analysis

Genomic regions from the 46-way placental mammal track were downloaded from the UCSC Genome Browser. Only regions of size 100 bp or more were used. Coding exon sequences were extracted from the RefSeq hg19 database. BedTools (Quinlan and Hall 2010) was used to subtract exonic coding sequences from conserved regions to find noncoding conserved regions. Regions were scanned for motifs, and TFBS scores were assigned to genes in the same way as open chromatin regions.

Discriminative motif finding using MEME

We used MEME to first calculate a position-specific prior using the psp-gen tool. For example, if we wanted to identify the motifs in the DHS sequences related to the UR genes when classifying against DR genes, the positive file was the UR genes and the negative file was the DR genes. This was then used to identify motifs from 1000 positive sequences, UR gene DHS sequences in the above example. The motif width was chosen to be between 8 and 12 bp, and motifs were allowed to have zero or one occurrence per sequence. Motifs with an E-value < 0.05 were then compared with either the top 10 nonredundantly expressed TFs or the full list of TFs using STAMP. A STAMP E-value < 0.05 was accepted as a good match of a motif to a TF.

ChIP analysis

ChIP-seq peak coordinates were obtained from the ENCODE webpage in BED (narrowPeaks) file format. To assess whether DHSs are really bound by the TF or not, we checked for overlap between the coordinates of DHS and ChIP-seq peaks for that TF. To calculate AuROC, the TFBS scores in the DHS were compared against ChIP-seq peaks. Overlaps with the 60-bp window around the TFBS site were considered true positives, and others were considered false positives for AuROC calculations.

We thank Terry Furey and Darin London for initial processing of DNase-seq data, Molly Megraw and Ashlee Benjamin for providing code for feature generation and classification, and Danielle Maatouk for help with conservation analysis. We also thank Alexander Hartemink for helpful discussions. This work is supported by NHGRI grant U54 HG004563 to G.E.C. and NHGRI award R01 HG004065 to U.O. N.C.S. is supported by an NSF Graduate Research Fellowship and the Norwegian Research Council.