Mercurial > repos > jcb-mpl > orbi_h5_to_dat
changeset 0:03c9a3b59377 draft
Uploaded
| author | jcb-mpl |
|---|---|
| date | Mon, 26 Apr 2021 16:28:21 +0000 |
| parents | |
| children | dcbde1b60b1b |
| files | 2018_07_mz_hdf5_to_dat.xml src/mz_1_correct_time.sci src/mz_1_read_txt.sci src/mz_1_read_txt_ms1_ms2.sci src/mz_hdf5nantestosci.sci src/mz_hdf5tosci.sci src/mz_sci2hdf5.sci test-data/._.DS_Store test-data/VI2016_AC_4Areduit_1.dat test-data/VI2016_AC_4Areduit_1.h5 test-data/VI2016_AC_4Areduit_2.h5 test-data/VI2016_AC_4Areduit_3.h5 |
| diffstat | 12 files changed, 660 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/2018_07_mz_hdf5_to_dat.xml Mon Apr 26 16:28:21 2021 +0000 @@ -0,0 +1,114 @@ +<tool id="2021_05_hdf5_to_dat" name="MS .h5 -> MS .dat " version="0.0.1"> + +<description> from an Orbitrap MS </description> + + +<requirements> + <requirement type="package" >scilab</requirement> + </requirements> + + +<stdio> + <exit_code range="1:" level="fatal" /> +</stdio> + + +<command> <![CDATA[ + if [ -d $__root_dir__/packages/scilab-6.1.0 ]; then $__root_dir__/packages/scilab-6.1.0/bin/scilab-cli -nb -quit -f $* < ${script_file}; else scilab-cli -nb -quit -f $* < ${script_file}; fi +]]> +</command> + + + <configfiles> + <configfile name="script_file"> + <![CDATA[ exec("$__tool_directory__/src/mz_hdf5tosci.sci",-1); ... + if ~isdef('pls') then ... + atomsSystemUpdate(); ... + atomsInstall('FACT'); ... + atomsLoad('FACT'); ... + end; ... + lasterror(); ... + res=mz_hdf5tosci("${filein}"); ... + save("${res}",'res'); ... + if ~isempty(lasterror(%f)) then ... + write(0,lasterror()); ... + end ]]> + </configfile> + </configfiles> + + + <inputs> + <param name="filein" format="h5" type="data" label="Mass spectrum, .h5 format" /> + </inputs> + + + <outputs> + + <data name="res" format="mat" label="${filein.name}.dat" /> + + </outputs> + + + <tests> + <test> + <param name="filein" value="VI2016_AC_4Areduit_1.h5"/> + <output name="res" file="VI2016_AC_4Areduit_1.dat" compare="sim_size" delta="500" /> + </test> + </tests> + + <help> + + +**Author** Jean-Claude Boulet (INRAE). + + +--------------------------------------------------- + +================= +MS .h5 -> MS .dat +================= + +----------- +Description +----------- + +This function converts an Orbitrap HDF5 mass spectrum into a Scilab mass spectrum + +----- +Input +----- + +**an .h5 hrms** + +Previously, Orbitrap devices yielded .raw files. + +MSconvert selected MS1 and performed the conversion .raw -> .md5 + +Then HDF5View performed the conversion .md5 -> .h5 +which is the format of the input file. + + +------- +Outputs +------- + +**a .dat mass spectrum** + +This is an HDF5 file whose structure is specific to Scilab. + +It contains two fields: + +- time: the retention times, a vector; + +- mzdata: a list containing a MS1 spectrum (m/z + signals) for each retention time. + + +</help> + +<citations> + +</citations> + + + +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_1_correct_time.sci Mon Apr 26 16:28:21 2021 +0000 @@ -0,0 +1,24 @@ +function xout=mz_1_correct_time(xin,pas) + + // entrée: xin.time + // xin.mzdata + // sortie: xout.time + // xout.mzdata + // but: corriger le dernier point d'acquision qui est la somme de touts les RT avec Orbitrap + + // pas: pour ne pas garder tous les temps // modifie le 23avril19 + + n=max(size(xin.time)); + + indexes=[1:pas:n-1]; + + n_i=max(size(indexes)); + + xout.time=xin.time(indexes); + xout.mzdata=list(); + for i=1:n_i; + xout.mzdata(i)=xin.mzdata(indexes(i)); + end + + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_1_read_txt.sci Mon Apr 26 16:28:21 2021 +0000 @@ -0,0 +1,120 @@ +function [res]= mz_1_read_txt(xfichier, pas ) + + // lecture d'un fichier texte, obtenu par exportation d'un fichier .mzML + + // xfichier: le nom du fichier texte ex: xfichier -> '151123-vinrose_acet3%_T5h_cpl-pos.txt' + // pas: le pas d'acquisition des données: pas=10 -> 1/10, pas = 1 -> 1/1 = toutes + // pas est un multiple de 10 + // res: une structure + + // ex: res=mz_read_txt('151123-vinrose_acet3%_T5h_cpl-pos.txt',2); + + // compatibilité version 5.5.2 et 6.0.0 + version=getversion(); + version=strsplit(version,'-'); + version=version(2); + version=part(version,1); + version=strtod(version); + if version < 6 then + stacksize('max') + end + + // pas par défaut: pas =1 / on garde toutes les données + if argn(2)==1 then + pas=1; + end + + // ouverture du fichier: + id=mopen(xfichier,'rt'); + + // initialisation: + i=1; + x_out=list(); + x_masses=[]; + x_time=[]; + end_file=0; + list_index=1; + + // début de la boucle while, pour collecter les données de masse + while end_file==0; + + line=mgetl(id,1); + line=stripblanks(line) //enlève les blancs de début et fin, pas les blancs intermédiaires + + if meof(id)==0 & line~=[] then // on continue... + + // temps départ de l'acquisition: + if regexp(line,'/scan start time/')~=[] then + [a]=strsplit(line,','); + time=strtod(a(2)); + end + + // arrêt de boucle à la fin des spectres de masse: + if (i<=10*pas) | (x_time($-1)<x_time($)) then // 10*pas pour avoir des valeurs dans x_time + + // extraction des données: + if regexp(line,'/binary:/')~=[] then + + a=strsplit(line,']'); + xtemp=stripblanks(a(2)); // extrait les données (1 x 1) et enlève les blancs de début et de fin + xtemp=strsplit(xtemp,' '); // met xtemp sous forme d'un vecteur vectical de chaines de caractères + xtemp=strtod(xtemp); // mise en alpha-numérique + + xtemp=xtemp(isnan(xtemp)==%F); // suppression des nan rajoutés accidentellement en fin de fichier xtemp, 16mars17 + + if x_masses==[] then + x_masses=xtemp; + else + if size(x_masses,1)==size(xtemp,1) then + x_masses=[x_masses xtemp]; + else + disp('erreur de dimensions ligne 73 de mz_read_txt' ) + disp(size(x_masses),'size(x_masses)=') + disp(size(xtemp),'size(x_temp)=') + pause + clear x_masses + end + + end + end + + + // sélection d'un spectre de masse et d'une unité de temps par unité de pas: + if size(x_masses,2)==2 then + if floor(i/pas)==ceil(i/pas) then // on en garde 1/pas + if x_time==[] then + x_time=time; + else + x_time=[x_time;time]; + end + x_out(list_index)=x_masses; + //disp(size(x_out(list_index)),'size(x_out(list_index))=') + list_index=list_index+1; + end + //disp(i,'i=') + //disp(list_index,'list_index=') + x_masses=[]; + i=i+1; + end + else // ligne 53 + end_file=1; + end // ligne 53 + + else // ligne 44 + //end_file=meof(id) // on arrête la boucle while + end_file=1; + end // ligne 44 + + end // ligne 39 + + mclose(id) + + res.time=x_time; + if size(res.time,1)==1 then + res.time=res.time'; + end + res.mzdata=x_out; + + + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_1_read_txt_ms1_ms2.sci Mon Apr 26 16:28:21 2021 +0000 @@ -0,0 +1,144 @@ +function [res]= mz_1_read_txt_ms1_ms2(xfichier, ms_level) + + // lecture d'un fichier texte, obtenu par exportation d'un fichier .mzML + + // xfichier: le nom du fichier texte ex: xfichier -> '151123-vinrose_acet3%_T5h_cpl-pos.txt' + // ms_level: 1 (MS1) ou 2 (MS2) + + // res: une structure + + // ex: res=mz_read_txt('151123-vinrose_acet3%_T5h_cpl-pos.txt',2); + + + + + // initialisation: + i=1; + x_out=list(); + x_masses=[]; + x_time=[]; + x_selected_ion=[]; + x_signal_sel_ion=[]; + + end_file=0; + list_index=1; + flag2=1; + + + mclose('all') // si fichier resté ouvert lorsd'une session précédente + + // ouverture du fichier: + id=mopen(xfichier,'rt'); + + // début de la boucle while, pour collecter les données de masse + while meof(id)==0 & flag2==1; + line=mgetl(id,1); + line=stripblanks(line) //enlève les blancs de début et fin, pas les blancs intermédiaires + + if line~=[] & flag2==1 then // on étudie la ligne + + if grep(line,'chromatogramList')~=[] | grep(line,'electromagnetic')~=[] | grep(line,'TIC')~=[] then // 7nov20: on a fini le MS2 - on arrête là + flag=0; + flag2=0; + end + + if regexp(line,'/ms level/')~=[] & flag2==1 then // on identifie un nouveau spectre + + [a]=strsplit(line,','); + mslevel=strtod(a(2)); + + if mslevel==ms_level then; // on va garder ce spectre + + // initialisation des sorties + time=-1; + selected_ion=-1; + signal_sel_ion=-1; + x_masses=[]; + + flag=1; + while flag==1 & flag2==1; //23nov + line=mgetl(id,1); + line=stripblanks(line); + if line ~=[] then + + if flag~=0 & flag2~=0 then + if regexp(line,'/scan start time/')~=[] then + a=strsplit(line,','); + time=strtod(a(2)); + end + if regexp(line,'/selected ion/')~=[] then + a=strsplit(line,','); + selected_ion=strtod(a(2)); + end + if regexp(line,'/peak intensity/')~=[] then + a=strsplit(line,','); + signal_sel_ion=strtod(a(2)); + end + + // extraction des données: + if regexp(line,'/binary:/')~=[] then + a=strsplit(line,']'); + xtemp=stripblanks(a(2)); // extrait les données (1 x 1) et enlève les blancs de début et de fin + xtemp=strsplit(xtemp,' '); // met xtemp sous forme d'un vecteur vectical de chaines de caractères + xtemp=strtod(xtemp); // mise en alpha-numérique + xtemp=xtemp(isnan(xtemp)==%F); // suppression des nan rajoutés accidentellement en fin de fichier xtemp, 16mars17 + // le fichier x_masses doit contenir 2 colonnes: m/z et signal + if x_masses==[] then + x_masses=xtemp; + else + if size(x_masses,1)==size(xtemp,1) then + x_masses=[x_masses xtemp]; + else + ('erreur de dimensions ligne 78 de mz_read_txt' ) + (size(x_masses),'size(x_masses)=') + (size(xtemp),'size(x_temp)=') + end + end + end + + if regexp(line,'/defaultArrayLength/')~=[] then // on passe au spectre suivant + if size(x_masses,2)==2 then + if x_time==[] then + x_time=time; + if ms_level==2 then + x_selected_ion=selected_ion; + x_signal_sel_ion=signal_sel_ion; + end + else + x_time=[x_time;time]; + if ms_level==2 then + x_selected_ion=[x_selected_ion;selected_ion]; + x_signal_sel_ion=[x_signal_sel_ion;signal_sel_ion]; + end + end; + x_out(list_index)=x_masses; + list_index=list_index+1; + else + error('error in the dimension of x_masses') + end + flag=0; // on a fini pour ce spectre + end + end + end + end + end + end + end + end + + + mclose(id) + + res.time=x_time; + if size(res.time,1)==1 then + res.time=res.time'; + end + res.mzdata=x_out; + if ms_level==2 then + res.mz_targetted=x_selected_ion; + res.mz_targetted_signal=x_signal_sel_ion; + end + + + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_hdf5nantestosci.sci Mon Apr 26 16:28:21 2021 +0000 @@ -0,0 +1,101 @@ +function x_mz=mz_hdf5nantestosci(file_in,cvparam_txt,max_scans) + + // entrees: ------------------------- + // file_in: un fichier au format md5 (HDF5) du format de BIBS-Nantes + // rappel: mzxml de Nantes -> md5 avec msconvert + // cvparam_txt: l'export en txt de CVParam avec HDFView; un fichier texte + // 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 + + + + // lecture du fichier HDF5 ---------------------------------- + a=h5open(file_in,'a'); + //h5dump(a.root.CVParam) + //spectrum_time=h5read(a, "/CVParam",[1 22],[1 1]); + spectrum_intensity=h5read(a.root.SpectrumIntensity); // intensite du signal; un vecteur (1 x 4E08) + spectrum_mz=h5read(a.root.SpectrumMZ); // m/z; un vecteur (1 x 4E08) + spectrum_index=h5read(a.root.SpectrumIndex); // index d'identification des spectres + // dans les variables précédentes + // ex: 112166 227841 -> 1° spectre de 1 à 112166 + // 2° spectre de 112167 à 227841 + // etc + h5close(a) + + if argn(2)<3 then + nbr_scans=max(size(spectrum_index)); + else + nbr_scans=max_scans; + end + + // 1- reconstruction du spectre -------------------------------------------- + time=spectrum_index(1:nbr_scans); + mzdata=list(); + + 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 + + // 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 + + + // 2- Reconstruction des temps de retention -------------------------------- + + // lecture du fichier + a=mopen(cvparam_txt); + txt=mgetl(a) + mclose(a) + + // tri des lignes + extraction des RT + n=size(txt,1); + x_rt=[]; + + for i=1:n; + x_temp=strsplit(txt(i,:),' 18'); // valeur 18 trouvee avec tab avant + if max(size(x_temp))>1 then + x_temp2=stripblanks(x_temp(1)); + x_temp2=strtod(x_temp2); + if x_rt==[] then + x_rt=x_temp2; + else + x_rt=[x_rt;x_temp2]; + end + end + end + + // 3- prepaaration de la sortie -------------------------------------------- + x_mz.time=x_rt'; + x_mz.mzdata=mzdata; + + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_hdf5tosci.sci Mon Apr 26 16:28: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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mz_sci2hdf5.sci Mon Apr 26 16:28:21 2021 +0000 @@ -0,0 +1,54 @@ +function res_out=mz_sci2hdf5(res,h5name) + + // entrees: ------------------------- + // res: une structure avec les champs .time et .mzdata + + // sorties: + // newname2: le nom du fichier de sortie, au format .h5 + // avec les champs: + // ChromatogramIndex + // ChomatogramTime + // SpectrumIntensity + // SpectrumMZ + // SpectrumIndex + + // preparer les donnes .dat au format .h5 + // load(filedat_in) // donne res.time et res.bary + + n=length(res.time); + + ChromatogramIndex=n; + + ChomatogramTime=res.time; + + SpectrumIntensity=[]; + SpectrumMZ=[]; + SpectrumIndex=[]; + + for i=1:n; + mz_i=res.mzdata(i); // une matrice de 2 colonnes + diff_i=mz_i(2:$,1)-mz_i(1:$-1,1); + if i==1 then + SpectrumIntensity=mz_i(:,2); + SpectrumMZ=[mz_i(1,1); diff_i]; + SpectrumIndex=size(mz_i,1); + else + SpectrumIntensity=[SpectrumIntensity; mz_i(:,2)]; + SpectrumMZ=[SpectrumMZ; mz_i(1,1); diff_i]; + SpectrumIndex=[SpectrumIndex;SpectrumIndex($)+size(mz_i,1)]; + end + end + + // on a reconstitue les champs au format h5 + + b=h5open(h5name,"w"); + h5write(b,"ChromatogramIndex",ChromatogramIndex); + h5write(b,"ChomatogramTime",ChomatogramTime'); + h5write(b,"SpectrumIntensity",SpectrumIntensity'); + h5write(b,"SpectrumMZ",SpectrumMZ'); + h5write(b,"SpectrumIndex",SpectrumIndex'); + h5close(b) + + res_out=1; + +endfunction