# HG changeset patch # User jcb-mpl # Date 1619531541 0 # Node ID 0d2d20ef6e727fe36be9ba7fe5a46f08927a0d6b Uploaded diff -r 000000000000 -r 0d2d20ef6e72 2018_07_centroids.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/2018_07_centroids.xml Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,207 @@ + + + + + + scilab + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +**Author** Jean-Claude Boulet (INRAE). + + +--------------------------------------------------- + +==================== +CENTROID CALCULATION +==================== + +----------- +Description +----------- + + +Centroids are calculated according to the reference cited below. + + +------ +Inputs +------ + +**MS format** + +Two formats are possible: HDF5 obtained from Scilab (.dat) or HDF5 obtained from MSconvert then HDFView (.h5) + +The Scilab format contains a structure with two fields: + +- time: the retention times, a vector; + +- mzdata: a list of same length as time; each element of the list is a matrix of two columns representing a MS1 mass spectrum: m/zvalues and associated signals. + +The HDF5 (.h5) file contains at least the following fields: + +- ChomatogramTime: a vector containing the retention times of each mass spectrum. + +- ChromatogramIndex: the indexes of the signals identified in ChomatogramTime; the last value of ChromatogramIndex is the size of ChomatogramTime. + +- ChromatogramIntensity: the sum of intensities of the mass spectra associated to the RTs of ChomatogramTime; ChomatogramTime and ChromatogramIntensity are vectors of same length. + +- SpectrumIndex: a vecteur of integers corresponding to the numbers of signals of each MS in SpectrumMZ et SpectrumIntensity. For example, if the two first values of SpectrumIndex are 1421 and 2887, then the first MS is between 0 à 1421, the second between 1422 and 2887. The last value of spectrumIndex is the number of values in SpectrumMZ and in SpectrumIntensity. + +- SpectrumMz: the m/z values. For each mass spectrum, the first value is the real m/z value, but the following are obtained by the sum of the previous values. For example, if the three first values are: 220.0501, 0.1171 and 0.9856, then the real m/z values to be calculated will be: 220.0501, 220.1672=220.0501+0.1171, 220.2528=220.0501+0.1171+0.9856. + +- SpectrumIntensity: a vector containing all the intensities recorded. + +Even if the MS spectrum has been obtained with the centroid option of MSconvert, this function should be applied to remove all the zeros kept by the MSconcert peak picking option. + + +**MS data** + +The mass spectrum + + +---------- +Parameters +---------- + + +**Baseline threshold** + +Signals under the threshold are dropped. Default value = 0. + +Most of the peaks should be separated by regions with signal=0. + + +**Mode during spectrum acquisition** + +Profile or centroid. + + + +------- +Outputs +------- + +**MS-centroid** + +For each RT, the signals associated to the same m/z value have been replaced by their centroid. + +The output is a structure with two fields: + +- time: retention times; + +- mzdata: m/z values and centroid signals. + + + + + + + +@article{sgcentroid, + title={High-resolution mass spectrometry (HRMS): focus on the m/z values estimated by the Savitzky-Golay first derivative}, + author={Boulet, J.C. and Meudec, E. and Vallverdu-Queralt, A. and Cheynier, V.}, + journal={Rapid Communications in Mass Spectrometry}, + year={2020}, + doi={10.1002/rcm.9036} + } + + + + + + diff -r 000000000000 -r 0d2d20ef6e72 src/mz_1bis_signal_threshold.sci --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_1bis_signal_threshold.sci Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,22 @@ +function res2=mz_1bis_signal_threshold(res,sigthresh); + + res2=res; + + n=max(size(res2.time)); + + for i=1:n; + resi=res2.mzdata(i); + + resi_col2=resi(:,2); + + // application du seuil + resi_col2(resi_col24 then + // calcul des dérivées 1° + d1=savgol_xl(xin(:,2)',5,1); + // Obtention des signes: + signe_d1=sign(d1.d); + + // recherche des valeurs de signal nulles //8dec20 + ind_null=find(xin(:,2)==0); + signe_d1(ind_null)=0; + + // pause + // correction des valeurs nulles + + // 8dec20 +// // correction des bords de Savgol: +// signe_d1(1)=1; +// signe_d1(2)=1; +// signe_d1($-1)=-1; +// signe_d1($)=-1; + + // suppression le 31 mai 18 -------------------------- -> pas de lissage si pics intercallés, ex: ex: -1 -1 1 -1 -1 + + // identification des pics par des n°: 1, 2, 3, etc. 27juin17 + no_pic=zeros(n,1); + no_pic(1)=1; + for i=2:n; + if signe_d1(i-1)<0 & signe_d1(i)>=0 then + no_pic(i)=no_pic(i-1)+1; + else + no_pic(i)=no_pic(i-1); + end + end + + + // identification des pics + all_pics=[]; + // pour chaque pic + i=1; + pic=zeros(1,3); + + while i<=n; + if pic(1)==0 then + pic(1)=xin(i,1); + end + if i signe_d1(i) & pic(2)==0 then + pic=zeros(1,3); // 3oct18 fin de la boucle sans rien faire + + elseif i signe_d1(i) & pic(2)~=0 then // une vallée + pic(3)=xin(i,1); + if all_pics==[] then + all_pics=pic; + else + all_pics=[all_pics;pic]; + end + pic=zeros(1,3); + elseif i==n & pic(1)~=0 & pic(2)~=0 then + pic(3)=xin(i,1); + if all_pics==[] then + all_pics=pic; + else + all_pics=[all_pics;pic]; + end + end + i=i+1; + end + + if size(all_pics,1)== max(no_pic) then // 3oct18 + break + end + end + + // sortie du résultat + xout=all_pics; + + xout2=xin; // 27juin17 + n_out=size(xout,1); + if no_pic($)~=0 & no_pic($) ~=[] then + n_pics=no_pic($); + //if n_out ~=n_pics then //3oct18 + //disp('ligne 154') + //pause + //error('n_out different de n_pics') + //end + for i=1:n_out; + liste_pics=find(no_pic==i); + xout2(liste_pics,1)=xout(i,2); + end + end + + // somme des signaux 9mai18 + xout2_unique=unique(xout2(:,1)); + n_unique=max(size(xout2_unique)); + xout3=zeros(n_unique,2); + xout3(:,1)=xout2_unique; + for i=1:n_unique; + tri=find(xout2(:,1)==xout3(i,1)); + sum_signals=sum(xout2(tri,:),'r'); + xout3(i,2)=sum_signals(2); + end + + + // tri de xout3 (enlever les valeurs de signal nulles) 7dec20 + tri2=find(xout3(:,2))~=0; + xout3=xout3(tri2,:); + + +endfunction diff -r 000000000000 -r 0d2d20ef6e72 src/mz_3_remove_zeros.sci --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_3_remove_zeros.sci Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,18 @@ +function x_out=mz_3_remove_zeros(x_in) + + // x_in: un fichier au format barycentrique avec 2 champs: .time et .mzdata + // si le calcul centroide a été fait par MSconvert, le fichier contient des m/z avec signal=0: on va les enlever + + // initialisation + x_out.time=x_in.time; + x_out.mzdata=list(); + + n=max(size(x_in.time)); + + for i=1:n; + mzdata_i=x_in.mzdata(i); + index=find(mzdata_i(:,2)>0); + x_out.mzdata(i)=mzdata_i(index,:); + end + +endfunction diff -r 000000000000 -r 0d2d20ef6e72 src/mz_eic5_step1.sci --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_eic5_step1.sci Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,45 @@ +function resbary=mz_eic5_step1(res,sigthresh) + + // cette fonction mobilise les fonctions suivantes: + // exec('mz_1bis_signal_threshold.sci',-1); + // exec('mz_2_barycenter_adjust2.sci',-1); + // exec(‘mz_2_savgol2.sci’,-1); + + // Arguments: + // res: une variable avec 2 champs: time et mzdata + // sigthresh: un seuil pour un seul signal sur un seul TR ; ex: sigthresh=10000 + + // resbary: un fichier barycentrique + + + res=mz_1bis_signal_threshold(res,sigthresh); + + resbary=mz_2_barycenter_adjust2(res); + + +endfunction + + + + + + + + + + + + + + + + + + + + + + + + + diff -r 000000000000 -r 0d2d20ef6e72 src/mz_hdf5tosci.sci --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_hdf5tosci.sci Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,103 @@ +function x_mz=mz_hdf5tosci(file_in,max_scans) + + // entrees: ------------------------- + // file_in: un fichier au format HDF5 + // max_scans: valeur introduite pour accelerer la fonction de test + + // sorties: --------------------------- + // x_mz: une variable avec deux champs: + // x_mz.time + // x_mz.mzdata une liste dont chaque élément est un vecteur de 2 colonnes: m/z et intensités de signal + + // JCB 2 aout18 + + h5close() + + // lecture du fichier HDF5 ---------------------------------- + a=h5open(file_in,'r','stdio'); // que lecture 6jan20:stdio + //pause + spectrum_chromatogram_index=h5read(a.root.ChromatogramIndex); + spectrum_chromatogram_index=spectrum_chromatogram_index(1); // une valeur + spectrum_time=h5read(a.root.ChomatogramTime); + spectrum_intensity=h5read(a.root.SpectrumIntensity); + spectrum_mz=h5read(a.root.SpectrumMZ); + spectrum_index=h5read(a.root.SpectrumIndex); + //spectrum_index=spectrum_index(1:spectrum_chromatogram_index); + h5close(a) + + // reconstruction du spectre ----------------- + if argn(2)<2 then + nbr_scans=max(size(spectrum_index)); + else + nbr_scans=max_scans; + end + + + // 1- extraction du temps + //25juin20 + time=spectrum_time; + if size(time,1)==1 then + time=time'; + end + //pas=(spectrum_time($)-spectrum_time(1))/(nbr_scans-1); + //time=[spectrum_time(1):pas:spectrum_time($)+pas]; // rajout de 1 par securite pour être sur d'arriver a nbr_scans + //time=time(1:nbr_scans); // nbr exact de RT + + mzdata=list(); + + //disp(nbr_scans,'nbr_scans=') + + for i=1:nbr_scans; + if i==1 then + plage_i=[1,double(spectrum_index(i))]; + else + plage_i=[double(spectrum_index(i-1))+1,double(spectrum_index(i))]; + end + //disp(plage_i,'plage_i=') + // 2- extraction et reconstruction des m/z + // dans le fichier HDF5, ne sont enregistrées que les différences de m/z entre 2 mesures + mz_i=spectrum_mz(plage_i(1):plage_i(2)); + + n=size(mz_i,2); + mz_i2=zeros(1,n); + + mz_i2(1)=mz_i(1); + for j=2:n; + mz_i2(j)=mz_i2(j-1)+mz_i(j); + end + + // 3- extraction des signaux + sig_i=spectrum_intensity(plage_i(1):plage_i(2)); + + // regroupement des m/z et signaux + sp_i=[mz_i2;sig_i]; + + mzdata(i)=sp_i'; + + + end + + // verif. que time est un vecteur colonne + if size(time,2)>1 then + time=time'; + end + + + // amelioration de la coherence des donnees + index=time(2:$)-time(1:$-1); + index2=find(index<0); + if index2 ~=[] then + index3=min(index2); + else; + index3=max(size(time)); + end + + x_mz.time=time(1:index3); + mzdata2=list(); + for i=1:min(index3,size(mzdata)); + mzdata2(i)=mzdata(i); + end + x_mz.mzdata=mzdata2; + + +endfunction diff -r 000000000000 -r 0d2d20ef6e72 src/savgol_xl.sci --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/savgol_xl.sci Tue Apr 27 13:52:21 2021 +0000 @@ -0,0 +1,192 @@ +function[res_savgol]=savgol_xl(saisir,window2,d012,degre_polynom) + // + // Algorithme basé sur Savitzsky-Golay (SG), NON modifiée cad projection + // sur trois axes non orthogonaux entre eux + // + // ------------------------------------------------------------------------ + // Entrées: + // saisir: matrice au format Saisir ou bien matrice classique 2D + // plage: nombre entier positif impair >=3 + // d012: dérivée (0, 1 ou 2) + // degre_polynom: degré du polynome; par défaut: 2 + // Sorties: + // res_savgol spectres après application de d012 + // filter_d012: un filtre F tel que: X après dérivée d012= X.F + // ------------------------------------------------------------------------ + // + // Avertissement: + // SG ne gère pas les demi-plages à droite et à gauche. Celles-ci + // sont approximées au mieux par régression aux moindres carrés, mais le + // résultat est toujours nettement moins bien que lorsque SG est applicable. + // En conséquence les extrémités droite et gauche sont perturbées, + // et la perturbation s'amplifie avec le degré de dérivation + // + // Avec les filtres, c'est encore pire pour d1 et d2 + // ------------------------------------------------------------------------ + + + //Vérification/transformation en matrice Saisir: + saisir=div(saisir); + + //vérification de l'absence de valeurs manquantes + test=sum(sum(isnan(saisir.d))); + if(test>0) then + xdisp('There are ',test,' Nan values in the data table'); + xdisp('It is necessary to resolve the problem before using savgol'); + error(' '); + end + + + // choix du degré du polynome + if argn(2)==4 then + degre=degre_polynom; + else + degre=2; // degré 2 + end + + if degre<2 | degre>5 then + error('the degree of the polynom should be between 2 and 5') + end + + + //Initialisation: + [n,p]=size(saisir.d); + filter_d0.d=zeros(p,p); //9mars18 + filter_d1.d=zeros(p,p); + filter_d2.d=zeros(p,p); + + rayon=floor((window2-1)/2); + plage=2*rayon+1; + + // verification de la fenêtre + if rayon<1 | plage ~= window2 then + error('the window should be an odd integer, equal to 5 or more '); + end + + d0_savgol=saisir.d; + d1_savgol=ones(n,p); + d2_savgol=ones(n,p); + m=zeros(6,plage); + + + // --------------------------------------------------------- + // calculs de Savitzsky-Golay sur le spectre, cotés exceptés + // --------------------------------------------------------- + + //Calcul de m + vecteur=[-rayon:1:rayon]; + m(1,:)=vecteur.^0; + m(2,:)=vecteur.^1; + m(3,:)=vecteur.^2; + m(4,:)=vecteur.^3; + m(5,:)=vecteur.^4; + m(6,:)=vecteur.^5; + + // ajustement de m à la dimension du polynome + m=m(1:degre+1,:); + + c=inv(m*m')*m; + + for a=rayon+1:p-rayon; // -1; 12mars18 + d0_savgol(:,a)=saisir.d(:,a-rayon:a+rayon)*c(1,:)'; + d1_savgol(:,a)=saisir.d(:,a-rayon:a+rayon)*c(2,:)'; + d2_savgol(:,a)=2*saisir.d(:,a-rayon:a+rayon)*c(3,:)'; + //la multiplication par 2 correspond à la dérivée seconde de: a+bx+cx**2 + + // calcul des filtres + filter_d0.d(a-rayon:a+rayon,a)=c(1,:)'; + filter_d1.d(a-rayon:a+rayon,a)=c(2,:)'; + filter_d2.d(a-rayon:a+rayon,a)=2*c(3,:)'; + end; + + + // --------------------------------------------------- + // calculs de Savitzsky-Golay sur les cotés du spectre + // --------------------------------------------------- + + + // 1 - calcul le plus précis + + //d0 à gauche + calcul=saisir.d(:,1:2*rayon+1)*c'*m; + d0_savgol(:,1:rayon)=calcul(:,1:rayon); + //d0 à droite + calcul=saisir.d(:,p-2*rayon:p)*c'*m; + d0_savgol(:,p-rayon:p)=calcul(:,rayon+1:2*rayon+1); + + //d1 à gauche + d1=d0_savgol(:,2:p)-d0_savgol(:,1:p-1); + d1_savgol(:,1:rayon)=d1(:,1:rayon); + calcul=d1_savgol(:,1:2*rayon+1)*m'*inv(m*m')*m; + d1_savgol(:,1:rayon)=calcul(:,1:rayon); + //d1 à droite + d1_savgol(:,p-rayon:p)=d1(:,p-rayon-1:p-1); //perte d'un indice à cause de la dérivée première + calcul=d1_savgol(:,p-2*rayon:p)*m'*inv(m*m')*m; + d1_savgol(:,p-rayon:p)=calcul(:,rayon+1:2*rayon+1); + + //d2 à gauche + d2=d1_savgol(:,2:p-1)-d1_savgol(:,1:p-2); + d2_savgol(:,1:rayon)=d2(:,1:rayon); + calcul=d2_savgol(:,1:2*rayon+1)*m'*inv(m*m')*m; + d2_savgol(:,1:rayon)=calcul(:,1:rayon); + //d2 à droite + d2_savgol(:,p-rayon:p)=d2(:,p-rayon-2:p-2); //perte de deux indices à cause de la dérivée seconde + calcul=d2_savgol(:,p-2*rayon:p)*m'*inv(m*m')*m; + d2_savgol(:,p-rayon:p)=calcul(:,rayon+1:2*rayon+1); + + +// // 2- calcul sous forme de filtres (+ généralisable mais moins précis) +// for i=1:rayon; +// mi=zeros(6,2*i+1); +// vi=[-i:1:i]; +// mi(1,:)=vi.^0; +// mi(2,:)=vi.^1; +// mi(3,:)=vi.^2; +// mi(4,:)=vi.^3; +// mi(5,:)=vi.^4; +// mi(6,:)=vi.^5; +// mi=mi(1:degre+1,:); +// ci=inv(mi*mi')*mi; //dimensions: (degre+1) x (2*rayon+1) +// filter_d0.d(1:2*i+1,i+1)= ci(1,:)'; // à gauche +// filter_d0.d(p-(2*i+1)+1:p,p-i)= ci(1,:)'; // à droite +// filter_d1.d(1:2*i+1,i+1)= ci(2,:)'; // à gauche +// filter_d1.d(p-(2*i+1)+1:p,p-i)= ci(2,:)'; // à droite +// filter_d2.d(1:2*i+1,i+1)= 2*ci(3,:)'; // à gauche +// filter_d2.d(p-(2*i+1)+1:p,p-i)= 2*ci(3,:)'; // à droite +// if i==1 then +// filter_d0.d(1,1)=1; // meilleur lissage pour 1 valeur = elle-même +// filter_d0.d(p,p)=1; +// filter_d1.d(1:2,1)=[-1;1]; // on soustrait le 1° au 2° +// filter_d1.d(p-1:p,p)=[-1;1]; +// filter_d2.d(1:2*i+1,1)= 2*ci(3,:)'; // on recopie les valeurs à coté; pas trouvé mieux +// filter_d2.d(p-(2*i+1)+1:p,p)= 2*ci(3,:)'; +// end +// end + + + //sélection des sorties + if d012==0 then + res_savgol.d=d0_savgol; + // filter_d012=filter_d0; + elseif d012==1 then + res_savgol.d=d1_savgol; + // filter_d012=filter_d1; + elseif d012==2 then + res_savgol.d=d2_savgol; + // filter_d012=filter_d2; + end + + + // div + res_savgol.i=saisir.i; + res_savgol.v=saisir.v; + res_savgol=div(res_savgol); + +// filter_d012.i=saisir.v; +// filter_d012.v=saisir.v; +// filter_d012=div(filter_d012); + + + +endfunction; + diff -r 000000000000 -r 0d2d20ef6e72 test-data/VI2016_AC_4Areduit_1.dat Binary file test-data/VI2016_AC_4Areduit_1.dat has changed diff -r 000000000000 -r 0d2d20ef6e72 test-data/VI2016_AC_4Areduit_1.h5 Binary file test-data/VI2016_AC_4Areduit_1.h5 has changed diff -r 000000000000 -r 0d2d20ef6e72 test-data/VI2016_AC_4Areduit_1centroide.mat Binary file test-data/VI2016_AC_4Areduit_1centroide.mat has changed