peakset_r.m

利用蚁群算法和PSO算法实现MALDI_TOF数据分析

function [x, y] = peakSet_R(peaksByMass, minimumSN, secondSN, tickSeparation, massSeparation)
%
% [X, Y] = peakSet(peaksByMass, minSN, secondSN, massSeparation)
% Coalesces (aligns) peaks found in multiple spectra to determine which ones 
% represent the same molecular stuff.
%
% INPUTS:
%	peaksByMass is a matrix with one row for each peak and five columns
%		that represent
%			spectrum ID
%			peak location in clock ticks
%			peak location in mass units
%			signal-to-noise ratio of the peak
%			normalized baseline-corrected intensity of the peak
%	minimumSN = signal-to-noise ratio required for a single peak to count; usually 10
%	secondSN  = signal-to-noise needed to be included on second pass; usually 2
%	tickSeparation = minimum distance between resolvable peaks, in clock
%		ticks, usually 7
%	massSeparation = minimum distance between resolvable peaks as a relative
%		fraction of mass; usually 0.003.
%
% OUTPUTS:
%	X = peak class information, a matrix with one row per peak class and
%	thirteen columns, representing
%		unique class ID
%		min tick
%		max tick
%		min mass
%		max mass
%		number of unique spectra having this peak
%		number of peaks collected in this class
%		min S/N
%		max S/N
%		mean S/N
%		median tick
%		mean of normalized intensity
%		standard deviation of normalized intensity
%	Y = mapping from peaks to classes, a matrix with one row per peak and
%	five columns representing
%		peak class ID (as in X)
%		spectrum ID
%		tick location
%		mass location
%		s/N

% Copyright (c) 2003, 2004 UT M.D. Anderson Cancer Center. All rights reserved.
% See the accompanying file "license.txt" for details.

% First version written by Kevin Coombes, December 2003

% We'd like to assume that the peaks have already been sorted in order of
% increasing mass. Since we're trying to write this as a reusable function,
% we'll resort ourselves to make certain.
[y i] = sort(peaksByMass(:, 3));
pbm = peaksByMass(i,:);

[prows, pcols] = size(pbm);	% how many rows are in the peak matrix
peakCount = 1;						% number of distinct peak classes
lastTick = pbm(1, 2);			% current tick location
lastMass = pbm(1, 3);			% current mass location
peakClasses = [0 0 0 0];		% start of storage for peak class information
peakIDS = repmat(0, [1 floor(prows/20)]);	% storage for which peaks are in this class
% collector = [1];					% Changed by H. Ressom (was as below)
collector = [];					% temporary accumulator

% In the first pass, we look for peaks with S/N greater than the minimum
% cutoff. We group them together if their distance in clock ticks or in
% relative mass is less than the specified values. In the peakClasses matrix,
% we store four items about each peak class:
%		minimum tick where it occurs
%		minimum mass
%		maximum tick
%		maximum mass
% In the peakIDS matrix, we store a list of indices that point to the rows
% in the pbm matrix that belong to this class.

for i = 1:prows,
   if pbm(i,4) < minimumSN && numel(collector)== 0   % added by Rency 01/31/05
      lastTick = pbm(i, 2);;                         % added by Rency 01/31/05
      lastMass = pbm(i, 3);                          % added by Rency 01/31/05
   end                                               % added by Rency 01/31/05
   if pbm(i, 4) >= minimumSN,	% ignore smaller S/N on first pass
      tick = pbm(i, 2);
      mass = pbm(i, 3);
      if (tick - lastTick) < tickSeparation | (mass - lastMass)/lastMass < massSeparation,
         collector(end+1) = i;
      else
         p = pbm(collector(1),:);
         q = pbm(collector(end), :);
         peakClasses(peakCount,:) = [p(2) p(3) q(2) q(3)];
         peakIDS(peakCount, 1:length(collector)) = collector;
         collector = [i];
         peakCount = peakCount+1;
      end
      lastTick = tick;
      lastMass=mass;
   end
end
% We haven't yet included the very last peak on the list
p = pbm(collector(1),:);
q = pbm(collector(end), :);
peakClasses(peakCount,:) = [p(2) p(3) q(2) q(3)];
peakIDS(peakCount, 1:length(collector)) = collector;

% The first pass on the 24 WCX2 QC spectra used for the wavelet denoising
% paper yields 174 peaks at this point. 

% the point of the seocnd pass is to go back and fill in with smaller peaks
% if they align with existing peaks. Here "smaller" means with S/N greater
% than the parameter secondSN. There is probably a better way to do this,
% since I repeat the same code three times, but the logical statement was
% too unwieldy otherwiase.
[crow, ccol] = size(peakClasses);
pointer = 1;					% which established peak class am I looking at?
for i = 1:prows,
     if pbm(i, 4) < minimumSN & pbm(i, 5) >= secondSN,
      tick = pbm(i, 2);
      mass = pbm(i, 3);
      while (tick - peakClasses(pointer, 3)) > tickSeparation & (mass - peakClasses(pointer, 4))/peakClasses(pointer, 4) > massSeparation,
         pointer = pointer + 1;
         if pointer > crow,
            break; % from inner loop
         end
      end %while
      if pointer > crow,
         break; % from outer loop
      end
      if tick < peakClasses(pointer, 1) & peakClasses(pointer, 1) - tick < tickSeparation,
         temp = peakIDS(pointer,:);	% which actual peaks belong to this class?
         lastID = sum(temp>0);
         peakIDS(pointer, lastID+1) = i;
      elseif tick >= peakClasses(pointer, 1) & tick <= peakClasses(pointer, 3),
         temp = peakIDS(pointer,:);	% which actual peaks belong to this class?
         lastID = sum(temp>0);
         peakIDS(pointer, lastID+1) = i;
      elseif tick > peakClasses(pointer, 3) & tick - peakClasses(pointer, 3) < tickSeparation,
         temp = peakIDS(pointer,:);	% which actual peaks belong to this class?
         lastID = sum(temp>0);
         peakIDS(pointer, lastID+1) = i;
      end
   end
end

% now we reassemble the data along with some summary statistics for
% each peak class in order to return something useful.
x = repmat(0, [crow 13]);	% storage for returning class information
y = repmat(0, [prows 7]);		% storage for returning peak information
ypoint = 0;
for i = 1:crow,	% loop over peak sets
   temp = peakIDS(i,:);
   temp = temp(temp>0);
   n = length(temp);
   nuni = length(unique(pbm(temp,1)));
   reup = [1 3 2 4];
   SNs = pbm(temp,6);
   signals = log2(pbm(temp, 6));
   x(i,:) = [i peakClasses(i,reup) nuni n min(SNs) max(SNs) mean(SNs) floor(median(pbm(temp,2))) mean(signals) std(signals)];
   for j = 1:length(temp),
      ypoint = ypoint + 1;
      y(ypoint,:) = [i pbm(temp(j), :)];
   end
end
ygood = y(:,1) > 0;	% only keep peaks that got matched to something with big S/N
y=y(ygood,:);