www.gusucode.com > bioinfo 案例源码程序 matlab代码 > bioinfo/EstimateSynonymousAndNonsynonymousSubRatesExample.m
%% Estimate synonymous and nonsynonymous substitution rates between two nucleotide sequences % This example shows how to estimate synonymous and nonsynonymous % substitution rates between two nucleotide sequences that are not % codon-aligned. % 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] = dnds(humanHEXA_aligned,mouseHEXA_aligned,'verbose',true)