www.gusucode.com > bioinfo 案例源码程序 matlab代码 > bioinfo/EstimateSynonymousAndNonsynonymousSubRatesMLExample.m

    %% Estimate synonymous and nonsynonymous substitution rates between two nucleotide sequences using maximum likelihood method
% This example shows how to estimate synonymous and nonsynonymous
% substitution rates between two nucleotide sequences that are not
% codon-aligned using maximum likelihood method.

% Copyright 2015 The MathWorks, Inc.


%%
% This example uses two nucleotide sequences representing the human HEXA gene
% (accession number: NM_000520) and mouse HEXA gene (accession number:
% AK080777). 

%%
% If you have live internet connection, you can use |getgenbank| function
% to retrieve the sequence information from the NCBI data repository and
% load the data into MATLAB(R).
%
%   humanHEXA = getgenbank('NM_000520');
%   mouseHEXA = getgenbank('AK080777');

%%
% For your convenience, MATLAB provides these two sequences in the
% following mat file. Note that data in public databases are frequently
% updated and curated, and the results in this example may slightly differ
% if you use the latest data.
load hexosaminidase.mat

%%
% Extract the coding regions from the two nucleotide sequences.
humanHEXA_cds = featureparse(humanHEXA,'feature','CDS','Sequence',true);
mouseHEXA_cds = featureparse(mouseHEXA,'feature','CDS','Sequence',true);

%%
% Align the amino acid sequences converted from the nucleotide sequences.
[sc,al] = nwalign(nt2aa(humanHEXA_cds),nt2aa(mouseHEXA_cds),'extendgap',1);


%%
% Use the |seqinsertgaps| function to copy the gaps from the aligned amino
% acid sequences to their corresponding nucleotide sequences, thus
% codon-aligning them.
humanHEXA_aligned = seqinsertgaps(humanHEXA_cds,al(1,:))
mouseHEXA_aligned = seqinsertgaps(mouseHEXA_cds,al(3,:))

%%
% Estimate the synonymous and nonsynonymous substitutions rates of the
% codon-aligned nucleotide sequences and also display the codons considered
% in the computations and their amino acid translations.
[nonsynSubRate,synSubRate] = dndsml(humanHEXA_aligned,mouseHEXA_aligned,'verbose',true)