250 lines
9.8 KiB
Matlab
250 lines
9.8 KiB
Matlab
function t = seqNeighJoin(D, method, names, varargin)
|
|
%SEQNEIGHJOIN neighbor-joining method for phylogenetic tree reconstruction.
|
|
%
|
|
% TREE = SEQNEIGHJOIN(DIST) computes a phylogenetic tree object from the
|
|
% pairwise distances DIST between the species or products applying the
|
|
% neighbor-joining method. The input DIST is a matrix (or vector) such as
|
|
% is generated by SEQPDIST.
|
|
%
|
|
% TREE = SEQNEIGHJOIN(DIST,METHOD) selects the method to compute the
|
|
% distances of the new nodes to all other nodes at every iteration. The
|
|
% general expression to calculate the distances between the new node (n)
|
|
% (after joining i and j) and all others nodes (k) is given by:
|
|
%
|
|
% D(n,k) = a*D(i,k) + (1-a)*D(j,k) - a*D(n,i) - (1-a)*D(n,j)
|
|
%
|
|
% This expression guarantees to find the correct tree with additive data
|
|
% (a.k.a. minimum variance reduction). The options for METHOD are:
|
|
%
|
|
% 'equivar' --- Assumes equal variance and independence of
|
|
% (default) evolutionary distance estimates (a = 1/2). Such as
|
|
% in Studier and Keppler, JMBE (1988).
|
|
% 'firstorder' --- Assumes a first order model of the variances and
|
|
% covariances of evolutionary distance estimates, 'a'
|
|
% is adjusted at every iteration to a value between 0
|
|
% and 1. Such as in Gascuel, JMBE (1997).
|
|
% 'average' --- New distances are the weighted average of previous
|
|
% distances, the branch distances are ignored:
|
|
% D(n,k) = [ D(i,k) + D(j,k) ] /2
|
|
% As in the original neighbor-joining algorithm by
|
|
% Saitou and Nei, JMBE (1987)
|
|
%
|
|
% TREE = SEQNEIGHJOIN(DIST,METHOD,NAMES) passes a list of names to label
|
|
% the leaf nodes (e.g. species or products) in the phylogenetic tree
|
|
% object. NAMES can be a vector of structures with the fields 'Header' or
|
|
% 'Name' or a cell array of strings. In both cases the number of elements
|
|
% provided must comply with the number of samples used to generate the
|
|
% pairwise distances in DIST.
|
|
%
|
|
% TREE = SEQNEIGHJOIN(...,'REROOT',false) excludes rerooting the resulting
|
|
% tree. This is useful to observe the original linkage order followed by
|
|
% the algorithm. By default SEQNEIGHJOIN reroots the resulting tree using
|
|
% the mid-point method.
|
|
%
|
|
% Example:
|
|
%
|
|
% % Load a multiple alignment of amino acids:
|
|
% seqs = fastaread('pf00002.fa');
|
|
%
|
|
% % Measure the 'Jukes-Cantor' pairwise distances:
|
|
% dist = seqpdist(seqs,'method','jukes-cantor','indels','pair');
|
|
%
|
|
% % Build the phylogenetic using the neighbor-joining algorithm
|
|
% tree = seqneighjoin(dist,'equivar',seqs)
|
|
% view(tree)
|
|
%
|
|
% See also MULTIALIGN, PHYTREE, PHYTREE/REROOT, PHYTREE/VIEW, SEQLINKAGE,
|
|
% SEQPDIST.
|
|
|
|
% References:
|
|
% [1] Saitou N, Nei M.The neighbor-joining method: a new method for
|
|
% reconstructing phylogenetic trees. Mol Biol Evol.(1987) 4(4):406-25
|
|
% [2] Gascuel O. BIONJ: An improved version of the NJ algorithm based on a
|
|
% simple model of sequence data. Mol. Biol. Evol. (1997) 14:685-695
|
|
% [3] Studier JA, Keppler KJ. A note on the neighbor-joining algorithm of
|
|
% Saitou and Nei. Mol Biol Evol. (1988) 5(6):729-31.
|
|
|
|
% Copyright 2003-2006 The MathWorks, Inc.
|
|
% $Revision: 1.1.8.2 $ $Date: 2006/05/17 20:48:43 $
|
|
|
|
rerootTree = true;
|
|
checkForNames = true;
|
|
|
|
% check the input distances
|
|
if isnumeric(D)
|
|
[m, n] = size(D);
|
|
if isvector(D)
|
|
n = numel(D);
|
|
m = (1+sqrt(1+8*n))/2; % number of leaves
|
|
if m ~= fix(m)
|
|
error('Bioinfo:seqneighjoin:DbadSize',...
|
|
'Size of DIST not compatible with the output of the SEQPDIST function.');
|
|
end
|
|
D = squareform(D);
|
|
elseif m~=n
|
|
error('Bioinfo:seqneighjoin:DnotSquare',...
|
|
'Size of DIST not compatible with the output of the SEQPDIST function.');
|
|
end
|
|
else
|
|
error('Bioinfo:seqneighjoin:DnotNumeric',...
|
|
'DIST must be a numeric vector compatible with the output of the SEQPDIST function.');
|
|
end
|
|
|
|
% Selects appropiate method
|
|
if nargin == 1
|
|
method = 'e'; % set default method
|
|
checkForNames = false;
|
|
else
|
|
okmethods = {'equivar','firstorder','average','reroot'};
|
|
methodkeys = {'e','f','a'};
|
|
k = find(strncmpi(method,okmethods,numel(method)));
|
|
if isempty(k)
|
|
error('Bioinfo:seqneighjoin:UnknownMethod',...
|
|
'Unknown method name: %s.',method);
|
|
elseif length(k)>1
|
|
error('Bioinfo:seqneighjoin:IncorrectMethod',...
|
|
'Ambiguous method name: %s.',method);
|
|
elseif k<4
|
|
method = methodkeys{k};
|
|
else % case that second and third input are optional paired inputs
|
|
if numel(varargin) || nargin==2
|
|
error('Bioinfo:seqneighjoin:IncorrectInputArguments',...
|
|
'Incorrect format of input arguments.')
|
|
end
|
|
varargin = {method,names};
|
|
checkForNames = false;
|
|
method = 'e'; % set default method
|
|
end
|
|
end
|
|
|
|
% detects the names
|
|
if checkForNames && (nargin >2) % names were supplied, check validity
|
|
if iscell(names) || isfield(names,'Header') || isfield(names,'Name') || isfield(names,'LocusName')
|
|
if isfield(names,'Header') % if struct put them in a cell
|
|
names = {names(:).Header};
|
|
elseif isfield(names,'Name') % if struct put them in a cell
|
|
names = {names(:).Name};
|
|
elseif isfield(names,'LocusName') % if struct put them in a cell
|
|
names = {names(:).LocusName};
|
|
end
|
|
names = names(:);
|
|
namesSupplied = true;
|
|
if numel(names)~=m
|
|
error('Bioinfo:seqneighjoin:IncorrectSize',...
|
|
'NAMES must have the same size as number of leaves in the tree')
|
|
end
|
|
elseif strncmpi(names,{'reroot'},numel(names))
|
|
varargin = {names varargin{:}};
|
|
namesSupplied = false;
|
|
else
|
|
error('Bioinfo:seqneighjoin:IncorrectInputType',...
|
|
'NAMES must be a cell with char arrays or a vector of structures.')
|
|
end
|
|
else
|
|
namesSupplied = false;
|
|
end
|
|
|
|
% check optional input arguments
|
|
nvarargin = numel(varargin);
|
|
if nvarargin
|
|
if rem(nvarargin,2)
|
|
error('Bioinfo:seqneighjoin:IncorrectNumberOfArguments',...
|
|
'Incorrect number of arguments to %s.',mfilename);
|
|
end
|
|
okargs = {'reroot',''};
|
|
for j=1:2:nvarargin
|
|
pname = varargin{j};
|
|
pval = varargin{j+1};
|
|
k = find(strncmpi(pname,okargs,numel(pname)));
|
|
if isempty(k)
|
|
error('Bioinfo:seqneighjoin:UnknownParameterName',...
|
|
'Unknown parameter name: %s.',pname);
|
|
elseif length(k)>1
|
|
error('Bioinfo:seqneighjoin:AmbiguousParameterName',...
|
|
'Ambiguous parameter name: %s.',pname);
|
|
else
|
|
rerootTree = opttf(pval);
|
|
if isempty(rerootTree)
|
|
error('Bioinfo:seqneighjoin:rerootInputOptionNotLogical',...
|
|
'%s must be a logical value, true or false.',...
|
|
upper(char(okargs(k))));
|
|
end
|
|
end
|
|
end
|
|
end
|
|
|
|
% ------------------------------------------------------------------------
|
|
% Algorithm starts here:
|
|
%
|
|
% D - pairwise distance matrix (full form)
|
|
% method - 'a' for Saitou & Nei, 'e' for Studier & Keppler, and 'f' for
|
|
% Gascuel
|
|
% names - labels for leaf nodes
|
|
|
|
N = size(D,1);
|
|
Z = zeros(N-1,2);
|
|
bd = zeros(N*2-1,1);
|
|
|
|
p = 1:N; % pointers in matrix
|
|
bc=1; % branch counter
|
|
if method=='f'
|
|
V=D; %initialize the variance matrix for Gascuel method
|
|
end
|
|
|
|
for n = N:-1:3
|
|
R = sum(D)/(n-2); % sums of columns
|
|
Q = D-repmat(R,n,1)-repmat(R',1,n)+diag(inf(n,1)); %Studier & Keppler optimized
|
|
% S = (sum(sum(D))/(n-2)+Q)/2 % original total sum in Saitou & Nei
|
|
[m,g] = min(Q(:)); %#ok
|
|
[i,j] = ind2sub(n,g); % find minimum
|
|
if i>j
|
|
k=i;i=j;j=k; % j>i always
|
|
end
|
|
pp = p([i j]); % pointers to join
|
|
bl = (R([i j])*[1 -1;-1 1] + D(j,i))/2; % branch lengths
|
|
bl = max(bl,0);
|
|
bd(pp) = bl; % save branch lengths
|
|
Z(bc,:) = pp; % save pointers
|
|
h = [1:i-1 i+1:j-1 j+1:n];
|
|
switch method % distances to new node
|
|
case 'a' % Saitou & Nei method
|
|
d = (sum(D(h,[i j]),2))/2;
|
|
D = [[D(h,h) d];[d' 0]]; % update distance matrix
|
|
case 'e' % Studier & Keppler method
|
|
d = (sum(D(h,[i j]),2)-D(j,i))/2;
|
|
if any(d<0)
|
|
d = max(d,0);
|
|
end
|
|
D = [[D(h,h) d];[d' 0]]; % update distance matrix
|
|
case 'f' % Gascuel method
|
|
if V(i,j)
|
|
lambda = max(0,min(1,(1+sum(V(h,i)-V(h,j))/(n-2)/V(i,j))/2));
|
|
else
|
|
lambda = 1/2;
|
|
end
|
|
d = D(h,[i j])*[lambda;1-lambda] - bl*[lambda;1-lambda];
|
|
if any(d<0)
|
|
d = max(d,0);
|
|
end
|
|
D = [[D(h,h) d];[d' 0]]; % update distance matrix
|
|
v = V(h,[i j])*[lambda;1-lambda] - V(i,j)*lambda*(1-lambda);
|
|
if any(v<0)
|
|
v = max(v,0);
|
|
end
|
|
V = [[V(h,h) v];[v' 0]]; % update variance matrix
|
|
end
|
|
p = [p(h),N+bc]; % update pointers
|
|
bc = bc+1; % update branch counter
|
|
end
|
|
Z(bc,:) = p; % pointers of last branch
|
|
bd(p) = D(2)/2; % lengths of last branch
|
|
|
|
% convert data to a phylogenetic tree object
|
|
if namesSupplied
|
|
t = phyTree(Z,bd,names);
|
|
else
|
|
t = phyTree(Z,bd);
|
|
end
|
|
if rerootTree
|
|
t = reRoot(t); % reroot with mid-point method
|
|
end
|