view src/mz_9_compares_peaks.sci @ 2:edc94aa0b5c4 draft default tip

Deleted selected files
author jcb-mpl
date Thu, 12 May 2022 12:15:16 +0000
parents 4660bf9c8059
children
line wrap: on
line source

function x_peaks_final=mz_9_compares_peaks(peaks_list,peaks_detail,x836,diff_mz,k)
    
    // peaks_list et peaks_detail: les pics obtenus avec: 2018_07_mz_eics_to_features_v2
    // x836: les pics de référence ; un div avec 2 colonnes: mz puis RT 
    // k:identification des k plus proches voisins; par défaut: 1 

    // depuis le répertoire: tampon2020/projet_publi_eics/jan20


    x836=div(x836);
    
    if argn(2) <4 then
        diff_mz=0.0050;
    end
    
    if argn(2)<5 then
        k=1;
    end

    diff_mz2=diff_mz;

    
    // tri selon m/z croissant annule le 11fev20
    //[nul,index]=gsort(x836.d(:,1),'g','i');
    //x836_mz=x836(index,:);
    x836_mz=x836;
    
    // bilan : ----------------
    // peaks_list = 4434 x 5 TRmin /TRmax/TRmedian/mz/signal
    // x836_mz: 836 x 2      mz/TR 
    
    // identification des proches voisins 
    //diff_mz2=0.0050; 
    n=size(x836_mz.d,1);
    xout.d=[];
    
    for i=1:n;
        peaks_list3=peaks_list.d(:,4);    // mz
        diff_temp=abs(peaks_list3-x836_mz.d(i,1));
        tri1=find(diff_temp<diff_mz2);
        tri2=find(diff_temp>=diff_mz2);
        //pause
        peaks_list3(tri1)=0;      // distance nulle 
        peaks_list3(tri2)=10;     // forte distance
        //pause
        d=(peaks_list.d(:,3)-x836_mz.d(i,2))**2 + (peaks_list3)**2;
        //pause
        d2=[[1:size(peaks_list.d,1)]' d];     // rajout des indices  
        [d_trie,tri]=gsort(d,'g','i');             // choix des k premiers 
        tri2=tri(1:k); 
        d_trie2=sqrt(d_trie(1:k));                  // distances  
        xout_d=[ones(k,1)*x836_mz.d(i,:) peaks_list.d(tri2,:)]; 
        xout_d=[xout_d(:,[1 6 2 4 7]) tri2];
        if xout.d==[] then           
            xout.d=xout_d;
        else    
            xout.d=[xout.d; xout_d]; 
        end
        //disp(i,'i=')
    end
       
    // mise en forme + rajout des différences 
    xout2.d=[xout.d(:,1:2) abs(xout.d(:,1)-xout.d(:,2)) xout.d(:,3:4)  abs(xout.d(:,3)-xout.d(:,4)) xout.d(:,5:6)];
    xout2.d(:,1:3)=0.0001*round(10000*xout2.d(:,1:3));
    xout2.d(:,4:6)=0.01*round(100*xout2.d(:,4:6));
    xout2.v=['mz-ref';'mz-obs';'diff_mz';'RT-ref';'RT-obs';'diff_RT';'signal';'n° in peak_list'];
    
    x_peaks_final=div(xout2);
    
    // même pic associé à plusieurs mz ref?
    nrepet=ones(k*n,1);
    for i=1:k*n;
        tri=find(xout2.d(:,8)==xout2.d(i,8));
        ntri=max(size(tri));
        nrepet(i)=ntri;
    end
 
    // sorties    
    label_temp=x836.i;
    label_temp=repmat(label_temp,[1,k]);
    label_temp=label_temp';
    label_temp=matrix(label_temp,[k*n,1]);
    x_peaks_final.d=[x_peaks_final.d nrepet]
    x_peaks_final.v=[x_peaks_final.v;'nbr of identifications']
    x_peaks_final.i=label_temp;
    x_peaks_final=div(x_peaks_final);
    
    // rajout des no ref 
    no_ref=[1:n]';
    no_ref=no_ref*ones(1,k);
    no_ref=no_ref';
    no_ref=matrix(no_ref,[k*n,1]);
    
    x_peaks_final.d=[no_ref x_peaks_final.d];
    x_peaks_final.v=['no_ref';x_peaks_final.v];

    // remise en ordre
    x_peaks_final=x_peaks_final(:,[1,9,10,2:8]);
    
    
endfunction