A Clustering Approach for Motif Discovery in ChIP-Seq Dataset.

Chromatin immunoprecipitation combined with next-generation sequencing (ChIP-Seq) technology has enabled the identification of transcription factor binding sites (TFBSs) on a genome-wide scale. To effectively and efficiently discover TFBSs in the thousand or more DNA sequences generated by a ChIP-Seq data set, we propose a new algorithm named AP-ChIP. First, we set two thresholds based on probabilistic analysis to construct and further filter the cluster subsets. Then, we use Affinity Propagation (AP) clustering on the candidate cluster subsets to find the potential motifs. Experimental results on simulated data show that the AP-ChIP algorithm is able to make an almost accurate prediction of TFBSs in a reasonable time. Also, the validity of the AP-ChIP algorithm is tested on a real ChIP-Seq data set.

1. Introduction

Transcription factor binding sites (TFBSs) [1] are short, degenerated nucleotide fragments (usually ≤30 bps) located in specific DNA regions. They play an important role in regulating gene expression. The Planted Motif Search (PMS) problem [2] is a popular motif model for the identification of TFBSs (i.e., motif discovery) in bioinformatics, and is formally defined as follows:

Given a set of DNA sequences

According to the different values of q representing the distribution of motif instances, there are three different motif discovery sequence models [3]: (i) Exactly one motif occurrence per sequence (the OOPS model), (ii) zero or one motif occurrences per sequence (the ZOOPS model), or (iii) zero or more motif occurrences per sequence (the TCM model). For the OOPS model, , for the ZOOPS model and TCM model, .

Generally, there are two kinds of algorithms for solving the PMS problem: exact algorithms and approximate algorithms. Exact algorithms [4,5,6,7,8,9,10] always use consensus sequences [11] to represent motifs and can find all motifs. Most exact algorithms are pattern-driven algorithms, which attempt to enumerate all possible l-mers (substring patterns of length l) to find the l-mer with the maximum number of approximate occurrences. Approximate algorithms [12,13,14,15,16,17,18] usually use a position weight matrix (PWM) [19] to describe the most likely occurring motifs and can report results in a short time, but do not always identify all motifs. Most approximate algorithms use probabilistic analysis to maximize the score function which describes how likely it is for an l-mer pattern to be a motif instance.

Recently, chromatin immunoprecipitation combined with next-generation sequencing (ChIP-Seq) technology has produced extremely valuable information for the genome-wide identification of transcription factor binding sites (TFBSs) and in the field of epigenetics, which mainly focus on DNA methylation, nucleosome localization, and histone modification. For transcription factors, ChIP-Seq is widely used to study the binding of transcription factors for the analysis of gene expression regulation on a genome-wide scale. For histones, ChIP-Seq performs high-throughput histone modification sequencing in the whole genome with sufficient sequencing depth and range, which not only improves the sensitivity and specificity of sequencing, but can also transform qualitative sequencing methods into quantitative detection.

In this paper, we focus on an algorithm to discover transcription factor binding sites in a ChIP-Seq data set. A ChIP-Seq data set is a set of peak regions containing TFBSs obtained through ChIP-Seq experiments, read mapping, and peak calling. It contains hundreds or more DNA sequences, which increases the difficulty of accurate and efficient identification of TFBSs.

Some algorithms have recently been proposed to discover TFBSs in ChIP-Seq data sets [20,21,22,23,24,25,26,27,28,29]. However, none of them has proven to be absolutely superior, compared to the rest. Some of these are tailored versions of previous motif discovery algorithms, specifically tailored towards ChIP-Seq data sets, such as MEME-ChIP [20] and HMS [21]. MEME-ChIP [20], which incorporates two complementary motif discovery algorithms, known as MEME and DREME [22], can identify motifs without restriction on the size or number of sequences, allowing very large ChIP-Seq data sets to be analyzed. HMS [21], which is an improved version of Gibbs Sampler, combines stochastic sampling and a deterministic greedy search step, which improves computation efficiency. DREME [22] is specifically designed to find short, core DNA-binding motifs of eukaryotic TFs, and is optimized to analyze very large ChIP-Seq data sets in just minutes. One may speed up the existing motif discovery algorithms by integrating some information, such as in the cases of STEME [23] and ChIP-Munk [24]. STEME [23] accelerates MEME by indexing sequences with suffix trees. ChIP-Munk [24] combines a greedy approach with an expectation-maximization (EM) algorithm to achieve a high efficiency. There are also exhaustive methods for determining exact motifs in ChIP-Seq data sets, such as FMotif [25] and Weeder [26]. FMotif [25] first constructs a mismatched suffix tree to scan and count all possible motif instances, and then implements a depth-first search to enumerate all possible motifs. However, the run time of FMotif becomes unrealistic with increasing values of l and d. Others use word enumeration methods to process full-size ChIP-Seq data sets, such as CisFinder [28] and MCES [29]. CisFinder [28] employs a word clustering method to group short l-mers ( 8, or 9), but struggles to find exact motifs with larger values of l and d in ChIP-Seq data sets. MCES [29], a new planted motif discovery algorithm, mines and combines substrings to rapidly identify exact motifs in full-size ChIP-Seq data sets.

In this paper, we propose a new motif discovery algorithm, named AP-ChIP, which is specially designed for better discovering TFBSs in ChIP-Seq data sets. The algorithm first constructs and then further filters cluster subsets using probabilistic analysis. Then, Affinity Propagation (AP) clustering [30] is applied to the candidate cluster subsets in order to discover optimal motifs. Experimental results show that the AP-ChIP Algorithm 1 can find TFBSs in a ChIP-Seq data set very efficiently and effectively.

2. Method

A ChIP-Seq data set has the following fundamental features: (i) Some of the sequences may contain no motifs at all; and (ii) thousands of sequences lead to huge amount of background l-mers. To cater to ChIP-Seq data sets, we design the AP-ChIP Algorithm 1 under the ZOOPS model and set some proper thresholds to filter redundant background l-mers. More specifically, the AP-ChIP Algorithm 1 consists of the following three steps:

2.1. Construct Cluster Subsets

We introduce the observation that any two motif instances and , each of which differs from the same motif x up to d positions, must have Hamming distance of no more than , denoted as . Consequently, if an l-mer in one sequence is a motif instance, all other motif instances in the remaining sequences will be gathered in the corresponding cluster subset. Under the ZOOPS model, we choose sequences as the reference sequences to ensure that at least one sequence among the h sequences contains a motif instance [31]. In general, we use the first h sequences as the reference sequences. As the l-mer which is the motif instance is not known in advance, we consider all l-mers in the first h sequences as the reference subsequences.

The ideal cluster subsets are expected to contain as few background l-mers as possible and, also, to include sufficient motif instances. Therefore, we set a threshold k, so that, for each reference subsequence , all l-mers , in the whole sequences, except the ith sequence such that , are selected to construct a cluster subset; that is where represents the set of the selected l-mers in the th sequence and if and only if is an mer of the sequence .

To set a proper threshold k, two probabilistic expressions are employed. The first is the probability of the Hamming distance between two random l-mers and being no more than k [4]:

The other is the probability of the Hamming distance between two selected motif instances and being no more than k: [18].

Now, we describe the method for calculating . Given two motif instances, and , of the same motif with up to d mutations, the distances between , and satisfy and , . Thus, the probability can be calculated as where represents the conditional probability of given , , such that and represents the probability of and ; that is

As and are independent, we have

Combining Equations (5)–(9), we have

Having calculated the two probabilities and , we now describe how to set the proper threshold k, in order to construct the cluster subsets which contain as few background l-mers as possible while including sufficient motif instances. A larger value of indicates more motif instances belong to the cluster subset; however, a smaller value of suggests that fewer background l-mers appear in the same cluster subset. Therefore, the threshold k should be set in a way that ensures that the value of is large enough, compared to the value of .

To demonstrate this issue, let us consider the (18, 5) problem instance as an example. The values of and are shown in Table 1. When , the value of is 0.7414, which allows us to obtain sufficient motif instances, whereas the value of is 0.0012 which, in turn, allows us to reduce the scale of background l-mers in the same cluster subset. Therefore, the optimal value of k is 7.

Table 1

Values of and under different values of k for (18, 5) problem instance.

k	5	6	7	8	9	10
p_k	3.42×10^−5	2.31×10^−4	0.0012	0.0054	0.0193	0.0569
p_kmotif	0.2303	0.4741	0.7414	0.9242	0.9915	1

2.2. Filter Cluster Subsets

As is known, the true motif instances must exist in one of these cluster subsets , which are constructed with the reference subsequences from the first h sequences. However, with a great number of total cluster subsets, the identification of the cluster subsets that contain the true motif instances is highly time-consuming as most of these cluster subsets have redundant background l-mers.

To filter the interference cluster subsets, a threshold (i.e., an occurrence frequency) [29] is employed with the purpose of analyzing the probability of a random motif instance occurring in a given sequence. where is the probability of a character mutation in one position. For each sequence , if the number of the selected l-mers in , denoted as , is greater than or equal to , then the sequence may contain a motif instance and is stored in a set . For each , if the number of the selected sequences in , denoted as , is not less than q, then there are at least q possible motif instances in the corresponding cluster subset , the reference subsequence is considered as a potential motif instance, and the corresponding cluster subset is a candidate cluster subset .

2.3. Refine Cluster Subsets

Due to the occurrence frequency and the relatively small Hamming distance k between the reference subsequence and the selected l-mers, only a limited number of cluster subsets, which contains a small number of the selected l-mers, are retained as candidate cluster subsets for further Affinity Propagation (AP) clustering. In order to quickly produce highly conserved cluster subsets, we apply AP clustering to each candidate cluster subset. Compared to other clustering approaches, AP clustering can cluster large-scale data sets efficiently by exchanging messages between data points.

2.3.1. Affinity Propagation (AP) Clustering

As demonstrated in our recent work [32], it is possible to speed up AP clustering and improve its accuracy with the adapted similarity , which is based on pair-wise constraints and a variable-similarity measure [33]. The adapted similarity, , between two l-mers and is defined as where

Note that and .

For each candidate cluster subset, based on the adapted similarity , AP clustering recursively calculates two types of messages. The first type is the responsibility , which reflects the suitability of point as the exemplar for point . The other type is the availability , which indicates how suitable it would be for a point to choose the point as its exemplar:

When the AP clustering converges, a set of l-mers in the produced cluster subset are selected as exemplars associated to the point :

2.3.2. Cluster Subset Refinement

To select an adequate number of desired cluster subsets with more motif instances but less background l-mers, we set an interval [min size, max size] = to further refine the cluster subsets , which is produced, by AP clustering, using a reference l-mer . Regarding the number of the l-mers in each cluster subset , we conclude that there are three cases:

(i) In the case where the number is less than (for example, , such a small number of l-mers is not enough to create a cluster subset which represents a true motif, so we consider it as an invalid cluster subset.

(ii) In the case where the number is between and t (for example , we consider it as a valid cluster subset.

(iii) In the case where the number is more than t (for example ), the size of the cluster subset is so large that it may include too many background l-mers and, so, we use a greedy strategy to select tl-mers from to form . First, is initialized with the AP clustering exemplar . Then, an l-mer from is repeatedly chosen, following Equations (21) and (22), and added to until : where is the length of the maximum intersection of and .

To appropriately sort the valid cluster subsets , we use the information content (IC) [34], and the cluster subset with the maximum value of IC is considered as the true motif model: where is the probability of the character at the position m of the l-mer , and is the corresponding background probability.

Based on the above described three steps, the whole AP-ChIP Algorithm 1 is described as follows:

Steps 2–17 describe the process of constructing the cluster subsets. Steps 18–21 describe the filtration of the cluster subsets. Steps 22–31 describe the refinement of the cluster subsets and the verification of the motif with the maximum IC score.

3. Results

3.1. Results on Simulated Data

Simulated data provide quantitative measures to test the performance of the AP-ChIP Algorithm 1. As in [29], we generate the simulated data as follows:

First, we generated t independent and identically distributed (i.i.d) sequences of length n and a motif m of length l. Second, we randomly generated q motif instances, each of which differed from the motif m at up to d positions. Third, the q motif instances were placed in a random position in a random selection of q sequences selected out of the t sequences. We, then, implemented the AP-ChIP Algorithm 1, using Matlab on a computer with a 2.6 GHZ processor and 4 Gbyte memory. The final experimental results consisted of the averages of five trials of simulated data experiments.

To evaluate the motif prediction accuracy, the nucleotide level performance coefficient , as defined by Pevezner and Sze [2], was used: where K is the set of nucleotide positions in the true motif and P is the corresponding set of nucleotide positions in the predicted motif. The value of is between 0 and 1; the larger the value of , the higher the accuracy of the predicted motif. We used the neighborhood probability [8] to select a group of motif instances. The larger the value of , the weaker the corresponding problem instance becomes:

In what follows, according to different values of (i.e., the ratio of the sequences containing motif instances to all the sequences), we designed two groups of experiments, both of which consisted of two sub-experiments, to test the performance of the AP-ChIP Algorithm 1 on simulated data sets.

In the first group of experiments, we set for both sub-experiments. We compared the performance of the AP-ChIP Algorithm 1 with that of the widely used motif finding algorithms MEME [13], VINE [14], and Projection [16].

In the first sub-experiment, we set the number of sequences as with sequence length . The running time and the predicted accuracy of these algorithms are shown in Table 2. For the instances (12, 2) and (15, 3) with , the AP-ChIP Algorithm 1 achieved near-optimal predicted accuracy within a relatively short time. For the instances (15, 4), (14, 4), (25, 8), and (21, 7) with , the predicted accuracy of the AP-ChIP Algorithm 1 was over and the running time of the AP-ChIP Algorithm 1 remained competitive, compared to the other algorithms.

Table 2

Comparisons on problem instances with , , and .

(l,d)	p_2d	MEME	VINE	Projection	AP-ChIP
(12, 2)	0.0028	0.68 (4 s)	1.00 (8 s)	0.86 (10 s)	0.98 (18 s)
(15, 3)	0.0042	0.73 (7 s)	1.00 (9 s)	0.82 (1.3 m)	1.00 (23 s)
(15, 4)	0.0566	0.87 (8 s)	0.96 (5.6 m)	0.89 (4.2 m)	0.97 (36 s)
(14, 4)	0.1117	0.84 (10 s)	0.95 (8.3 m)	0.80 (27.4 m)	0.96 (47 s)
(25, 8)	0.1494	0.91 (12 s)	0.93 (9.8 m)	0.78 (32.6 m)	0.94 (1.1 m)
(21, 7)	0.2564	0.87 (28 s)	0.92 (11.2 m)	0.76 (48.7 m)	0.91 (58 s)

In general, it is easy to find the true motif by increasing the sequence number and decreasing its length. Therefore, in the second sub-experiment, we set the sequence number as with sequence length . Table 3 shows that, for all problem instances with , the predicted accuracy of the AP-ChIP Algorithm 1 is over with the computational costs being satisfactory.

Table 3

Comparisons on problem instances with , , .

(l,d)	p_2d	MEME	VINE	Projection	AP-ChIP
(12, 3)	0.0540	0.94 (8 s)	0.96 (2.4 m)	0.91 (3.1 m)	0.97 (21 s)
(11, 3)	0.1146	0.86 (10 s)	0.95 (5.2 m)	0.84 (4.3 m)	0.96 (33 s)
(13, 4)	0.2060	0.83 (10 s)	0.93 (8.1 m)	0.78 (36.4 m)	0.95 (36 s)
(15, 5)	0.3135	0.78 (11 s)	0.84 (9.6 m)	0.74 (46.7 m)	0.93 (34 s)
(17, 6)	0.4261	0.70 (13 s)	0.83 (18.6 m)	0.72 (53.6 m)	0.92 (38 s)
(19, 7)	0.5346	0.68 (17 s)	0.75 (24.5 m)	0.70 (1.2 h)	0.90 (40 s)

Next, in the second group of experiments, to simulate a real ChIP-Seq data set, we set for both sub-experiments. This is because in real ChIP-Seq data set, most but not all of the sequences contain motif instances.

In the first sub-experiment, we test the validity of the AP-ChIP Algorithm 1 on the simulated ChIP-Seq data set for the identification of motifs with and . We choose to select a group of problem instances. The reason for this choice is that is approximately the same as the value of the (15, 4) problem instance, which is one of the most popular benchmarks for problem instance. The running time and the predicted motif by the AP-ChIP Algorithm 1 are shown in Table 4. For each problem instance, the AP-ChIP Algorithm 1 finds almost the same motif as the published one and also runs quite efficiently.

Table 4

The results on problem instances with , .

(l,d)	Time	Predicted Motif	Published Motif
(9, 2)	43 s	TTATCCCTC	TTATCCCTC
(12, 3)	34 s	TTTCCCGTCTGC	CTTTCCCGTCTG
(15, 4)	42 s	GGTTGRAGCTTAGGG	GGTTGGAGCTTAGGG
(18, 5)	38 s	CTTTGCCATATCCATAGG	TTTGCCATATCCATAGGC
(21, 6)	36 s	CAGGTAAACCATATTAAATTA	AGGTAAACCATATTAAATTAC

R: A,G.

In the second sub-experiment, in order to further demonstrate the performance of the AP-ChIP Algorithm 1 on the simulated ChIP-Seq data set, we compared the AP-ChIP Algorithm 1 against some established motif-finding algorithms in the following two aspects: (i) Different values of , ranging from 0.05 to 0.5, with fixed , and (ii) different values of floating from 0.7 to 1 with fixed . Although a genome-wide ChIP-Seq data set typically has thousands to tens of thousands of sequences, using to of the ChIP-Seq data set is usually adequate for obtaining a good estimate of the true motifs. MEME-ChIP, a well-known algorithm for discovering motifs in ChIP-Seq data sets, is able to well identify motifs with only 600 sequences. Thus, it was reasonable to set the sequence number as and sequence length as for motif discovery in our experiments.

First, we compared the running time and prediction accuracy of the AP-ChIP Algorithm 1 with those of the three compared algorithms, MEME-ChIP [20], ChIP-Munk [24], and FMotif [25], on different values of and with fixed . As shown in Table 5, for each problem instance, the AP-ChIP Algorithm 1 could solve it in a relatively short time, and its prediction accuracy was better than those of the three compared algorithms. Specifically, with increasing values of l and d, FMotif found it difficult to find exact motifs.

Table 5

Comparison of problem instances with , , and .

(l,d)	p_2d	MEME-ChIP	ChIP-Munk	FMotif	AP-ChIP
(9, 2)	0.049	0.96 (12 s)	0.96 (1.8 m)	1.00 (47 s)	1.00 (43 s)
(11, 3)	0.114	0.94 (24 s)	0.92 (2.0 m)	0.99 (7.9 m)	0.98 (46 s)
(13, 4)	0.205	0.90 (38 s)	0.83 (2.4 m)	0.98 (1.45 h)	0.93 (58 s)
(15, 5)	0.319	0.85 (42 s)	0.80 (8.2 m)	–	0.92 (1.1 m)
(17, 6)	0.426	0.80 (45 s)	0.78 (9.6 m)	–	0.89 (1.3 m)
(19, 7)	0.534	0.78 (48 s)	0.76 (10.7 m)	–	0.87 (1.6 m)

Next, we compared the prediction accuracy of AP-ChIP Algorithm 1 with those of the algorithms MEME [13], MEME-ChIP [20], DREME [22], and FMotif [25], on different values of floating from 0.7 to 1 and fixed . As the ratio of a motif instance has a strong effect on the prediction accuracy, it was necessary for us to test the prediction accuracy of AP-ChIP Algorithm 1 on different values of . It is rather cumbersome to identify the true motif when the value of is small. FMotif is a powerful, exhaustive algorithm for finding exact short motifs contained in ChIP-Seq data sets. As shown in Figure 1, the prediction accuracy of AP-ChIP Algorithm 1 was nearly the same as that of FMotif, and was higher than that of MEME-ChIP, MEME, and DREME.

Figure 1

Prediction accuracy for different values of .

3.2. Results on Real Data

First, we tested the validity of the AP-ChIP Algorithm 1 for the identification of real motifs using a diverse set of real ChIP-Seq data sets; specifically, on 12 differently sized mESC data sets (Nanog, Oct4, Sox2, Esrrb, Zfx, Klf4, c-Myc, n-Myc, Tcfcp21l, Smad1, STAT3, and CTCF) [35]. We compared the motifs detected by the AP-ChIP Algorithm 1 with the motifs published by Chen et al. [35], and presented them in the form of sequence logos [36], which graphically represent the degree of motif conservation, as measured by relative entropy. Table 6 shows the running times and the predicted motifs of the AP-ChIP Algorithm 1. For each data set, the AP-ChIP Algorithm 1 was capable of finding motifs highly similar to the published ones within a reasonable time.

Table 6

Results on the mESC data set.

Data Set (Seq #)	Time	Predicted Motif	Published Motif
c-Myc (3422)	125 s
CTCF (39609)	19 s
Esrrb (21647)	10 s
Klf4 (10875)	138 s
Nanog (10343)	12 s
n-Myc (7182)	36 s
Oct4 (3761)	48 s
Smad1 (1126)	12 s
Sox2 (4525)	15 s
STAT3 (2546)	27 s
Tcfcp211 (26910)	11 s
Zfx (10338)	68 s

Moreover, to better show the results, we compared AP-ChIP 1 with MEME-ChIP on nine differently sized ENCODE data sets (Nfyb, Hnf4, Elf1, Ets, Egr1, Yy1, Six5, Srf, and Tal1) [37], where the TFBSs were referenced in the JASPAR database [38]. Table 7 shows the published motifs and the motifs predicted by the two algorithms. The motifs are also shown in the form of sequence logos. AP-ChIP 1 could successfully find a motif similar to the published motif for each data set, while, for some data sets MEME-ChIP failed to accurately predict the motif (e.g., in the Elf1 data set), or lost information on individual bases (e.g., in the Tal dataset).

Table 7

Results on the ENCODE dataset.

Data Set (Seq #)	AP-ChIP Predicted Motif	MEME-ChIP Predicted Motif	Published Motif
Nfyb (10096)
Hnf4 (11045)
Elf1 (8611)
Ets  (5525)
Egr1  (15400)
Yy1  (2077)
Six5  (4664)
Srf  (4903)
Tal1  (25507)

4. Concluding Remarks

In this paper, our goal was to find a method providing balance between time performance and prediction accuracy for TFBS discovery in ChIP-Seq data sets. Consequently, we aimed to obtain these results with high prediction accuracy in a relatively short time. To do so, we proposed a novel clustering-based algorithm named AP-ChIP 1. Firstly, to achieve high prediction accuracy, we set a threshold k to restrict the number of the selected l-mers in the candidate cluster subsets. Next, to obtain good time performance, we set the threshold in terms of probabilistic analysis, in order to filter the interferential candidate cluster subsets. Furthermore, a powerful data clustering method, AP clustering, was used to obtain the almost accurate motifs. Experimental results on both simulated and real ChIP-Seq datasets showed that the AP-ChIP Algorithm 1 not only discovers the motifs as consistently as the published ones, but also does so quite efficiently. This demonstrates that the AP-ChIP Algorithm 1 is a powerful new approach for ChIP-Seq data set analysis which provides a good trade-off between time performance and prediction accuracy.