BACKGROUND
matlab integrates programming, visualization and computation in an easy-to-use environment and is widely used in scientific and engineering studies.
one of the most attractive features of matlab is that the basic data element of the system is a matrix that does not require dimensioning.
this allows users to solve many technical computing problems, especially those with matrix and vector formulations, in a very effective way.
the matlab environment itself offers a comprehensive set of built-in functions and many toolboxes have been developed, and are often freely available, for more specialized needs.
however, to our knowledge, these advantages in the matlab environment have not been fully utilized in the area of molecular biology and evolution.
only a few matlab toolboxes or functions are freely available for data analysis, exploration, and visualization of nucleotide and protein sequences.
mathworks has recently provided a bioinformatics toolbox, however this toolbox has relatively limited functions for molecular evolutionary studies.
mbetoolbox, is presented here to fulfil the most obvious needs in sequence manipulation, genetic distance estimation and phylogeny inference under the matlab environment.
moreover, this toolbox provides an extensible, functional framework to formulate and solve problems in evolutionary data analysis.
it facilitates the rapid construction of both general applications, as well as special-purpose tools for evolutionary biologists, in a fraction of the time it would take to write a program in a scalar, noninteractive language such as c or fortran.
implementation
mbetoolbox is written in the matlab language and has been tested on the windows platform with matlab version  <dig> .
<dig>  the main functions implemented are: sequence manipulation, computation of evolutionary distances derived from nucleotide-, amino acid- or codon-based substitution models, phylogenetic tree construction, sequence statistics and graphics functions to visualize the results of analyses.
although it implements only a small fraction of the multiplicity of existing methods used in molecular evolutionary analyses, interested users can easily extend the toolbox.
input data and formats
mbetoolbox requires a single ascii file containing the nucleotide or amino acid sequence alignment in either phylip  <cit> , clustalw  <cit>  or fasta format.
the toolbox does provide a built-in clustalw  <cit>  interface if an unaligned sequence file is provided.
protein-coding dna sequences can be automatically aligned based on the corresponding protein alignment with the command alignseqfile.
after input, in common with the mathworks bioinformatics toolbox, mbetoolbox represents the alignment as a numeric matrix with every element standing for a nucleic or amino acid character.
nucleotides a, c, g and t are converted to integers  <dig> to  <dig>  and the  <dig> amino acids are converted to integers  <dig> to  <dig>  a header, containing information about the names and type of the sequences as well as the relevant genetic code for protein-coding nucleotides, is attached to the alignment matrix to form a matlab structure.
an example alignment structure, aln, in matlab code follows:
aln =
seqtype: 2
geneticcode: 1
seqnames: { <dig> × n cell}
seq:
where n is the number of sequences and m is the length of the aligned sequences.
the type of sequence is denoted by  <dig>   <dig> or  <dig> for sequences of non-coding nucleotides, protein coding nucleotides and amino acids, respectively.
sequence manipulation and statistics
the alignment structure, aln, can be manipulated using the matlab language.
for example, aln.seq will extract the xth sequence from the alignment, while aln.seq will extract columns i to j from the alignment.
users may easily extract more specific positions by using functions developed in the toolbox, such as extractpos or extractdegeneratesites to obtain the third codon positions or fourfold degenerate sites, respectively.
for each sequence, some basic statistics such as the nucleotide composition  and gc content, can be reported.
other functions include the calculation of the relative synonymous codon usage  and the codon adaptation index , counts of segregating sites, taking the reverse complement or translating a sequence, and determining the sequence complexity.
evolutionary distances
the evolutionary distance is one of the important measures in molecular evolutionary studies.
it is required to measure the diversity among sequences and to infer distance-based phylogenies.
mbetoolbox contains a number of functions to calculate evolutionary distances based on the observed number of differences.
the formulae used in these functions are analytical solutions of a variety of markov substitution models, such as jc <dig>  <cit> , k2p  <cit> , f <dig>  <cit> , hky  <cit>  .
given the stationarity condition, the most general form of markov substitution models is the general time reversible  model  <cit> .
there is no analytical formula to calculate the gtr distance directly.
a general method, described by rodriguez et al.
<cit> , has been implemented here.
in this method a matrix f, where fij denotes the proportion of sites for which sequence  <dig>  has an i and sequence  <dig>  has a j, is formed.
the gtr distance between s <dig> and s <dig> is then given by
where ∏ denotes the diagonal matrix with values of nucleotide equilibrium frequencies on the diagonal, and trdenotes the trace of matrix a. the above formula can be expressed in matlab syntax directly as:
>> d = -trace*f))
mbetoolbox also calculates the gamma distribution distance and the logdet distance  <cit>  .
for alignments of codons, the toolbox provides calculation or estimation of the synonymous  and non-synonymous  substitution rates by the counting method of nei and gojobori  <cit> , the degenerate methods of li, wu and luo  <cit>  and the method of li or pamilo and bianchi  <cit> , as well as the maximum likelihood method through paml  <cit> .
all these methods for calculating ks and ka require that the input sequences are aligned in the appropriate reading frame, which can be performed by the function alignseqfile.
unresolved codon sites will be removed automatically.
in addition, several quantities, including the number of substitutions per site at only synonymous sites, at only non-synonymous sites, at only four-fold-degenerate sites, or at only zero-fold-degenerate sites can be calculated.
the output from these calculations are distance matrices which can be exported into text or excel files, or used directly in further operations.
phylogeny inference
two distance-based tree creation algorithms, unweighted pair group method with arithmetic mean  and neighbor-joining   <cit>  are provided and trees from these methods can be displayed or exported.
maximum parsimony and maximum likelihood algorithms can be applied to nucleotide or amino acid alignments through an interface to the phylip package  <cit> .
as properly implemented maximum likelihood methods are the best vehicles for statistical inference of evolutionary relationships among species from sequence data, several maximum likelihood functions have been explicitly implemented in mbetoolbox.
these functions allow users to incorporate various evolutionary models, estimate parameters and compare different evolutionary trees.
the simplest case of estimation of the evolutionary distance between two sequences, s <dig> and s <dig>  can be considered as the estimation of the branch length  separating ancestor and descendent nodes.
branch lengths, relative to a calibrated molecular clock, can reveal the time interval for this separation.
a continuous time markov process is generally used to model evolution along the branch from s <dig> to s <dig>  a transition rate matrix, q, is used to indicate the rate of changing from one state to another.
for a specified time interval or distance, t, the transition probability matrix is calculated from p = eqt.
if there are n sites, the full likelihood is
in this equation,  and  are the ith bases of sequences  <dig> and  <dig> respectively;  is the expected frequency of base .
in mbetoolbox, to calculate the likelihood, l, at a given time interval  t, we have to specify a substitution model by using an appropriate model defining function, such as modeljc, modelk2p or modelgtr for non-coding nucleotides, modeljtt or modeldayhoff for amino acids, or modelgy <dig> for codons.
these functions return a model structure composed of an instantaneous rate matrix, r, and an equilibrium frequency vector, pi which give q, ).
once the model is specified, the function likelidist can calculate the log likelihood of the alignment of the two sequences, s <dig> and s <dig>  with respect to the time or distance, t, under the substitution model, model.
in most cases we wish to estimate t instead of calculating l as a function of t, so the function optimlikelidist  will search for the t that maximises the likelihood by using the nelder-mead simplex  method, while holding the other parameters in the model at fixed values.
this constraint can be relaxed by allowing every parameter in the model to be estimated by functions, such as optimlikelidistk2p, that can estimate both t and the model's parameters.
figure  illustrates the estimation of the evolutionary distance between two ribonuclease genes through the fixed- and free-parameter k2p models, respectively.
when the k2p model's parameter, kappa, is fixed, the result and trace of the optimisation process is illustrated by the graph of l and t .
when kappa is a free parameter, a surface shows the result and trace of the optimisation process .
when calculating the likelihood of a phylogenetic tree, where s <dig> and s <dig> are two  nodes in a tree joined to an internal  node, sa, we must sum over all possible assignments of nucleotides to sa to get the likelihood of the distance between s <dig> and s <dig>  consequently, the number of possible combinations of nucleotides becomes too large to be enumerated for even moderately sized trees.
the pruning algorithm introduced by felsenstein  <cit>  takes advantage of the tree topology to evaluate the summation in a computationally efficient  manner.
this and a simple and elegant mapping from a 'parentheses' encoding of a tree to the matrix equation for calculating the likelihood of a tree, developed in the matlab software, phyllab  <cit> , have been adopted in likelitree.
combination of functions
basic operations can be combined to give more complicated functions.
a simple combination of the function to extract the fourfold degenerate sites with the function to calculate gc content produces a new function  that determines the gc content at 4-fold degenerate sites .
a subfunction for calculating synonymous and nonsynonymous differences between two codons, getsynnonsyndiff, can be converted into a program for calculating codon volatility  <cit>  with trivial effort.
similarly, karlinsig which returns karlin's genomic signature  for a given sequence can be easily re-formulated to estimate relative di-codon frequencies, which may be a new index of biological signals in a coding sequence.
in addition, the menu-driven user interface, mbegui, is also a good example illustrating the power of combination of basic mbetoolbox functions.
graphics and gui
good visualisation is essential for successful numerical model building.
leveraging the rich graphics functionality of matlab, mbetoolbox provides a number of functions that can be used to create graphic output, such as scatterplots of ks vs ka, plots of the number of transitions and transversions against genetic distance, sliding window analyses on a nucleotide sequence and the z-curve .
a simple menu-driven graphical user interface  has been developed by using guide in matlab.
the top menu includes file, sequences, distances, phylogeny, graph, polymorphism and help submenus .
it aids the usage of the most frequently required functions so that users do not have to run any scripts or functions from the matlab command line in most cases.
RESULTS
vectorization simplifies programming
matlab is a matrix language, which means it is designed for vector and matrix operations.
programming can be simplified and made more efficient by using algorithms that take advantage of vectorization .
the matlab compiler in version  <dig>  will automatically recognize and vectorize loops without recursion.
an example of vectorization is the calculation of z-scores  <cit>  for smith-waterman alignments  <cit>  to give a measure of the significance of an alignment score against a background of scores from randomly generated sequences with the same composition and length.
hence, z-scores are designed to overcome the bias due to the composition of the alignment and are usually calculated by comparing an actual alignment score with the scores obtained on a set of random sequences generated by a monte-carlo process.
the z-score is defined as:
z =  - mean)/standard deviation
where s is the smith-waterman  score between two sequences a and b. the mean and standard deviation are taken from realignments of the permuted sequences.
the algorithm is implemented as follows in matlab with as few as  <dig> lines of code:
function  = zscores
ml = length;
m <dig> = length;
% initialise two vectors holding z-score of
% s1_rep and s2_rep, i.e., replicate samples
% of sequences s <dig> and s <dig>
v_z <dig> = zeros;
v_z <dig> = zeros;
z_raw = smithwaterman;
for ,
s1_rep = s1);
v_z <dig> = smithwaterman ;
s2_rep = s2);
v_z <dig> = smithwaterman;
end
z <dig> = )./std;
z <dig> = )./std;
z = min;
where randperm is a vector function returning a random permutation of the integers from  <dig> to n and smithwaterman performs local alignment by the standard dynamic programming technique.
extensibility
an important distinction between compiled languages with subroutine libraries and interactive environments like matlab is the ease with which problems can be specified and solved in the latter.
moreover, matlab toolboxes are traditionally organised in a less object-oriented mode and, consequently, functions are more independent of each other and easier to combine and extend.
several examples were given in the implementation section.
comparison with other toolboxes
some other toolboxes have been developed in matlab for bioinformatics related analyses.
these include phyllab  <cit>  and matarray  <cit>  as well as the bioinformatics toolbox developed by mathworks.
other examples can be found at the link and file exchange maintained at matlab central  <cit> .
phyllab is a molecular phylogeny toolbox which also provides some functions for sequence and tree input and manipulation.
its main focus is on creating a maximum likelihood tree based on bayesian principles using a markov chain monte carlo method to compute posterior parameter distributions.
matarray is focussed on the analysis of gene expression data from microarrays and provides normalization and clustering functions but does not address molecular evolution.
the bioinformatics toolbox from mathworks provides a range of bioinformatics functions, including some related to molecular evolution.
mbetoolbox provides a much broader range of molecular evolution related functions and phylogenetic methods than either the more specialized phyllab project or the more general bioinformatics toolbox from mathworks.
these extra functions include io in phylip format, statistical and sequence manipulation functions relevant to molecular evolution , evolutionary distance calculation for nucleic and amino acid sequences, phylogeny inference functions and graphic plots relevant to molecular evolution .
as such it makes an important contribution to the bioinformatics analyses that can be performed in the matlab environment.
a novel enhanced window analysis
to test for the selective pressures in the different lineages of a phylogenetic tree, the nonsynonymous to synonymous rate ratio  is normally estimated  <cit> .
values of ka/ks =  <dig>  >  <dig>  or <  <dig> indicate neutrality, positive selection, or purifying selection, respectively.
however, ks and ka are measurements of average synonymous and nonsynonymous substitutions per site along the whole length of the sequences.
average ks and kavalues give neither the pattern of intragenic fluctuation of selective constraints, nor region- or site-specific information.
a sliding window method is usually adopted to examine the intragenic pattern of the substitution rates and to test for the occurrence of significant clusters of variant regions  <cit> .
significant heterogeneity in ks would indicate that the neutral substitution rate varies across the gene, whereas heterogeneity in kamay indicate that selective constraints vary along the gene.
the results and accuracy of sliding window methods, either overlapping or non-overlapping, depend on both the size of the window and the moving distance adopted.
large window lengths may obliterate the details of patterns in ks or ka, whereas small window lengths usually result in larger statistical fluctuations.
hence, the resolution of a sliding window is usually limited.
a mathematical formalism, similar to the z'-curve  <cit> , is introduced here to solve this problem.
consider a subsequence based analysis of ks or ka.
in the n-th step, count the cumulative numbers of ks or ka occurring from the first to the n-th nucleotide position in the gene sequences being inspected.
let  denote either ks or ka and  denote the cumulative  at the n-th sequence position.
is usually an approximately mono-increasing linear function of n. the points , n =  <dig>   <dig>  ..., n are fit by a least square method to a linear function, f() = βn, to give a straight line with β being its slope.
we define
the two-dimensional curve of  gives an alternative representation of the normal sliding window curve.
to compare these two curve representations, the example dataset of suzuki and gojobori  <cit> , which contains the coding regions of two hepatitis c virus strains , was used.
the entire coding sequence is divided into eight regions .
some of the coding regions have been combined as these short orfs are unlikely to yield meaningful ks and ka values.
the reduction of ks in the c, el and ns5b regions, as well as its elevation in ns <dig>  which have been shown in previous studies  <cit> , are not clear in a standard sliding window representation .
in contrast a sharp increase in the  curve , indicates an increase in , while a drop in the curve indicates a decrease in .
this new method has been implemented in the function plotslidingkaks.
since it is derived from the sliding window method, it is called the enhanced sliding window method.
limitations
the current version of this toolbox lacks novel algorithms yet it implements a variety of existing algorithms.
there are some limitations in the practical use of mbetoolbox.
first, though the toolbox provides many methods to infer and handle sequence and evolutionary analyses, the full range of these features can only be accessed through the matlab command line interface, as in the majority of matlab packages.
second, some of the functions cannot handle ambiguous nucleotide or amino acid codes in the sequences.
the future development of mbetoolbox will overcome these present limitations.
CONCLUSIONS
the mbetoolbox project is an ongoing effort to provide an easy-to-use, yet powerful, analysis environment for molecular biology and evolution.
currently, it offers a substantial set of frequently used functions to manipulate sequences, to calculate genetic distances, to infer phylogenetic trees and related analyses.
mbetoolbox is a useful tool which should inspire evolutionary biologists to take advantage of the matlab environment.
availability and requirements
project name: mbetoolbox
project web page:
operating system: windows 95/98/2000/xp
programming language: matlab  <dig>  or higher
other requirements: statistics toolbox
license: gpl
any restrictions on use by non-academics: license needed
authors' contributions
jjc designed and implemented the software and wrote the draft of the manuscript.
dks participated in the design and revised the manuscript.
xx participated in the design and provided suggestions for future development.
kyy supervised and participated in the design of the study.
all authors read and approved the final version of the manuscript.

 matlab integrates programming, visualization and computation in an easy-to-use environment and is widely used in scientific and engineering studies. one of the most attractive features of matlab is that the basic data element of the system is a matrix that does not require dimensioning. this allows users to solve many technical computing problems, especially those with matrix and vector formulations, in a very effective way. the matlab environment itself offers a comprehensive set of built-in functions and many toolboxes have been developed, and are often freely available, for more specialized needs. however, to our knowledge, these advantages in the matlab environment have not been fully utilized in the area of molecular biology and evolution. only a few matlab toolboxes or functions are freely available for data analysis, exploration, and visualization of nucleotide and protein sequences. mathworks has recently provided a bioinformatics toolbox, however this toolbox has relatively limited functions for molecular evolutionary studies. mbetoolbox, is presented here to fulfil the most obvious needs in sequence manipulation, genetic distance estimation and phylogeny inference under the matlab environment. moreover, this toolbox provides an extensible, functional framework to formulate and solve problems in evolutionary data analysis. it facilitates the rapid construction of both general applications, as well as special-purpose tools for evolutionary biologists, in a fraction of the time it would take to write a program in a scalar, noninteractive language such as c or fortran. implementation mbetoolbox is written in the matlab language and has been tested on the windows platform with matlab version. the main functions implemented are: sequence manipulation, computation of evolutionary distances derived from nucleotide-, amino acid- or codon-based substitution models, phylogenetic tree construction, sequence statistics and graphics functions to visualize the results of analyses. although it implements only a small fraction of the multiplicity of existing methods used in molecular evolutionary analyses, interested users can easily extend the toolbox. input data and formats mbetoolbox requires a single ascii file containing the nucleotide or amino acid sequence alignment in either phylip, clustalw or fasta format. the toolbox does provide a built-in clustalw interface if an unaligned sequence file is provided. protein-coding dna sequences can be automatically aligned based on the corresponding protein alignment with the command alignseqfile. after input, in common with the mathworks bioinformatics toolbox, mbetoolbox represents the alignment as a numeric matrix with every element standing for a nucleic or amino acid character. nucleotides a, c, g and t are converted to integers to and the amino acids are converted to integers to a header, containing information about the names and type of the sequences as well as the relevant genetic code for protein-coding nucleotides, is attached to the alignment matrix to form a matlab structure. an example alignment structure, aln, in matlab code follows: aln = seqtype: 2 geneticcode: 1 seqnames: { × n cell} seq: where n is the number of sequences and m is the length of the aligned sequences. the type of sequence is denoted by or for sequences of non-coding nucleotides, protein coding nucleotides and amino acids, respectively. sequence manipulation and statistics the alignment structure, aln, can be manipulated using the matlab language. for example, aln.seq will extract the xth sequence from the alignment, while aln.seq will extract columns i to j from the alignment. users may easily extract more specific positions by using functions developed in the toolbox, such as extractpos or extractdegeneratesites to obtain the third codon positions or fourfold degenerate sites, respectively. for each sequence, some basic statistics such as the nucleotide composition and gc content, can be reported. other functions include the calculation of the relative synonymous codon usage and the codon adaptation index, counts of segregating sites, taking the reverse complement or translating a sequence, and determining the sequence complexity. evolutionary distances the evolutionary distance is one of the important measures in molecular evolutionary studies. it is required to measure the diversity among sequences and to infer distance-based phylogenies. mbetoolbox contains a number of functions to calculate evolutionary distances based on the observed number of differences. the formulae used in these functions are analytical solutions of a variety of markov substitution models, such as jc, k2p, f, hky. given the stationarity condition, the most general form of markov substitution models is the general time reversible model. there is no analytical formula to calculate the gtr distance directly. a general method, described by rodriguez et al., has been implemented here. in this method a matrix f, where fij denotes the proportion of sites for which sequence has an i and sequence has a j, is formed. the gtr distance between s and s is then given by where ∏ denotes the diagonal matrix with values of nucleotide equilibrium frequencies on the diagonal, and trdenotes the trace of matrix a. the above formula can be expressed in matlab syntax directly as: >> d = -trace*f)) mbetoolbox also calculates the gamma distribution distance and the logdet distance. for alignments of codons, the toolbox provides calculation or estimation of the synonymous and non-synonymous substitution rates by the counting method of nei and gojobori, the degenerate methods of li, wu and luo and the method of li or pamilo and bianchi, as well as the maximum likelihood method through paml. all these methods for calculating ks and ka require that the input sequences are aligned in the appropriate reading frame, which can be performed by the function alignseqfile. unresolved codon sites will be removed automatically. in addition, several quantities, including the number of substitutions per site at only synonymous sites, at only non-synonymous sites, at only four-fold-degenerate sites, or at only zero-fold-degenerate sites can be calculated. the output from these calculations are distance matrices which can be exported into text or excel files, or used directly in further operations. phylogeny inference two distance-based tree creation algorithms, unweighted pair group method with arithmetic mean and neighbor-joining are provided and trees from these methods can be displayed or exported. maximum parsimony and maximum likelihood algorithms can be applied to nucleotide or amino acid alignments through an interface to the phylip package. as properly implemented maximum likelihood methods are the best vehicles for statistical inference of evolutionary relationships among species from sequence data, several maximum likelihood functions have been explicitly implemented in mbetoolbox. these functions allow users to incorporate various evolutionary models, estimate parameters and compare different evolutionary trees. the simplest case of estimation of the evolutionary distance between two sequences, s and s can be considered as the estimation of the branch length separating ancestor and descendent nodes. branch lengths, relative to a calibrated molecular clock, can reveal the time interval for this separation. a continuous time markov process is generally used to model evolution along the branch from s to s a transition rate matrix, q, is used to indicate the rate of changing from one state to another. for a specified time interval or distance, t, the transition probability matrix is calculated from p = eqt. if there are n sites, the full likelihood is in this equation, and are the ith bases of sequences and respectively; is the expected frequency of base. in mbetoolbox, to calculate the likelihood, l, at a given time interval t, we have to specify a substitution model by using an appropriate model defining function, such as modeljc, modelk2p or modelgtr for non-coding nucleotides, modeljtt or modeldayhoff for amino acids, or modelgy for codons. these functions return a model structure composed of an instantaneous rate matrix, r, and an equilibrium frequency vector, pi which give q, ). once the model is specified, the function likelidist can calculate the log likelihood of the alignment of the two sequences, s and s with respect to the time or distance, t, under the substitution model, model. in most cases we wish to estimate t instead of calculating l as a function of t, so the function optimlikelidist will search for the t that maximises the likelihood by using the nelder-mead simplex method, while holding the other parameters in the model at fixed values. this constraint can be relaxed by allowing every parameter in the model to be estimated by functions, such as optimlikelidistk2p, that can estimate both t and the model's parameters. figure illustrates the estimation of the evolutionary distance between two ribonuclease genes through the fixed- and free-parameter k2p models, respectively. when the k2p model's parameter, kappa, is fixed, the result and trace of the optimisation process is illustrated by the graph of l and t. when kappa is a free parameter, a surface shows the result and trace of the optimisation process. when calculating the likelihood of a phylogenetic tree, where s and s are two nodes in a tree joined to an internal node, sa, we must sum over all possible assignments of nucleotides to sa to get the likelihood of the distance between s and s consequently, the number of possible combinations of nucleotides becomes too large to be enumerated for even moderately sized trees. the pruning algorithm introduced by felsenstein takes advantage of the tree topology to evaluate the summation in a computationally efficient manner. this and a simple and elegant mapping from a'parentheses' encoding of a tree to the matrix equation for calculating the likelihood of a tree, developed in the matlab software, phyllab, have been adopted in likelitree. combination of functions basic operations can be combined to give more complicated functions. a simple combination of the function to extract the fourfold degenerate sites with the function to calculate gc content produces a new function that determines the gc content at 4-fold degenerate sites. a subfunction for calculating synonymous and nonsynonymous differences between two codons, getsynnonsyndiff, can be converted into a program for calculating codon volatility with trivial effort. similarly, karlinsig which returns karlin's genomic signature for a given sequence can be easily re-formulated to estimate relative di-codon frequencies, which may be a new index of biological signals in a coding sequence. in addition, the menu-driven user interface, mbegui, is also a good example illustrating the power of combination of basic mbetoolbox functions. graphics and gui good visualisation is essential for successful numerical model building. leveraging the rich graphics functionality of matlab, mbetoolbox provides a number of functions that can be used to create graphic output, such as scatterplots of ks vs ka, plots of the number of transitions and transversions against genetic distance, sliding window analyses on a nucleotide sequence and the z-curve. a simple menu-driven graphical user interface has been developed by using guide in matlab. the top menu includes file, sequences, distances, phylogeny, graph, polymorphism and help submenus. it aids the usage of the most frequently required functions so that users do not have to run any scripts or functions from the matlab command line in most cases.  vectorization simplifies programming matlab is a matrix language, which means it is designed for vector and matrix operations. programming can be simplified and made more efficient by using algorithms that take advantage of vectorization. the matlab compiler in version will automatically recognize and vectorize loops without recursion. an example of vectorization is the calculation of z-scores for smith-waterman alignments to give a measure of the significance of an alignment score against a background of scores from randomly generated sequences with the same composition and length. hence, z-scores are designed to overcome the bias due to the composition of the alignment and are usually calculated by comparing an actual alignment score with the scores obtained on a set of random sequences generated by a monte-carlo process. the z-score is defined as: z = - mean)/standard deviation where s is the smith-waterman score between two sequences a and b. the mean and standard deviation are taken from realignments of the permuted sequences. the algorithm is implemented as follows in matlab with as few as lines of code: function = zscores ml = length; m = length; % initialise two vectors holding z-score of % s1_rep and s2_rep, i.e., replicate samples % of sequences s and s v_z = zeros; v_z = zeros; z_raw = smithwaterman; for,    s1_rep = s1);    v_z = smithwaterman;    s2_rep = s2);    v_z = smithwaterman; end z = )./std; z = )./std; z = min; where randperm is a vector function returning a random permutation of the integers from to n and smithwaterman performs local alignment by the standard dynamic programming technique. extensibility an important distinction between compiled languages with subroutine libraries and interactive environments like matlab is the ease with which problems can be specified and solved in the latter. moreover, matlab toolboxes are traditionally organised in a less object-oriented mode and, consequently, functions are more independent of each other and easier to combine and extend. several examples were given in the implementation section. comparison with other toolboxes some other toolboxes have been developed in matlab for bioinformatics related analyses. these include phyllab and matarray as well as the bioinformatics toolbox developed by mathworks. other examples can be found at the link and file exchange maintained at matlab central. phyllab is a molecular phylogeny toolbox which also provides some functions for sequence and tree input and manipulation. its main focus is on creating a maximum likelihood tree based on bayesian principles using a markov chain monte carlo method to compute posterior parameter distributions. matarray is focussed on the analysis of gene expression data from microarrays and provides normalization and clustering functions but does not address molecular evolution. the bioinformatics toolbox from mathworks provides a range of bioinformatics functions, including some related to molecular evolution. mbetoolbox provides a much broader range of molecular evolution related functions and phylogenetic methods than either the more specialized phyllab project or the more general bioinformatics toolbox from mathworks. these extra functions include io in phylip format, statistical and sequence manipulation functions relevant to molecular evolution, evolutionary distance calculation for nucleic and amino acid sequences, phylogeny inference functions and graphic plots relevant to molecular evolution. as such it makes an important contribution to the bioinformatics analyses that can be performed in the matlab environment. a novel enhanced window analysis to test for the selective pressures in the different lineages of a phylogenetic tree, the nonsynonymous to synonymous rate ratio is normally estimated. values of ka/ks = > or < indicate neutrality, positive selection, or purifying selection, respectively. however, ks and ka are measurements of average synonymous and nonsynonymous substitutions per site along the whole length of the sequences. average ks and kavalues give neither the pattern of intragenic fluctuation of selective constraints, nor region- or site-specific information. a sliding window method is usually adopted to examine the intragenic pattern of the substitution rates and to test for the occurrence of significant clusters of variant regions. significant heterogeneity in ks would indicate that the neutral substitution rate varies across the gene, whereas heterogeneity in kamay indicate that selective constraints vary along the gene. the results and accuracy of sliding window methods, either overlapping or non-overlapping, depend on both the size of the window and the moving distance adopted. large window lengths may obliterate the details of patterns in ks or ka, whereas small window lengths usually result in larger statistical fluctuations. hence, the resolution of a sliding window is usually limited. a mathematical formalism, similar to the z'-curve, is introduced here to solve this problem. consider a subsequence based analysis of ks or ka. in the n-th step, count the cumulative numbers of ks or ka occurring from the first to the n-th nucleotide position in the gene sequences being inspected. let denote either ks or ka and denote the cumulative at the n-th sequence position. is usually an approximately mono-increasing linear function of n. the points, n =..., n are fit by a least square method to a linear function, f() = βn, to give a straight line with β being its slope. we define the two-dimensional curve of gives an alternative representation of the normal sliding window curve. to compare these two curve representations, the example dataset of suzuki and gojobori, which contains the coding regions of two hepatitis c virus strains, was used. the entire coding sequence is divided into eight regions. some of the coding regions have been combined as these short orfs are unlikely to yield meaningful ks and ka values. the reduction of ks in the c, el and ns5b regions, as well as its elevation in ns which have been shown in previous studies, are not clear in a standard sliding window representation. in contrast a sharp increase in the curve, indicates an increase in, while a drop in the curve indicates a decrease in. this new method has been implemented in the function plotslidingkaks. since it is derived from the sliding window method, it is called the enhanced sliding window method. limitations the current version of this toolbox lacks novel algorithms yet it implements a variety of existing algorithms. there are some limitations in the practical use of mbetoolbox. first, though the toolbox provides many methods to infer and handle sequence and evolutionary analyses, the full range of these features can only be accessed through the matlab command line interface, as in the majority of matlab packages. second, some of the functions cannot handle ambiguous nucleotide or amino acid codes in the sequences. the future development of mbetoolbox will overcome these present limitations.  the mbetoolbox project is an ongoing effort to provide an easy-to-use, yet powerful, analysis environment for molecular biology and evolution. currently, it offers a substantial set of frequently used functions to manipulate sequences, to calculate genetic distances, to infer phylogenetic trees and related analyses. mbetoolbox is a useful tool which should inspire evolutionary biologists to take advantage of the matlab environment. availability and requirements project name: mbetoolbox project web page: operating system: windows 95/98/2000/xp programming language: matlab or higher other requirements: statistics toolbox license: gpl any restrictions on use by non-academics: license needed authors' contributions jjc designed and implemented the software and wrote the draft of the manuscript. dks participated in the design and revised the manuscript. xx participated in the design and provided suggestions for future development. kyy supervised and participated in the design of the study. all authors read and approved the final version of the manuscript. 
it can be used to manipulate aligned sequences, calculate evolutionary distances, estimate synonymous and nonsynonymous substitution rates, and infer phylogenetic trees. however, few functions are freely available in matlab to perform the sequence data analyses specifically required for molecular biology and evolution. we have developed a matlab toolbox, called mbetoolbox, aimed at filling this gap by offering efficient implementations of the most needed functions in molecular biology and evolution. mbetoolbox is a useful tool that can aid in the exploration, interpretation and visualization of data in molecular biology and evolution. moreover, it provides an extensible, functional framework for users with more specialized requirements to explore and analyze aligned nucleotide or protein sequences from an evolutionary perspective. matlab is a high-performance language for technical computing, integrating computation, visualization, and programming in an easy-to-use environment.
 matlab is a high-performance language for technical computing, integrating computation, visualization, and programming in an easy-to-use environment. it has been widely used in many areas, such as mathematics and computation, algorithm development, data acquisition, modeling, simulation, and scientific and engineering graphics. however, few functions are freely available in matlab to perform the sequence data analyses specifically required for molecular biology and evolution.  we have developed a matlab toolbox, called mbetoolbox, aimed at filling this gap by offering efficient implementations of the most needed functions in molecular biology and evolution. it can be used to manipulate aligned sequences, calculate evolutionary distances, estimate synonymous and nonsynonymous substitution rates, and infer phylogenetic trees. moreover, it provides an extensible, functional framework for users with more specialized requirements to explore and analyze aligned nucleotide or protein sequences from an evolutionary perspective. the full functions in the toolbox are accessible through the command-line for seasoned matlab users. a graphical user interface, that may be especially useful for non-specialist end users, is also provided.  mbetoolbox is a useful tool that can aid in the exploration, interpretation and visualization of data in molecular biology and evolution. the software is publicly available at and. 