home *** CD-ROM | disk | FTP | other *** search
/ ftp.disi.unige.it / 2015-02-11.ftp.disi.unige.it.tar / ftp.disi.unige.it / pub / .person / BarlaA / sw / matlab / Cromwell / peakSet.m < prev    next >
Text File  |  2007-12-04  |  6KB  |  161 lines
  1. function [x, y] = peakSet(peaksByMass, minimumSN, secondSN, tickSeparation, massSeparation)
  2. %
  3. % [X, Y] = peakSet(peaksByMass, minSN, secondSN, massSeparation)
  4. % Coalesces (aligns) peaks found in multiple spectra to determine which ones 
  5. % represent the same molecular stuff.
  6. %
  7. % INPUTS:
  8. %    peaksByMass is a matrix with one row for each peak and five columns
  9. %        that represent
  10. %            spectrum ID
  11. %            peak location in clock ticks
  12. %            peak location in mass units
  13. %            signal-to-noise ratio of the peak
  14. %            normalized baseline-corrected intensity of the peak
  15. %    minimumSN = signal-to-noise ratio required for a single peak to count; usually 10
  16. %    secondSN  = signal-to-noise needed to be included on second pass; usually 2
  17. %    tickSeparation = minimum distance between resolvable peaks, in clock
  18. %        ticks, usually 7
  19. %    massSeparation = minimum distance between resolvable peaks as a relative
  20. %        fraction of mass; usually 0.003.
  21. %
  22. % OUTPUTS:
  23. %    X = peak class information, a matrix with one row per peak class and
  24. %    thirteen columns, representing
  25. %        unique class ID
  26. %        min tick
  27. %        max tick
  28. %        min mass
  29. %        max mass
  30. %        number of unique spectra having this peak
  31. %        number of peaks collected in this class
  32. %        min S/N
  33. %        max S/N
  34. %        mean S/N
  35. %        median tick
  36. %        mean of normalized intensity
  37. %        standard deviation of normalized intensity
  38. %    Y = mapping from peaks to classes, a matrix with one row per peak and
  39. %    five columns representing
  40. %        peak class ID (as in X)
  41. %        spectrum ID
  42. %        tick location
  43. %        mass location
  44. %        s/N
  45.  
  46. % Copyright (c) 2003, 2004 UT M.D. Anderson Cancer Center. All rights reserved.
  47. % See the accompanying file "license.txt" for details.
  48.  
  49. % First version written by Kevin Coombes, December 2003
  50.  
  51. % We'd like to assume that the peaks have already been sorted in order of
  52. % increasing mass. Since we're trying to write this as a reusable function,
  53. % we'll resort ourselves to make certain.
  54. [y i] = sort(peaksByMass(:, 3));
  55. pbm = peaksByMass(i,:);
  56.  
  57. [prows, pcols] = size(pbm);    % how many rows are in the peak matrix
  58. peakCount = 1;                        % number of distinct peak classes
  59. lastTick = pbm(1, 2);            % current tick location
  60. lastMass = pbm(1, 3);            % current mass location
  61. peakClasses = [0 0 0 0];        % start of storage for peak class information
  62. peakIDS = repmat(0, [1 floor(prows/20)]);    % storage for which peaks are in this class
  63. collector = [];                    % temporary accumulator
  64.  
  65. % In the first pass, we look for peaks with S/N greater than the minimum
  66. % cutoff. We group them together if their distance in clock ticks or in
  67. % relative mass is less than the specified values. In the peakClasses matrix,
  68. % we store four items about each peak class:
  69. %        minimum tick where it occurs
  70. %        minimum mass
  71. %        maximum tick
  72. %        maximum mass
  73. % In the peakIDS matrix, we store a list of indices that point to the rows
  74. % in the pbm matrix that belong to this class.
  75.  
  76. for i = 1:prows,
  77.    if pbm(i, 4) >= minimumSN,    % ignore smaller S/N on first pass
  78.       tick = pbm(i, 2);
  79.       mass = pbm(i, 3);
  80.       if (tick - lastTick) < tickSeparation | (mass - lastMass)/lastMass < massSeparation,
  81.          collector(end+1) = i;
  82.       else
  83.          p = pbm(collector(1),:);
  84.          q = pbm(collector(end), :);
  85.          peakClasses(peakCount,:) = [p(2) p(3) q(2) q(3)];
  86.          peakIDS(peakCount, 1:length(collector)) = collector;
  87.          collector = [i];
  88.          peakCount = peakCount+1;
  89.       end
  90.       lastTick = tick;
  91.       lastMass=mass;
  92.    end
  93. end
  94. % We haven't yet included the very last peak on the list
  95. p = pbm(collector(1),:);
  96. q = pbm(collector(end), :);
  97. peakClasses(peakCount,:) = [p(2) p(3) q(2) q(3)];
  98. peakIDS(peakCount, 1:length(collector)) = collector;
  99.  
  100. % The first pass on the 24 WCX2 QC spectra used for the wavelet denoising
  101. % paper yields 174 peaks at this point. 
  102.  
  103. % the point of the seocnd pass is to go back and fill in with smaller peaks
  104. % if they align with existing peaks. Here "smaller" means with S/N greater
  105. % than the parameter secondSN. There is probably a better way to do this,
  106. % since I repeat the same code three times, but the logical statement was
  107. % too unwieldy otherwiase.
  108. [crow, ccol] = size(peakClasses);
  109. pointer = 1;                    % which established peak class am I looking at?
  110. for i = 1:prows,
  111.    if pbm(i, 4) < minimumSN & pbm(i, 4) >= secondSN,
  112.       tick = pbm(i, 2);
  113.       mass = pbm(i, 3);
  114.       while (tick - peakClasses(pointer, 3)) > tickSeparation & (mass - peakClasses(pointer, 4))/peakClasses(pointer, 4) > massSeparation,
  115.          pointer = pointer + 1;
  116.          if pointer > crow,
  117.             break; % from inner loop
  118.          end
  119.       end %while
  120.       if pointer > crow,
  121.          break; % from outer loop
  122.       end
  123.       if tick < peakClasses(pointer, 1) & peakClasses(pointer, 1) - tick < tickSeparation,
  124.          temp = peakIDS(pointer,:);    % which actual peaks belong to this class?
  125.          lastID = sum(temp>0);
  126.          peakIDS(pointer, lastID+1) = i;
  127.       elseif tick >= peakClasses(pointer, 1) & tick <= peakClasses(pointer, 3),
  128.          temp = peakIDS(pointer,:);    % which actual peaks belong to this class?
  129.          lastID = sum(temp>0);
  130.          peakIDS(pointer, lastID+1) = i;
  131.       elseif tick > peakClasses(pointer, 3) & tick - peakClasses(pointer, 3) < tickSeparation,
  132.          temp = peakIDS(pointer,:);    % which actual peaks belong to this class?
  133.          lastID = sum(temp>0);
  134.          peakIDS(pointer, lastID+1) = i;
  135.       end
  136.    end
  137. end
  138.  
  139. % now we reassemble the data along with some summary statistics for
  140. % each peak class in order to return something useful.
  141. x = repmat(0, [crow 13]);    % storage for returning class information
  142. y = repmat(0, [prows 6]);        % storage for returning peak information
  143. ypoint = 0;
  144. for i = 1:crow,    % loop over peak sets
  145.    temp = peakIDS(i,:);
  146.    temp = temp(temp>0);
  147.    n = length(temp);
  148.    nuni = length(unique(pbm(temp,1)));
  149.    reup = [1 3 2 4];
  150.    SNs = pbm(temp,4);
  151.    signals = log2(pbm(temp, 5));
  152.    x(i,:) = [i peakClasses(i,reup) nuni n min(SNs) max(SNs) mean(SNs) floor(median(pbm(temp,2))) mean(signals) std(signals)];
  153.    for j = 1:length(temp),
  154.       ypoint = ypoint + 1;
  155.       y(ypoint,:) = [i pbm(temp(j), :)];
  156.    end
  157. end
  158. ygood = y(:,1) > 0;    % only keep peaks that got matched to something with big S/N
  159. y=y(ygood,:);
  160.  
  161.