view src/mz_0_publi_features4_galaxy.sci @ 3:b82aec727e31 draft default tip

Deleted selected files

function [x_bilan,x_synthese,x_tr_mz_trouves]=mz_0_publi_features4_galaxy(mz_ref,x_tr_mz4,electron,threshold)

    // entrées: -------
    // mz_ref:      un fichier tabulation avec 1° colonne = identifiants, 2° colonne = valeurs de m/z; 1° ligne = ids des colonnes
    //              ce fichier contient les valeurs de référence  
    // x_tr_mz:     un fichier Scilab contenant une matrice d'EICs nommée x_tr_mz4
    // electron:    enlever l'électron (1) ou pas (0)
    
    // sorties: -------
    // x_bilan:     les résultats principaux avec en colonne les pas choisis (0,0005, 0,0010 etc)
    //              et en ligne les indices: recovery / unicity / rmsep / dw 
    // x_synthese.x_trouves: une liste de divs; pour chacun: .i = indice trouvé; .v = m/z exp. / m/z ref. / Durbin-Watson
    // x_synthese.x_reference: un div contenant les références trouvées dans le fichier mz_ref
    // x_tr_mz_trouve: la liste de matrices de x_tr_mz4 après extraction des m/z de référence trouvés
    // x_bilan: 
  
    
    if argn(2)<4 | threshold<0 then
        threshold=[0.0005,0.0010,0.0030,0.0050,0.0100];
    end

// 22avril21
    // charger le  fichier des EICs obtenu par Galaxy
//    load(x_tr_mz);       // la variable s'appelle: x_tr_mz4 

    // 24nov19
//    list_var=who_user(%f);
//    list_var=list_var(1);
//    for i=1:size(list_var,1);
//        execstr(msprintf("xtemp= %s ;",list_var(i)));
//        if type(xtemp)==16 then 
//            x_tr_mz4=xtemp;
//            //disp(sel_var,'sel_var=')
//        end
//    end

//    // patch pour corriger le mauvais fonctionnement du code ci-dessus
//    if isdef('x_tr_mz3') then
//        x_tr_mz4=x_tr_mz3;
//    end

    // si un M devant les valeurs de m/z: le supprimer 
    //n_v=size(x_tr_mz4.v(i),1);
    //for i=1:n_v;
    x_tr_mz4.v=strsubst(x_tr_mz4.v,'M','');
    //end
 
    // obtenir les valeurs de m/z
    x_mz=strtod(x_tr_mz4.v);


    // charger les donnees de reference en csv classique 
    x_anna0=glx_tab2div(mz_ref);

    //30juin20 on ne garde que 1° colonne
    x_anna0=x_anna0(:,1);
 
    // correction de la masse de l’electron
    // 1 -> on corrige; 0 -> on ne corrige pas 
    x_anna0.d=x_anna0.d-electron*0.00055;
    
    // arrondis et suppression des valeurs en double
    x_anna0.d=0.0001*round(10000*x_anna0.d);                    // arrondis à E-4  (au lieu de 10E-5 avant le mai18 car ca faisait des doubles)
    x_anna0.d=gsort(x_anna0.d,'g','i');                         // tri par ordre croissant
    [temp,indice_unique]=unique(x_anna0.d);                     // valeurs uniques a 10E-4 pres 
    x_anna=x_anna0;
    x_anna=x_anna(indice_unique,:);                             // on garde la precision a 10E-5
    //disp(size(x_anna.d),'size(x_anna.d)=')
    
    // recherche des composés d'anna:
    n_anna=size(x_anna.d,1);
    x_verif.d=[x_anna.d zeros(n_anna,1)];
    x_verif.i=x_anna.i;
    x_verif.v=[x_anna.v;'nbr of peacks identified'];
    x_verif=div(x_verif);

    // gestion de la dimension de threshold 
    n_thresh=max(size(threshold));
    
    // initialisation
    recovery=zeros(1,n_thresh);
    unicity=zeros(1,n_thresh);
    rmsep=zeros(1,n_thresh);
    dw=zeros(1,n_thresh);
    xtrouves_all=[];
    xtrouves_all2=list(); 
  
    
    for k=1:n_thresh;          // pour chacune des valeurs de thresh
    //disp(k,'k=')
        x_trouves=[];
        xtrouves_all=[];
    
        for i=1:n_anna;
            temp_diff=abs(x_mz-x_anna.d(i));
            trouve=find(temp_diff<threshold(k) & temp_diff==min(temp_diff));
            trouve2=find(temp_diff<threshold(k));
            n_trouve2=max(size(trouve2));
            if trouve==[] then
                x_verif.d(i,2)=0;
            elseif max(size(trouve))==1 then
                x_rajoute=[x_mz(trouve) x_anna.d(i) n_trouve2]; 
                        // dans l'ordre: mz calculee la plus proche de mz anna / mz Anna / nbre de mz plus proches de mz Anna que threshold
                if x_trouves==[] then
                    x_trouves=x_rajoute;
                else
                    x_trouves=[x_trouves;x_rajoute];
                end
            end
            
            // construction de xtrouves_all qui rassemble tous les resultats trouves
            if trouve2 ~=[] then
                for j=1:n_trouve2;
                    x_add=[x_mz(trouve2(j)) x_anna.d(i) ];
                    if xtrouves_all==[] then
                        xtrouves_all=x_add;
                        xtrouves_indices=trouve2(j);
                    else
                        xtrouves_all=[xtrouves_all;x_add];
                        xtrouves_indices=[xtrouves_indices;trouve2(j)];
                    end
                end
            end
            
            
        end    
  
        if isdef('xtrouves_indices') then   // cas ou on a trouve qq chose  13sept18
  
            // calcul du rmsep
            diff2=x_trouves(:,1)-x_trouves(:,2);
            recovery(k)=size(x_trouves,1);
            rmsep(k)=sqrt((diff2'*diff2)/recovery(k));
            unicity(k)=sum(x_trouves(:,3))/recovery(k);
            if unicity(k)==0 then 
                pause
            end 
      
            // calcul de Durbin-Watson:
            //pause
            x_dw=x_tr_mz4(:,xtrouves_indices);
            res_dw=mz_6bis_durbin_watson(x_dw);
            dw(k)=mean(res_dw) /unicity(k);
//pause
            res_dw=ones(size(xtrouves_all,1),1);
            // rajout de DW a x_trouves_all
            xtrouves_all2_temp.d=[xtrouves_all res_dw];
            xtrouves_all2_temp.i=string(xtrouves_indices);
            xtrouves_all2_temp.v=['m/z exp.';'m/z ref.';'Durbin-Watson'];
            xtrouves_all2_temp=div(xtrouves_all2_temp);
         
            xtrouves_all2(k)=xtrouves_all2_temp;

        else
            xtrouves_all2(k)=[];
            recovery(k)=-1;
            unicity(k)=-1;
            rmsep(k)=-1;
            dw(k)=-1;
            
        end

    end

    x_bilan.d=[recovery;unicity;rmsep;dw];
    x_bilan.i=['recovery';'unicity';'rmsep';'dw'];
    if size(threshold,2)==1 then
        threshold=threshold';
    end
    x_bilan.v='thresh_'+string(threshold);

    x_bilan=div(x_bilan); 
    
    
    
    // 27juin: recalcul de unicity pour corriger un bug 
    n_xtrouves_all2=size(xtrouves_all2);
    
    nbr_trouves=zeros(1,n_xtrouves_all2);
    
    for i=1:n_xtrouves_all2;
        if xtrouves_all2(i)~=[] then
            nbr_trouves(i)=max(size(xtrouves_all2(i).d));
        else
            nbr_trouves(i)=1;     // 1 car on divise par -1 qq lignes + bas -> on obtient -1 au final 
        end
    end
    
    
    // sorties 
    x_bilan.d(2,:)=nbr_trouves ./ x_bilan.d(1,:);
    
    x_synthese.x_trouves=xtrouves_all2;
    x_synthese.x_reference=x_anna;
    

    // matrice triee ne comprenant que les EICs trouvees
    x_tr_mz_trouves=list();
    
    for i=1:n_thresh;
        x_tr_mz_trouves(i)=x_tr_mz4(:,strtod(x_synthese.x_trouves(i).i));
    end

    
endfunction