Gapped k-mer SVM classifier outperforms k-mer SVM classifier

Because of the difficulty of reliably estimating long k-mer counts, we hypothesized that gkm-SVM would perform better than kmer-SVM, and that gapped k-mers would be most advantageous as features, when long TFBSs are important sequence elements in a given data set. To directly test this idea, we compared the classification performance of gkm-SVM to kmer-SVM in predicting the binding sites of CTCF [13] in the human genome, a TF whose binding specificity has been well-characterized [14]. As shown in Figure S1, CTCF recognizes very long DNA sequences (the full PWM is 19 bp), and the genomic CTCF bound regions are almost perfectly predicted by matches to the CTCF PWM (Figure S2): in the PWM analysis, we used as a predictor the best matching log-odd score to the PWM model in the region, and achieved area under the ROC curve (AUC) of 0.983. It is very rare for a single PWM to perform this well, and in our experience CTCF is unique in this regard. The CTCF dataset therefore provides an excellent opportunity to test our gapped k-mer classifier. We used the top 2,500 CTCF ChIP-seq signal enriched regions in the GM12878 cell line available at Gene Expression Omnibus (GSE19622) [13] as a positive dataset, and equal numbers of random genomic sequences (1×) as a negative dataset. We generated these negative sequences by matching length, GC and repeat fraction of the positive set [6].

We compared the performance of gkm-SVM and kmer-SVM on the CTCF data set for a range of oligomer lengths by varying either k (for kmer-SVM) or l (for gkm-SVM) from 6 to 20. We fixed the parameter k=6 for gkm-SVM. We then quantified the classification performance of each by calculating test-set AUC with standard five-fold cross validation (CV) (see Methods). Figure 1A shows a summary of the comparisons. As anticipated, gkm-SVM performs consistently better than kmer-SVM for all lengths. More significantly, while kmer-SVM suffers severely from overfitting when k is greater than 10, gkm-SVM is virtually unaffected by l. In fact, gkm-SVM achieves the best result (AUC=0.967) when l=14 and k=6, which is significantly better than the kmer-SVM (AUC=0.912 when k=10); the best ROC curve is shown in Figure 1C. It should be noted, however, that the PWM classification result (Figure S2) is still the best (AUC=0.983) among the three methods we tested in this analysis. A complicating factor is that while both kmer-SVM and gkm-SVM use entire sequences (average length is 316 bp) to calculate the prediction scores, the PWM scores are from the best matching 19 bp sub-sequence in the region. It may be that the extra ~300 bp sequences contribute noise in the SVM prediction scores, which slightly impairs the overall classification accuracy. In any event, the gkm-SVM is a significant improvement in accuracy over the kmer-SVM, and both gkm-SVM and the PWM are excellent predictors on this dataset.

Open in a separate window

Figure 1

gkm-SVM outperforms kmer-SVM over a wide range of k-mer length.

Both gkm-SVM and kmer-SVM were trained on (A) CTCF bound and (B) EP300 bound genomic regions using different word lengths (k for kmer-SVM and l for gkm-SVM). The parameter k for gkm-SVM was fixed at 6. While AUCs of the kmer-SVMs show significant overfitting in both cases as k gets larger (dotted), gkm-SVMs accuracy is higher for a broad range of larger l (solid). Results using the truncated Gkm-SVM with m _max=3 are shown as dashed lines and AUCs of these faster approximations are comparable when the difference between m_max and l – k are relatively small. ROC for the optimal k or l for each case are shown in (C) and (D). Gkm-SVMs (solid) consistently outperform kmer-SVMs (dashed) on both data sets. Error bars here and below represent 5-fold CV standard deviation.

Interestingly, gkm-SVM shows consistently better performance than kmer-SVM even if l is relatively small (l<10) (Figure 1A). This suggests that gkm-SVM may also be better at modeling diverse combinations of TFBSs than kmer-SVM. To test this hypothesis, we analyzed a mouse enhancer dataset of more varied sequence composition: genomic EP300 bound regions in embryonic mouse forebrain [15]. We have previously shown that our original kmer-SVM classifiers can accurately predict EP300 binding when mediated by sets of active TFBSs [3]. This EP300 data set thus provides a direct test of the effectiveness of using gapped k-mer features to detect more complex regulatory features. For this analysis, we defined a new set of the 1,693 400 bp sites that maximize the EP300 ChIP-seq signal within each of the peaks determined by MACS [16] after removing any regions which were more than 70% repeats. We repeated the k and l scaling with the EP300 data set and a 1× negative set, and again found that gkm-SVM consistently outperforms kmer-SVM for all feature lengths (Figure 1B). Analogous to the observations modeling CTCF binding, gkm-SVM AUC is high and does not degrade with large l. In contrast, the kmer-SVM accuracy drops rapidly as k increases. Moreover, although the difference in performance is smaller than found for the CTCF data set, the gkm-SVM achieves the best AUC (0.947) with l=9 and k=6, while the kmer-SVM achieves 0.932 with k=7, suggesting that longer k-mers with some flexibility do contain more complete information about TF binding (Figure 1D). At the same time, the gapped k-mer features are more robustly estimated (having more counts) and for this reason make more reliable predictors. The consequences of these improvements in AUC are significant when considering the genome-scale precision of the improved gkm-SVM classifiers. The rate of false positive predictions is dominated by the large neutral fraction of the genome, so the precision of a genome-scale classifier is best assessed by a Precision-Recall curve in combination with a much larger negative set, as discussed in Ref. [3]. The Precision-Recall curves for the gkm-SVM and kmer-SVM classifiers on a 100× negative set are shown in Figure S3. For CTCF, at a recall of 50%, the precision increases from 36% to 59%. These ranges of precision and recall are in the relevant range of experiments aiming to discover and test novel enhancers, and we therefore expect that predictions based on gkm-SVM will have up to a two-fold higher successful validation rate.

One further modification can substantially reduce the computational cost of using gapped k-mers with little degradation in performance. The algorithm using the k-mer tree data structure produces identical results to the direct evaluation of Equation (3), but typically is much faster when the number of mismatches, l – k, is smaller than four, and the number of training sequences is large. The k-mer tree algorithm can be made even more computationally efficient, if we prune the traversal of the tree, by ignoring any k-mer pairs that have more mismatches than a predetermined parameter, m_max. This provides an approximation to the exact kernel calculation, but the approximation error is usually negligible given that the coefficient h_m for large numbers of mismatches are generally much smaller compared to those with small m. This approximation significantly reduces the total number of calculations and allows the user to control the running time of the algorithm by setting the parameter m_max, and makes the use of longer word lengths l feasible for any given k. To systematically investigate the classification performance of this approximation, we applied the same analysis above using both CTCF and EP300 data sets (Figure 1A, B), and found that AUCs from the approximate method are virtually identical to the exact method when the difference between m_max and l – k are small. Interestingly, the approximation method achieved even higher AUC with CTCF data set in some cases.

Encouraged by the analyses of CTCF and EP300 data sets above, we systematically compared gkm-SVM to kmer-SVM using a very broad range of human data sets generated by the ENCODE project [17],[18]. We used 467 sets of ChIP-seq peaks produced by the ENCODE uniform processing pipeline containing at least 500 regions (see Methods). We truncated any data set with greater than 5,000 regions by random sampling. We then trained both kmer-SVM and gkm-SVM on each set against an equal size (1×) negative set of random genomic regions and calculated AUCs with five-fold cross validation. We used k=6 for kmer-SVM, and l=10 and k=6 for gkm-SVM, but as shown in Figure 1 the improvements are generally insensitive to these parameter choices. Strikingly, we find that gkm-SVM almost always outperforms kmer-SVM (Figure 2A). We also find that variances of AUCs from test CV sets are generally reduced, suggesting that gkm-SVM is more robust than kmer-SVM (Figure S4). More significantly, gkm-SVM performs much better especially for TFs with long binding sites. In this dataset, most of these long binding sites arise in ChIP-seq data sets for CTCF and members of the cohesin complex (RAD21, SMC3) known to be physically associated with CTCF [19]. On these CTCF associated factors gkm-SVM exhibits much higher AUC than kmer-SVM, as highlighted by the cluster of purple circles in Figure 2A. We have also compared gkm-SVM to the best single PWM AUC as shown in Ref. [6] (Figure 2B). As expected, gkm-SVM outperforms all datasets except CTCF, for which gkm-SVM performance is only marginally reduced. For a consistent analysis of this dataset, we used l=10 and k=6, although for CTCF the gkm-SVM performance is optimal at larger l, as seen in Figure 1A.

Open in a separate window

Figure 2

gkm-SVM consistently outperforms kmer-SVM and the best known PWM on human ENCODE ChIP-seq data sets.

(A) We trained gkm-SVM and kmer-SVM on the complete set of 467 ENCODE ChIP-seq data sets (with k=6 for kmer-SVM, and l=10 and k=6 for gkm-SVM). gkm-SVM AUC is consistently higher than kmer-SVM with only a few very minor exceptions. The gkm-SVM method specially outperforms the kmer-SVM for the data sets bound by members of the CTCF complex, highlighted as purple circles. (B) We also compared gkm-SVM and the best known PWM on the same data sets, and gkm-SVM AUCs are significantly higher than the PWM AUC in almost all cases. (C) The ENCODE data sets were divided into four groups: (1) no PWM, (2) only one PWM, (3) two PWMs, and (4) three or more PWMs identified by Wang et al. Then, for each group except the first one, we calculated the number of PWMs recovered by our method. At least one PWM was recovered for more than ~90% of the data sets.

Motif analysis of the ENCODE ChIP-seq data sets

The predictive sequence features that allow gkm-SVM to outperform the single best PWM imply that cooperative binding is the underlying molecular mechanisms that targets TFs to these regulatory regions. Previously we have typically focused on a handful of the highest SVM weight k-mers (say top ten positive and top ten negative weight k-mers) to interpret the classification results [3],[6],[20]. This simple method becomes unwieldy when applied to the gkm-SVM results because of the large number of very similar significant features (when l and/or k are large). Although the k-mers at the extreme top and bottom tails of the k-mer weight distribution are still important and biologically meaningful, those k-mers usually cover only a fraction of the significant feature set, and many more important features are included in the larger tails of the k-mer weight distribution. Therefore, more sophisticated algorithms are needed to extract the biologically relevant features from the classification results.

To directly address this issue, we developed a new method to combine multiple similar k-mers into more compact and interpretable PWMs and analyzed the 467 ENCODE data sets [18]. In this approach, we used a larger number of predictive k-mers to build de novo PWMs (see Methods). We used the top 1% of 10-mers from each of the gkm-SVMs trained on the ENCODE data sets and identified up to three distinct PWMs (Figure S5) from k-mers in this set. We then compared our results with the previous PWMs found in the same data sets using a conventional tool (MEME-ChIP) [18] [21]. Similar to our approach, Wang et al. analyzed 457 ENCODE ChIP-seq data sets (440 sets are in common with those we analyzed above) and identified five PWMs from each data set. Collectively, Wang et al. found 79 distinct PWMs enriched, of which our method recovered 74. Comparing each ChIP-seq data set individually, we recovered most of the PWMs reported by Wang et al. using our method (Figure 2C). Interestingly, while Wang et al. largely failed to identify biologically meaningful PWMs from most of the POL2 ChIP-seq data sets (47 out of 58 sets returned no meaningful PWMs), our methods frequently identified cell-specific TFs as well as promoter specific TFs (Figure S5). For example, the GATA1 TF identified from POL2 ChIP-seq in the erythroleukemic cell line K562 is known to play central roles in erythroid differentiation [22]. The ETS1 TF from HUVEC is another extensively studied TF, known to be important for angiogenesis [23]. A major difference between the two methods is the number of training sequences. While the previous study was limited to the top 500 of ChIP-seq peaks (ranked by ChIP-seq signal), we were able to use 10× larger numbers of ChIP-seq peaks (5,000 regions), and the large training sizes enabled us to robustly identify diverse combinatorial sequence features.

Comparison to previous kernels

Since the early development of k-mer based supervised machine learning techniques [24], there have been a number of improvements. Some of these extend the feature set to include imperfect matches, similar in spirit to our gkm-SVM. The mismatch string kernel [8] is one such method, originally motivated by the fact that homologous protein sequences are not usually identical and have many frequently mutated positions. The mismatch kernel also uses k-mers as features, but allows some mismatches when counting k-mers and building feature vectors. The wildcard kernel [9] is another variant of the original string kernel, which introduces a wildcard character that matches any single letter in the given alphabet. More recently, an alternative di-mismatch kernel [10] has been proposed to directly model TFBSs, and has been successfully applied to protein binding microarray (PBM) data sets [25] and several other ChIP-seq data sets [10],[11]. The di-mismatch method tries to overcome the limitation of the mismatch kernel by favoring k-mers with consecutive mismatches. However, in a recent comparison of methods for modeling transcription factor sequence specificity, full k-mer methods outperformed the di-nucleotide approaches when applied to PBM data [26].

To further evaluate our proposed method, we directly compared the gkm-kernel with the aforementioned three alternative methods, Mismatch kernel [8], Wildcard kernel [9], and Di-mismatch kernel [10],[11], using the mouse forebrain EP300 data set. As shown in Figure 3, gkm-kernel outperforms the other three existing methods both in terms of the classification accuracy and running time. The best AUC we achieved for gkm-kernel is 0.947 as compared to 0.937, 0.935, and 0.944 for the wildcard kernel, mismatch kernel, and di-mismatch kernel, respectively (Figure 3A). Although the wildcard kernel and gkm-kernel are quite similar, the systematic improvement in gkm-kernel AUCs is primarily due to the incorporation of reverse complement sequences. We directly tested this by adding reverse complement sequences to the feature set for the previously published methods, and indeed found that with this modification, these methods were also able to achieve comparable AUCs (Figure S6).

Open in a separate window

Figure 3

Comparison of gkm-SVM and existing methods on the mouse forebrain EP300 data set.

(A) For each method, averages of 5-CV AUCs are shown as a function of the word length with the optimal number of mismatches, m, held fixed. Also shown are gkm-SVM results using fixed k=6 and varying m _max. (B) Running time for each of the kernel computations shown in (A). Gkm-kernels show better classification performance and significantly more efficient computation at peak AUC.

To further evaluate our proposed method, we directly compared the gkm-kernel with the aforementioned three alternative methods, Mismatch kernel [8], Wildcard kernel [9], and Di-mismatch kernel [10],[11], using the mouse forebrain EP300 data set. As shown in Figure 3, gkm-kernel outperforms the other three existing methods both in terms of the classification accuracy and running time. The best AUC we achieved for gkm-kernel is 0.947 as compared to 0.937, 0.935, and 0.944 for the wildcard kernel, mismatch kernel, and di-mismatch kernel, respectively (Figure 3A). Although the wildcard kernel and gkm-kernel are quite similar, the systematic improvement in gkm-kernel AUCs is primarily due to the incorporation of reverse complement sequences. We directly tested this by adding reverse complement sequences to the feature set for the previously published methods, and indeed found that with this modification, these methods were also able to achieve comparable AUCs (Figure S6).

More significantly, when we compare running times at parameters which maximize AUC for each method, our gkm-SVM implementation (l=9,l-k=3) is roughly two orders of magnitude faster than di-mismatch (10,3), and slightly more than one order of magnitude faster than mismatch (l=10,m=2) and wildcard (l=8,m=3) on the EP300 data set (Figure 3B and Figure S7). Also, by fixing k=6 and the parameter m_max in our algorithm, the AUC becomes less sensitive to the feature length l, compared to a scan at fixed m, varying k (Figure 3A). Direct running time comparisons using our tree structure in the mismatch and wildcard kernels (described below) are shown in Figure S7B and S7C. We should note that we were only able to test the di-mismatch kernel up to l=10, because it required more than 128 GB of memory and did not finish within 2000 minutes when using l=11.

Interestingly, we also note that both Mismatch kernel and Wildcard kernel are special cases of the more general class of kernels, defined by Equation (3). This unification allows direct application of the methods we developed for mismatch profile computation and therefore gives more efficient methods for computation of these existing methods (see Methods).

Calculation and performance of estimated l-mer frequencies for gkm-SVM

As an alternative to the gapped k-mer feature set, we also developed an alternative kernel by replacing the k-mer counts with robust l-mer count estimates [7] in our original kmer-SVM framework. We have developed efficient methods to compute this new kernel (see Methods). In Ref. [7], we considered the mapping from l-mers to gapped k-mers. Among all possible sets of l-mer frequencies that could produce the same gapped k-mer frequency distribution, we developed a method to estimate the “most likely” l-mer frequency set. Full details of this method are described in the Ref. [7]. In brief, we first define a gapped k-mer count vector similar to the definition of the gapped k-mer feature vector for gkm-SVM as shown above. Then, the count estimate, , for l-mer u is given by

(5)

The weight in Equation (5) was shown to only depend on the number of mismatches, m, between the gapped k-mer corresponding to and u, and takes the following form:

(6)

where b is the alphabet size and is equal to four in case of DNA sequences (A,C,G and T). Since the above equation is applied to every l-mer, it would provide a non-zero frequency even for an l-mer that does not have any exact match appearing in any training set sequence.

Direct calculation of Equation (5), however, requires actual counting of all of the M gapped k-mers, which becomes computationally intractable for large l and k in a way similar to Equation (1). Besides, summing up a large set of floating point numbers may result in poor numerical precision. To overcome these issues, we developed a simple method, referred to as the gkm-filter, to more efficiently calculate the robust l-mer count estimates, , without calculating the intermediate gapped k-mer counts (see Methods). In summary, in the calculation of the robust l-mer count estimates, we give a non-zero weight to l-mers with few numbers of mismatches. The k-mer frequency estimation method is not constrained to produce non-negative frequencies and may occasionally generate negative count estimates. To obtain strictly positive frequencies, we used a revised version of the gkm-filter method, which we call the truncated gkm-filter. Finally, we developed a method to directly calculate the kernels using these feature sets (see Methods). An important result here is that the evaluation of the gkm-kernel (the inner product of the l-mer count estimates vectors) is still given by Equation (3), but with a new set of weights c_lk(m) given by Equation (14), below, replacing h_lk(m). Therefore, efficient algorithms for pairwise mismatch profiles that we developed for the gkm-kernel can be directly used for this new feature set without any modification. Because of this symmetry, we also refer to this method as gkm-kernel with (full or truncated) filter. A numerical example using count estimates on two short sequences is provided in Text S1.

To systematically compare the classification performance of these new methods with the original gapped k-mers, we repeated the previous analysis with the ENCODE ChIP-seq data sets. Using the truncated gkm-filter yields results highly comparable to the original gkm-SVM for most datasets with modestly but consistently better relative performance when AUC is greater than 0.9 (shown as purple circles in Figure S8A). Any improvement in the range of high AUC (>0.9) typically strongly reduces the classifier's False Prediction Rate [27], therefore, we generally recommend the truncated filter method as the method of choice for most analyses. Compared to the original gkm-SVM, using the gkm-SVM with full filter yields lower AUCs (Figure S8B) although it is still significantly higher compared to the kmer-SVM method.

Direct computation of Gkm-kernel

For direct computation of the gkm-SVM kernel matrix, we represent each training sequence with a list of l-mers and corresponding count for each l-mer. Then for each pair of sequences, we compute the number of mismatches for all pairs of l-mers and use the corresponding coefficient h_m to obtain the inner product of Equation (3). As the number of unique l-mers in each sequence is L and the number of sequences is N, this algorithm would require O(N²L ²) comparisons. In addition, a naive algorithm for counting the number of mismatches between two l-mers (i.e. the hamming distance) would be O(l). Our implementation employs bitwise operators, providing a constant-factor speedup. Briefly, using two bits to represent each base (A,C,G and T), we used an integer variable to represent non-overlapping substrings of t base pairs of the l-mer, therefore using total integers to represent each l-mer, where is the ceiling function. For counting the number of mismatches, we take the bitwise XOR (exclusive OR) of the integer representations of the two l-mers and use a precomputed look-up table to obtain the total number of mismatches using the XOR result. This method requires a look-up table of size 2^2t. The optimal value of t depends on the processor architecture and amount of cache memory. We used t=6 for our analysis.

Gkm-kernel with k-mer tree data structure

As depicted in Figure 5, we use a k-mer tree to hold all the l-mers in the collection of all of the sequences. We construct the tree by adding a path for every l-mer observed in a training sequence. Each node t_i at depth d represents a sub-sequence of length d, denoted by s(t_i), which is determined by the path from the root of the tree to the node t_i. Each terminal leaf node of the tree represents an l-mer, and holds the list of training sequence labels in which that l-mer appeared and the number of times that l-mer appeared in each sequence. As an example, Figure 5 shows the tree that stores all the substrings of length l=3 in three sequences S₁=AAACCC, S₂=ACC, and S₃=AAAAA. Then, to evaluate the mismatch profile we traverse the tree in a depth-first search (DFS) [35] order. In contrast to the mismatch tree used in Ref. [8], here for each node t_i, at depth d, we store the list of pointers to all the nodes t_j at depth d for which s(t_i) and s(t_j) have at most l – k number of mismatches. We also store the number of mismatches between s(t_i) and s(t_j). Similar to the mismatch tree [8], we do not need to store these values for all the nodes in the tree, but we compute them recursively as we traverse the tree. When reaching a leaf node, we increment the corresponding mismatch profile N_m(S_i, S_j) for each pair of sequences S_i in that leaf node's sequence list, and all the S_j's in the list of sequences in the pointer list for that leaf node. At the end of one DFS traversal of the tree, the mismatch profiles for all pairs of sequences are completely determined.

Open in a separate window

Figure 5

Fast computation of mismatch profiles using k-mer tree structure.

As an example, we use l=3 and three sequences S₁=AAACCC, S₂=AAAAA, and S₃=ACC to build the k-mer tree. The leaves (nodes at depth d=l=3) correspond to 3-mers AAA, AAC, ACC, and CCC. The sequence ID and the number of times each 3-mer appeared in each sequence are stored for each leaf. Each node t_i at depth d represents a sequence of length d, denoted by s(t_i), which is determined by the path from the root of the tree to t_i. For example, s(t₂)=C and s(t₄)=AC. DFS is started at the root node, t₀. When visiting each node t_i, at depth d, we compute the list of all the nodes t_j at depth d for which s(t_i) and s(t_j) have at most m_max mismatches. We also compute the number of mismatches between s(t_i) and s(t_j). When reaching a leaf, we increment the corresponding mismatch profile N_m(S_i, S_j) for each pair of sequences S_i in that leaf and S_j in the list.

To increase the speed further, we also introduce an optional parameter m_max, which limits the maximum number of mismatches. By setting m_max smaller than l – k, we only consider l-mer pairs that have at most m_max number of mismatches. This can reduce calculation significantly by ignoring l-mer pairs which potentially contribute less to the overall similarity scores. This method provides fast and efficient approximations of the exact solution. In addition, we only compute the lower triangle of the matrix because of the symmetry in the kernel matrix. Hence, at each node t_i, we exclude the nodes t_j in the list that have maxID(t_i)<minID(t_j), where minID(t_i) and maxID(t_j) are the maximum and minimum sequence ID in the subtrees of t_i and t_j respectively and are computed and stored for each node at the time we build the tree.

Analysis of de novo PWMs from gkm-SVM

We developed a new method for building de novo PWMs by systematically merging the most predictive k-mers from a trained gkm-SVM. We first determined a set of predictive k-mers by scoring all possible 10-mers and selecting the top 1% of the high-scoring 10-mers. We then found a set of distinct PWM models from these predictive 10-mers using a heuristic iterated greedy algorithm. Specifically, we first built an initial PWM model from the highest scoring 10-mer. Then, for each of the remaining predictive 10-mers, we calculated the log-odd ratios of all possible alignments of the 10-mer to the PWM model, and identified the best alignment (i.e. the position and the orientation that give rise to the highest log-odd ratio value). Since multiple distinct classes of TFBSs are expected to be identified in most cases, we only considered 10-mers with good alignments (i.e. we used threshold of 5.0 for log-odd ratio scores relative to a genomic GC=0.42 background). After each of the 10-mers was aligned, we updated the PWM model only with successfully aligned 10-mers. To further refine the PWM, we repeated this by iterating through all of the top 1% 10-mers until no changes were made. When updating the PWM model, we assumed that the contribution of each k-mer is exponentially weighted proportional to its SVM score, using exp(α w_i), with α=3.0. The 10-mers used for creating the 1^st PWM were then removed from the list, and the process was repeated on the remaining predictive k-mers, to find up to three PWMs. Lastly, we matched our PWMs to the previously identified PWMs [18] using TOMTOM [36] software. Each of the PWMs identified by our method were associated with Ref. [18] PWMs if the q-value (false discovery rate) <0.05.

Implementation of mismatch and wildcard kernels using gkm-kernel framework

In the gkm-kernel, we define the feature vector to consist of the frequency of all the l-mers with exactly k known bases and l – k gaps. In contrast, the wildcard kernel [9] also includes all the l-mers with l – k wildcards, where l – k ranges from 0 to the maximum number of wildcards allowed, M. Thus in the wildcard kernel, the parameter M replaces k in our gkm-kernel. In the sum, these are weighted by λ^{l - k} to penalize sequences with more wildcards [9]. We derived an equation to directly compute the inner products from the mismatch profiles without the need to calculate the actual gapped k-mer counts. Here we show that a similar approach can be used to calculate the wildcard kernel. We derive a new set of coefficients that can substitute h_m, in Equation (3). To evaluate we only need to consider the contribution of each pair of l-mers with m mismatches in the inner product of the corresponding feature vectors of the two sequences. Equation (8) gives those weights:

(8)

Using the above form allows us to directly use the fast algorithms we have developed for calculation of the mismatch profiles to calculate the wildcard kernels. Although there are similarities between our tree algorithm and the tree algorithm described in Ref. [9], there are some key differences. In the Ref. [9], the algorithm literally transverses all the possible gapped l-mers (with maximum M number of gaps) while our algorithm takes advantage of the fact that the final inner product will only depend on the number of pairwise mismatches and hence only traverses all the l-mers that are present in the input data. Another difference is that Ref. [9] uses a list of all partially matching l-mers at each node of the tree, while we use a list of pointers to tree nodes instead. So, for example, at the beginning of the algorithm (at depth d=0) they start with a large list consisting of all the possible l-mers in the input data, while in our algorithm the list at depth d=0 consists of only one node (the root of the tree). Using this representation of all the partially matching l-mers, we can more efficiently perform the comparisons at each step of the algorithm when the tree is dense.

In the mismatch string kernel described in Ref. [8] and [9], the feature vectors consist of the counts for all the l-mers with maximum distance M from the l-mers in the sequence. Here we show that the approach above can be used to implement the mismatch kernel. Again, the only difference is in the set of weights used in Equation (3). To calculate the mismatch string kernel value for two sequences we replace h_lk(m) in Equation (3) by :

(9)

where b is the alphabet size (b=4 for DNA sequences) and r=m ₁+m ₂−m−2t. Given two l-mers x ₁ and x ₂ where x ₁ and x ₂ differ in exactly m places, the term inside the summations counts the number of all possible l-mers that exactly differ x ₁ in m ₁ places and x ₂ in m ₂ places t of which fall in the common l-m bases of x ₁ and x ₂. (See Figure S9). So the result of the summation is the number of all l-mers that differ x ₁ and x ₂ in at most M places. This form for the mismatch string kernel has the advantage that we can directly use equation (3) to compute the kernels by only having the mismatch profiles that can be computed more efficiently.

Gkm-filters for computation of the robust l-mer count estimates

To compute the l-mer count estimates by using Equation (5), one should first calculate the gapped k-mer counts, y_i, and then use Equation (5) to combine the y_i with a weight corresponding to the number of mismatches, given by Equation (6). This is shown schematically in Figure S10. The mapping from observed l-mer counts to gapped k-mer counts is performed by the matrix A, whose elements are a_ij. If the gapped k-mer v_i matches l-mer u_j, then a_ij=1, otherwise a_ij=0. There is a second matrix W, which performs the mapping from gapped k-mer counts to estimated l-mer counts, and whose elements are w_ij. In a previous work we showed that matrix W is the Penrose-Moore pseudo-inverse of A [7]. The element w_ij only depends on the number of mismatches between the l-mer u_i and the gapped k-mer v_j, and is given by Equation (6). Here we show that, for efficient computation, we can combine the two mapping matrices, A and W, and directly calculate the minimum norm l-mer count estimates from actual l-mer counts in a sequence. We refer to this combined mapping as the gkm-filter. The combined mapping matrix G=WA, has elements g_ij, shown on the bottom of Figure S10. As shown below, g_ij also only depends on the number of mismatches, m, between the l-mers u_i and u_j. We denote these values by g_lk(m) and refer to this as the gkm-filter since the domain and range of this mapping is the same.

To obtain the element g_lk(m), that gives the weight for the contribution of an observed l-mer u_i in the training set to the minimum norm l-mer count estimate u_j that has exactly m mismatches with u_i, we sum over the contribution of all the gapped k-mers v_τ that match u_i. Note that a_ij=0 for all other gapped k-mers. There exist different gapped k-mers that match u_i and have exactly m mismatches with u_j. Figure S11 shows how we enumerate all these gapped k-mers. The black solid circles denote the m mismatch positions of u_i and u_j, the gray circles denote the l – m match positions and the empty dotted circles denote the l – k gap positions. For a gapped k-mer to have exactly t mismatches with u_j, there are ways to select the t mismatch positions and ways to select the k – t match positions. Now considering the weight w(t) for the gapped k-mers with t mismatches, the gapped k-mer filter elements, g_lk(m) can be obtained as follows:

(10)

In other words, there are different ways we can construct a gapped k-mer that matches u_i, and has exactly t mismatches with u_j, by selecting t positions from the m mismatch positions and k – t positions from the l – m match positions as explained above (Figure S11). It can be easily shown that g_lk(m) is a polynomial of degree k in m. Now using the weights given in Equation (10), for any given l-mer, u we finally obtain the minimum norm l-mer count estimate as follows:

(11)

where N _tr(u, m) is the number of l-mers with exactly m mismatches with u in the training set. For large values of l and k, the number of all possible gapped k-mers gets exponentially large and since this method avoids evaluating the gapped k-mer counts, it significantly reduces the cost of calculating the l-mer count estimates compared to the original method we developed in Ref. [7].

In summary, we defined a generalized k-mer count (referred to as the robust l-mer count estimates) by giving a non-zero weight to l-mers with few number of mismatches (In the conventional k-mer count only perfectly matching k-mers are counted). These weights are given by g_lk(m). Figure S12A shows the plots for g_lk(m) for l=20 and various values of k. Each plot is normalized so that weight corresponding to zero mismatches is equal to one. The case with l=k is equivalent to the conventional k-mer count. Also Figure S12B shows g_lk(m) for and various values of l. With a fixed length l, higher values of k result in smaller coefficients for larger mismatches, and therefore less smoothing of the estimated counts (Figure S12). Moreover, g_lk(m) can become slightly negative for large numbers of mismatches. This is because in our estimation of the frequencies we did not restrict the frequencies to be positive, and doing so would yield a more complicated expression. The assumed Gaussian distribution allows non-physical negative frequencies to have non-zero probability. A beta-distribution would not have this problem but would introduce offsetting complications. In cases where the estimated counts are required to be strictly positive, such as when we need to calculate the logarithm or ratios of the estimated frequencies, we truncate the gkm-filter g_lk(m) by setting g_lk(m)=0 for every m≥m ₀, where m ₀ is the smallest number of mismatches for which g_lk(m ₀)<0. This will give an approximation to the value of in Equation (5), so it will no longer strictly be the minimum norm estimate, but it will guarantee that all the count estimates are non-negative.

Gkm-kernel with l-mer count estimates

Given a sequence S, we define an l-mer count estimate vector where N is the number of all l-mers (4^l in case of DNA sequences), and is the estimated count of the i ^th l-mer appearing in sequence S using Equation (11). Then, we can calculate a standard linear kernel simply by using this vector in Equation (1). Similar to the gkm-kernel method, we can further simplify this equation using the same technique introduced in Equation (2) which does not involve the computation of individual l-mer estimates. We show that the inner product of the two l-mer count estimate vectors can be obtained as follows:

(12)

where n ₁ and n ₂ are the number of l-mers in S ₁ and S ₂, and is the i'th l-mer in S ₁ and is the j'th l-mer in S₂. If u ₁ and u ₂ have exactly m mismatches then c(u ₁, u ₂)=c_m. Grouping all the l-mer pairs with m mismatches, we can rewrite Equation (12) as follows:

(13)

where N_m(S ₁, S ₂) is the mismatch profile of S ₁ and S ₂ as previously defined in Equation (3). We show that the weight c_lk(m), denoted in short by c_m, can be obtained as:

(14)

where r=m ₁+m ₂−2t−m, b is the alphabet size. The summations are taken over the range 0 to l. Figure S13 shows how we obtained the equation for c_m, similar to the previous development shown in Figure S11. Given two l-mers u ₁ and u ₂, with mismatches and l – m matched positions, we want to enumerate the number of all possible l-mers, u, that have m ₁ mismatches with u ₁ and m ₂ mismatches with u ₂. For this, we assume that t of the m ₁ mismatches are among the l – m match positions and m ₁ – t of them are among the m mismatch positions. There are ways to choose these m ₁ positions and (b – 1)^t choices for the values of the t mismatches. These t mismatches plus the m−(m ₁−t) unselected mismatch positions also do not match u ₂. Then, for the remaining r=m ₂−(t+m−(m ₁−t)) mismatches for u ₂ there are ways to select the positions and (b – 2)^r ways to select the values. Hence the total number of l-mers, u with m ₁ mismatches with u ₁ and m ₂ mismatches with u ₂, where t of the mismatches of u ₁ and u are among the l – m match positions of u ₁ and u ₂ is given by .

Using matrix notation, we can further show that c_m=g_m if we use the full filter g_lk(m). To see this, note that where and are the l-mer count vectors for S ₁ and S ₂. Given G=WA, we have . Hence, . Here A is the binary incidence matrix that maps l-mer counts to gapped k-mer counts as defined in Ref. [7] and W is the Moore-Penrose pseudo-inverse of A. Note that this result does not hold for the truncated filter g_m. In that case, we directly use Equation (14) to obtain c_m coefficients.

ROC curves

To compare the performance of different classification methods, we calculated the area under the receiver operating characteristic (ROC) curve for each classifier. To plot the ROC curves and calculate area under the curves (AUCs) we used the ROCR package [37] in R.

Cross validation

Following standard five-fold cross validation procedures, we divided the positive and negative sets into five segments, left one segment out as the test set and used the other four segments for training. We repeated for all of the five segments and calculated the mean and standard error of the prediction accuracy on the test set elements.

ENCODE ChIP-seq datasets

The ENCODE ChIP-seq datasets were downloaded from ftp://ftp.ebi.ac.uk/pub/databases/ensembl/encode/integration_data_jan2011/byDataType/peaks/jan2011/spp/optimal/hub/.

We gratefully acknowledge contribution of source code for the wildcard and di-mismatch approaches from Rui Kuang and Christina Leslie for the direct comparisons of these methods.