SOMEA: self-organizing map based extraction algorithm for DNA motif identification with heterogeneous model

Performance on real datasets

The eight test datasets used in this experiments are composed of seven datasets used in [7] and a dataset collected from the Promoter Database of S. cerevisiae (SCPD) [18]. Each sequence contains at least one true binding site. These datasets consist of motifs from Escherichia coli(CRP), homo sapiens(ERE, MEF2, SRF, CREB, E2F, MYOD) and S. cerevisiae(GCN4).

SOMEA was run with map sizes that were arbitrarily selected between 10 × 10 to 20 × 20 depending on the size of the dataset. In each case, SOMEA was trained for 100 epochs with a motif length value in [l – 3, l + 3], where l is the known motif consensus length. The top 10 highest ranked motifs according to their MAP score [19] were saved for evaluation purpose. A 3rd order markov chain model [15] was used to compute MAP score. The learning rate parameter was fixed at 0.005 in all the experiments. Whereas, the neighborhood function parameter value, σ was set at 3.0.

For WEEDER, we used the online tool [20] with the following options: sites might appear more than once, both strands, and normal or complete scan. The interesting motifs and their instances that scored at least 90 were used in the evaluation. SOMBRERO was run with the default map sizes and random initialization method. The standalone tool was downloaded from [21]. We evaluated all the “best-motifs” returned by the tools. MEME was run with the “any number” model option and minimum and maximum length value as discussed above. AlignACE was run online [22] with default arguments in most cases.

It can be noticed that, SOMEA performance is comparable or better than ALIGNACE, MEME and WEEDER. For example, in terms of F-measure rates, SOMEA produces the best results for five of the eight datasets due to its higher precision rates (note that, both SOMEA and ALIGNACE achieve the same F-measure value for the CRP dataset). SOMEA’s average F-measure value for all datasets (i.e. 0.72) is found better than MEME (0.65), ALIGNACE(0.69) and SOMBRERO(0.55) and equally good as WEEDER(0.72).

It should be noted that, the comparison results between programs cannot be completely fair as every program has its own strengths and weaknesses. For example, some programs might perform rather well for strong motifs; whereas some are designed to discover motifs with certain characteristics (e.g. gapped motifs). The nature of the datasets can be an influential factor to the success of each program. Therefore, the results reported here should serve only as reference.

Performance on artificial datasets with multiple planted motifs

Practically, we can often find multiple motifs in upstream region of a set of co-regulated genes. These motifs often work as cis-regulatory module to regulate gene expressions. Motif discovery programs should be able to return all of these potential motifs. Local search algorithms, such as MEME, perform a search for single motif at one time; whereas SOMEA and SOMBRERO search for all motifs simultaneously. It is interesting to compare these two strategies.

We have prepared five artificial DNA datasets generated from Annotated regulatory Binding Sites (ABS, v1.0) database [23]. Every DNA dataset has twenty(20) sequences (each with 500bp in length) with three planted real motifs. We run MEME, WEEDER, SOMEA and SOMBRERO five times on each dataset. We asked SOMEA and MEME to return the top 20 motifs for the evaluation purposes. Again, we evaluate all best motifs returned by WEEDER and SOMBRERO.

It should be noted, there are some biases in the comparisons for two reasons. Firstly, both SOMEA and SOMBRERO are rather sensitive to the motif length parameter. As the motif consensuses in a dataset have different lengths, a single run with a fixed length value might not be suited for all motifs. On the contrary, MEME is able to find a length value that suits better each motif. Consequently, in some of SOMEA/SOMBRERO runs, some motifs might appear to be performed better in the experiments. Secondly, the lower precision rates for SOMEA and SOMBRERO could be explained by the fact that the optimal map sizes are not known. Improper map sizes can, to some extend, affect the results for multiple motif datasets.

Robustness analysis

We have conducted some analysis on the robustness of SOMEA with respect to different map sizes. We have computed recall, precision and F-measure analysis on SOMEA using the eight real datasets for map sizes 10 × 10, 15 × 15, and 20 × 20. Each dataset is run for five times and their average recall, precision, and F-measure is computed.

The computational time of SOMEA is mainly imposed by three operations: a) finding winner node for each kmer; b) updating winner and its neighboring nodes models; and c) updating node model at the end of an epoch. The time complexity of SOMEA with map size R × C is O((L × R × C) + (P × L) + (R × C)), where L is the total length of DNA sequences and P is the size of neighborhood during k-mers assignment. Here, the O(L × R × C) term is due to the computation of finding winning node for every k-mer; (P × L) operations are needed for the computation of updating the temporary model variables during the k-mers assignment stage; and (R × C) operations for updating the node models at the end of an epoch. Self-organizing map based algorithm is known to suffer from heavy computational time due to the global search to simultaneously discover all clusters. We have recorded the execution time of SOMEA for the eight real datasets and found that it has the highest average computational time of 1364s as compared to WEEDER (825s), SOMBRERO (326s), MEME (126.7s) and ALIGNACE (101s). The slower computational time of SOMEA compared to SOMBRERO is due to the fact that we have to update the ∆M _pssm and ∆M _mc parameters (see Methods) for the winner and its neighborhood during k-mers assignment (i.e. see Eqs (9) and (10)). In SOMBRERO, the update of node models only occurs at the end of an epoch. Also, some heuristic optimizations are included in it to reduce computational time. It can be observed that, current version of SOMEA requires slightly larger computational time, however, its better sensitivity and specificity performances can offer a good trade-off.

Overview

The main idea of our SOMEA algorithm is to use a hybrid node model, where a model is heterogeneously composed of PSSM and MC. We assume each node on the map is a fuzzy composition between a motif signal and background noise. Since we do not have prior knowledge on the type of each node, we use a soft competitive weighting scheme for the two components (i.e., PSSM and MC) of each node model. We refer it as intra-node competition. Our framework design is inspired by the fact that, the two sequence classes (i.e. motif and background noise) in the DNA dataset have distinctive properties. Subsequently, it is necessary to represent them using appropriate signal’s models.

The two steps above define a recursive regression process [9], where the optimal models parameters are estimated by iteratively applying the k-mers to the system. After some training epochs, similar k-mers from supposing motif or background class are projected onto the same or adjacent nodes on the 2D grid map. The k-mers projected in the vicinity on the map, generally forming clusters. This implies the similarity of their respective features. Once the nodes’ models have been stabilized, we can identify candidate motifs using a motif model evaluation metric.

Basic concepts and problem formulation

We first give some notations used in this paper, and then describe the SOMEA algorithm. Denoted by D = {S₁, S₂,…, S _N }, a DNA dataset with N sequences. Let a k-mer K _i = (b₁b₂…b _k ) be a continuous subsequence of length k in a sequence, and i = 1,…, Z, with Z is the total number of k-mers in that sequence. For a sequence with length L, there are L – k + 1 number of k-mers can be produced using k length window shifting process.

We can represent a k-mer K as a 4 × k matrix [27]. Let the matrix representation be e(K) = [a _ij ]_4×k, where (a_{1 j} , a_{2 j} , a_{3 j} , a_{4 j} ) = (A, C, G, T) and j = 1,…, k. The matrix has a column j representing certain nucleotide i at that position j in the k-mer.

A 2D SOM map is a lattice of R × C nodes, where R, C is the number of rows and columns respectively. Each node V _ij , i = 1,…, R and j = 1,…, C, has a subset of k-mers assigned to it. For convenience, we use the notation V _l to represent a node, where 1 ≤ l ≤ (R × C). The coordinate of a node V _l in the lattice is expressed as z _l = (i, j). Then, each node V _l has a parameterized model Θ _l associated with it.

Let us formulate the clustering based motif discovery task. Clustering on the k-mers dataset aims to partition the dataset into a set of non-overlapping clusters {C₁, C₂,…, C _U }, where each cluster C _i holds a subset of k-mers. In our study, each node in the SOM 2D-lattice represents a cluster (i.e. U = R × C). After forming the clusters, each cluster C _i ’s potential is evaluated as true motif using motif model evaluation metric and rank the clusters based on their obtained scores. In SOMEA, we used Maximum A Posteriori score (MAP score) as the model evaluation metric. Then, top H highest ranked clusters are selected as putative motifs, and k-mers from those clusters indicate the motif locations in the sequences.

PSSM based motif model M _pssm

We use the Position-Specific-Scoring-Matrix (PSSM) [2] to model the motif signals. The PSSM based motif model, let it denoted by M _pssm , is a matrix, i.e., M _pssm = [f(b _i , i)]_4×k, where b _i ∈ {A, C, G, T} and i = 1,…, k. Here, each entry f(b _i , i) represents the probability of nucleotide b _i in position i and . In our SOMEA, the M _pssm for a node V _l can be calculated from the k-mers in a node using the maximum likelihood principle, with a pseudo-count value added as under sample correction to the probabilistic model. We follow the Bayesian estimation method for this purpose [28]. The PSSM entries are computed as follows:

f(b _i , i) = (c(b _i , i) + g(b _i )) /(N + 1), (1)

where N is the number of k-mers, c(b _i , i) is the frequency of nucleotide b _i at position i of a set of k-mers in a node, g(b _i ) = [n(b _i ) + 0.25]/(N × k + 1) and .

Markov chain based background model M _mc

In our approach, the background signal is modeled by using the markov chain (MC) model [15]. The MC is a commonly used background signal model to distinguish over-represented motifs from background signals (e.g. in [16, 19]). The stochastic and temporal nature of this model can effectively model the complex relationship of the background bases. The MC model assumes that, the probability of occurrence of a nucleotide b _i at position i in a DNA sequence is dependent only on the occurrences of m previous nucleotides. This relationship can be expressed by the conditional probability p(b _i |b _{i – m} …b _{i –1}), where b _{i – m} …b _{i –1} are bases that precede base b _i , and m is the markov order. In our approach, the first order MC (i.e. m = 1) is used because higher order model usually requires more input data to avoid over-fitting. The maximum likelihood estimation of the conditional probability p(b _i |b _{i – m} …b _{i –1}) is given by [15] as:

(2)

where c′(x) is the number of times sub-sequence x found in a set of k-mers in a node.

Let us denote π (a, a′) to represent the conditional probability p(a′|a) of the first order MC, where a, a′ ∈ {A, C, G, T}. Then the MC transition matrix gives the background model M _mc to be used in SOMEA, i.e., M _mc = [π(a, a′)]_4×4 , where .

Similarity score

A similarity metric is needed for k-mers assignment to the nodes during the learning. The score of a k-mer K _j = (b₁b₂ …b _k ) in respect with the PSSM based model assigned to node V _l , can be computed as,

(3)

Here, k is the length of k-mer, and f(b _i , i) represents the probability of nucleotide b _i in position i. Then, the score of a k-mer K _j to the MC [15] based model M _mc of node V _l is computed as:

(4)

Here, p(b₁) is the independent and identically distribution (i.i.d) probability of nucleotide b₁ in current node, which is estimated from the k-mers of node V _l .

Hybrid model

In practice, we are unable to certainly deduce if a SOMEA’s node is a motif or background at any stage of the learning process. Also, before the system converged, the members of a node are likely to be composed of mixed signals. Therefore, neither PSSM or MC based models (i.e.M _pssm and M _mc ) alone would satisfactorily model such composition. However, we can weigh the fitness of MC and PSSM models with respect to the k-mers in a node. In other words, when a set of k-mers fit with a certain model, (i.e., either motif model given by M _pssm or background model given by M _mc ), it is more likely that those k-mers represent that class. Note that both signal models, can represent signal features from opposite class to some extent.

In this work, we aimed to combine the expression abilities of both of the models (i.e., i.e.M _pssm and M _mc ) in an unified mechanism to improve the distinguishing ability of the system, since each node given by SOMEA (or any clustering based approach) contains a fuzzy mixture of motif signals and background signals.

In implementation, we adopted a simple linear weighting scheme to combine these two models for a node V _l as follows:

(5)

Equation (5) gives a linear combination of the two models to produce a heterogeneous model for a node V _l . Here, λ is a scaling factor, for simplicity default value of λ is set as 0.5. If a k-mer K _j gets a higher score by this heterogeneous model based scoring Θ _l (K _j ), that indicates K _j has a better fit to the combined model of node V _l .

Motif ranking

Once the SOMEA is stabilized after training, we have to perform an evaluation on the nodes in order to identify the most prominent candidate motifs. The candidate motifs can be identified using either motif evaluation metric or statistical significance value. These metrics usually require the use of background sequences model for computation.

In this work, we adopt the Maximum A Posteriori score (MAP score) [19] for motif ranking. The MAP score measures the conservation property of a motif with respect to the species background sequences [19]. Since, rare motifs in the background can achieve a higher MAP score, this measure can be used to distinguish a true motif from false ones based on their scores ranking. The background sequences can be modeled by using the markov chain model generated from the intergenic sequences of a species under study. This model can be used to assign a probability of a K, namely , under the background model given by . The MAP score of a node V _l can be calculated as follows:

(6)

where N _l is the number of k-mers in node V _l and refer to background probability of a k-mer K in respect with background model . can be written as,

(7)

Here, m is the Markov chain order, k is the length of k-mers, p(b₁, b₂,…, b_m) is the probability of subsequence b₁,b₂, …, b _m , and p(b _i |b _{i – m} , b _{i – m +1},…, b _{i –1}) is the conditional probability of the subsequence b _i under b _i–m , b _{i–m +1,…,}b _{i –1} occurrence constraints. For instance, using the 3rd order model, the probability of the sequence ATGCG can be calculated as: . This background probability is usually pre-computed on the sequences of interest. In Eq (6), E(V _l ) is the Shannon’s entropy, that can be written as,

(8)

Here, f(b_i, i) is the probability of nucleotide base b _i ∈ {A, C, G, T} to occur in i-th position of the PSSM.

Adaptation process

We opted for the more speed efficient batch training scheme to update the nodes’ model parameters. This method delays the update of the model parameters at the end of an epoch. Heuristic rules are proposed to update each node’s PSSM and MC model parameters. We associate each node with three computing components including: two matrices ∆M _pssm , ∆M _mc and a counter r. Let be BMU of an input k-mer K = (b₁,b₂,…,b _k ). Denoted ∆M _pssm = [∆f(b _i , i)]_{4× k} for b _i ∈ {A, C, G, T} and i = 1,… ,k. Similarly, let ∆M _mc = [∆π(a, a′)]_4×4 for a, a′ ∈ {A, C, G, T}. We initialize all entries in both matrices ∆M _pssm and ∆M _mc as 0. Also let r = 0. Once a winning node for a k-mer K is found, the matrices of a node are updated as follows.

(9)

(10)

where is an entry of the binary matrix e(K), count(a, a′) is the frequency of di-nucleotide (aa’) in k-mer K and h is a neighborhood function. The neighborhood function h is defined as

(11)

where σ is the variance whose value is fixed throughout the learning stage. We also update r = r + 1. Upon completion of an epoch, all nodes’ model parameters will be updated as follows:

(12)

(13)

where η is the learning rate and f(b _i , i) and π (a, a′) is defined in Eq (1) and Eq (2) respectively. Note that, in the computation of Eq (12) and Eq (13), we first compute f(b _i , i) and π (a, a′) using the current set of kmers assigned to a node.

It is also necessary to update the weighting parameters α and β. Assuming a set of N _l k-mers {K₁ …, K _Nl } is assigned to a node V _l at the end of an epoch, the weighting parameters update equations are

(14)

and

β _new = 1 – α _new . (15)

Training

Let the BMU index for a k-mer K is q(K).

for epoch=1 to max_epoch do

for each K ∈ X do

end for

Update model parameters of all nodes using Eqs (12) and (13).