How to convert MATLAB code to C++ code

10 次查看（过去 30 天）
Israa Amjad
Israa Amjad 2015-12-30
回答： Walter Roberson 2015-12-30
seq1='TACGGGTAT';
seq2='GGACGTACG';
[s,a]=nwalign(seq1,seq2,'scoringmatrix','blosum62','gapopen',11,'extendgap',1)
function [score, alignment, startat, matrices] = nwalign(seq1,seq2,varargin)
match = 5;
mismatch= -2;
gab=-6;
isAminoAcid = true;
scale=1;
% If the input is a structure then extract the Sequence data.
if isstruct(seq1)
seq1 = seqfromstruct(seq1);
end
if isstruct(seq2)
seq2 = seqfromstruct(seq2);
end
if nargin > 2
if rem(nargin,2) == 1
error('Bioinfo:IncorrectNumberOfArguments',...
'Incorrect number of arguments to %s.',mfilename);
end
okargs = {'scoringmatrix','gapopen','extendgap','alphabet','scale','showscore','glocal'};
for j=1:2:nargin-2
pname = varargin{j};
pval = varargin{j+1};
k = find(strncmpi(pname, okargs,numel(pname)));
if isempty(k)
error('Bioinfo:UnknownParameterName',...
'Unknown parameter name: %s.',pname);
elseif length(k)>1
error('Bioinfo:AmbiguousParameterName',...
'Ambiguous parameter name: %s.',pname);
else
switch(k)
case 1 % scoring matrix
if isnumeric(pval)
ScoringMatrix = pval;
else
if ischar(pval)
pval = lower(pval);
end
try
[ScoringMatrix,ScoringMatrixInfo] = feval(pval);
catch allExceptions
error('Bioinfo:InvalidScoringMatrix','Invalid scoring matrix.');
end
end
case 2 %gap open penalty
gapopen = -pval;
case 3 %gap extend penalty
gapextend = -pval;
setGapExtend = true;
case 4 %if sequence is nucleotide
isAminoAcid = bioinfoprivate.optAlphabet(pval,okargs{k}, mfilename);
case 5 % scale
scale=pval;
case 6 % showscore
showscore = bioinfoprivate.opttf(pval, okargs{k}, mfilename);
case 7 % glocal
glocal = bioinfoprivate.opttf(pval, okargs{k}, mfilename);
end
end
end
end
% setting the default scoring matrix
if ~exist('ScoringMatrix','var')
if isAminoAcid
[ScoringMatrix,ScoringMatrixInfo] = blosum50;
else
[ScoringMatrix,ScoringMatrixInfo] = nuc44;
end
end
% getting the scale from ScoringMatrixInfo, if it exists
if exist('ScoringMatrixInfo','var') && isfield(ScoringMatrixInfo,'Scale')
scale=scale*ScoringMatrixInfo.Scale;
end
% handle properly "?" characters typically found in pdb files
if isAminoAcid
if ischar(seq1)
seq1 = strrep(seq1,'?','X');
else
seq1(seq1 == 26) = 23;
end
if ischar(seq2)
seq2 = strrep(seq2,'?','X');
else
seq2(seq2 == 26) = 23;
end
end
% check input sequences
if isAminoAcid && ~(isaa(seq1) && isaa(seq2))
error('Bioinfo:InvalidAminoAcidSequences',...
'Both sequences must be amino acids, use ALPHABET = ''NT'' for aligning nucleotides.');
elseif ~isAminoAcid && ~(isnt(seq1) && isnt(seq2))
error('Bioinfo:InvalidNucleotideSequences',...
'When ALPHABET = ''NT'', both sequences must be nucleotides.');
end
% use numerical arrays for easy indexing
if ischar(seq1)
seq1=upper(seq1); %the output alignment will be all uppercase
if isAminoAcid
intseq1 = aa2int(seq1);
else
intseq1 = nt2int(seq1);
end
else
intseq1 = uint8(seq1);
if isAminoAcid
seq1 = int2aa(intseq1);
else
seq1 = int2nt(intseq1);
end
end
if ischar(seq2)
seq2 = upper(seq2); %the output alignment will be all uppercase
if isAminoAcid
intseq2 = aa2int(seq2);
else
intseq2 = nt2int(seq2);
end
else
intseq2 = uint8(seq2);
if isAminoAcid
seq2 = int2aa(intseq2);
else
seq2 = int2nt(intseq2);
end
end
m = length(seq1);
n = length(seq2);
if ~n||~m
error('Bioinfo:InvalidLengthSequences','Length of input sequences must be greater than 0');
end
scoringMatrixSize = size(ScoringMatrix,1);
highestVal = max([intseq1, intseq2]);
if highestVal > scoringMatrixSize
% if the matrix contains the 'Any' we map to that
if isAminoAcid
anyVal = aa2int('X');
else
anyVal = nt2int('N');
end
if scoringMatrixSize >= anyVal
intseq1(intseq1>scoringMatrixSize) = anyVal;
intseq2(intseq2>scoringMatrixSize) = anyVal;
else
error('Bioinfo:InvalidSymbolsInInputSequences',...
'Sequences contain symbols that cannot be handled by the given scoring matrix.');
end
end
if glocal
algorithm = 3;
else
algorithm = 1;
end
if setGapExtend
% flip order of input sequences for consistency with older versions
if showscore % return the score matrices
[score, path(:,2), path(:,1), F(:,:,1), F(:,:,2), F(:,:,3)] = ...
affinegapmex(intseq2, intseq1, gapopen, gapextend, ScoringMatrix, algorithm);
elseif nargout == 4 % return score matrices and pointer matrices
[score, path(:,2), path(:,1), F(:,:,1), F(:,:,2), F(:,:,3), pointer] = ...
affinegapmex(intseq2, intseq1, gapopen, gapextend, ScoringMatrix, algorithm);
pointer = shiftdim(pointer,1); % for backward compatibility
else % return only score and alignment
[score, path(:,2), path(:,1)] = affinegapmex(intseq2, intseq1, ...
gapopen, gapextend, ScoringMatrix, algorithm);
end
else
% flip order of input sequences for consistency with older versions
if showscore % return the score matrices
[score, path(:,2), path(:,1), F] = ...
simplegapmex(intseq2, intseq1, gapopen, ScoringMatrix, algorithm);
elseif nargout == 4
[score, path(:,2), path(:,1), F, pointer] = ...
simplegapmex(intseq2, intseq1, gapopen, ScoringMatrix, algorithm);
else
[score, path(:,2), path(:,1)] = ...
simplegapmex(intseq2, intseq1, gapopen, ScoringMatrix, algorithm);
end
end
path = path(sum(path,2)>0,:);
path = flipud(path);
% re-scaling the output score
score = scale * score;
if nargout<=1 && ~showscore
return
end
% setting the size of the alignment
alignment = repmat(('- -')',1,size(path,1));
% adding sequence to alignment
alignment(1,path(:,1)>0) = seq1;
alignment(3,path(:,2)>0) = seq2;
% find locations where there are no gaps
h=find(all(path>0,2));
if isAminoAcid
noGaps1=aa2int(alignment(1,h));
noGaps2=aa2int(alignment(3,h));
else
noGaps1=nt2int(alignment(1,h));
noGaps2=nt2int(alignment(3,h));
end
% erasing symbols that cannot be scored
htodel=max([noGaps1;noGaps2])>scoringMatrixSize;
h(htodel)=[];
noGaps1(htodel)=[];
noGaps2(htodel)=[];
% score pairs with no gap
value = ScoringMatrix(sub2ind(size(ScoringMatrix),double(noGaps1),double(noGaps2)));
% insert symbols of the match string into the alignment
alignment(2,h(value>=0)) = ':';
alignment(2,h(noGaps1==noGaps2)) = '|';
startat = [1;1];
% undocumented fourth output -- score and pointer matrices
if nargout > 3
matrices.Scores = F;
matrices.Pointers = pointer;
end
if showscore
figure
F=scale.*max(F(2:end,2:end,:),[],3);
clim=max(max(max(abs(F(~isinf(F))))),eps);
imagesc(F,[-clim clim]);
colormap(privateColorMap(1));
set(colorbar,'YLim',[min([F(:);-eps]) max([F(:);eps])])
title('Scoring Space and Winning Path')
xlabel('Sequence 1')
ylabel('Sequence 2')
hold on
plot(path(all(path>0,2),1),path(all(path>0,2),2),'k.')
end
function pcmap = privateColorMap(selection)
%PRIVATECOLORMAP returns a custom color map
switch selection
case 1, pts = [0 0 .3 20;
0 .1 .8 25;
0 .9 .5 15;
.9 1 .9 8;
1 1 0 26;
1 0 0 26;
.4 0 0 0];
otherwise, pts = [0 0 0 128; 1 1 1 0];
end
xcl=1;
for i=1:size(pts,1)-1
xcl=[xcl,i+1/pts(i,4):1/pts(i,4):i+1]; %#ok<AGROW>
end
pcmap = interp1(pts(:,1:3),xcl);

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2015-12-30
The graphics part of that cannot be converted automatically.
The tool to convert MATLAB code to C++ code is MATLAB Coder, which costs about $US7600 for the commercial version. It is also available for Academic license, but not for Student Version or Home license.

类别

在 Help Center 和 File Exchange 中查找有关 Phylogenetic Analysis 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by