view src/mz_5_tr_matrix6bis.sci @ 0:245b2c2b3d75 draft

Uploaded

    function [xout,calculs_stats_tr_mz]=mz_5_tr_matrix6bis(xbary,allmz,thresh,signal_thresh)
        
        // entrées:------------

        // xbary:  une structure de type barycentrique avec les champs x.time et x.mzdata
        // all_mz:  les valeurs de m/z identifiées avec mz_scan_mz_raw2; une matrice de deux colonnes, m/z et signaux 
        // thresh:  le seuil de regroupement, ex: 0,0024
        // signal_thresh: un seuil mini de signal, ex: 30000

        // xout: une matrice de type (TR x m/z)
        //      xout.d:  une matrice contenant les TR en ligne et les valeurs m/z en colonne 
        //      xout.i:  les valeurs de TR
        //      xout.v:  les valeurs de m/z       
        // calculs_stats_tr_mz: une matrice avec nbre lignes=nbre m/z et 7 colonnes contenant differents calculs, pour des statistiques ulterieures
        //      une matrice avec les m/z en ligne et 7 colonnes:
	    //      nbre de TR / somme des TR / somme des TR**2 / somme des mz / somme des mz**2 / somme des mz*signal / somme des signaux 
        
        
        // application d'un filtre sur le niveau minimum de signal 
        tri=find(allmz(:,2)>signal_thresh);
        all_mz=allmz(tri,1);

           
        n=max(size(xbary.time));                                                // nombre de TR
        n_mz=max(size(all_mz));                                                 // nombre de m/z 
     
        // initialisation de xout:
        xout=div(0); 
        xout.d=zeros(n,n_mz);
        xout.i=string(xbary.time);
        xout.v=string(all_mz);
        
        // initialisation de calcul_stats_tr_mz:
        calculs_stats_tr_mz=zeros(n_mz,7);
        
        //calcul_stats_tr_mz contient: 1=n; 2=somme des tr; 3=somme des carres des tri; 4=somme des m/z; 5=somme des carrés des mzi;
        // 6=somme des m/z*signaux; 7=somme des signaux
        
        
        // remplissage de xout et de calculs_stats_tr_mz:
        for i=1:n;
            //xi=xbary.mzdata(i);
            xbaryi=xbary.mzdata(i);             //tous m/z remplacés par leurs barycentres avec mz_barycenter_adjust
            n_xi=size(xbaryi,1);                // nbre de m/z
            for j=1:n_xi;
                abs_differ=abs(all_mz-xbaryi(j,1));
                if min(abs_differ)<thresh then    // 18avril18   on ne remplit que si différence faible 
                    indice=find(abs_differ==min(abs_differ));
                    imax=max(size(indice));
                    for k=1:imax;                                                          //9jan20 boucle rajoutee car il peut y avoir 2 solutions
                        // remplissage de xout:
                        xout.d(i,indice(k))=xout.d(i,indice(k))+ xbaryi(j,2);             //26juil; possibilité de rajouter 2 valeurs dans 1 case
                        // remplissage de calculs_stats_tr_mz:
                        calculs_stats_tr_mz(indice(k),1)=calculs_stats_tr_mz(indice(k),1)+ 1;                           // rajout de +1
                        calculs_stats_tr_mz(indice(k),2)=calculs_stats_tr_mz(indice(k),2)+ xbary.time(i);               // rajout de TR
                        calculs_stats_tr_mz(indice(k),3)=calculs_stats_tr_mz(indice(k),3)+ xbary.time(i)**2;            // rajout de TR**2
                        calculs_stats_tr_mz(indice(k),4)=calculs_stats_tr_mz(indice(k),4)+ xbaryi(j,1);                 // somme des m/z
                        calculs_stats_tr_mz(indice(k),5)=calculs_stats_tr_mz(indice(k),5)+ xbaryi(j,1)**2;              // somme des m/z**2
                        calculs_stats_tr_mz(indice(k),6)=calculs_stats_tr_mz(indice(k),6)+ xbaryi(j,1)*xbaryi(j,2);     // somme des m/z * signal
                        calculs_stats_tr_mz(indice(k),7)=calculs_stats_tr_mz(indice(k),7)+ xbaryi(j,2);                 // somme des signaux
                    end
                end
            end
        
        end
    
        // ajustement des m/z  17mai18
        stats_tr_mz_5= calculs_stats_tr_mz(:,6)./calculs_stats_tr_mz(:,7);
        xout.v= string(stats_tr_mz_5);

         // suppression des signaux trop faibles
        tri=find(calculs_stats_tr_mz(:,7)>signal_thresh);   
        xout=xout(:,tri);                             // 31jan20 verif que ce tri marche 1000X+vite qu'enlever les colonnes 
        calculs_stats_tr_mz= calculs_stats_tr_mz(tri,:);         
        
     
    
    endfunction