diff src/mz_8_peakpicking_dw.sci @ 0:c0a2e5e78725 draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/mz_8_peakpicking_dw.sci	Tue Apr 27 14:48:05 2021 +0000
@@ -0,0 +1,140 @@
+function  [peaks_list,peaks_detail,no_eic]=mz_8_peakpicking_dw(x_tr_mz3_trouves,list_mz,plage_tr,plage_min_tr,dw_max,signal_min,option_plot)
+    
+    // x_tr_mz3_trouves:    un div (n x p)
+    // list_mz:             les n° des m/z à traiter; valeurs entre 1 et p
+    // plage tr:            une valeur; par defaut: +/- 0,2 minutes
+    // plage_min_tr:        une valeur; par defaut: 0,05
+    // dw_max:              une valeur entre 0 et 2; par défaut: 0.5
+    // signal_min:          une valeur; par défaut: 100000
+    // option_plot:         1=figures; 0= non (par défaut)
+    
+    if argn(2)<7 then
+        option_plot=0;
+    end
+    
+    if argn(2)<6 then
+        signal_min=100000;
+    end
+    
+    if argn(2)<5 then
+        dw_max=0.5;
+    end
+    
+    if argn(2)<4 then
+        plage_min_tr=0.05;
+    end
+    
+    if argn(2)<3 then
+        plage_tr=0.2; 
+    end
+    
+    
+    no_eic=[];
+
+    n_mz=size(x_tr_mz3_trouves.d,2);
+    list_mz_possible=[1:n_mz]';
+    list_mz=list_mz_possible(list_mz);  // pour accepter le $
+   
+    n=max(size(list_mz));
+    
+    peaks_list.d=[]; 
+    peaks_list.i=[];
+    peaks_list.v=['TRmin';'TRmax';'TRmedian';'m/z';'signal'];
+    
+    peaks_detail=list();
+    index_peaks_detail=1;
+    
+    // transformation de plage_tr (des valeurs de temps) en n_tr (un nbre de TR)   
+    ntemp=100;  // pour ne pas commencer au début 
+    n_tr=0;
+    val_tr=0;
+    
+    while val_tr<=plage_tr;
+        n_tr=n_tr+1;
+        val_tr=strtod(x_tr_mz3_trouves.i(ntemp+n_tr))-strtod(x_tr_mz3_trouves.i(ntemp));
+    end
+    
+    nbr_total_tr=size(x_tr_mz3_trouves.i,1);             
+    
+    for i=1:n;
+        xi=x_tr_mz3_trouves(:,list_mz(i));
+        
+        signal_max=max(xi.d)                // 1 seule colonne 
+        
+        index_selected=0*strtod((xi.i));  // un vecteur de 0 de même taille que xi.i
+        
+        while signal_max > signal_min;
+            // identification de la plage du pic
+            index=find(xi.d==signal_max & index_selected==0);
+            index_min=max(1,index-n_tr);
+            index_max=min(nbr_total_tr,index+n_tr); 
+            //ajustement du min pour ne pas retomber sur une zone deja selectionnee
+            index_selected_min=index_selected(index_min);
+            while index_selected_min>0;
+                index_min=index_min+1;
+                index_selected_min=index_selected(index_min);
+            end 
+             //ajustement du max pour ne pas retomber sur une zone deja selectionnee
+            index_selected_max=index_selected(index_max);
+            while index_selected_max>0;
+                index_max=index_max-1;
+                index_selected_max=index_selected(index_max);
+            end
+           
+            // options de tri: DW et plage_min_tr
+            xi_sel=xi(index_min:index_max,:);
+            
+            if (strtod(xi.i(index_max))-strtod(xi.i(index_min)))> plage_min_tr then
+            
+                dw=mz_6bis_durbin_watson(detrending(xi_sel')');       // car si ldb élevée et régulière: DW= faible 
+           
+                if dw<dw_max then 
+           
+                    // calcul du signal de la plage 
+                    vect_tr=[index_min:index_max]
+                    signal=sum(xi.d(vect_tr));  
+           
+                    // enregistrement du resultat
+                    if peaks_list.d==[] then
+                        peaks_list.d=[strtod(xi.i(index_min))  strtod(xi.i(index_max)) strtod(xi.i(index)) strtod(xi.v) signal];
+                        peaks_list.i=i;
+                    else 
+                        peaks_list.d=[peaks_list.d; [strtod(xi.i(index_min)) strtod(xi.i(index_max)) strtod(xi.i(index)) strtod(xi.v) signal]];
+                        peaks_list.i=[peaks_list.i;i];
+                    end
+                
+                    peaks_detail(index_peaks_detail)=xi([index_min:index_max],:);
+                    index_peaks_detail=index_peaks_detail+1;
+        
+                    // edition de figure (selon option)
+                    if option_plot==1 then 
+                        figure;
+                        curves(xi([index_min:index_max],:)); 
+                    end          
+            
+                    // rajout des n° d'EICs de x_tr_mz3_trouves ou x_tr_mz4
+                    if no_eic==[] then
+                        no_eic=i;
+                    else
+                        no_eic=[no_eic;i];
+                    end
+                end
+            
+            end 
+            // neutralisation de la plage precedente
+            index_selected(index_min:index_max)=1;      // on supprime
+            xi.d(index_min:index_max)=0;
+            
+            signal_max=max(xi.d);
+        end
+        
+    end
+    
+    
+    // préparation des sorties
+    peaks_list=div(peaks_list);
+    
+    
+
+    
+endfunction