www.pudn.com > maldicode.zip > bPeakDetectCalibrate.m

 
clear all; 
close all;clc 
 
load casecirr_3030_prepro 
 
% Peak Finding 
maxi = 1; 
mini = .5; 
y1 = min(MZ); 
y2 = max(MZ); 
p    = (mini-maxi)/(y2-y1); 
q    = maxi-y1*p;  
 
YN = sample_tr; 
 
[m,n]=size(YN); 
peaksByMass = []; 
 
for i = 1: n 
    slopeSign = diff(YN(:,i))>0; 
    slopeSignChange = diff(slopeSign)<0; 
     
    h = find(slopeSignChange) + 1; % finds all peaks 
    h( YN(h,i) <  p*MZ(h)+q )= []; % eliminate small peaks 
 
for j = 1: length(h) 
           a    =   h(j); 
           b    =   MZ(a); 
           c    =   YN(a,i); 
           f    =   [i a b c c c]; 
           peaksByMass  =   [peaksByMass; f]; 
    end 
end 
 
% save PeakMass2 peaksByMass 
 
% %% 
% clear all; 
% close all;clc 
%  
% load casecirr_3030_prepro 
% load PeakMass2.mat 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%     PEAKSET     %%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
% PARAMETERS 
minimumSN = 4; 
minSN_end = 1; 
 
secondSN = 1; 
secSN_end = .5; 
 
y1 = min(MZ); 
y2 = max(MZ); 
 
p1  = (minSN_end-minimumSN)/(y2-y1); 
q1  = minimumSN-y1*p1; 
 
p2  = (secSN_end-secondSN)/(y2-y1); 
q2  = secondSN-y1*p2; 
 
% peaksByMass  
for i = 1: length(peaksByMass) 
    peaksByMass(i,4) = peaksByMass(i,4) + (minimumSN - (p1*peaksByMass(i,3)+q1)); 
    peaksByMass(i,5) = peaksByMass(i,5) + (secondSN - (p2*peaksByMass(i,3)+q2)); 
end 
 
tickSeparation = 7; 
massSeparation = 0.0009; 
tic 
[x,y] = peakSet_R(peaksByMass, minimumSN, secondSN, tickSeparation, massSeparation); 
toc 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Remove windows with less that 5 spectra having a peak in that window 
indr =find(x(:,6) <= 5); 
tempy = []; 
tempx = []; 
for i = 1: length(indr) 
    ty = find(y(:,1) == indr(i)); 
    tempy =[tempy; ty]; 
     
    tx = find(x(:,1) == indr(i)); 
    tempx =[tempx; tx]; 
end 
y(tempy,:)=[]; 
x(tempx,:)=[]; 
 
% Extract windows 
    m        =   length(unique(y(:,1))); 
    kk       =   unique(y(:,1)); 
    windows  =   []; 
    win_tick =   []; 
    m_in     =   []; 
    for k = 1: m 
        ind = find(y(:,1)==kk(k)); 
        v   = y(ind,4); 
        vv  = y(ind,3); 
        windows  = [windows; min(v) max(v)]; 
        win_tick = [win_tick;min(vv) max(vv)]; 
    end 
 
save Windows x y windows win_tick ; 
 
%% 
  
% PLOT WINDOWS 
YN = sample_tr; 
figure(); 
hold on; 
for i=1:m; 
    maa = windows(i,1); 
    plot([maa,maa],[0 99],'LineWidth',1.5,'Color',[.9 0.5 0.5]); 
    while maa < windows(i,2) 
     maa = maa+.1; 
     plot([maa,maa],[0 99],'LineWidth',1.5,'Color',[.9 0.5 0.5]);   
    end 
end 
  
% Plot spectrum 
plot(MZ,YN,'LineWidth',1,'Color',[0.2 0.2 0.2]); hold on; 
title('Peak Alignment','FontSize',26); 
xlabel('m/z','FontSize',24); 
ylabel('Normalized Intensity','FontSize',24); 
axis([919 9980 -2 100]); 
 
ss  = YN; 
lab = label_tr; 
ca = mean(ss(:,find(lab==1)),2); 
co = mean(ss(:,find(lab==2)),2); 
cr = mean(ss(:,find(lab==3)),2); 
 
plot(MZ,ca,'r'); hold on; 
plot(MZ,co,'g'); hold on; 
plot(MZ,cr,'b');  
axis([919 9980 -2 100]);