BACKGROUND
eukaryotic primary mrnas consist of exons and introns.
the mature transcript as the substrate for translation is produced by removing introns in a process called splicing.
splicing can be either constitutive, always producing the same mrna, or alternative, by skipping of variable parts of the primary transcript.
alternative splicing is a mechanism for producing transcript and protein diversity  <cit> .
it is particularly widespread in higher eukaryotes, especially mammals.
various studies have estimated that up to 74% of all human genes are alternatively spliced.
large scale detection of alternative splicing is usually done using expressed sequence tags   <cit>  or microarrays .
since alternative splicing can be highly specific for tissues or developmental stages, these methods can only detect splice events that occur in the underlying probe samples with sufficient frequencies and/or are limited to by the microarray design.
furthermore, nowadays whole genome shotgun  sequencing projects are churning out genomic data at a higher rate than corresponding transcriptome data – the number of ests in genbank release  <dig> had increased by 19% in one year, compared to a gain of 39% in the number of contigs in the wgs genbank division  <cit> .
thus it can be expected that in the foreseeable future, we shall have several genomes without the level of corresponding extensive transcript coverage required to reveal the extent of alternative splicing, and hence transcriptomic and proteomic variability.
accordingly, there is a need for in silico methods of detecting alternative splicing.
moreover, such methods can provide further insights into the mechanisms of alternative splicing.
exon skipping, whereby a given exon in its entirety is either included in, or excluded from the mature transcript, is the most prevalent form of alternative splicing in humans  <cit> .
it has been shown that sequence-based features, derived from the exon and its flanking introns, can be used to predict skipping of exons that are conserved between human and mouse and alternatively spliced in both species; denoted conserved exon skipping events  <cit> .
previous studies have used such features with state-of-the-art classifiers such as support vector machines   <cit>  and regularized least-squares classifier  <cit> , and achieved success in predicting exon skipping.
other approaches use protein domain information  <cit>  and evolutionary conservation  <cit>  to detect alternative splice events.
here, we use bayesian networks , a state-of-the-art machine learning method, to predict conserved exon skipping events.
bns are an increasingly popular machine learning approach to data modeling and classification  <cit> .
the ability of bns to cope with features of various value ranges and to learn dependencies between features makes them especially versatile and suited to a large variety of applications.
bns allow multiple dependencies between variables, impose no fixed ordering of variables, allow integration of arbitrary features, and the network structure can be automatically learned.
this makes bns a flexible choice for biological sequence data analysis  <cit> .
we introduce several novel features that distinguish alternative exons from constitutive exons, including features based on the single-strandedness of exonic splicing enhancers and silencers , and features involving intronic splicing regulatory elements .
by validating our classifiers on various datasets, we identify features which are discriminative irrespective of dataset-specific biases, and provide independent measures of the predictive power of the bns.
even though conservation based features have proved to be among the most discriminative features for predicting exon skipping, it is desirable to be able to predict alternative splicing using only information from a single genome.
we show that our approach can still predict exon skipping without using conservation-based features.
methods
datasets and genome browser
we used the dataset of  <cit> , henceforth called dataset d <dig>  consisting of  <dig> alternative and  <dig>  constitutive exons, kindly provided by gideon dror.
in this dataset, constitutive exons are supported by at least four ests each in human and mouse with no est evidence for exon skipping, whereas alternative exons are skipped in both species.
the second dataset is the acescan training set  <cit> , henceforth called dataset d <dig>  which comprises  <dig>  constitutive and  <dig> alternative exons.
for validation purposes we use the genome builds hg <dig> for human and mm <dig> for mouse of the ucsc genome browser  <cit>  were used.
features for machine learning
in total, we used  <dig> features in this study .
thereof,  <dig> were previously used by  <cit> :  exon length,  symmetry, that is, divisibility of exon length by 3;  percent identity of the alignment between the exon and its mouse ortholog;  length of and percent identity of the best local alignment between the up- and downstream  <dig> nt intronic flanks and their mouse orthologs, which are in total four features;  trimer counts for the exon and the  <dig> nt flanking intronic regions, which are a total of  <dig> ×  <dig> =  <dig> features;  intensity of the poly-pyrimidine tract  as the number of pyrimidines in the window - <dig> to - <dig> from the 3'ss; and  nucleotides at the 5'ss positions - <dig> to - <dig> and + <dig> to + <dig>  which are a total of  <dig> ×  <dig> =  <dig> features.
the features , used for the first time in  <cit> , were kindly provided by gideon dror, along with the dataset d <dig>
note that the total number of features used is  <dig> whereas the sum of the entries here is  <dig>  because some features have been counted in more than one category .
we also used six features from  <cit> : the percent identity of the global alignments between up- and downstream  <dig> nt intronic flanks and their mouse orthologs, lengths of the upstream and downstream flanking introns, and the strength of the 3' and 5' splice sites  computed by maxentscan  <cit> .
we used the programs "needle" and "water" from the emboss software suite  <cit>  for aligning the exons and the intronic flanks with their mouse orthologs and computing the conservation based features.
among the new features we added were dinucleotide counts for the exon and the  <dig> nt intronic flanks, a total of  <dig> ×  <dig> =  <dig> features.
as it has been shown that exon skipping is more prevalent in regions of low gc content  <cit> , we used the gc content of the exon and the intronic flanks as three additional features.
to use features based on eses and esss, we applied neighbourhood inference  scores  <cit> .
briefly, each hexamer has an ni score between - <dig> and  <dig>  with negative scores indicating a tendency towards acting as an ess, and a positive score, a tendency to act like an ese.
hexamers with a score of  <dig> or - <dig> are considered "trusted" eses and esss, respectively, and those with a score of greater than  <dig>  or smaller than - <dig>  are considered to have "strong" ese or ess activity.
we used the density of ni scores, defined as the number of hexamers with ni scores  <dig>  ≥  <dig> , >  <dig>  <  <dig>  ≤ - <dig> , - <dig>  normalized by the number of hexamers in the exon .
additionally, the distribution of eses and esss may have a bearing on splicing as well.
therefore, we used the variance of ni scores for "trusted" and "strong" eses and esss .
since the density of eses and esss near splice junctions has been suggested to be important in determining splicing outcome  <cit> , we also measured the densities in the first and last  <dig> nucleotides of the exon .
we also designed features using very recently published datasets of conserved isres enriched in the upstream and downstream intron flanks of all exons, as well isres enriched in upstream and downstream introns flanking alternative exons  <cit> .
we used the density of isres from these four lists in both upstream and downstream  <dig> nt flanking intronic regions, giving us eight novel features.
secondary structure can influence alternative splicing  <cit> .
the single-strandedness of ese, ess or isre motifs was characterized using pu  values  <cit> , which represent the probability that all the bases in the given motif are unpaired.
since local rna folding is influenced by the length of the sequence context  <cit> , we minimized potential biases by using  <dig> to  <dig> nt symmetrical context lengths up- and downstream of a given hexamer, and computing the average of the  <dig> pu values thus obtained  <cit> .
we pre-computed pu values in this manner for all the hexamers in the exons, and combined the ni scores with pu values.
various thresholds were used for absolute ni score value  and a pu value of  <dig> .
two kinds of combinations were used:  a "boolean" combination, that is, counting the number of hexamers with ni and pu values both above the thresholds; and  the product of ni and pu values .
similarly, we used pu values in conjunction with isre information to characterize the single-strandedness of intronic splicing regulatory elements .
mutations around the splice junctions can effect splicing.
therefore, we designed a feature to measure how well the immediate neighbourhood of the splice junctions was conserved.
we formed two 12-mers consisting of the bases from positions - <dig> to + <dig> around the 3'ss and the 5'ss.
the number of identical nucleotides between the human and mouse 12-mers result in two new features.
we also used several motifs from a recent study characterizing conserved motifs associated with constitutive and alternative splicing  <cit> .
however, depending on the partition, these features were either not discriminative or weakly so, indicating that they are important only for a small minority of the alternative exons.
to count the number of binding sites for the polypyrimidine-tract-binding protein , a well-studied repressive regulator of alternative splicing  <cit> , we counted the simplest known motifs for its binding sites – ucuu and cucucu, as well as the sum.
the density of ptb binding sites in the  <dig> nt intron flanks and the exon gives nine features.
lastly, we used novel features derived from features already known to be discriminative.
for example, while it is known that skipped exons tend, on the average, to be shorter than constitutive exons, it has been shown that long exons can be skipped if flanked by very long introns  <cit> .
furthermore, it is possible that the shorter the exon is with respect to the flanking introns, the harder it is for the spliceosome to reliably recognize it.
consequently, we used the ratio of upstream and downstream intron length to exon length, as well of the intron lengths, as three features.
we also used the pairwise products of human-mouse identity of the exon and each  <dig> nt intron flank as well as of the exon and both flanks, in order to capture information about simultaneous conservation of the exon and the intronic flanks .
information gain and information gain ratio
to compare the information content of the features, we used information gain, and information gain ratio, which are established measures of the usefulness of features in the field of machine learning  <cit> .
the formula for information gain is:
ig = h - h
where h is the entropy of the class variable, and h is the conditional entropy of the class variable, given the feature.
we used the weka package  <cit>  for computing information gain and information gain ratio, in order to rank the features according to how informative they were.
bayesian networks
we used the algorithms for feature selection, model learning and classification as described in  <cit> , and made available via the public webserver biobayesnet  <cit> .
feature subset selection
given a training set, we selected features in a three step procedure.
first, we use an entropy based method developed by  <cit>  to find partitions of the feature ranges which best separate the given classes .
features for which the entire feature range is partitioned into at least two intervals, such that the distributions of the two classes differ significantly in these intervals, are called "discriminative" and they are the basis of further analysis.
on the other hand, those features for which no such intervals are found are essentially non-informative, or "non-discriminative" features for our purposes.
once the discretization algorithm has chosen the set of discriminative features, an optimal  subset can be selected using the sequential floating feature selection  method  <cit> .
briefly, this algorithm starts with an initially empty feature subset, and at each step, adds the feature which most improves a specific quality measure.
after this addition, all previously added features are deleted from the subset, unless doing so worsens the quality measure.
this is done in order to avoid getting trapped in local minima.
the algorithm stops when neither inserting new features nor deleting existing ones improves the quality measure provided by the subset.
thirdly, one can enforce inclusion or exclusion of any given feature manually.
the manual feature selection consists only of removing a few "weak" features  as they are unlikely to generalize well to unseen data; and addition of a few "strong" features , which were selected by the discretizer but not by the sffs algorithm.
learning the bayesian network
we restrict the structure of the bns by using the so-called tree-augmented naïve bayes  structure  <cit> .
in a naïve bayes classifier/network, the attributes are assumed to be independent, given the class, that is, the node representing the class variable is a parent of all other nodes, and there are no other edges in the network.
a tan classifier augments the underlying naïve bayes classifier by allowing at most one additional parent per node, that is, each node is the child of the class attribute node, and of at most one more node.
we use tan classifiers because while learning the best bn structure, given some training data, is in general an np-hard problem  <cit> , for tan networks there exist efficient structure-learning algorithms that reduce the problem of determining the optimal tree structure to finding a maximum-weighted spanning tree  <cit> .
once the structure of the network has been learned, the  probability distributions over the feature values of each feature  are estimated in a straightforward manner from count statistics derived from learning data.
finally, bayesian inference of marginal probabilities can be approximately calculated by the efficient technique of variable elimination  <cit> .
data partition
given a dataset , we partitioned the data into three equal parts as carried out in  <cit> .
then, in turn, we used two-thirds of the data to train the bns, and the remaining one-third was used for testing.
the test set remained untouched while the training set was used for discretization, feature selection, and learning the bn  <cit> .
finally, the bn which had been learned on the training set was used to classify the samples in the test set.
this procedure was repeated twice for the other two one-thirds, and the average of the three runs was taken as the final performance.
for comparing 2-fold, 3-fold, 5-fold and 10-fold cross-validation, we used weka  <cit> .
RESULTS
improved prediction of conserved cassette exons by bayesian networks
as pointed out by  <cit> , good performance at low false positive rates is especially important for the task of distinguishing alternative exons from constitutive exons on a genome-wide scale, since the latter comprise the majority of exons.
furthermore, a low number of false positives is critical in case of experimental verification of predictions.
to this end, we measure the true positive rate  at false positive rate  of  <dig> %, and call it tp <dig> .
we also compute the receiver operating curve  and measure the area under the roc curve , which is a standard measure of the quality of a classifier  <cit> .
we used the dataset and the cross-validation scheme described in  <cit> .
this dataset contains  <dig> alternative and  <dig>  constitutive exons and is called d <dig> in the following.
the overall performance obtained, using novel features in addition to those described in the literature , was tp <dig>  = 61%, and auc =  <dig>  , compared to tp <dig>  = 50%, and auc =  <dig>  reported in  <cit>  using svms.
this substantial improvement demonstrates that many of the novel features are informative and discriminative for conserved exon skipping events.
feature selection
the number of features studied in machine learning tasks is often very high, and many  of them might be irrelevant, or redundant  <cit> .
therefore, it is customary to preprocess the data in order to select a useful subset of features – this is called "feature selection".
feature selection can be carried out in three stages within the biobayesnet framework  <cit> .
firstly, a "discretizer" applying the algorithm of fayyad and irani  <cit>  discards features for which no suitable discriminative intervals are found.
secondly, the sequential feature subset selection  algorithm  <cit>  can be applied to select a subset of the remaining features.
thirdly, one can enforce inclusion or exclusion of any given feature manually.
the manual feature selection  typically involved the addition and removal of  <dig> or fewer features each, given a feature subset of 20– <dig> features obtained using the two automated approaches.
the performance on d <dig> using only feature selection using the "discretizer" was tp <dig>  = 39%, and auc =  <dig> .
using the sffs algorithm for further feature selection resulted in tp <dig>  = 47%, and auc =  <dig> , whereas the manual inclusion/exclusion of features gave the final performance of tp <dig>  = 61%, and auc =  <dig> .
this illustrates that the overall quality of classification, as measure by the auc, is quite robust, and we get good performance even when only the "discretizer" is used, but the performance at low false positive rates is quite sensitive to small changes in the feature subset, so the other two methods of feature selection result in significant improvement.
we note that manual feature selection is only needed to improve tp <dig>  – if we consider more global measures of classification performance such the auc or balanced sensitivity and specificity, the automated feature selection methods suffice.
using only automated feature selection, we routinely achieve auc values in the  <dig> – <dig>  range, and balanced sensitivity and specificity in the 87%–91% range.
discriminative features
eses and esss are motifs bound by proteins which either enhance or suppress splicing.
it has been shown that alternative and constitutive exons differ in the density of eses and esss  <cit> .
we used neighborhood inference  scores to infer ese and ess activity for all hexamers  <cit> .
we used the density of eses and esss, with various thresholds for the ni scores.
the constitutive exons have a slightly higher density of eses than do alternative exons , as well as esss  this was also confirmed using other ese/ess datasets  <cit>  and is in agreement with previous studies  <cit> .
depending on the split, the density of eses and esss was either not discriminative, or weakly so.
varying the threshold of the ni score did not change this.
on the other hand, some of the novel features using ni scores were discriminative on most splits – for instance, the average of all positive ni scores, as well as the average of all negative ni scores.
similarly, the average of all "strong eses"  and "strong esss"  were discriminative features.
however, the density of eses and esss near the splice sites was not found to be discriminative.
splicing regulatory elements are found in introns as well  <cit> .
consequently, we also designed features using very recently published datasets of conserved intronic splicing regulatory elements   <cit> .
similar to ese and ess based features, we also used the density of isres in the upstream and downstream intronic flanks.
four sets of isres  are given in  <cit> .
we found three of these sets to be discriminative .
for these three discriminative features, alternative exons have a higher density of isres than constitutive exons, which agrees well with the finding that the set of isres has an overlap with esss, and thus many of them may have silencing tendencies  <cit> .
secondary structure can influence alternative splicing  <cit> .
to the best of our knowledge, existing methods to predict alternative splice events do not use secondary structure related properties.
previously, we found that functionally important splicing motifs are preferentially located in single-stranded mrna secondary structures  <cit>  and that ab initio motif finding benefits from taking the single-strandedness of motif occurrences into account  <cit> .
thus, we use features based on a measure of the single-strandedness of ese, ess or isre occurrences, reasoning that single-strandedness of enhancer and silencer occurrences influence the propensity of proteins to bind them.
interestingly, we found that the single-strandedness of ese motifs is informative.
the density of single-stranded eses is higher in constitutive than in alternative exons .
moreover, the information gain of this feature was more than that of the ese density feature .
as single-stranded motifs are expected binding sites for splicing regulatory proteins, this observation adds to previous evidence that mrna secondary structures influence alternative splicing  <cit> .
the density of ptb binding sites in the exon and the upstream  <dig> nt intronic flank were weakly discriminative, indicating that they are important only for a small minority of the alternative exons.
the density of ptb binding sites in the downstream  <dig> nt intronic flank was not discriminative.
the conservation around the splice site, as measured by the number of human-mouse identical positions in a window of  <dig> nt  around the exon boundaries, is a highly discriminative feature, despite other features already capturing both conservation information as well as splice site strength.
it is interesting that while alternative exons have weaker splice sites, they have stronger conservation around the splice junctions.
while only  <dig> % of the constitutive exons have identical matches from positions - <dig> to + <dig> at the 3'ss, the corresponding figure for alternative exons is as high as  <dig> %.
at the 5'ss, the corresponding numbers are  <dig> % and  <dig> %, respectively.
this is consistent with a previous study  <cit> .
the gc content of the upstream intronic flank was found to be a useful discriminative feature, and was lower for flanks of alternative exons than of constitutive exons , in agreement with previous studies  <cit> .
however, neither the gc content of the exon nor of the downstream intronic flank was found to be discriminative.
we tested if the di- and trimer  frequency in exons and intron flanks is different for alternative and constitutive exons.
we found that the frequency of di- and trimers in exons is often much more discriminative than the intronic di- and trimer frequencies.
this suggests that splice regulatory elements governing exon skipping are more common in alternative exons than in introns flanking them.
apart from introducing novel features, we also used features derived from known features.
these combinations were often more informative than the individual features.
for example, the ratios of intron lengths to exon length were more informative than the lengths themselves.
the ratio of the length of the downstream intron to the exon length was an especially useful feature, suggesting that exon skipping may occur when the spliceosome finds it difficult to accurately "spot" an exon upstream of a relatively much longer intron.
splice site strength, first used by  <cit> , was also found be a discriminative feature, with alternative exons having both weaker 3'ss as well as 5'ss than constitutive exons .
most informative features
next we asked which were the most informative features using the information gain, a well established measure in machine learning.
information gain is the reduction in the entropy of the class variable, given the feature.
while information gain is a good measure of the quality of features, it tends to prefer features with a large number of possible values  <cit> .
a measure which avoid this is the information gain ratio, which divides the information by the information of the feature itself, thus penalizing features with a high inherent information.
the top ten features according to the information gain and the information gain ratio criteria are given in table  <dig>  two features, exon identity and length of best alignment in the upstream intron flank, appear in both lists.
a bayesian network with an optimized set of  <dig> features
all three methods of feature selection available in the biobayesnet framework were used to arrive at an optimized subset of  <dig> features.
a performance of tp <dig>  = 61% , and auc =  <dig>   was achieved using the same subset of  <dig> features with each fold.
the bn learned on the entire dataset with the same features, with  <dig> nodes and  <dig> edges, can be seen in fig.
<dig>  we would like to point out some interesting edges in this network which confirm and may extend our biological knowledge of the splicing process:
- "length of the best local alignment of the upstream intron flank and its mouse ortholog"  and "density of intronic splice regulatory elements  enriched in introns flanking as exons, in the upstream intron flank" : since alternative conserved exons  tend to have longer conserved regions and a higher density of isres in their intron flanks, this is a biologically meaningful dependence.
- "ratio of the lengths of the downstream intron and the exon"  and "sum of the maxentscan scores of the 3' and 5' splice sites" : aces tend to have high ratios of intron to exon lengths, and weak splice sites, when compared to constitutive exons  <cit> .
- "density of single-stranded eses"  and "density of tctt in exon" : aces are enriched in exons with multiple occurrences of tctt, which is a binding site of the splice-repressor, ptb  <cit> , and tend to have a lower density of single-stranded eses when compared to constitutive exons.
- "maxentscan score of the 5'ss"  and nodes representing positions in the 5'ss region : the node representing the strength of the 5'ss has an edge to the node representing the binary variable "gat5ssplus4", which indicates whether the nucleotide at the position + <dig> is a g or not, and this node has an edge to the node representing the variable for a t at position + <dig>  both of these nucleotides are different from the splice site consensus at their respective positions, and thus contribute to lowering the splice site score.
furthermore, it is also known that there are dependencies among the 5'ss positions  <cit> .
- "density of intronic splice regulatory elements  enriched in intronic flanks of as exons, in the upstream intron flank"  and "abundance of cgg in the exon" : isres which are enriched in the flanks of alternative exons tend to be cg-rich  <cit> , so the link to cgg motifs in the exon might indicate a subclass of alternative exons found in cg-rich regions.
- "abundance of cca in the exon"  and "average of negative ni scores in exon" : these nodes correspond to features representing the density of the trimer cca in the exon, and the average ni score of all hexamers with negative ni scores, i.e.
ess-like tendencies.
the cca motif is ~35-fold less frequent in esss than eses , so the bn captures the association of cca abundance with the average of negative ni scores.
- "abundance of cgg and gca in the exon"  and "ratio of eses to esss in the exon" : these nodes represent the density of the trimers cgg and gca in the exon, and the ratio of "trusted" eses and esss .
the motif cgg occurs in  <dig> %  of the trusted eses, but in only 1%  of the trusted esss.
similarly, the motif cgg occurs in  <dig> %  of the trusted eses, but in only  <dig> %  of the trusted esss.
thus there is a correlation between the abundance of the motifs cgg and gca in the exon, and the ratio of eses to esss.
some of the other edges can be explained in a trivial manner, for instance those involving the density of overlapping motifs .
we note that one must be careful in interpreting the edges, as not all of them may lend themselves to meaningful biological interpretation.
while not all edges can be interpreted with biological knowledge, they definitely help in our classification since a classifier omitting all edges  performs worse  <cit> .
comparison of 2-fold, 3-fold, 5-fold, and 10-fold cross-validation
we used 3-fold cross-validation in order to compare our results to  <cit> , who did the same.
however, since it is common in machine learning to use 2-fold, 5-fold, 10-fold, or "leave one out"  cross-validation, we compared the performance of these different approaches on the dataset d <dig>  using the weka package  <cit> , and the optimized set of  <dig> features described above.
the results for the 2/3/5/10/loo cross-validations were: tp <dig>  = 57%/60%/57%/58%/59%, and auc =  <dig> / <dig> / <dig> / <dig> / <dig> .
performance using the same features as the svm
to assess the factors behind the improved performance of bns, we used the same  <dig> features as reported in  <cit> , and obtained the same overall quality of prediction  and slightly improved tp <dig>  .
this indicates that accurate classification of conserved exon skipping depends more on the features used rather than the choice of classifier.
performance of bayesian networks on a second dataset
next, we tested our approach on a different dataset of conserved exon skipping events, the acescan training set  <cit>  henceforth called dataset d <dig>  which comprises  <dig>  constitutive and  <dig> alternative exons.
using the basic set of  <dig> features  <cit> , the bn achieved values of tp <dig>  = 52%, and auc =  <dig> .
after incorporating the novel features, and performing feature selection as described above, the best performance achieved on d <dig> was: tp <dig>  = 59%, and auc =  <dig>  .
thus, we achieve a good performance, similar to that on d <dig>  on the dataset d <dig> as well.
however, the number of discriminative features is smaller than for d <dig>  this trend continues with the addition of novel features – of all the  <dig> features, typically 110– <dig> are discriminative on a 2/ <dig> split of d <dig>  whereas only 65– <dig> are discriminative on a 2/ <dig> split of d <dig>  a possible reason for this could be the different criteria used in the construction of the two datasets, resulting in possibly different extents of corruption of the sets of constitutive exons by alternative exons, because the dataset d <dig> requires  <dig> identical ests for an exon to be considered constitutive, whereas the dataset d <dig> does not.
furthermore, d <dig> has more exaggerated differences among the two classes for several features – for example, while 74% of alternative exons preserve the reading frame compared to 37% of the constitutive exons, the corresponding numbers for d <dig> are 67% and 39%.
thus, the subset of conserved exon skipping events in d <dig> seems to be characterized by more strongly discriminative features.
cross-validation by training and testing on the two independently constructed datasets
it is usual in machine learning to divide the available data into training and testing partitions, and optimize the classification using these.
it is then assumed that similar performance can be achieved on other datasets of a similar nature.
however, given that there are often differences in the way independent datasets are prepared by different groups of scientists, it may be optimistic to presume this.
we suggest that testing on an independent dataset is likely to give a better indication of the level of performance that can be expected when scaling to a genome-wide prediction.
to use d <dig> and d <dig> for this purpose, we removed from d <dig> the exons already present in d <dig> – leaving  <dig> alternative and  <dig>  constitutive exons in d <dig>  to minimize any biases introduced by different ratios of the numbers of samples in each class, we then randomly sampled constitutive exons from d <dig> to have the same ratio  of constitutive to alternative exons, leaving  <dig> alternative and  <dig>  constitutive exons in d <dig>  we then used the optimal feature subsets obtained on d <dig> and d <dig> earlier to train bns on the respective entire datasets.
when we used the bn trained on d <dig> to test d <dig>  the performance achieved was tp <dig>  = 27%, and auc =  <dig> .
the corresponding performance achieved with training on d <dig> and testing on d <dig> was tp <dig>  = 26%, and auc =  <dig> .
while an auc value of  <dig>   indicates good overall classification, this is less than the  <dig>  achieved when tested on unseen data from the same source.
the effect on tp <dig>  is quite dramatic.
we think that these figures might be a more accurate estimate of what to expect when a classifier is used to classify independently produced data.
performance will tend to be  slightly worse on independently produced data than on unseen data from the same source, something which is true of all classifiers in general.
assessing over-fitting
to assess whether our increased performance is due to over-fitting, we randomly permutated the labels 'alternative' and 'constitutive' between the data points and trained the bn on the relabelled datasets d <dig> and d <dig>  in case of overfitting, we would still expect a good performance, while the auc value of a random classifier should be close to  <dig>   <cit> .
after relabelling, most features are no longer discriminative.
in fact, only  <dig> features remained discriminative, and these were the same for both datasets – symmetry, and the  <dig> features describing the positional biases in the 5'ss region.
the auc achieved was  <dig>  on dataset d <dig>  and  <dig>  on d <dig>  this shows that our approach has no problems with over-fitting.
to further rule out overfitting, we used a random three way split: 60% of the data for training, 20% for validation and optimization, and 20% for testing.
we obtained tp <dig>  = 63% and auc =  <dig>  on the validation set; using the same set of features, the performance on the test set was tp <dig>  = 59% and auc =  <dig> .
using this "train-validate-test" approach on d <dig>  we obtained tp <dig>  = 58% and auc =  <dig>  on the validation set, and tp <dig>  = 60% and auc =  <dig>  on the test set.
since the performance on both datasets is very similar to the performance achieved using our three-stage feature selection approach, we conclude that the improvement is not mainly due to manual feature selection.
however, manual selection is not ideal, and an automated feature selection algorithm designed to optimize performance in the low false-positive region would be more satisfying.
this is one of the possible future directions of work.
as a first approach to entirely automated feature selection, we performed the following experiment: we randomly chose 75% of d <dig> for training, and 25% for testing.
feature selection was done using only the training part, and the test part was touched only once at the very end of the procedure.
the feature selection was as follows: starting with the full set of features, we iteratively discarded one feature at a time, and performed 10-fold cross-validated classification using a bn  with the remaining features.
features were discarded in order of increasing information gain, that is, the least informative features were discarded first.
we re-inserted a feature only if at least one of tp <dig>  or auc decreased as a result of omitting it.
this was done only in one pass, and features once discarded, were not considered again.
this is clearly not an optimal strategy, and leads to bigger feature subsets than the approach used before, but still yields good results.
using the subset of  <dig> features thus obtained on d <dig> led to performance of tp <dig>  = 54% and auc =  <dig>  on the training set, and tp <dig>  = 57% and auc =  <dig>  on the test set.
on d <dig>  this approach yielded a subset of  <dig> features and a performance of tp <dig>  = 56% and auc =  <dig>  on the training set, and tp <dig>  = 52% and auc =  <dig>  on the test set.
thus, we can also obtain good performance on unseen data using a feature selection strategy which, though suboptimal, is easy to automate.
moreover, we also used the feature sets obtained in the "train-validate-test" setting with a naïve bayes classifier  and obtained tp <dig>  = 47% and auc =  <dig>  for a 10-fold cross validation on d <dig>  and tp <dig>  = 43% and auc =  <dig>  for a 10-fold cross validation on d <dig>  which are both better than the performance using nbc reported in  <cit>  .
compared to the bns, nbcs achieve a higher sensitivity but lower specificity.
this indicates that the novel features help in improving classification performance, and similar improvements should be possible using other classifiers like svms, neural networks and so on.
false positives with high posterior probability are likely true alternative exons
next, we carefully looked at exons that are labelled constitutive but obtained a high posterior probability of being alternative exons from the bn.
since they seemed to be more similar to aces than to other constitutive exons, we hypothesized that newer est/cdna data provides evidence for exon skipping, or any other kind of alternative splicing at these exons.
out of  <dig>  exons in d <dig> that were labelled constitutive,  <dig> were assigned a p – posterior probability of being an ace – of  <dig>  or more.
a detailed inspection using the ucsc genome browser  <cit>  revealed that seven have est and/or mrna evidence of alternative splicing in at least one of human and mouse  and that two of them are alternatively skipped in both species, that is, have evidence of being aces.
of the remaining seven exons, one has evidence of being a cassette exon in orangutan .
the results on d <dig> are even more impressive – there are  <dig> exons labelled constitutive and with p ≥  <dig> , of which  <dig> have evidence of exon skipping or another alternative splicing event .
thus, most fp predictions with high posterior probabilities of being cassette exons in both d <dig> and d <dig> datasets are actually alternative despite being labelled constitutive at the time the datasets were prepared.
this further demonstrates the good performance of the bn.
predicting exon skipping without using conservation based features
while conservation based features have proved to be the most discriminative, it is desirable to be able to predict alternative splice events using only features that are available to the spliceosome.
the performance on this test is also indicative of our understanding of the process of exon skipping.
hence, we should also aim to predict splicing using only information available from a single genome.
we predicted exon skipping omitting all conservation based features – the best performance achieved was tp <dig>  = 29%, and auc =  <dig>  on dataset d <dig> and tp <dig>  = 26%, and auc =  <dig>  on dataset d <dig>
while this performance is noticeably poorer than that achieved using conservation based features, we would like to note that the datasets d <dig> and d <dig> consist of exons that are either constitutively spliced in both human and mouse, or cassette exons in both.
thus, we are still distinguishing only between conserved constitutive splicing and conserved exon skipping, leaving out cases of species-specific splicing, as well as of alternative splicing of species-specific exons, which form the majority of alternative exons  <cit> .
CONCLUSIONS
using bayesian networks  and several novel features that emerged from recent studies of alternative splicing, we have achieved considerably improved classification of conserved cassette exons.
we were able to improve the performance described in  <cit>  due to the incorporation of novel features.
to the best of our knowledge, this is the first time that features involving secondary structure and intronic splice regulatory elements have been employed for distinguishing alternative exons from constitutive ones.
we also compared our performance on two datasets, and showed that the bn is able to produce accurate classification on both.
however, it is worth noting that these datasets differ with respect to discriminative properties.
one direction of future work would be to consolidate various datasets of constitutive and alternative exons, and compile sets of features, which are discriminative over each of them, and the intersection of these sets, which is discriminative over all datasets.
another interesting line to pursue is to predict other kinds of alternative splicing.
here, we focused on exon skipping, which is the most prevalent form of alternative splicing in human and higher vertebrates.
however, other major forms of alternative splicing such as alternative donor and acceptor sites  <cit>  are also of biological importance, and it would be worthwhile to develop similarly accurate classifiers for these events.
ideally, we should be able to predict splicing outcomes without conservation based information, as the information required by the spliceosome is present in the given genome.
we report our performance at this task, while it is a promising beginning, clearly there is much work to be done.
it should be noted that we have ignored two prominent subclasses of alternative exons – namely orthologous exons which are alternatively spliced in a species-specific manner, and species-specific exons which are alternatively spliced.
both these classes are potentially quite important: as up to 50% of all human alternative exons may be human-specific, and up to 60% of all conserved exons which are alternatively spliced may be alternative in a species-specific manner  <cit> .
classifiers for these tasks are yet to be developed.
authors' contributions
rs designed new features, performed feature extraction, used biobayesnet to perform the classification, did the analyses using weka, analyzed results, and drafted the manuscript.
mh helped with feature design, extraction of secondary structure related features, analysis of the results, and manuscript development.
rp implemented the biobayesnet webserver, which was the framework for using bns, and helped by providing various options for feature discretization, feature selection, learning bns and evaluating performance.
ug and mp participated in the analysis of the results, and helped draft the manuscript.
rb participated in feature design and helped draft the manuscript.
rb and mp supervised the study.
all authors have read and approved the final manuscript.
supplementary material
additional file 1
the top false positives in the datasets d <dig> and d <dig>  results of investigating the top false positives  =  <dig> ) in the datasets d <dig> and d <dig>  using the ucsc genome browser .
click here for file
acknowledgements
we thank gideon dror for providing the dataset d <dig> as well as a few related files, and helpful discussions.
we would like to acknowledge juan valcarcel for suggesting the use of ptb binding sites as a feature, klaus huse for a critical reading of the manuscript, karol szafranski and stefanie schindler for stimulating discussions, and swetlana nikolaweja for help regarding the biobayesnet webserver.
this work was supported by grants from the deutsche forschungsgemeinschaft  and the german ministry of education and research .

 eukaryotic primary mrnas consist of exons and introns. the mature transcript as the substrate for translation is produced by removing introns in a process called splicing. splicing can be either constitutive, always producing the same mrna, or alternative, by skipping of variable parts of the primary transcript. alternative splicing is a mechanism for producing transcript and protein diversity. it is particularly widespread in higher eukaryotes, especially mammals. various studies have estimated that up to 74% of all human genes are alternatively spliced. large scale detection of alternative splicing is usually done using expressed sequence tags or microarrays. since alternative splicing can be highly specific for tissues or developmental stages, these methods can only detect splice events that occur in the underlying probe samples with sufficient frequencies and/or are limited to by the microarray design. furthermore, nowadays whole genome shotgun sequencing projects are churning out genomic data at a higher rate than corresponding transcriptome data – the number of ests in genbank release had increased by 19% in one year, compared to a gain of 39% in the number of contigs in the wgs genbank division. thus it can be expected that in the foreseeable future, we shall have several genomes without the level of corresponding extensive transcript coverage required to reveal the extent of alternative splicing, and hence transcriptomic and proteomic variability. accordingly, there is a need for in silico methods of detecting alternative splicing. moreover, such methods can provide further insights into the mechanisms of alternative splicing. exon skipping, whereby a given exon in its entirety is either included in, or excluded from the mature transcript, is the most prevalent form of alternative splicing in humans. it has been shown that sequence-based features, derived from the exon and its flanking introns, can be used to predict skipping of exons that are conserved between human and mouse and alternatively spliced in both species; denoted conserved exon skipping events. previous studies have used such features with state-of-the-art classifiers such as support vector machines and regularized least-squares classifier, and achieved success in predicting exon skipping. other approaches use protein domain information and evolutionary conservation to detect alternative splice events. here, we use bayesian networks, a state-of-the-art machine learning method, to predict conserved exon skipping events. bns are an increasingly popular machine learning approach to data modeling and classification. the ability of bns to cope with features of various value ranges and to learn dependencies between features makes them especially versatile and suited to a large variety of applications. bns allow multiple dependencies between variables, impose no fixed ordering of variables, allow integration of arbitrary features, and the network structure can be automatically learned. this makes bns a flexible choice for biological sequence data analysis. we introduce several novel features that distinguish alternative exons from constitutive exons, including features based on the single-strandedness of exonic splicing enhancers and silencers, and features involving intronic splicing regulatory elements. by validating our classifiers on various datasets, we identify features which are discriminative irrespective of dataset-specific biases, and provide independent measures of the predictive power of the bns. even though conservation based features have proved to be among the most discriminative features for predicting exon skipping, it is desirable to be able to predict alternative splicing using only information from a single genome. we show that our approach can still predict exon skipping without using conservation-based features. methods datasets and genome browser we used the dataset of, henceforth called dataset d consisting of alternative and constitutive exons, kindly provided by gideon dror. in this dataset, constitutive exons are supported by at least four ests each in human and mouse with no est evidence for exon skipping, whereas alternative exons are skipped in both species. the second dataset is the acescan training set, henceforth called dataset d which comprises constitutive and alternative exons. for validation purposes we use the genome builds hg for human and mm for mouse of the ucsc genome browser were used. features for machine learning in total, we used features in this study. thereof, were previously used by: exon length, symmetry, that is, divisibility of exon length by 3; percent identity of the alignment between the exon and its mouse ortholog; length of and percent identity of the best local alignment between the up- and downstream nt intronic flanks and their mouse orthologs, which are in total four features; trimer counts for the exon and the nt flanking intronic regions, which are a total of × = features; intensity of the poly-pyrimidine tract as the number of pyrimidines in the window - to - from the 3'ss; and nucleotides at the 5'ss positions - to - and + to + which are a total of × = features. the features, used for the first time in, were kindly provided by gideon dror, along with the dataset d note that the total number of features used is whereas the sum of the entries here is because some features have been counted in more than one category. we also used six features from: the percent identity of the global alignments between up- and downstream nt intronic flanks and their mouse orthologs, lengths of the upstream and downstream flanking introns, and the strength of the 3' and 5' splice sites computed by maxentscan. we used the programs "needle" and "water" from the emboss software suite for aligning the exons and the intronic flanks with their mouse orthologs and computing the conservation based features. among the new features we added were dinucleotide counts for the exon and the nt intronic flanks, a total of × = features. as it has been shown that exon skipping is more prevalent in regions of low gc content, we used the gc content of the exon and the intronic flanks as three additional features. to use features based on eses and esss, we applied neighbourhood inference scores. briefly, each hexamer has an ni score between - and with negative scores indicating a tendency towards acting as an ess, and a positive score, a tendency to act like an ese. hexamers with a score of or - are considered "trusted" eses and esss, respectively, and those with a score of greater than or smaller than - are considered to have "strong" ese or ess activity. we used the density of ni scores, defined as the number of hexamers with ni scores ≥, > < ≤ -, - normalized by the number of hexamers in the exon. additionally, the distribution of eses and esss may have a bearing on splicing as well. therefore, we used the variance of ni scores for "trusted" and "strong" eses and esss. since the density of eses and esss near splice junctions has been suggested to be important in determining splicing outcome, we also measured the densities in the first and last nucleotides of the exon. we also designed features using very recently published datasets of conserved isres enriched in the upstream and downstream intron flanks of all exons, as well isres enriched in upstream and downstream introns flanking alternative exons. we used the density of isres from these four lists in both upstream and downstream nt flanking intronic regions, giving us eight novel features. secondary structure can influence alternative splicing. the single-strandedness of ese, ess or isre motifs was characterized using pu values, which represent the probability that all the bases in the given motif are unpaired. since local rna folding is influenced by the length of the sequence context, we minimized potential biases by using to nt symmetrical context lengths up- and downstream of a given hexamer, and computing the average of the pu values thus obtained. we pre-computed pu values in this manner for all the hexamers in the exons, and combined the ni scores with pu values. various thresholds were used for absolute ni score value and a pu value of. two kinds of combinations were used: a "boolean" combination, that is, counting the number of hexamers with ni and pu values both above the thresholds; and the product of ni and pu values. similarly, we used pu values in conjunction with isre information to characterize the single-strandedness of intronic splicing regulatory elements. mutations around the splice junctions can effect splicing. therefore, we designed a feature to measure how well the immediate neighbourhood of the splice junctions was conserved. we formed two 12-mers consisting of the bases from positions - to + around the 3'ss and the 5'ss. the number of identical nucleotides between the human and mouse 12-mers result in two new features. we also used several motifs from a recent study characterizing conserved motifs associated with constitutive and alternative splicing. however, depending on the partition, these features were either not discriminative or weakly so, indicating that they are important only for a small minority of the alternative exons. to count the number of binding sites for the polypyrimidine-tract-binding protein, a well-studied repressive regulator of alternative splicing, we counted the simplest known motifs for its binding sites – ucuu and cucucu, as well as the sum. the density of ptb binding sites in the nt intron flanks and the exon gives nine features. lastly, we used novel features derived from features already known to be discriminative. for example, while it is known that skipped exons tend, on the average, to be shorter than constitutive exons, it has been shown that long exons can be skipped if flanked by very long introns. furthermore, it is possible that the shorter the exon is with respect to the flanking introns, the harder it is for the spliceosome to reliably recognize it. consequently, we used the ratio of upstream and downstream intron length to exon length, as well of the intron lengths, as three features. we also used the pairwise products of human-mouse identity of the exon and each nt intron flank as well as of the exon and both flanks, in order to capture information about simultaneous conservation of the exon and the intronic flanks. information gain and information gain ratio to compare the information content of the features, we used information gain, and information gain ratio, which are established measures of the usefulness of features in the field of machine learning. the formula for information gain is: ig = h - h where h is the entropy of the class variable, and h is the conditional entropy of the class variable, given the feature. we used the weka package for computing information gain and information gain ratio, in order to rank the features according to how informative they were. bayesian networks we used the algorithms for feature selection, model learning and classification as described in, and made available via the public webserver biobayesnet. feature subset selection given a training set, we selected features in a three step procedure. first, we use an entropy based method developed by to find partitions of the feature ranges which best separate the given classes. features for which the entire feature range is partitioned into at least two intervals, such that the distributions of the two classes differ significantly in these intervals, are called "discriminative" and they are the basis of further analysis. on the other hand, those features for which no such intervals are found are essentially non-informative, or "non-discriminative" features for our purposes. once the discretization algorithm has chosen the set of discriminative features, an optimal subset can be selected using the sequential floating feature selection method. briefly, this algorithm starts with an initially empty feature subset, and at each step, adds the feature which most improves a specific quality measure. after this addition, all previously added features are deleted from the subset, unless doing so worsens the quality measure. this is done in order to avoid getting trapped in local minima. the algorithm stops when neither inserting new features nor deleting existing ones improves the quality measure provided by the subset. thirdly, one can enforce inclusion or exclusion of any given feature manually. the manual feature selection consists only of removing a few "weak" features as they are unlikely to generalize well to unseen data; and addition of a few "strong" features, which were selected by the discretizer but not by the sffs algorithm. learning the bayesian network we restrict the structure of the bns by using the so-called tree-augmented naïve bayes structure. in a naïve bayes classifier/network, the attributes are assumed to be independent, given the class, that is, the node representing the class variable is a parent of all other nodes, and there are no other edges in the network. a tan classifier augments the underlying naïve bayes classifier by allowing at most one additional parent per node, that is, each node is the child of the class attribute node, and of at most one more node. we use tan classifiers because while learning the best bn structure, given some training data, is in general an np-hard problem, for tan networks there exist efficient structure-learning algorithms that reduce the problem of determining the optimal tree structure to finding a maximum-weighted spanning tree. once the structure of the network has been learned, the probability distributions over the feature values of each feature are estimated in a straightforward manner from count statistics derived from learning data. finally, bayesian inference of marginal probabilities can be approximately calculated by the efficient technique of variable elimination. data partition given a dataset, we partitioned the data into three equal parts as carried out in. then, in turn, we used two-thirds of the data to train the bns, and the remaining one-third was used for testing. the test set remained untouched while the training set was used for discretization, feature selection, and learning the bn. finally, the bn which had been learned on the training set was used to classify the samples in the test set. this procedure was repeated twice for the other two one-thirds, and the average of the three runs was taken as the final performance. for comparing 2-fold, 3-fold, 5-fold and 10-fold cross-validation, we used weka.  improved prediction of conserved cassette exons by bayesian networks as pointed out by, good performance at low false positive rates is especially important for the task of distinguishing alternative exons from constitutive exons on a genome-wide scale, since the latter comprise the majority of exons. furthermore, a low number of false positives is critical in case of experimental verification of predictions. to this end, we measure the true positive rate at false positive rate of %, and call it tp. we also compute the receiver operating curve and measure the area under the roc curve, which is a standard measure of the quality of a classifier. we used the dataset and the cross-validation scheme described in. this dataset contains alternative and constitutive exons and is called d in the following. the overall performance obtained, using novel features in addition to those described in the literature, was tp = 61%, and auc =, compared to tp = 50%, and auc = reported in using svms. this substantial improvement demonstrates that many of the novel features are informative and discriminative for conserved exon skipping events. feature selection the number of features studied in machine learning tasks is often very high, and many of them might be irrelevant, or redundant. therefore, it is customary to preprocess the data in order to select a useful subset of features – this is called "feature selection". feature selection can be carried out in three stages within the biobayesnet framework. firstly, a "discretizer" applying the algorithm of fayyad and irani discards features for which no suitable discriminative intervals are found. secondly, the sequential feature subset selection algorithm can be applied to select a subset of the remaining features. thirdly, one can enforce inclusion or exclusion of any given feature manually. the manual feature selection typically involved the addition and removal of or fewer features each, given a feature subset of 20– features obtained using the two automated approaches. the performance on d using only feature selection using the "discretizer" was tp = 39%, and auc =. using the sffs algorithm for further feature selection resulted in tp = 47%, and auc =, whereas the manual inclusion/exclusion of features gave the final performance of tp = 61%, and auc =. this illustrates that the overall quality of classification, as measure by the auc, is quite robust, and we get good performance even when only the "discretizer" is used, but the performance at low false positive rates is quite sensitive to small changes in the feature subset, so the other two methods of feature selection result in significant improvement. we note that manual feature selection is only needed to improve tp – if we consider more global measures of classification performance such the auc or balanced sensitivity and specificity, the automated feature selection methods suffice. using only automated feature selection, we routinely achieve auc values in the – range, and balanced sensitivity and specificity in the 87%–91% range. discriminative features eses and esss are motifs bound by proteins which either enhance or suppress splicing. it has been shown that alternative and constitutive exons differ in the density of eses and esss. we used neighborhood inference scores to infer ese and ess activity for all hexamers. we used the density of eses and esss, with various thresholds for the ni scores. the constitutive exons have a slightly higher density of eses than do alternative exons, as well as esss this was also confirmed using other ese/ess datasets and is in agreement with previous studies. depending on the split, the density of eses and esss was either not discriminative, or weakly so. varying the threshold of the ni score did not change this. on the other hand, some of the novel features using ni scores were discriminative on most splits – for instance, the average of all positive ni scores, as well as the average of all negative ni scores. similarly, the average of all "strong eses" and "strong esss" were discriminative features. however, the density of eses and esss near the splice sites was not found to be discriminative. splicing regulatory elements are found in introns as well. consequently, we also designed features using very recently published datasets of conserved intronic splicing regulatory elements. similar to ese and ess based features, we also used the density of isres in the upstream and downstream intronic flanks. four sets of isres are given in. we found three of these sets to be discriminative. for these three discriminative features, alternative exons have a higher density of isres than constitutive exons, which agrees well with the finding that the set of isres has an overlap with esss, and thus many of them may have silencing tendencies. secondary structure can influence alternative splicing. to the best of our knowledge, existing methods to predict alternative splice events do not use secondary structure related properties. previously, we found that functionally important splicing motifs are preferentially located in single-stranded mrna secondary structures and that ab initio motif finding benefits from taking the single-strandedness of motif occurrences into account. thus, we use features based on a measure of the single-strandedness of ese, ess or isre occurrences, reasoning that single-strandedness of enhancer and silencer occurrences influence the propensity of proteins to bind them. interestingly, we found that the single-strandedness of ese motifs is informative. the density of single-stranded eses is higher in constitutive than in alternative exons. moreover, the information gain of this feature was more than that of the ese density feature. as single-stranded motifs are expected binding sites for splicing regulatory proteins, this observation adds to previous evidence that mrna secondary structures influence alternative splicing. the density of ptb binding sites in the exon and the upstream nt intronic flank were weakly discriminative, indicating that they are important only for a small minority of the alternative exons. the density of ptb binding sites in the downstream nt intronic flank was not discriminative. the conservation around the splice site, as measured by the number of human-mouse identical positions in a window of nt around the exon boundaries, is a highly discriminative feature, despite other features already capturing both conservation information as well as splice site strength. it is interesting that while alternative exons have weaker splice sites, they have stronger conservation around the splice junctions. while only % of the constitutive exons have identical matches from positions - to + at the 3'ss, the corresponding figure for alternative exons is as high as %. at the 5'ss, the corresponding numbers are % and %, respectively. this is consistent with a previous study. the gc content of the upstream intronic flank was found to be a useful discriminative feature, and was lower for flanks of alternative exons than of constitutive exons, in agreement with previous studies. however, neither the gc content of the exon nor of the downstream intronic flank was found to be discriminative. we tested if the di- and trimer frequency in exons and intron flanks is different for alternative and constitutive exons. we found that the frequency of di- and trimers in exons is often much more discriminative than the intronic di- and trimer frequencies. this suggests that splice regulatory elements governing exon skipping are more common in alternative exons than in introns flanking them. apart from introducing novel features, we also used features derived from known features. these combinations were often more informative than the individual features. for example, the ratios of intron lengths to exon length were more informative than the lengths themselves. the ratio of the length of the downstream intron to the exon length was an especially useful feature, suggesting that exon skipping may occur when the spliceosome finds it difficult to accurately "spot" an exon upstream of a relatively much longer intron. splice site strength, first used by, was also found be a discriminative feature, with alternative exons having both weaker 3'ss as well as 5'ss than constitutive exons. most informative features next we asked which were the most informative features using the information gain, a well established measure in machine learning. information gain is the reduction in the entropy of the class variable, given the feature. while information gain is a good measure of the quality of features, it tends to prefer features with a large number of possible values. a measure which avoid this is the information gain ratio, which divides the information by the information of the feature itself, thus penalizing features with a high inherent information. the top ten features according to the information gain and the information gain ratio criteria are given in table two features, exon identity and length of best alignment in the upstream intron flank, appear in both lists. a bayesian network with an optimized set of features all three methods of feature selection available in the biobayesnet framework were used to arrive at an optimized subset of features. a performance of tp = 61%, and auc = was achieved using the same subset of features with each fold. the bn learned on the entire dataset with the same features, with nodes and edges, can be seen in fig. we would like to point out some interesting edges in this network which confirm and may extend our biological knowledge of the splicing process: - "length of the best local alignment of the upstream intron flank and its mouse ortholog" and "density of intronic splice regulatory elements enriched in introns flanking as exons, in the upstream intron flank": since alternative conserved exons tend to have longer conserved regions and a higher density of isres in their intron flanks, this is a biologically meaningful dependence. - "ratio of the lengths of the downstream intron and the exon" and "sum of the maxentscan scores of the 3' and 5' splice sites": aces tend to have high ratios of intron to exon lengths, and weak splice sites, when compared to constitutive exons. - "density of single-stranded eses" and "density of tctt in exon": aces are enriched in exons with multiple occurrences of tctt, which is a binding site of the splice-repressor, ptb, and tend to have a lower density of single-stranded eses when compared to constitutive exons. - "maxentscan score of the 5'ss" and nodes representing positions in the 5'ss region: the node representing the strength of the 5'ss has an edge to the node representing the binary variable "gat5ssplus4", which indicates whether the nucleotide at the position + is a g or not, and this node has an edge to the node representing the variable for a t at position + both of these nucleotides are different from the splice site consensus at their respective positions, and thus contribute to lowering the splice site score. furthermore, it is also known that there are dependencies among the 5'ss positions. - "density of intronic splice regulatory elements enriched in intronic flanks of as exons, in the upstream intron flank" and "abundance of cgg in the exon": isres which are enriched in the flanks of alternative exons tend to be cg-rich, so the link to cgg motifs in the exon might indicate a subclass of alternative exons found in cg-rich regions. - "abundance of cca in the exon" and "average of negative ni scores in exon": these nodes correspond to features representing the density of the trimer cca in the exon, and the average ni score of all hexamers with negative ni scores, i.e. ess-like tendencies. the cca motif is ~35-fold less frequent in esss than eses, so the bn captures the association of cca abundance with the average of negative ni scores. - "abundance of cgg and gca in the exon" and "ratio of eses to esss in the exon": these nodes represent the density of the trimers cgg and gca in the exon, and the ratio of "trusted" eses and esss. the motif cgg occurs in % of the trusted eses, but in only 1% of the trusted esss. similarly, the motif cgg occurs in % of the trusted eses, but in only % of the trusted esss. thus there is a correlation between the abundance of the motifs cgg and gca in the exon, and the ratio of eses to esss. some of the other edges can be explained in a trivial manner, for instance those involving the density of overlapping motifs. we note that one must be careful in interpreting the edges, as not all of them may lend themselves to meaningful biological interpretation. while not all edges can be interpreted with biological knowledge, they definitely help in our classification since a classifier omitting all edges performs worse. comparison of 2-fold, 3-fold, 5-fold, and 10-fold cross-validation we used 3-fold cross-validation in order to compare our results to, who did the same. however, since it is common in machine learning to use 2-fold, 5-fold, 10-fold, or "leave one out" cross-validation, we compared the performance of these different approaches on the dataset d using the weka package, and the optimized set of features described above. the results for the 2/3/5/10/loo cross-validations were: tp = 57%/60%/57%/58%/59%, and auc = / / / /. performance using the same features as the svm to assess the factors behind the improved performance of bns, we used the same features as reported in, and obtained the same overall quality of prediction and slightly improved tp. this indicates that accurate classification of conserved exon skipping depends more on the features used rather than the choice of classifier. performance of bayesian networks on a second dataset next, we tested our approach on a different dataset of conserved exon skipping events, the acescan training set henceforth called dataset d which comprises constitutive and alternative exons. using the basic set of features, the bn achieved values of tp = 52%, and auc =. after incorporating the novel features, and performing feature selection as described above, the best performance achieved on d was: tp = 59%, and auc =. thus, we achieve a good performance, similar to that on d on the dataset d as well. however, the number of discriminative features is smaller than for d this trend continues with the addition of novel features – of all the features, typically 110– are discriminative on a 2/ split of d whereas only 65– are discriminative on a 2/ split of d a possible reason for this could be the different criteria used in the construction of the two datasets, resulting in possibly different extents of corruption of the sets of constitutive exons by alternative exons, because the dataset d requires identical ests for an exon to be considered constitutive, whereas the dataset d does not. furthermore, d has more exaggerated differences among the two classes for several features – for example, while 74% of alternative exons preserve the reading frame compared to 37% of the constitutive exons, the corresponding numbers for d are 67% and 39%. thus, the subset of conserved exon skipping events in d seems to be characterized by more strongly discriminative features. cross-validation by training and testing on the two independently constructed datasets it is usual in machine learning to divide the available data into training and testing partitions, and optimize the classification using these. it is then assumed that similar performance can be achieved on other datasets of a similar nature. however, given that there are often differences in the way independent datasets are prepared by different groups of scientists, it may be optimistic to presume this. we suggest that testing on an independent dataset is likely to give a better indication of the level of performance that can be expected when scaling to a genome-wide prediction. to use d and d for this purpose, we removed from d the exons already present in d – leaving alternative and constitutive exons in d to minimize any biases introduced by different ratios of the numbers of samples in each class, we then randomly sampled constitutive exons from d to have the same ratio of constitutive to alternative exons, leaving alternative and constitutive exons in d we then used the optimal feature subsets obtained on d and d earlier to train bns on the respective entire datasets. when we used the bn trained on d to test d the performance achieved was tp = 27%, and auc =. the corresponding performance achieved with training on d and testing on d was tp = 26%, and auc =. while an auc value of indicates good overall classification, this is less than the achieved when tested on unseen data from the same source. the effect on tp is quite dramatic. we think that these figures might be a more accurate estimate of what to expect when a classifier is used to classify independently produced data. performance will tend to be slightly worse on independently produced data than on unseen data from the same source, something which is true of all classifiers in general. assessing over-fitting to assess whether our increased performance is due to over-fitting, we randomly permutated the labels'alternative' and'constitutive' between the data points and trained the bn on the relabelled datasets d and d in case of overfitting, we would still expect a good performance, while the auc value of a random classifier should be close to. after relabelling, most features are no longer discriminative. in fact, only features remained discriminative, and these were the same for both datasets – symmetry, and the features describing the positional biases in the 5'ss region. the auc achieved was on dataset d and on d this shows that our approach has no problems with over-fitting. to further rule out overfitting, we used a random three way split: 60% of the data for training, 20% for validation and optimization, and 20% for testing. we obtained tp = 63% and auc = on the validation set; using the same set of features, the performance on the test set was tp = 59% and auc =. using this "train-validate-test" approach on d we obtained tp = 58% and auc = on the validation set, and tp = 60% and auc = on the test set. since the performance on both datasets is very similar to the performance achieved using our three-stage feature selection approach, we conclude that the improvement is not mainly due to manual feature selection. however, manual selection is not ideal, and an automated feature selection algorithm designed to optimize performance in the low false-positive region would be more satisfying. this is one of the possible future directions of work. as a first approach to entirely automated feature selection, we performed the following experiment: we randomly chose 75% of d for training, and 25% for testing. feature selection was done using only the training part, and the test part was touched only once at the very end of the procedure. the feature selection was as follows: starting with the full set of features, we iteratively discarded one feature at a time, and performed 10-fold cross-validated classification using a bn with the remaining features. features were discarded in order of increasing information gain, that is, the least informative features were discarded first. we re-inserted a feature only if at least one of tp or auc decreased as a result of omitting it. this was done only in one pass, and features once discarded, were not considered again. this is clearly not an optimal strategy, and leads to bigger feature subsets than the approach used before, but still yields good results. using the subset of features thus obtained on d led to performance of tp = 54% and auc = on the training set, and tp = 57% and auc = on the test set. on d this approach yielded a subset of features and a performance of tp = 56% and auc = on the training set, and tp = 52% and auc = on the test set. thus, we can also obtain good performance on unseen data using a feature selection strategy which, though suboptimal, is easy to automate. moreover, we also used the feature sets obtained in the "train-validate-test" setting with a naïve bayes classifier and obtained tp = 47% and auc = for a 10-fold cross validation on d and tp = 43% and auc = for a 10-fold cross validation on d which are both better than the performance using nbc reported in. compared to the bns, nbcs achieve a higher sensitivity but lower specificity. this indicates that the novel features help in improving classification performance, and similar improvements should be possible using other classifiers like svms, neural networks and so on. false positives with high posterior probability are likely true alternative exons next, we carefully looked at exons that are labelled constitutive but obtained a high posterior probability of being alternative exons from the bn. since they seemed to be more similar to aces than to other constitutive exons, we hypothesized that newer est/cdna data provides evidence for exon skipping, or any other kind of alternative splicing at these exons. out of exons in d that were labelled constitutive, were assigned a p – posterior probability of being an ace – of or more. a detailed inspection using the ucsc genome browser revealed that seven have est and/or mrna evidence of alternative splicing in at least one of human and mouse and that two of them are alternatively skipped in both species, that is, have evidence of being aces. of the remaining seven exons, one has evidence of being a cassette exon in orangutan. the results on d are even more impressive – there are exons labelled constitutive and with p ≥, of which have evidence of exon skipping or another alternative splicing event. thus, most fp predictions with high posterior probabilities of being cassette exons in both d and d datasets are actually alternative despite being labelled constitutive at the time the datasets were prepared. this further demonstrates the good performance of the bn. predicting exon skipping without using conservation based features while conservation based features have proved to be the most discriminative, it is desirable to be able to predict alternative splice events using only features that are available to the spliceosome. the performance on this test is also indicative of our understanding of the process of exon skipping. hence, we should also aim to predict splicing using only information available from a single genome. we predicted exon skipping omitting all conservation based features – the best performance achieved was tp = 29%, and auc = on dataset d and tp = 26%, and auc = on dataset d while this performance is noticeably poorer than that achieved using conservation based features, we would like to note that the datasets d and d consist of exons that are either constitutively spliced in both human and mouse, or cassette exons in both. thus, we are still distinguishing only between conserved constitutive splicing and conserved exon skipping, leaving out cases of species-specific splicing, as well as of alternative splicing of species-specific exons, which form the majority of alternative exons.  using bayesian networks and several novel features that emerged from recent studies of alternative splicing, we have achieved considerably improved classification of conserved cassette exons. we were able to improve the performance described in due to the incorporation of novel features. to the best of our knowledge, this is the first time that features involving secondary structure and intronic splice regulatory elements have been employed for distinguishing alternative exons from constitutive ones. we also compared our performance on two datasets, and showed that the bn is able to produce accurate classification on both. however, it is worth noting that these datasets differ with respect to discriminative properties. one direction of future work would be to consolidate various datasets of constitutive and alternative exons, and compile sets of features, which are discriminative over each of them, and the intersection of these sets, which is discriminative over all datasets. another interesting line to pursue is to predict other kinds of alternative splicing. here, we focused on exon skipping, which is the most prevalent form of alternative splicing in human and higher vertebrates. however, other major forms of alternative splicing such as alternative donor and acceptor sites are also of biological importance, and it would be worthwhile to develop similarly accurate classifiers for these events. ideally, we should be able to predict splicing outcomes without conservation based information, as the information required by the spliceosome is present in the given genome. we report our performance at this task, while it is a promising beginning, clearly there is much work to be done. it should be noted that we have ignored two prominent subclasses of alternative exons – namely orthologous exons which are alternatively spliced in a species-specific manner, and species-specific exons which are alternatively spliced. both these classes are potentially quite important: as up to 50% of all human alternative exons may be human-specific, and up to 60% of all conserved exons which are alternatively spliced may be alternative in a species-specific manner. classifiers for these tasks are yet to be developed. authors' contributions rs designed new features, performed feature extraction, used biobayesnet to perform the classification, did the analyses using weka, analyzed results, and drafted the manuscript. mh helped with feature design, extraction of secondary structure related features, analysis of the results, and manuscript development. rp implemented the biobayesnet webserver, which was the framework for using bns, and helped by providing various options for feature discretization, feature selection, learning bns and evaluating performance. ug and mp participated in the analysis of the results, and helped draft the manuscript. rb participated in feature design and helped draft the manuscript. rb and mp supervised the study. all authors have read and approved the final manuscript. supplementary material additional file 1 the top false positives in the datasets d and d results of investigating the top false positives = ) in the datasets d and d using the ucsc genome browser. click here for file acknowledgements we thank gideon dror for providing the dataset d as well as a few related files, and helpful discussions. we would like to acknowledge juan valcarcel for suggesting the use of ptb binding sites as a feature, klaus huse for a critical reading of the manuscript, karol szafranski and stefanie schindler for stimulating discussions, and swetlana nikolaweja for help regarding the biobayesnet webserver. this work was supported by grants from the deutsche forschungsgemeinschaft and the german ministry of education and research.
at a stringent false positive rate of %, our bn achieves an improved true positive rate of 61%, compared to a previously reported 50% on the same dataset using support vector machines. finally, we predict exon skipping without using conservation-based features, and achieve a true positive rate of 29% at a false positive rate of %. incorporating several novel discriminative features such as intronic splicing regulatory elements leads to the improvement. remarkably, we could show that about half of the exons which are labelled constitutive but receive a high probability of being alternative by the bn, are in fact alternative exons according to the latest est data. bns can be used to achieve accurate identification of alternative exons and provide clues about possible dependencies between relevant features. conservation based features continue to be the most informative, and hence distinguishing alternative exons from constitutive ones without using conservation based features remains a challenging problem. here, we show that the use of bayesian networks allows accurate prediction of evolutionary conserved exon skipping events. features related to mrna secondary structure increase the prediction performance, corroborating previous findings that secondary structures are important for exon recognition. currently, large scale detection of alternative splicing using expressed sequence tags or microarrays does not capture all alternative splicing events. moreover, for many species genomic data is being produced at a far greater rate than corresponding transcript data, hence in silico methods of predicting alternative splicing have to be improved.
 alternative splicing is a major contributor to the diversity of eukaryotic transcriptomes and proteomes. currently, large scale detection of alternative splicing using expressed sequence tags or microarrays does not capture all alternative splicing events. moreover, for many species genomic data is being produced at a far greater rate than corresponding transcript data, hence in silico methods of predicting alternative splicing have to be improved.  here, we show that the use of bayesian networks allows accurate prediction of evolutionary conserved exon skipping events. at a stringent false positive rate of %, our bn achieves an improved true positive rate of 61%, compared to a previously reported 50% on the same dataset using support vector machines. incorporating several novel discriminative features such as intronic splicing regulatory elements leads to the improvement. features related to mrna secondary structure increase the prediction performance, corroborating previous findings that secondary structures are important for exon recognition. random labelling tests rule out overfitting. cross-validation on another dataset confirms the increased performance. when using the same dataset and the same set of features, the bn matches the performance of an svm in earlier literature. remarkably, we could show that about half of the exons which are labelled constitutive but receive a high probability of being alternative by the bn, are in fact alternative exons according to the latest est data. finally, we predict exon skipping without using conservation-based features, and achieve a true positive rate of 29% at a false positive rate of %.  bns can be used to achieve accurate identification of alternative exons and provide clues about possible dependencies between relevant features. the near-identical performance of the bn and svm when using the same features shows that good classification depends more on features than on the choice of classifier. conservation based features continue to be the most informative, and hence distinguishing alternative exons from constitutive ones without using conservation based features remains a challenging problem. 