BACKGROUND
alternative splicing is widespread during gene expression and has been studied thoroughly in mammals and insect species in recent years.
the proportion of dna that codes for proteins is greatly reduced in eukaryotes, but the number of ways in which the coding blocks called exons are combined to form new proteins is surprisingly large  <cit> .
during transcription, intervening sequences  in pre-mrna molecules are spliced out by the spliceosome machinery and then exons are joined together to form mature mrna molecules.
alternative splicing occurs when the exons of the rna molecules produced by transcription of a gene are reconnected in multiple ways.
little is known about why one exon is chosen over another under specific circumstances.
pairs of exons that are never included together in the final mrna transcript are called mutually exclusive  exons.
during the processing of pre-mrna transcripts, accurate discrimination of exons and introns requires additional regulatory elements in addition to splice site  signals at the 5'- and 3'-ends of exons.
these conserved motifs have been termed exonic and intronic splicing regulatory sequences  because they occur in the exons or introns of a gene sequence and can be subdivided into exonic and intronic splicing enhancers  and exonic and intronic splicing silencers  that activate or repress splicing, respectively  <cit> .
since these regulatory sequences are relatively short, usually 4- <dig> nucleotides, most esr studies have focused on hexamers  <cit> .
a recent study suggested that the sequence specificity of rna-binding proteins that target splicing regulatory sequences is conserved from insects to mammals  <cit> .
several groups have used microarrays in conjunction with manipulation of splicing regulator expression or crosslinking immunoprecipitation  of splicing regulators to identify their indirect or direct targets  <cit> .
such studies provide the most valuable data for dissecting alternative splicing regulation centered on one splicing regulator of interest.
several systematic computational approaches combined with in vivo or in vitro selection methods have been employed to identify motifs in genomic sequences  <cit> .
bioinformatics approaches such as comparative genomics analysis have provided significant sequence and functional insights into the regulatory sequences that occur within exonic regions of a transcript.
in combination with motif analysis, one can further study motif enrichment in a group of tissue-specific alternative exons  <cit> .
me exons have been shown to be expressed in different tissues and cell lines or at different stages during development  <cit> .
a recent study showed that most identified exonic regulatory elements were found to contribute to the alternative splicing between two tissues, while some were important in multiple tissues  <cit> .
a systematic analysis of complete alternative splicing events in a recent study focused on the identification of exon skipping and me splicing events  <cit> .
recently, several alternative splicing databases have been built from abundant expressed sequence tags   <cit> .
one of these databases, the drosophila melanogaster exon database , uses de-bruijn graphs constructed by matching together similar est segments  <cit> .
dedb contains entire gene sequences in splicing graph format, which are used to represent different splicing events  <cit> .
to generate the splicing graphs, est transcripts were clustered on the basis of overlapping genomic positions they occupy.
exons and introns with identical start and end positions were merged into nodes and connections respectively.
the nodes were then linked via connections to form the complete splicing graph for a gene and stored as xml files  <cit> .
methods
we propose a new method, ests-to-esrs , that uses splicing graphs from the drosophila melanogaster exon database  to identify mutually exclusive  exons and subsequently finds putative esr motifs by looking for conserved fragments in a multiple alignment of these me exons.
our focus was on me exons as we were interested in identifying tissue-specific cis-regulatory splicing sequences, which are thought to regulate the inclusion of one exon over another depending on the tissue in which they reside.
we applied e2e to  <dig> species of insects, including  <dig> drosophila species, apis mellifera, tribolium castaneum and anopheles gambiae.
we compared 100% conserved motifs with those obtained using the e2e sliding window spikes method and we compared the insect exonic motifs with those published for mouse and human  <cit> .
finally, we also looked at some mutually exclusive exons in the drosophila melanogaster and compared our findings with those from a recent study  <cit> .
the e2e method uses an in silico analysis pipeline that includes three steps:  identification of mutually exclusive exons in the drosophila splicing graphs,  acquisition and pre-processing of the multiple alignment of mutually exclusive exons and of the phylogenetic tree of  <dig> species of insects, and  a sliding window comparative analysis that identifies putative esrs.
identification of me exons
a recent study  <cit>  has demonstrated the high degree of conservation for splicing regulatory sequences in different species.
this high degree of conservation further supports similar regulatory mechanisms of rna splicing and the biological importance of protein products containing mutually exclusive exons.
we focused on mutually exclusive exons because these have been shown to have evolved from exon duplications in 60% of the cases and are believed to be a major route for generating functional diversity during the evolution of multicellular eukaryotes  <cit> .
exon duplications and therefore me exons represent the origin of alternative splicing and deserve special attention.
in order to extract evolutionarily conserved me exon pairs in insect species, we started with the identification of me exon pairs from the transcriptome of drosophila melanogaster.
the drosophila melanogaster exon database  contains splicing graphs for gene sequences constructed from est transcripts  <cit> .
splicing graphs were downloaded as xml files and used to extract additional information, such as the type of splicing events: intron retention, cassette exon splicing, alternative termination, among others.
we focused on the identification of cassette exons.
a cassette exon event involves an exon that is sometimes included in a transcript and sometimes absent.
in particular, splicing graphs in xml format were used to identify mutually exclusive exons.
the dedb online database visually displays a splicing graph as exon and intron connections for alternative gene transcripts .
each splicing graph represents all observed contiguous mrna transcripts for one single gene.
as a first step in the identification of me exons, splicing graphs containing the phrase "<cassette_exon>" were downloaded from dedb in xml format.
the xml files were parsed and further classified by the flybase gene symbol, node id, start and end position of all exons and introns, and the entire sequences of all exons  and introns  were collected.
the sample splicing graph in figure 1a shows several types of alternative splicing events.
in the graph, individual transcripts are represented as a left to right traversal through a set of connected nodes.
cassette exons involved in mutually exclusive splicing are represented as blue squares .
a series of cassette exons  that overlap with exons  <dig> and  <dig> which are connected to either exon  <dig> or  <dig> are excluded as me exon candidates .
in order to be considered mutually exclusive, the cassette exon must appear in one transcript and be absent in another, and a different exon must be chosen instead.
exons  <dig> and  <dig>  for example, are mutually exclusive exons.
their sequences appear at different non-overlapping positions in the genome and they appear in different transcripts, never together.
the horizontal lines represent connections to another exon corresponding to the dna sequence that separates them in the gene sequence, including introns or other exons not present in a transcript.
identifying me exons requires more than just searching for exons that are not connected.
exons  <dig> and  <dig> in figure 1a are overlapping exons and therefore not mutually exclusive.
they share the same sequence for some of their length, but exon  <dig> starts later than exon  <dig>  exons  <dig> and  <dig> are not mutually exclusive even though they are not connected  because exons  <dig> and  <dig> appear on the same transcript and exon  <dig> overlaps with exon  <dig>
algorithm
the algorithm for identification of me exon pairs uses as input a splicing graph for a gene containing cassette exons.
the terms node  and exon  will be used interchangeably in the following paragraphs.
the graph is stored as a directed acyclic graph  to allow for efficient depth first searches .
figure  <dig> does not show the direction of the edges but there is an implied left to right traversal order for all paths through the graph.
the first step in the algorithm is to construct a list of all cassette exon pairs by identifying the dedb exons with the 'cassette exon' annotation and initialize them as being candidates for me exon pairs.
then one by one the pairs of me exon candidates are excluded from the list if they do not satisfy a set of rules.
a dfs traversal is used to exclude connected exons.
this approach discards exon pairs that could potentially occur in the same transcript, based on the observed est evidence.
for example, exons  <dig> and  <dig> in figure 1b are not considered mutually exclusive because they reside in separate splicing subgraphs connected by regular exons .
in the algorithm described below, each pair will be flagged with a status: "m" for mutually exclusive, "o" for overlap, "c" for connected, "t" for rule three, "f" for rule four.
algorithm for identifying me exons
input: splicing graph s with n cassette exon nodes  from a single gene
output: list of mutually exclusive exon node pairs
<dig>  initialize a list of cassette pairs with
•  "o" if they overlap or "m" as me candidates
<dig>  loop through the cassette exon list
•  run a depth-first search  for each cassette exon a and for each cassette exon b encountered, flag the pair  with "c"
<dig>  for all "m" exon pairs
•  if exon a in a pair  overlaps with an exon c in another pair  which has been marked as "c", mark  as "t"  are not me exons because  <dig> overlaps with  <dig> and  <dig> is connected to 9)
•  if an exon a in a pair  overlaps with an exon c in another pair  marked as "m", then if length<length mark pair  as "f", else mark  as "f"  and  are me exon pairs, but for our analysis purposes because  <dig> and  <dig> overlap we only pick )
retrieval of multiple alignment of me exons and a phylogenetic tree of  <dig> insect species
to obtain the multiple sequence alignment of me exons, the drosophila melanogaster genome sequence was downloaded from the ucsc ftp server: http://hgdownload.cse.ucsc.edu/downloads.html.
the cassette exons of  <dig> genes were mapped onto the downloaded drosophila melanogaster genome using blast.
based on the scaffolds created by the downloaded drosophila melanogaster genome, reading frames were created using the node id, chromosome, start- and end-positions of the cassette nodes in the splicing graphs.
these reading frames created starting points, as well as information regarding the length of the nodes.
ests that were used to generate the splicing graphs were used as guidelines to verify correctness of exon start, end positions and orientation.
manual inspections of the data lead to the removal of incongruent cassette exons.
multiple alignments of the  <dig> cassette exon containing  <dig> different insects, incuding  <dig> drosophila species, apis mellifera, tribolium castaneum and anopheles gambiae, were downloaded from the ucsc ftp server in maf format.
the maf alignments were converted into fasta alignments using a perl script which generates the reverse complement for the exons transcribed from the other strand and also verifies that the downloaded sequences cover the exons in question.
the multiple alignments from ucsc were chosen over other alignment options  to ensure that the sequences used to extract esrs were the homologous sequences of the genes from individual species as there are several repetitive extracellular domains for some genes .
a phylogenetic tree of the  <dig> insect species was also obtained from the ucsc drosophila database.
the ucsc tree was generated with phylofit and phastcons to estimate conserved and non-conserved branch lengths  <cit> .
the ucsc tree was the preferred choice as it represents the consensus based on the whole genome sequences of all  <dig> species and a self-generated tree would have been based on a limited number of sequences.
identification of esrs
the aligned me exon sequences of  <dig> insects were searched for conserved words of lengths 5- <dig> bp as esrs are relatively short and recent studies focused mainly on hexamers  <cit> .
two different approaches were used to identify these conserved words as putative esrs:
<dig>  words that were within 100% conserved regions regardless of flanking sequences
<dig>  words that were significantly more conserved than flanking sequences in phylogenetic clusters using the welch's t test
of special interest were motifs that appeared in one me exon but not the other within a me exon pair.
100% conserved regions
we used a simple sliding window algorithm to scan all columns of an exon multiple alignment with an initial window length of  <dig> bp.
if no gaps were present in that window and the conservation was 100%, the size of the window was increased by one until a non-conserved column was found.
the drosophila melanogaster nucleotide sequences corresponding to the conserved segments in me exons were saved to a file.
in a second step, all words of length n, with n =  <dig> to  <dig>  were extracted from the conserved segments.
the frequency of the words and their location within an exon were subsequently determined.
conservation spikes in phylogenetic clusters
this method implemented a sliding window approach along the multiple alignments of exons.
a conserved segment is not required to be 100% conserved but has to be significantly more conserved than the neighboring segments and above a predefined phylogenetic conservation threshold q+.
an appropriate value for q+ can be found in the evolutionary conservation scores computed from a phylogenetic tree of the  <dig> insect species.
additionally, phylogenetically related insect species were grouped into overlapping clusters.
the tree topology was used to separate insect species into  <dig> phylogenetic clusters as shown in figure  <dig>  q+ is therefore calculated once for each cluster in the phylogenetic tree and compared to the conservation score of hexamers in the sequence alignments.
phylogenetic tree distances were used to calculate the tree conservation scores.
the "branch lengths"  were replaced with tree conservation scores  in the phylogenetic tree.
the tree conservation score q  was defined as the probability that any given base which is not under selection has not mutated in the time separating the ancestor and the descendant.
if the tree specifies the number of synonymous substitutions per site s, then as suggested in  <cit> , the proximity may be reasonably approximated as
q=e-s
let q+ = q + a, be the global conservation score threshold for each cluster, with a = 1-q.
for each window, we compared the conservation score, which is in the interval  <dig> <cit> , with '1' indicating 100% conservation, to the global conservation score q+.
the dro <dig> cluster was given a threshold of  <dig> because of the close evolutionary relatedness of the species in this cluster.
in our calculated conservation scores , a =  <dig> .
q
q+
a putative esr was identified if a particular window showed a spike in the conservation score.
that is, when the window had a conservation score no less than the q+ threshold and was flanked by columns with a lower average conservation score.
the conservation score at position × in the multiple alignment for a window of size w is given by:
c=1w ∑i=1wmis
where mi is the number of matches in column i and s the number of species,  <dig> in our case.
a spike in conservation scores for several clusters and most prominently the dro_anogam cluster can be observed in figure  <dig>  between positions 13- <dig>  between 31- <dig>  64- <dig>  and 79- <dig>  the statistical significance of such spikes was assessed with the welch's t test for each cluster.
below, we present the spike identification rules:
<dig>  for a window's conservation score at position x, c, the following should be true:
• c ≥ q+ for a particular cluster
• c < c > c with w as the window length
<dig>  a spike is deemed significant if px-w <  <dig>  and px+w <  <dig> , with px-w, the left p-value for a t-test comparing c with c and px+w, the right p-value for a t-test comparing c with c
in order to conduct the t-test, the proportions were transformed into arcsine values to comply with the assumption of normality for t-tests as proportions have fixed limits 0- <dig>  the arcsine values and standard deviations were used to calculate the welch's t statistic.
spikes for windows of size  <dig> were identified by using a sliding window approach with a step size of  <dig>
for spikes at the start and end of sequences  we increased the q+ threshold by  <dig> .
the p-value for the difference in spike magnitude between windows of length w starting at positions x <dig> and x <dig> with standard deviations s <dig> and s <dig> was computed using the welch's t-test:
t=c-cs12w+s22w
multiple testing correction
the false discovery rate  multiple comparison correction developed by benjamini and hochberg  <cit>  was applied to the data to control the false positive rate.
unlike earlier analysis, which only considered windows with p-values <  <dig> , all p-values for every comparison at every positions of the alignment across all  <dig> cassette exons were considered for the multiple comparison correction.
adjusted p-values and false discovery rates were computed separately for each left and right t-test window comparison regardless of the phylogenetic cluster.
RESULTS
100% conserved words
when the 100% conservation rule was applied without consideration of surrounding bases,  <dig> words of length 5-  <dig> bp were found in the multiple alignments of  <dig> cassette exons from  <dig> genes of  <dig> insect species.
words of  <dig> bp in length or more were not observed more than  <dig> to  <dig> times.
the most frequent word of length  <dig> bp appeared  <dig> times in different exons: atgga.
atga was a frequent motif occurring  <dig>   <dig>   <dig> times, respectively.
the reverse  of the atgga sequence often appears in the 5' splice site.
however, this particular conserved sequence is not found near the splice site in any of those  <dig> occurrences.
none of the most frequent 100% conserved words matched the words found using the conservation spikes method described below.
spikes in sequence conservation
a list of  <dig> statistically significant "spike" hexamers was found by the e2e method.
about half or  <dig> elements had a left or right p-value <  <dig> .
only  <dig> elements had both left and right side spike magnitude p-values <  <dig> .
we were interested in identifying elements in the last list that appear in one mutually exclusive exon but not in the other, investigating how frequently this occurs and if possible, deducing their impact on tissue-specific splicing.
of these highly significant elements only five genes had me exons where the elements in one me exon differed from those in the other: snf4aγ , doa , zasp <dig> , imp , pdp <dig>  and l82fd .
by applying the fdr multiple comparison correction, we found  <dig> hexamer,  <dig> heptamer and  <dig> octamer motifs with significant left or right corrected p-value <  <dig>  .
a single motif had both a left and right corrected p-value <  <dig>  as the fdr significance threshold amounted to an uncorrected p-value threshold of < <dig> .
the motif was found in gene zasp <dig> .
among the list of fdr significant hexamers,  <dig> genes had significant spikes in one me exon that did not appear in the other.
five of them appeared in the previous list: pdp <dig>  snf4aγ, doa, zasp <dig>  l82fd, and six new genes were added to the list: heph , par- <dig> , bbc , sdt , cg <dig>  and cg <dig>  even though alternatively spliced transcripts have been annotated for all these genes in the flybase, no information on the tissue-specific expression patterns of alternatively spliced transcripts is available.
in order to assess the relevance of the predicted cis-acting elements, we compared the sequence similarities between the predicted elements with predicted binding sites from splicing regulator prediction databases of mouse and human  <cit> .
among the  <dig> statistically significant "spike" hexamers,  <dig> of the predicted elements could be validated with the wang  <dig> database ,  <dig> elements matched with predicted exonic splicing regulators from the ast lab,   <cit> .
only  <dig> elements matched with all three databases, including elements predicted by the burge lab, which were available only within the wang data set.
a majority of these  <dig> elements have left  or right  side spike magnitude p-values <  <dig>  .
four of those appear close to a splice site and four are in the dro_anogam cluster.
<dig> of the  <dig> hexamers in the fdr significant spikes list were also found in the wang list.
sequence logos for the consensus sequences of the significant hexamers were drawn using weblogo.
no discernable frequency pattern was found in the complete list of  <dig> fdr significant hexamers nor among the  <dig> hexamers in common between the ast lab and the fdr spikes list.
surprisingly, hexamers that were also found in the wang list have a sequence logo  that shows a marked preference for the a base at the first two positions.
*dm <dig> stands for drosophila melanogaster, l and r stand for left/right spike increase/decrease respectively.
esrs in the drosophila n-cadherin gene
a recent study looked into the spatiotemporal differential expression of the n-cadherin me exons 7a and 7b, 13a and 13b, 18a and 18b of drosophila and tribolium   <cit> .
this study found that transcripts containing me-13a are expressed only in the cns, while those containing me-13b are only expressed in the early mesoderm.
furthermore, the non-neuronal expression of me-13b drops sharply before synapses begin to form in the embryos.
based on the significant spikes approach, we found a putative esr motif, gatgcag, in 13b for the complete alignment, close to the 3' splice site.
a second esr, aaattg, was found in the dro_apimel cluster, close to the 5' splice site.
no significantly conserved spikes were found in exon 13a.
exons 7a and 7b have less divergent protein sequences and our analysis found the conserved hexamer tgggat in the 7a exon and aaagccag as a significant conservation spike near the 3' splice site for the dro_apimel cluster in 7b.
paralogous alternative exons 18a and 18b exhibit great sequence diversity from each other.
no conservation spikes were detected in 18b.
a conserved word of length  <dig>  tggggcga, appears in exon 18a in the dro_anogam cluster but the l and r p-values are just below  <dig> .
CONCLUSIONS
we presented a bioinformatics workflow protocol to extract exonic splicing regulatory sequences  from mutually exclusive exons of drosophila melanogaster.
we identified mutually exclusive exons in est splicing graphs and used a phylogenetic conservation threshold to identify spikes in conservation through sliding window analysis of exon multiple alignments of  <dig> insect species, as a higher degree of nucleotide sequence conservation is frequently observed in the alternatively spliced exons and/or flanking introns than in the constitutive exons  <cit> .
we looked at sequences with 100% nucleotide conservation, but found no hexamers among the cassette exons that appeared at high frequency and many conserved sequences were found to be part of highly conserved regions within the exons.
incorporating measures of phylogenetic relationships and variability of sequence conservation is a more qualitative approach for the identification of esrs as it is possible to calculate the statistical significance of the identified motifs.
the e2e spikes method found  <dig> putative esrs of which  <dig> elements had a left or right p-value <  <dig>  and  <dig> elements had both left and right side spike magnitude p-values <  <dig> .
the fdr correction process resulted in  <dig> motifs with either a left side or right side corrected p-value <  <dig> , but only  <dig> motif with both left and right corrected p-values <  <dig> .
the reason for this is that t-tests for a sample size of  <dig> will produce higher p-values leading to fewer significant motifs when a correction is applied.
among the  <dig> elements in the uncorrected list, we found five genes and among the  <dig> motifs in the corrected list we found  <dig> genes whose me exons had putative esrs that differed from esrs in the other paired me exons.
it is believed that recognition of esrs in the me exons drives expression of mutually exclusive exons in a tissue-specific manner  <cit> .
unfortunately there is little information in the literature about the tissue-specific expression patterns of alternatively spliced isoforms of drosophila genes, but we found information that linked certain genes to differing spatiotemporal expression.
the gene zasp functions in the formation of integrin adhesion sites and is active during different stages of development.
at late embryonic stages, zasp expression is particularly strong in mesodermal tissues such as visceral, pharyngeal, and somatic muscles  <cit> .
doa is expressed in ectoderm and mesoderm during early embryonic stages and in nervous system during late embryonic stages.
sr and sr-like proteins can be phosphorylated by doa  <cit> .
doa's ability to regulate activities of other splicing factors along with its many spliced isoforms may add more levels of complexity in gene expression regulation.
alternatively spliced sdt transcripts were shown to determine a developmental switch in mrna localization in which apical transcripts were only found during early stages of epithelial development  <cit> .
in the n-cadherin gene, we observed that the dro_anogam cluster which includes dro <dig>  has a spike in 18a but not dro <dig>  itself, which may be due to divergent evolution of the flanking regions around the conserved word in more distant species.
in vivo results from a recent study of the n-cadherin gene, showed no differential subcellular localization between transmembrane domain isoforms containing exons 18a or 18b in drosophila and showed the conserved in situ spatiotemporal expression patterns in drosophila and tribolium  <cit> .
although extensive knowledge of annotated spliced isoforms is available, more information on the in situ isoform-specific expression patterns would be required to validate this approach at a functional level.
a recent study identified and classified hundreds of alternative splicing events that are affected by one splicing regulator and showed that binding sites for the factor are conserved from insects to mammals  <cit> .
another study looked at alternative splicing patterns regulated by four drosophila homologues of the mammalian hnrnp a/b family  <cit> .
both demonstrated conservation of splicing regulatory mechanism between insects and mammals.
thus, we compared esrs predicted by the e2e method with splicing regulatory sequences from three other databases.
<dig> of the  <dig> predicted elements can be validated with the wang  <dig> database of human splicing regulatory sequences and  <dig> elements with the ast lab's list of human and mouse esrs  <cit> .
the regulation of rna splicing is a complicated process as it involves the interplay of many splicing factors and their target sequences.
some exons are constitutively expressed while others are expressed in a tissue- and stage- specific manner.
high-throughput sequencing of transcriptomes has recently demonstrated a high degree of alternative splicing, in up to 92~94% of protein coding genes in the human genome  <cit> .
with a genome-wide approach, e2e method identified several esrs within me exon pairs based on the nucleotide sequence conservation among insect species.
some of these esrs have been previously identified in mammalian splicing databases.
although the tissue-specific expression data of drosophila genes are limited, a few genes containing these identified esrs exhibit differential spatiotemporal expression patterns.
this demonstrates that the e2e method is a powerful tool to help identify conserved splicing signals which might be of high importance in biological functions and shed light on the tissue- and stage-specific splicing regulations.
future work will include more in depth mining of the literature on tissue-specific splicing to validate esrs found with the e2e method.
we will also expand the e2e method to other alternatively spliced exons and introns and apply it to multiple alignments of alternatively spliced genes from human tumor tissue of different patients as it is known that aberrations in alternative splicing occur in many cancers  <cit> .
competing interests
the authors declare that they have no competing interests.
authors' contributions
snh collected and formatted the data and provided the biological background knowledge.
snh and pb designed the study.
pb carried out the statistical analysis and drafted the manuscript.
rl developed the splicing graph parsing algorithm.
jt developed the esrs comparative analysis method.
all authors helped edit the manuscript and approved the final version.
supplementary material
additional file 1
- fdr corrected list of motifs and matches to other human and mouse regulatory element databases.
the table shows the list of motifs with an fdr corrected p-value >  <dig> .
motifs near the splice site were excluded.
a star in l_sig or r_sig indicates a left or right side significant element.
dm stands for drosophila melanogaster, a  <dig> or  <dig> in the wang lab and goren ast lab columns indicates a match with an element in those data sets.
click here for file
acknowledgements
we thank dr. akira chiba for his involvement in the initial phase of this project.
we also thank dr. xin wang's team for providing us with their list of regulatory elements.
this article has been published as part of bmc genomics volume  <dig> supplement  <dig>  2012: selected articles from the first ieee international conference on computational advances in bio and medical sciences : genomics.
the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/13/s <dig>

 alternative splicing is widespread during gene expression and has been studied thoroughly in mammals and insect species in recent years. the proportion of dna that codes for proteins is greatly reduced in eukaryotes, but the number of ways in which the coding blocks called exons are combined to form new proteins is surprisingly large. during transcription, intervening sequences in pre-mrna molecules are spliced out by the spliceosome machinery and then exons are joined together to form mature mrna molecules. alternative splicing occurs when the exons of the rna molecules produced by transcription of a gene are reconnected in multiple ways. little is known about why one exon is chosen over another under specific circumstances. pairs of exons that are never included together in the final mrna transcript are called mutually exclusive exons. during the processing of pre-mrna transcripts, accurate discrimination of exons and introns requires additional regulatory elements in addition to splice site signals at the 5'- and 3'-ends of exons. these conserved motifs have been termed exonic and intronic splicing regulatory sequences because they occur in the exons or introns of a gene sequence and can be subdivided into exonic and intronic splicing enhancers and exonic and intronic splicing silencers that activate or repress splicing, respectively. since these regulatory sequences are relatively short, usually 4- nucleotides, most esr studies have focused on hexamers. a recent study suggested that the sequence specificity of rna-binding proteins that target splicing regulatory sequences is conserved from insects to mammals. several groups have used microarrays in conjunction with manipulation of splicing regulator expression or crosslinking immunoprecipitation of splicing regulators to identify their indirect or direct targets. such studies provide the most valuable data for dissecting alternative splicing regulation centered on one splicing regulator of interest. several systematic computational approaches combined with in vivo or in vitro selection methods have been employed to identify motifs in genomic sequences. bioinformatics approaches such as comparative genomics analysis have provided significant sequence and functional insights into the regulatory sequences that occur within exonic regions of a transcript. in combination with motif analysis, one can further study motif enrichment in a group of tissue-specific alternative exons. me exons have been shown to be expressed in different tissues and cell lines or at different stages during development. a recent study showed that most identified exonic regulatory elements were found to contribute to the alternative splicing between two tissues, while some were important in multiple tissues. a systematic analysis of complete alternative splicing events in a recent study focused on the identification of exon skipping and me splicing events. recently, several alternative splicing databases have been built from abundant expressed sequence tags. one of these databases, the drosophila melanogaster exon database, uses de-bruijn graphs constructed by matching together similar est segments. dedb contains entire gene sequences in splicing graph format, which are used to represent different splicing events. to generate the splicing graphs, est transcripts were clustered on the basis of overlapping genomic positions they occupy. exons and introns with identical start and end positions were merged into nodes and connections respectively. the nodes were then linked via connections to form the complete splicing graph for a gene and stored as xml files. methods we propose a new method, ests-to-esrs, that uses splicing graphs from the drosophila melanogaster exon database to identify mutually exclusive exons and subsequently finds putative esr motifs by looking for conserved fragments in a multiple alignment of these me exons. our focus was on me exons as we were interested in identifying tissue-specific cis-regulatory splicing sequences, which are thought to regulate the inclusion of one exon over another depending on the tissue in which they reside. we applied e2e to species of insects, including drosophila species, apis mellifera, tribolium castaneum and anopheles gambiae. we compared 100% conserved motifs with those obtained using the e2e sliding window spikes method and we compared the insect exonic motifs with those published for mouse and human. finally, we also looked at some mutually exclusive exons in the drosophila melanogaster and compared our findings with those from a recent study. the e2e method uses an in silico analysis pipeline that includes three steps: identification of mutually exclusive exons in the drosophila splicing graphs, acquisition and pre-processing of the multiple alignment of mutually exclusive exons and of the phylogenetic tree of species of insects, and a sliding window comparative analysis that identifies putative esrs. identification of me exons a recent study has demonstrated the high degree of conservation for splicing regulatory sequences in different species. this high degree of conservation further supports similar regulatory mechanisms of rna splicing and the biological importance of protein products containing mutually exclusive exons. we focused on mutually exclusive exons because these have been shown to have evolved from exon duplications in 60% of the cases and are believed to be a major route for generating functional diversity during the evolution of multicellular eukaryotes. exon duplications and therefore me exons represent the origin of alternative splicing and deserve special attention. in order to extract evolutionarily conserved me exon pairs in insect species, we started with the identification of me exon pairs from the transcriptome of drosophila melanogaster. the drosophila melanogaster exon database contains splicing graphs for gene sequences constructed from est transcripts. splicing graphs were downloaded as xml files and used to extract additional information, such as the type of splicing events: intron retention, cassette exon splicing, alternative termination, among others. we focused on the identification of cassette exons. a cassette exon event involves an exon that is sometimes included in a transcript and sometimes absent. in particular, splicing graphs in xml format were used to identify mutually exclusive exons. the dedb online database visually displays a splicing graph as exon and intron connections for alternative gene transcripts. each splicing graph represents all observed contiguous mrna transcripts for one single gene. as a first step in the identification of me exons, splicing graphs containing the phrase "<cassette_exon>" were downloaded from dedb in xml format. the xml files were parsed and further classified by the flybase gene symbol, node id, start and end position of all exons and introns, and the entire sequences of all exons and introns were collected. the sample splicing graph in figure 1a shows several types of alternative splicing events. in the graph, individual transcripts are represented as a left to right traversal through a set of connected nodes. cassette exons involved in mutually exclusive splicing are represented as blue squares. a series of cassette exons that overlap with exons and which are connected to either exon or are excluded as me exon candidates. in order to be considered mutually exclusive, the cassette exon must appear in one transcript and be absent in another, and a different exon must be chosen instead. exons and for example, are mutually exclusive exons. their sequences appear at different non-overlapping positions in the genome and they appear in different transcripts, never together. the horizontal lines represent connections to another exon corresponding to the dna sequence that separates them in the gene sequence, including introns or other exons not present in a transcript. identifying me exons requires more than just searching for exons that are not connected. exons and in figure 1a are overlapping exons and therefore not mutually exclusive. they share the same sequence for some of their length, but exon starts later than exon exons and are not mutually exclusive even though they are not connected because exons and appear on the same transcript and exon overlaps with exon algorithm the algorithm for identification of me exon pairs uses as input a splicing graph for a gene containing cassette exons. the terms node and exon will be used interchangeably in the following paragraphs. the graph is stored as a directed acyclic graph to allow for efficient depth first searches. figure does not show the direction of the edges but there is an implied left to right traversal order for all paths through the graph. the first step in the algorithm is to construct a list of all cassette exon pairs by identifying the dedb exons with the'cassette exon' annotation and initialize them as being candidates for me exon pairs. then one by one the pairs of me exon candidates are excluded from the list if they do not satisfy a set of rules. a dfs traversal is used to exclude connected exons. this approach discards exon pairs that could potentially occur in the same transcript, based on the observed est evidence. for example, exons and in figure 1b are not considered mutually exclusive because they reside in separate splicing subgraphs connected by regular exons. in the algorithm described below, each pair will be flagged with a status: "m" for mutually exclusive, "o" for overlap, "c" for connected, "t" for rule three, "f" for rule four. algorithm for identifying me exons input: splicing graph s with n cassette exon nodes from a single gene output: list of mutually exclusive exon node pairs initialize a list of cassette pairs with • "o" if they overlap or "m" as me candidates loop through the cassette exon list • run a depth-first search for each cassette exon a and for each cassette exon b encountered, flag the pair with "c" for all "m" exon pairs • if exon a in a pair overlaps with an exon c in another pair which has been marked as "c", mark as "t" are not me exons because overlaps with and is connected to 9) • if an exon a in a pair overlaps with an exon c in another pair marked as "m", then if length<length mark pair as "f", else mark as "f" and are me exon pairs, but for our analysis purposes because and overlap we only pick ) retrieval of multiple alignment of me exons and a phylogenetic tree of insect species to obtain the multiple sequence alignment of me exons, the drosophila melanogaster genome sequence was downloaded from the ucsc ftp server: http://hgdownload.cse.ucsc.edu/downloads.html. the cassette exons of genes were mapped onto the downloaded drosophila melanogaster genome using blast. based on the scaffolds created by the downloaded drosophila melanogaster genome, reading frames were created using the node id, chromosome, start- and end-positions of the cassette nodes in the splicing graphs. these reading frames created starting points, as well as information regarding the length of the nodes. ests that were used to generate the splicing graphs were used as guidelines to verify correctness of exon start, end positions and orientation. manual inspections of the data lead to the removal of incongruent cassette exons. multiple alignments of the cassette exon containing different insects, incuding drosophila species, apis mellifera, tribolium castaneum and anopheles gambiae, were downloaded from the ucsc ftp server in maf format. the maf alignments were converted into fasta alignments using a perl script which generates the reverse complement for the exons transcribed from the other strand and also verifies that the downloaded sequences cover the exons in question. the multiple alignments from ucsc were chosen over other alignment options to ensure that the sequences used to extract esrs were the homologous sequences of the genes from individual species as there are several repetitive extracellular domains for some genes. a phylogenetic tree of the insect species was also obtained from the ucsc drosophila database. the ucsc tree was generated with phylofit and phastcons to estimate conserved and non-conserved branch lengths. the ucsc tree was the preferred choice as it represents the consensus based on the whole genome sequences of all species and a self-generated tree would have been based on a limited number of sequences. identification of esrs the aligned me exon sequences of insects were searched for conserved words of lengths 5- bp as esrs are relatively short and recent studies focused mainly on hexamers. two different approaches were used to identify these conserved words as putative esrs: words that were within 100% conserved regions regardless of flanking sequences words that were significantly more conserved than flanking sequences in phylogenetic clusters using the welch's t test of special interest were motifs that appeared in one me exon but not the other within a me exon pair. 100% conserved regions we used a simple sliding window algorithm to scan all columns of an exon multiple alignment with an initial window length of bp. if no gaps were present in that window and the conservation was 100%, the size of the window was increased by one until a non-conserved column was found. the drosophila melanogaster nucleotide sequences corresponding to the conserved segments in me exons were saved to a file. in a second step, all words of length n, with n = to were extracted from the conserved segments. the frequency of the words and their location within an exon were subsequently determined. conservation spikes in phylogenetic clusters this method implemented a sliding window approach along the multiple alignments of exons. a conserved segment is not required to be 100% conserved but has to be significantly more conserved than the neighboring segments and above a predefined phylogenetic conservation threshold q+. an appropriate value for q+ can be found in the evolutionary conservation scores computed from a phylogenetic tree of the insect species. additionally, phylogenetically related insect species were grouped into overlapping clusters. the tree topology was used to separate insect species into phylogenetic clusters as shown in figure q+ is therefore calculated once for each cluster in the phylogenetic tree and compared to the conservation score of hexamers in the sequence alignments. phylogenetic tree distances were used to calculate the tree conservation scores. the "branch lengths" were replaced with tree conservation scores in the phylogenetic tree. the tree conservation score q was defined as the probability that any given base which is not under selection has not mutated in the time separating the ancestor and the descendant. if the tree specifies the number of synonymous substitutions per site s, then as suggested in, the proximity may be reasonably approximated as q=e-s let q+ = q + a, be the global conservation score threshold for each cluster, with a = 1-q. for each window, we compared the conservation score, which is in the interval, with'1' indicating 100% conservation, to the global conservation score q+. the dro cluster was given a threshold of because of the close evolutionary relatedness of the species in this cluster. in our calculated conservation scores, a =. q q+ a putative esr was identified if a particular window showed a spike in the conservation score. that is, when the window had a conservation score no less than the q+ threshold and was flanked by columns with a lower average conservation score. the conservation score at position × in the multiple alignment for a window of size w is given by: c=1w ∑i=1wmis where mi is the number of matches in column i and s the number of species, in our case. a spike in conservation scores for several clusters and most prominently the dro_anogam cluster can be observed in figure between positions 13- between 31- 64- and 79- the statistical significance of such spikes was assessed with the welch's t test for each cluster. below, we present the spike identification rules: for a window's conservation score at position x, c, the following should be true: • c ≥ q+ for a particular cluster • c < c > c with w as the window length a spike is deemed significant if px-w < and px+w <, with px-w, the left p-value for a t-test comparing c with c and px+w, the right p-value for a t-test comparing c with c in order to conduct the t-test, the proportions were transformed into arcsine values to comply with the assumption of normality for t-tests as proportions have fixed limits 0- the arcsine values and standard deviations were used to calculate the welch's t statistic. spikes for windows of size were identified by using a sliding window approach with a step size of for spikes at the start and end of sequences we increased the q+ threshold by. the p-value for the difference in spike magnitude between windows of length w starting at positions x and x with standard deviations s and s was computed using the welch's t-test: t=c-cs12w+s22w multiple testing correction the false discovery rate multiple comparison correction developed by benjamini and hochberg was applied to the data to control the false positive rate. unlike earlier analysis, which only considered windows with p-values <, all p-values for every comparison at every positions of the alignment across all cassette exons were considered for the multiple comparison correction. adjusted p-values and false discovery rates were computed separately for each left and right t-test window comparison regardless of the phylogenetic cluster.  100% conserved words when the 100% conservation rule was applied without consideration of surrounding bases, words of length 5- bp were found in the multiple alignments of cassette exons from genes of insect species. words of bp in length or more were not observed more than to times. the most frequent word of length bp appeared times in different exons: atgga. atga was a frequent motif occurring times, respectively. the reverse of the atgga sequence often appears in the 5' splice site. however, this particular conserved sequence is not found near the splice site in any of those occurrences. none of the most frequent 100% conserved words matched the words found using the conservation spikes method described below. spikes in sequence conservation a list of statistically significant "spike" hexamers was found by the e2e method. about half or elements had a left or right p-value <. only elements had both left and right side spike magnitude p-values <. we were interested in identifying elements in the last list that appear in one mutually exclusive exon but not in the other, investigating how frequently this occurs and if possible, deducing their impact on tissue-specific splicing. of these highly significant elements only five genes had me exons where the elements in one me exon differed from those in the other: snf4aγ, doa, zasp, imp, pdp and l82fd. by applying the fdr multiple comparison correction, we found hexamer, heptamer and octamer motifs with significant left or right corrected p-value <. a single motif had both a left and right corrected p-value < as the fdr significance threshold amounted to an uncorrected p-value threshold of <. the motif was found in gene zasp. among the list of fdr significant hexamers, genes had significant spikes in one me exon that did not appear in the other. five of them appeared in the previous list: pdp snf4aγ, doa, zasp l82fd, and six new genes were added to the list: heph, par-, bbc, sdt, cg and cg even though alternatively spliced transcripts have been annotated for all these genes in the flybase, no information on the tissue-specific expression patterns of alternatively spliced transcripts is available. in order to assess the relevance of the predicted cis-acting elements, we compared the sequence similarities between the predicted elements with predicted binding sites from splicing regulator prediction databases of mouse and human. among the statistically significant "spike" hexamers, of the predicted elements could be validated with the wang database, elements matched with predicted exonic splicing regulators from the ast lab,. only elements matched with all three databases, including elements predicted by the burge lab, which were available only within the wang data set. a majority of these elements have left or right side spike magnitude p-values <. four of those appear close to a splice site and four are in the dro_anogam cluster. of the hexamers in the fdr significant spikes list were also found in the wang list. sequence logos for the consensus sequences of the significant hexamers were drawn using weblogo. no discernable frequency pattern was found in the complete list of fdr significant hexamers nor among the hexamers in common between the ast lab and the fdr spikes list. surprisingly, hexamers that were also found in the wang list have a sequence logo that shows a marked preference for the a base at the first two positions. *dm stands for drosophila melanogaster, l and r stand for left/right spike increase/decrease respectively. esrs in the drosophila n-cadherin gene a recent study looked into the spatiotemporal differential expression of the n-cadherin me exons 7a and 7b, 13a and 13b, 18a and 18b of drosophila and tribolium. this study found that transcripts containing me-13a are expressed only in the cns, while those containing me-13b are only expressed in the early mesoderm. furthermore, the non-neuronal expression of me-13b drops sharply before synapses begin to form in the embryos. based on the significant spikes approach, we found a putative esr motif, gatgcag, in 13b for the complete alignment, close to the 3' splice site. a second esr, aaattg, was found in the dro_apimel cluster, close to the 5' splice site. no significantly conserved spikes were found in exon 13a. exons 7a and 7b have less divergent protein sequences and our analysis found the conserved hexamer tgggat in the 7a exon and aaagccag as a significant conservation spike near the 3' splice site for the dro_apimel cluster in 7b. paralogous alternative exons 18a and 18b exhibit great sequence diversity from each other. no conservation spikes were detected in 18b. a conserved word of length tggggcga, appears in exon 18a in the dro_anogam cluster but the l and r p-values are just below.  we presented a bioinformatics workflow protocol to extract exonic splicing regulatory sequences from mutually exclusive exons of drosophila melanogaster. we identified mutually exclusive exons in est splicing graphs and used a phylogenetic conservation threshold to identify spikes in conservation through sliding window analysis of exon multiple alignments of insect species, as a higher degree of nucleotide sequence conservation is frequently observed in the alternatively spliced exons and/or flanking introns than in the constitutive exons. we looked at sequences with 100% nucleotide conservation, but found no hexamers among the cassette exons that appeared at high frequency and many conserved sequences were found to be part of highly conserved regions within the exons. incorporating measures of phylogenetic relationships and variability of sequence conservation is a more qualitative approach for the identification of esrs as it is possible to calculate the statistical significance of the identified motifs. the e2e spikes method found putative esrs of which elements had a left or right p-value < and elements had both left and right side spike magnitude p-values <. the fdr correction process resulted in motifs with either a left side or right side corrected p-value <, but only motif with both left and right corrected p-values <. the reason for this is that t-tests for a sample size of will produce higher p-values leading to fewer significant motifs when a correction is applied. among the elements in the uncorrected list, we found five genes and among the motifs in the corrected list we found genes whose me exons had putative esrs that differed from esrs in the other paired me exons. it is believed that recognition of esrs in the me exons drives expression of mutually exclusive exons in a tissue-specific manner. unfortunately there is little information in the literature about the tissue-specific expression patterns of alternatively spliced isoforms of drosophila genes, but we found information that linked certain genes to differing spatiotemporal expression. the gene zasp functions in the formation of integrin adhesion sites and is active during different stages of development. at late embryonic stages, zasp expression is particularly strong in mesodermal tissues such as visceral, pharyngeal, and somatic muscles. doa is expressed in ectoderm and mesoderm during early embryonic stages and in nervous system during late embryonic stages. sr and sr-like proteins can be phosphorylated by doa. doa's ability to regulate activities of other splicing factors along with its many spliced isoforms may add more levels of complexity in gene expression regulation. alternatively spliced sdt transcripts were shown to determine a developmental switch in mrna localization in which apical transcripts were only found during early stages of epithelial development. in the n-cadherin gene, we observed that the dro_anogam cluster which includes dro has a spike in 18a but not dro itself, which may be due to divergent evolution of the flanking regions around the conserved word in more distant species. in vivo results from a recent study of the n-cadherin gene, showed no differential subcellular localization between transmembrane domain isoforms containing exons 18a or 18b in drosophila and showed the conserved in situ spatiotemporal expression patterns in drosophila and tribolium. although extensive knowledge of annotated spliced isoforms is available, more information on the in situ isoform-specific expression patterns would be required to validate this approach at a functional level. a recent study identified and classified hundreds of alternative splicing events that are affected by one splicing regulator and showed that binding sites for the factor are conserved from insects to mammals. another study looked at alternative splicing patterns regulated by four drosophila homologues of the mammalian hnrnp a/b family. both demonstrated conservation of splicing regulatory mechanism between insects and mammals. thus, we compared esrs predicted by the e2e method with splicing regulatory sequences from three other databases. of the predicted elements can be validated with the wang database of human splicing regulatory sequences and elements with the ast lab's list of human and mouse esrs. the regulation of rna splicing is a complicated process as it involves the interplay of many splicing factors and their target sequences. some exons are constitutively expressed while others are expressed in a tissue- and stage- specific manner. high-throughput sequencing of transcriptomes has recently demonstrated a high degree of alternative splicing, in up to 92~94% of protein coding genes in the human genome. with a genome-wide approach, e2e method identified several esrs within me exon pairs based on the nucleotide sequence conservation among insect species. some of these esrs have been previously identified in mammalian splicing databases. although the tissue-specific expression data of drosophila genes are limited, a few genes containing these identified esrs exhibit differential spatiotemporal expression patterns. this demonstrates that the e2e method is a powerful tool to help identify conserved splicing signals which might be of high importance in biological functions and shed light on the tissue- and stage-specific splicing regulations. future work will include more in depth mining of the literature on tissue-specific splicing to validate esrs found with the e2e method. we will also expand the e2e method to other alternatively spliced exons and introns and apply it to multiple alignments of alternatively spliced genes from human tumor tissue of different patients as it is known that aberrations in alternative splicing occur in many cancers. competing interests the authors declare that they have no competing interests. authors' contributions snh collected and formatted the data and provided the biological background knowledge. snh and pb designed the study. pb carried out the statistical analysis and drafted the manuscript. rl developed the splicing graph parsing algorithm. jt developed the esrs comparative analysis method. all authors helped edit the manuscript and approved the final version. supplementary material additional file 1 - fdr corrected list of motifs and matches to other human and mouse regulatory element databases. the table shows the list of motifs with an fdr corrected p-value >. motifs near the splice site were excluded. a star in l_sig or r_sig indicates a left or right side significant element. dm stands for drosophila melanogaster, a or in the wang lab and goren ast lab columns indicates a match with an element in those data sets. click here for file acknowledgements we thank dr. akira chiba for his involvement in the initial phase of this project. we also thank dr. xin wang's team for providing us with their list of regulatory elements. this article has been published as part of bmc genomics volume supplement 2012: selected articles from the first ieee international conference on computational advances in bio and medical sciences: genomics. the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/13/s 
analysis of phylogenetic relationships and variability of sequence conservation as implemented in the e2e spikes method may lead to improved identification of esrs. the majority of these genes have been shown to have regulated spatiotemporal expression. a putative esr motif, gatgcag, was identified in the me-13b but not in the me-13a of drosophila n-cadherin, a gene that has been shown to have a distinct spatiotemporal expression pattern of spliced isoforms in a recent study. several drosophila genes with spatiotemporal expression patterns were identified to contain putative esrs located in one exon of the me exon pairs but not in the other. genes were identified with highly significant motifs in one me exon but not in the other, suggesting regulation of me exon splicing through these highly conserved hexamers. we found that approximately half of the putative esrs in common between insects and mammals have a high statistical support. our method, ests-to-esrs, uses graph analysis of est splicing graphs to identify mutually exclusive exons and combines phylogenetic measures, a sliding window approach along the multiple alignment and the welch's t statistic to extract conserved esr motifs. in this paper, we describe an in silico analysis pipeline to extract putative exonic splicing regulatory sequences from a multiple alignment of species of insects. 3- february first ieee international conference on computational advances in bio and medical sciences orlando, fl, usa a recent study suggested that the regulatory codes of individual rna-binding proteins may be nearly immutable between very diverse species such as mammals and insects. during alternative splicing, the inclusion of an exon in the final mrna molecule is determined by nuclear proteins that bind cis-regulatory sequences in a target pre-mrna molecule.
 during alternative splicing, the inclusion of an exon in the final mrna molecule is determined by nuclear proteins that bind cis-regulatory sequences in a target pre-mrna molecule. a recent study suggested that the regulatory codes of individual rna-binding proteins may be nearly immutable between very diverse species such as mammals and insects. the model system drosophila melanogaster therefore presents an excellent opportunity for the study of alternative splicing due to the availability of quality est annotations in flybase. methods in this paper, we describe an in silico analysis pipeline to extract putative exonic splicing regulatory sequences from a multiple alignment of species of insects. our method, ests-to-esrs, uses graph analysis of est splicing graphs to identify mutually exclusive exons and combines phylogenetic measures, a sliding window approach along the multiple alignment and the welch's t statistic to extract conserved esr motifs.  the most frequent 100% conserved word of length bp in different insect exons was "atgga". we identified statistically significant "spike" hexamers, motifs with either a left or right fdr corrected spike magnitude p-value < and with both left and right uncorrected p <. genes were identified with highly significant motifs in one me exon but not in the other, suggesting regulation of me exon splicing through these highly conserved hexamers. the majority of these genes have been shown to have regulated spatiotemporal expression. elements were found to match three mammalian splicing regulator databases. a putative esr motif, gatgcag, was identified in the me-13b but not in the me-13a of drosophila n-cadherin, a gene that has been shown to have a distinct spatiotemporal expression pattern of spliced isoforms in a recent study.  analysis of phylogenetic relationships and variability of sequence conservation as implemented in the e2e spikes method may lead to improved identification of esrs. we found that approximately half of the putative esrs in common between insects and mammals have a high statistical support. several drosophila genes with spatiotemporal expression patterns were identified to contain putative esrs located in one exon of the me exon pairs but not in the other. 3- february first ieee international conference on computational advances in bio and medical sciences orlando, fl, usa 