Mercurial > repos > mmonsoor > probmetab
changeset 3:abcfa1648b66 draft
planemo upload commit c897279aa8cae0a4ad889bb05b143f32d2b6d712
| author | lecorguille |
|---|---|
| date | Fri, 07 Apr 2017 07:14:12 -0400 |
| parents | 3a9914b37f2f |
| children | 52b222a626b0 |
| files | ProbMetab.xml README.txt lib.r macros.xml probmetab.r static/images/probmetab_workflow.png test-data/faahOK-single.xset.merged.group.retcor.group.fillPeaks.annotate.negative.Rdata test-data/ko15.CDF test-data/ko16.CDF test-data/wt15.CDF test-data/wt16.CDF |
| diffstat | 11 files changed, 753 insertions(+), 564 deletions(-) [+] |
line wrap: on
line diff
--- a/ProbMetab.xml Mon Jul 04 11:58:10 2016 -0400 +++ b/ProbMetab.xml Fri Apr 07 07:14:12 2017 -0400 @@ -1,4 +1,4 @@ -<tool id="Probmetab" name="ProbMetab Tool" version="1.0.1"> +<tool id="Probmetab" name="ProbMetab Tool" version="1.1.0"> <description>Wrapper function for ProbMetab R package.</description> @@ -12,201 +12,243 @@ <command> @COMMAND_CAMERA_SCRIPT@ #if $acquisition_options.mode == "one": - mode_acquisition $acquisition_options.mode xa $acquisition_options.xa + mode_acquisition $acquisition_options.mode + image '$acquisition_options.image' ##if $acquisition_options.xsetnofill_options.option == "show": - ##xsetnofill $acquisition_options.xsetnofill_options.xsetnofill - ##end if + ##xsetnofill $acquisition_options.xsetnofill_options.xsetnofill + ##end if + + @COMMAND_FILE_LOAD_ONE@ + #else - mode_acquisition $acquisition_options.mode inputs_mode $acquisition_options.input_mode.option + mode_acquisition $acquisition_options.mode + inputs_mode $acquisition_options.input_mode.option #if $acquisition_options.input_mode.option== "two": - image_pos $acquisition_options.input_mode.image_pos image_neg $acquisition_options.input_mode.image_neg - ##if $acquisition_options.input_mode.xsetnofill_options.option == "show": - ##xsetPnofill $acquisition_options.input_mode.xsetnofill_options.xsetPnofill xsetNnofill $acquisition_options.input_mode.xsetnofill_options.xsetNnofill - ##end if + image_pos '$acquisition_options.input_mode.image_pos' + image_neg '$acquisition_options.input_mode.image_neg' + ##if $acquisition_options.input_mode.xsetnofill_options.option == "show": + ##xsetPnofill $acquisition_options.input_mode.xsetnofill_options.xsetPnofill + ##xsetNnofill $acquisition_options.input_mode.xsetnofill_options.xsetNnofill + ##end if + + @COMMAND_FILE_LOAD_POSITIVE@ + @COMMAND_FILE_LOAD_NEGATIVE@ ##else - ##image_combinexsannos $acquisition_options.input_mode.image_combinexsannos image_pos $acquisition_options.input_mode.image_pos + ##image_combinexsannos $acquisition_options.input_mode.image_combinexsannos + ##image_pos $acquisition_options.input_mode.image_pos #end if #end if - #if $option_toexclude.option == "show": - toexclude $option_toexclude.toexclude + ## Extraction of CAMERA annotation [get.annot] + allowMiss $getannot.allowMiss + #if $getannot.option_toexclude.option == "show": + toexclude $getannot.option_toexclude.toexclude #end if - allowMiss $allowMiss html $html kegg_db $kegg_db ppm_tol $ppm_tol - opt $opt corths $corths corprob $corprob pcorprob $pcorprob prob $prob - - @COMMAND_ZIPFILE_LOAD@ + + ## Database matching [create.reactionM] + kegg_db $db.kegg_db + ppm_tol $db.ppm_tol + + ## Probability calculations matrix export [export.class.table] + prob $export.prob + html $export.html + + ## Calculate the correlations and partial correlations and cross reference then with reactions [reac2cor] + opt $reac2cor.opt + corprob $reac2cor.corprob + pcorprob $reac2cor.pcorprob + corths $reac2cor.corths + + @COMMAND_LOG_EXIT@ </command> <inputs> <conditional name="acquisition_options"> - <param name="mode" type="select" label="Choose your acquisition mode" > - <option value="one" selected="true" >One acquisition charge mode</option> - <option value="two" >Two acquisition charge mode (positif and negatif)</option> - </param> - - <!-- One acquisition mode--> - <when value="one"> - <param name="xa" type="data" label="Annotate RData" format="rdata.camera.positive,rdata.camera.negative,rdata" help="Output file from annotate step " /> - <!-- - <conditional name="xsetnofill_options"> - <param name="option" type="select" label="RData group step" help="xcmsSet xcms object after missing data replacement, to retrieve SNR to isotopic peaks." > - <option value="show">show</option> - <option value="hide" selected="true">hide</option> - </param> - <when value="show"> - <param name="xsetnofill" type="data" label="Positive or Negative RData from group step before fillpeaks " format="rdata" help=" output from group step" /> - </when> - - </conditional> - --> - </when> - <!-- Two acquisition modes--> - <when value="two"> + <param name="mode" type="select" label="Choose your acquisition mode" > + <option value="one" selected="true" >One acquisition charge mode</option> + <option value="two" >Two acquisition charge mode (positif and negatif)</option> + </param> - - <conditional name="input_mode"> - <param name="option" type="select" label="Choose your input type method:" > - <!-- Bug combinexsannos TODO <option value="one">Input from combinexsAnnos step</option> --> - <option value="two" selected="true">Rdata inputs from annotate</option> - </param> - <!-- - <when value="one"> - <param name="image_combinexsannos" type="data" label="RData output from combinexsAnnos step" format="rdata" help="output file from combinexAnnos step " /> - <param name="image_pos" type="data" label="Positive RData ion mode from annotatediffreport step" format="rdata" help="output file from annotatediffreport step " /> - </when> - --> - <when value="two"> - <param name="image_pos" type="data" label="Positive annotate RData" format="rdata.camera.positive,rdata" help="output file from annotate step " /> - <param name="image_neg" type="data" label="Negative annotate RData" format="rdata.camera.negative,rdata" help="output file from annotate step" /> - <!-- - <conditional name="xsetnofill_options"> - <param name="option" type="select" label="Two RData group step (positive and negative)" help="xcmsSet xcms objects after missing data replacement from your two acquisition modes, to retrieve SNR to isotopic peaks." > - <option value="show">show</option> - <option value="hide" selected="true">hide</option> - </param> - <when value="show"> - <param name="xsetPnofill" type="data" label="Positive RData from group step before fillpeaks " format="rdata.xcms.group,rdata" help="" /> - <param name="xsetNnofill" type="data" label="Negative RData from group step before fillpeaks" format="rdata.xcms.group,rdata" help="" /> - </when> - </conditional> - --> - </when> - </conditional> - </when> - + <!-- One acquisition mode--> + <when value="one"> + <param name="image" type="data" label="Annotate RData" format="rdata.camera.positive,rdata.camera.negative,rdata" help="Output file from annotate step " /> + <!-- + <conditional name="xsetnofill_options"> + <param name="option" type="select" label="RData group step" help="xcmsSet xcms object after missing data replacement, to retrieve SNR to isotopic peaks." > + <option value="show">show</option> + <option value="hide" selected="true">hide</option> + </param> + <when value="show"> + <param name="xsetnofill" type="data" label="Positive or Negative RData from group step before fillpeaks " format="rdata" help=" output from group step" /> + </when> + </conditional> + --> + <expand macro="input_file_load"/> + </when> + <!-- Two acquisition modes--> + <when value="two"> + <conditional name="input_mode"> + <param name="option" type="select" label="Choose your input type method:" > + <!-- Bug combinexsannos TODO <option value="one">Input from combinexsAnnos step</option> --> + <option value="two" selected="true">Rdata inputs from annotate</option> + </param> + <!-- + <when value="one"> + <param name="image_combinexsannos" type="data" label="RData output from combinexsAnnos step" format="rdata" help="output file from combinexAnnos step " /> + <param name="image_pos" type="data" label="Positive RData ion mode from annotatediffreport step" format="rdata" help="output file from annotatediffreport step " /> + </when> + --> + <when value="two"> + <param name="image_pos" type="data" label="Positive annotate RData" format="rdata.camera.positive,rdata" help="output file from annotate step " /> + <param name="image_neg" type="data" label="Negative annotate RData" format="rdata.camera.negative,rdata" help="output file from annotate step" /> + <!-- + <conditional name="xsetnofill_options"> + <param name="option" type="select" label="Two RData group step (positive and negative)" help="xcmsSet xcms objects after missing data replacement from your two acquisition modes, to retrieve SNR to isotopic peaks." > + <option value="show">show</option> + <option value="hide" selected="true">hide</option> + </param> + <when value="show"> + <param name="xsetPnofill" type="data" label="Positive RData from group step before fillpeaks " format="rdata.xcms.group,rdata" help="" /> + <param name="xsetNnofill" type="data" label="Negative RData from group step before fillpeaks" format="rdata.xcms.group,rdata" help="" /> + </when> + </conditional> + --> + </when> + </conditional> + <expand macro="input_file_load" polarity="Positive"/> + <expand macro="input_file_load" polarity="Negative"/> + </when> </conditional> - - <param name="allowMiss" type="boolean" checked="true" truevalue="TRUE" falsevalue="FALSE" label="Retrieves peaks with no eviendence of adduct or isotope" help=" [allowMiss] (ionAnnot function) Logical, annotate also the peaks as single charged molecules [M+/-H]." /> - <conditional name="option_toexclude"> - <param name="option" type="select" label="Exclude samples" > - <option value="show">show</option> - <option value="hide" selected="true">hide</option> - </param> - <when value="show"> - <param name="toexclude" type="text" value="blank,medium,QC" label="samples to be excluded of peak counting to non-annotated peak selection." help="" /> - </when> - <when value="hide" /> - </conditional> - <!-- - <conditional name="useIso_options"> - <param name="option" type="select" label="Calculates the relative isotopic abundance ratio (Carbon 13)" > - <option value="show">Yes</option> - <option value="hide" selected="true">No</option> - </param> - <when value="show"> - <param name="var" type="select" label="var (incorporate.isotopes)" help="1 to use standard mean/sd estimators to carbon number prediction, 2 for median/mad estimators." > - <option value="1">1</option> - <option value="2" selected="true">2</option> - </param> - </when> - </conditional> - --> - <param name="html" type="boolean" checked="true" truevalue="TRUE" falsevalue="FALSE" label="Logical, check if you want to generate a HTML ProbMetab report" help=" [html] (export.class.table function).This parameter uses the raw data to plot EICs and may be time consuming." /> + <section name="getannot" title="Extraction of CAMERA annotation [get.annot]" expanded="True"> + <param name="allowMiss" type="boolean" checked="true" truevalue="TRUE" falsevalue="FALSE" label="Retrieves peaks with no eviendence of adduct or isotope" help=" [allowMiss] (ionAnnot function) Logical, annotate also the peaks as single charged molecules [M+/-H]." /> + <conditional name="option_toexclude"> + <param name="option" type="select" label="Exclude samples" > + <option value="show">show</option> + <option value="hide" selected="true">hide</option> + </param> + <when value="show"> + <param name="toexclude" type="text" value="blank,medium,QC" label="Samples to be excluded of peak counting to non-annotated peak selection." help="[toexclude]" /> + </when> + <when value="hide" /> + </conditional> + </section> - <param name="kegg_db" type="text" size="40" label="Search on KEGG database or multiple organisms " help="Search on all KEGG organisms or multiple organisms (id1,id2,id3,...).By default,the value is KEGG which means searching on all KEGG organism. The list of KEGG IDs are available at http://rest.kegg.jp/list/organism" value="KEGG" > - <validator type="empty_field"/> - </param> - <param name="ppm_tol" type="integer" value="8" label="Parts per million mass tolerance allowed in the mass search" help="[ppm.tol] (create.reactionMfunction) " /> - - - <param name="opt" type="select" label="Correlation option" help="[opt] (reac2cor function) cor for correlation, and pcor for partial correlation." > - <option value="cor" selected="true">cor</option> - <option value="pcor">pcor</option> - </param> - - - <param name="corprob" type="float" value="0.8" label="Probability that the correlation is considered significant" help="[corprob] (reac2cor function) " /> - + <section name="db" title="Database matching [create.reactionM]" expanded="True"> + <param name="kegg_db" type="text" size="40" label="Search on KEGG database or multiple organisms " help="Search on all KEGG organisms or multiple organisms (id1,id2,id3,...).By default,the value is KEGG which means searching on all KEGG organism. The list of KEGG IDs are available at http://rest.kegg.jp/list/organism" value="KEGG" > + <validator type="empty_field"/> + </param> + <param name="ppm_tol" type="integer" value="8" label="Parts per million mass tolerance allowed in the mass search" help="[ppm.tol]" /> + <!-- + <conditional name="useIso_options"> + <param name="option" type="select" label="Calculates the relative isotopic abundance ratio (Carbon 13)" > + <option value="show">Yes</option> + <option value="hide" selected="true">No</option> + </param> + <when value="show"> + <param name="var" type="select" label="var (incorporate.isotopes)" help="1 to use standard mean/sd estimators to carbon number prediction, 2 for median/mad estimators." > + <option value="1">1</option> + <option value="2" selected="true">2</option> + </param> + </when> + </conditional> + --> + </section> - <param name="pcorprob" type="float" value="0.8" label="Probability that the partial correlation is considered significant." help="[pcorprob](reac2cor function)" /> - <param name="corths" type="float" value="0.75" label="Correlation intensity threshold" help="[corths] (reac2cor function)" /> + <section name="export" title="Probability calculations matrix export [export.class.table]" expanded="True"> + <param name="prob" type="select" label=" Calculation of the probability to attribute a mass to a compound" help="[prob] Default is 'count'. See the tool help for more details." > + <option value="count" selected="true">Count</option> + <option value="mean">Mean</option> + </param> + + <param name="html" type="boolean" checked="true" truevalue="TRUE" falsevalue="FALSE" label="Logical, check if you want to generate a HTML ProbMetab report" help="[html] This parameter uses the raw data to plot EICs and may be time consuming." /> + </section> - <param name="prob" type="select" label=" Calculation of the probability to attribute a mass to a compound" help="[prob] (export.class.table function). Default is 'count'. See the tool help for more details." > - <option value="count" selected="true">Count</option> - <option value="mean">Mean</option> - </param> + <section name="reac2cor" title="Calculate the correlations and partial correlations and cross reference then with reactions [reac2cor]" expanded="True"> + <param name="opt" type="select" label="Correlation option" help="[opt] cor for correlation, and pcor for partial correlation." > + <option value="cor" selected="true">cor</option> + <option value="pcor">pcor</option> + </param> + + <param name="corprob" type="float" value="0.8" label="Probability that the correlation is considered significant" help="[corprob]" /> + + <param name="pcorprob" type="float" value="0.8" label="Probability that the partial correlation is considered significant." help="[pcorprob]" /> + <param name="corths" type="float" value="0.75" label="Correlation intensity threshold" help="[corths]" /> + </section> + <!-- + <section name="cytoscape" title="CytoScape options"> <param name="organismId" type="text" size="40" value="NULL" label="organismIdorganismId" help="(create.pathway.node.attributes function) KEGG organism id (http://www.kegg.jp/kegg/catalog/org_list.html) to filter possible pathways for known pathways for that organism. Only works for KEGG database for now. Default is NULL (all KEGG organisms). -" /> + " /> + </section> + --> - --> - - <expand macro="zipfile_load"/> </inputs> <outputs> <!-- <data name="output_image" format="rdata" from_work_dir="probmetab.RData" label="Probmetab.RData" /> --> <data name="html_output" format="html" from_work_dir="AnalysisExample.html" label="Probmetab.Analysis_Report_html" > - <filter>(html)</filter> + <filter>(export['html'])</filter> </data> <data name="tsv_output" format="tabular" from_work_dir="Analysis_Report.tsv" label="Probmetab.CytoScape_output_Attribute_List.tsv" /> <data name="eics" format="zip" from_work_dir="Analysis_Report.zip" label="Probmetab.Analysis_Report_EICs_plots.zip" > - <filter>(html)</filter> + <filter>(export['html'])</filter> </data> <data name="sif_output" format="tabular" from_work_dir="sif.tsv" label="Probmetab.CytoScape_output.sif" /> - <data name="log" format="txt" from_work_dir="probmetab.log" label="Probmetab.log" /> <data name="variableMetadata" format="tabular" from_work_dir="variableMetadata.tsv" label="variableMetadata.tsv" > - <filter>(acquisition_options['mode'] == 'one')</filter> + <filter>(acquisition_options['mode'] == 'one')</filter> </data> <data name="CombineMolIon" format="tabular" from_work_dir="CombineMolIon.tsv" label="CombineMolIon.tsv" > - <filter>(acquisition_options['mode'] == 'two')</filter> + <filter>(acquisition_options['mode'] == 'two')</filter> </data> <data name="variableMetadata_Positive" format="tabular" from_work_dir="variableMetadata_Positive.tsv" label="variableMetadata_Positive.tsv" > - <filter>(acquisition_options['mode'] == 'two')</filter> + <filter>(acquisition_options['mode'] == 'two')</filter> </data> <data name="variableMetadata_Negative" format="tabular" from_work_dir="variableMetadata_Negative.tsv" label="variableMetadata_Negative.tsv" > - <filter>(acquisition_options['mode'] == 'two')</filter> + <filter>(acquisition_options['mode'] == 'two')</filter> </data> </outputs> <tests> <test> - <param name="acquisition_options|mode" value="one" /> - <param name="acquisition_options|xa" value="faahOK.xset.group.retcor.group.fillPeaks.annotate.negative.Rdata" /> - <param name="zipfile_load_conditional|zipfile_load_select" value="yes" /> - <param name="zipfile_load_conditional|zip_file" value="faahKO_reduce.zip" ftype="zip" /> - <output name="log"> - <assert_contents> - <has_text text="ko15 ko16 wt15 wt16" /> - <has_text text="Step 1... determine cutoff point" /> - <has_text text="Step 2... estimate parameters of null distribution and eta0" /> - <has_text text="Step 3... compute p-values and estimate empirical PDF/CDF" /> - <has_text text="Step 4... compute q-values and local fdr" /> - </assert_contents> - </output> + <conditional name="acquisition_options"> + <param name="mode" value="one" /> + <param name="image" value="faahOK.xset.group.retcor.group.fillPeaks.annotate.negative.Rdata" /> + </conditional> + <expand macro="test_commun"/> + <expand macro="test_file_load_zip"/> + <assert_stdout> + <has_text text="Step 1... determine cutoff point" /> + <has_text text="Step 2... estimate parameters of null distribution and eta0" /> + <has_text text="Step 3... compute p-values and estimate empirical PDF/CDF" /> + <has_text text="Step 4... compute q-values and local fdr" /> + </assert_stdout> + </test> + <test> + <conditional name="acquisition_options"> + <param name="mode" value="one" /> + <param name="image" value="faahOK-single.xset.merged.group.retcor.group.fillPeaks.annotate.negative.Rdata" /> + </conditional> + <expand macro="test_commun"/> + <expand macro="test_file_load_single"/> + <assert_stdout> + <has_text text="Step 1... determine cutoff point" /> + <has_text text="Step 2... estimate parameters of null distribution and eta0" /> + <has_text text="Step 3... compute p-values and estimate empirical PDF/CDF" /> + <has_text text="Step 4... compute q-values and local fdr" /> + </assert_stdout> </test> </tests> <help> - + @HELP_AUTHORS@ ========= @@ -225,7 +267,7 @@ **Details** ProbMetab assumes peak detection, retention time correction and peak grouping [4, 5] in order to -perform mass peak to compound assignment. +perform mass peak to compound assignment. Once the initial annotation for different forms of the same ion (adducts and isotopes), is defined, one can seek for a non-redundant set of putative molecules (after charge and possible adduct @@ -235,7 +277,7 @@ experimental condition. In order to address this issue, a flexible workflow, which allows users to integrate different methods, would improve true molecular ions recovery. -The ion annotation table has the following core information: exact mass of putative molecule with experimental error; isotopic pattern associated; adduct form associated, and the original reference to raw data. +The ion annotation table has the following core information: exact mass of putative molecule with experimental error; isotopic pattern associated; adduct form associated, and the original reference to raw data. @@ -250,9 +292,9 @@ ========================= ========================================== ======= ========== Name Output file Format Parameter ========================= ========================================== ======= ========== -xcms.annotate xset.annotate_POS (or NEG).RData RData RData file +xcms.annotate xset.annotate_POS (or NEG).RData RData RData file ========================= ========================================== ======= ========== - + **General schema of the metabolomic workflow** @@ -295,7 +337,7 @@ **Calculate** -**intervals** +**intervals** A vector of SNR numerical intervals, to which different carbon offset should be added to predicted C-number. **offset** @@ -309,7 +351,7 @@ Which noise model to use, "erfc" to complementary error function, or "gaussian" to standard gaussian with two sd corresponding to the given p.p.m precision. -**precision** +**precision** Equipment mass accuracy, usually the same used in exact mass search. @@ -384,6 +426,12 @@ Changelog/News -------------- +**Version 1.1.0 - 06/04/2017** + +- IMPROVEMENT: add some sections within to separate the different parts of the process + +- IMPROVEMENT: Probmetab is now compatible with merged individual data from xcms.xcmsSet + **Version 1.0.1 - 16/05/2016** - TEST: refactoring to pass planemo test using conda dependencies @@ -393,10 +441,8 @@ - NEW: ProbMetab first version - + </help> <expand macro="citation" /> </tool> - -
--- a/README.txt Mon Jul 04 11:58:10 2016 -0400 +++ b/README.txt Fri Apr 07 07:14:12 2017 -0400 @@ -1,19 +1,17 @@ Changelog/News -------------- +**Version 1.1.0 - 06/04/2017** + +- IMPROVEMENT: add some sections within to separate the different parts of the process + +- IMPROVEMENT: Probmetab is now compatible with merged individual data from xcms.xcmsSet + **Version 1.0.1 - 16/05/2016** - TEST: refactoring to pass planemo test using conda dependencies - **Version 1.0.0 - 10/06/2015** - NEW: ProbMetab first version - -Test Status ------------ - -Planemo test using conda: passed - -Planemo shed_test : passed
--- a/lib.r Mon Jul 04 11:58:10 2016 -0400 +++ b/lib.r Fri Apr 07 07:14:12 2017 -0400 @@ -1,323 +1,400 @@ -# lib.r ProbMetab version="1.0.0" -# Author: Misharl Monsoor ABIMS TEAM mmonsoor@sb-roscoff.fr -# Contributors: Yann Guitton and Jean-francois Martin - - -##Main probmetab function launch by the Galaxy ProbMetab wrapper -probmetab = function(xa, xaP, xaN, variableMetadata, variableMetadataP, variableMetadataN, listArguments){ - ##ONE MODE ACQUISITION## - if(listArguments[["mode_acquisition"]]=="one") { - comb=NULL - - #Get the polarity from xa object - polarity=xa@polarity - #SNR option - if ("xsetnofill" %in% names(listArguments)) { - load(listArguments[["xsetnofill"]]) - xsetnofill=xset - } - else{ - xsetnofill=NULL - } - #Exclude samples - if ("toexclude" %in% names(listArguments)) { - toexclude=listArguments[["toexclude"]] - } - else { - toexclude=NULL - } - ionAnnot=get.annot(xa, polarity=polarity, allowMiss=listArguments[["allowMiss"]],xset=xsetnofill,toexclude=toexclude) - comb=NULL - } - - ##TWO MODES ACQUISITION## - #Mode annotatediffreport - else if(listArguments[["inputs_mode"]]=="two"){ - ##Prepare the objects that will be used for the get.annot function - comb=1 - - - xsetPnofill=NULL - xsetNnofill=NULL - # TODO: a reactiver - #if ("xsetPnofill" %in% names(listArguments)) { - # load(listArguments[["xsetPnofill"]]) - # xsetPnofill=xset - #} - #if ("xsetNnofill" %in% names(listArguments)) { - # load(listArguments[["xsetNnofill"]]) - # xsetNnofill=xset - #} - # include CAMERA non-annotated compounds, and snr retrieval - # comb 2+ - used on Table 1 - ionAnnotP2plus = get.annot(axP, allowMiss=listArguments[["allowMiss"]], xset=xsetPnofill,toexclude=listArguments[["toexclude"]]) - ionAnnotN2plus = get.annot(axN, polarity="negative", allowMiss=listArguments[["allowMiss"]], xset=xsetNnofill,toexclude=listArguments[["toexclude"]]) - ionAnnot = combineMolIon(ionAnnotP2plus, ionAnnotN2plus) - print(sum(ionAnnot$molIon[,3]==1)) - print(sum(ionAnnot$molIon[,3]==0)) - write.table(ionAnnot[1], sep="\t", quote=FALSE, row.names=FALSE, file="CombineMolIon.tsv") - #Merge variableMetadata Negative and positive acquisitions mode - - - #Mode combinexsannos TODO bug avec tableau issus de combinexsannos - #else { - #load(listArguments[["image_combinexsannos"]]) - #image_combinexsannos=cAnnot - ##Prepare the objects that will be used for the combineMolIon function - #load(listArguments[["image_pos"]]) - #image_pos=xa - #ionAnnot=combineMolIon(peaklist=cAnnot, cameraobj=image_pos, polarity="pos") - #} - - } - - ##DATABASE MATCHING## - if (listArguments[["kegg_db"]]=="KEGG"){ - DB=build.database.kegg(orgID = NULL) - } - else{ - table_list <<- NULL - ids=strsplit(listArguments[["kegg_db"]],",") - ids=ids[[1]] - if(length(ids)>1){ - for(i in 1:length(ids)){ - table_list[[i]] <- build.database.kegg(ids[i]) - } - db_table=do.call("rbind",table_list) - DB=unique(db_table) - } - else{ - DB=build.database.kegg(listArguments[["kegg_db"]]) - } - } - #Matching des mass exactes mesurees avec les masses des compounds KEGG (pas M+H ou M-H) - reactionM = create.reactionM(DB, ionAnnot, ppm.tol=listArguments[["ppm_tol"]]) - ##PROBABILITY RANKING## - # number of masses with candidates inside the fixed mass window - # and masses with more than one candidate - length(unique(reactionM[reactionM[,"id"]!="unknown",1])) - sum(table(reactionM[reactionM[,"id"]!="unknown",1])>1) - #if (listArguments[["useIso"]]){ - #BUG TODO - # Calculate the ratio between observed and theoretical isotopic patterns. - # If you don't have an assessment of carbon offset to carbon number prediction - # skip this step and use the reactionM as input to weigthM function. - #isoPatt < incorporate.isotopes(comb2plus, reactionM, , samp=12:23, DB=DB) - # calculate the likelihood of each mass to compound assignment using mass accuracy,and isotopic pattern, when present - #wl < weightM(isoPatt,intervals=seq(0,1000,by=500), offset=c(3.115712, 3.434146, 2.350798)) - - #isoPatt=incorporate.isotopes(ionAnnot, reactionM,comb=comb,var=listArguments[["var"]],DB=DB) - - #wl = weightM(reactionM, useIso=true) - #} - #else { - #wl = weightM(reactionM, useIso=FALSE) - #} - wl =weightM(reactionM, useIso=FALSE) - w = design.connection(reactionM) - # Probability calculations - x = 1:ncol(wl$wm) - y = 1:nrow(wl$wm) - conn = gibbs.samp(x, y, 5000, w, wl$wm) - ansConn = export.class.table(conn, reactionM, ionAnnot, DB=DB,html=listArguments[["html"]],filename="AnalysisExample",prob=listArguments[["prob"]]) - if(listArguments[["html"]]){ - #Zip the EICS plot - system(paste('zip -r "Analysis_Report.zip" "AnalysisExample_fig"')) - } - - # calculate the correlations and partial correlations and cross reference then with reactions - mw=which(w==1,arr.ind=TRUE) - #reac2cor function : Use the intensity of putative molecules in repeated samples to calculate correlations and partial - #correlation in a user defined threshold of false discovery rate for significance testing. After the - #correlation test the function also overlay significant correlations with all putative reactions between - #two masses. - #It generates a list of estimated correlations and reactions. - corList=reac2cor(mw,ansConn$classTable,listArguments[["opt"]],listArguments[["corths"]],listArguments[["corprob"]],listArguments[["pcorprob"]]) - ans=list("ansConn"=ansConn,"corList"=corList) - #Generate the siff table for CytoScape - cytoscape_output(corList,ansConn) +# lib.r ProbMetab version="1.0.0" +# Author: Misharl Monsoor ABIMS TEAM mmonsoor@sb-roscoff.fr +# Contributors: Yann Guitton and Jean-francois Martin + + +##Main probmetab function launch by the Galaxy ProbMetab wrapper +probmetab = function(xa, xaP, xaN, variableMetadata, variableMetadataP, variableMetadataN, listArguments){ + ##ONE MODE ACQUISITION## + if(listArguments[["mode_acquisition"]]=="one") { + comb=NULL + + #Get the polarity from xa object + polarity=xa@polarity + #SNR option + if ("xsetnofill" %in% names(listArguments)) { + load(listArguments[["xsetnofill"]]) + xsetnofill=xset + } + else{ + xsetnofill=NULL + } + #Exclude samples + if ("toexclude" %in% names(listArguments)) { + toexclude=listArguments[["toexclude"]] + } + else { + toexclude=NULL + } + ionAnnot=get.annot(xa, polarity=polarity, allowMiss=listArguments[["allowMiss"]],xset=xsetnofill,toexclude=toexclude) + comb=NULL + } + + ##TWO MODES ACQUISITION## + #Mode annotatediffreport + else if(listArguments[["inputs_mode"]]=="two"){ + ##Prepare the objects that will be used for the get.annot function + comb=1 + + + xsetPnofill=NULL + xsetNnofill=NULL + # TODO: a reactiver + #if ("xsetPnofill" %in% names(listArguments)) { + # load(listArguments[["xsetPnofill"]]) + # xsetPnofill=xset + #} + #if ("xsetNnofill" %in% names(listArguments)) { + # load(listArguments[["xsetNnofill"]]) + # xsetNnofill=xset + #} + # include CAMERA non-annotated compounds, and snr retrieval + # comb 2+ - used on Table 1 + ionAnnotP2plus = get.annot(xaP, allowMiss=listArguments[["allowMiss"]], xset=xsetPnofill,toexclude=listArguments[["toexclude"]]) + ionAnnotN2plus = get.annot(xaN, polarity="negative", allowMiss=listArguments[["allowMiss"]], xset=xsetNnofill,toexclude=listArguments[["toexclude"]]) + ionAnnot = combineMolIon(ionAnnotP2plus, ionAnnotN2plus) + print(sum(ionAnnot$molIon[,3]==1)) + print(sum(ionAnnot$molIon[,3]==0)) + write.table(ionAnnot[1], sep="\t", quote=FALSE, row.names=FALSE, file="CombineMolIon.tsv") + #Merge variableMetadata Negative and positive acquisitions mode + + + #Mode combinexsannos TODO bug avec tableau issus de combinexsannos + #else { + #load(listArguments[["image_combinexsannos"]]) + #image_combinexsannos=cAnnot + ##Prepare the objects that will be used for the combineMolIon function + #load(listArguments[["image_pos"]]) + #image_pos=xa + #ionAnnot=combineMolIon(peaklist=cAnnot, cameraobj=image_pos, polarity="pos") + #} + + } + + ##DATABASE MATCHING## + if (listArguments[["kegg_db"]]=="KEGG"){ + DB=build.database.kegg(orgID = NULL) + } + else{ + table_list <<- NULL + ids=strsplit(listArguments[["kegg_db"]],",") + ids=ids[[1]] + if(length(ids)>1){ + for(i in 1:length(ids)){ + table_list[[i]] <- build.database.kegg(ids[i]) + } + db_table=do.call("rbind",table_list) + DB=unique(db_table) + } + else{ + DB=build.database.kegg(listArguments[["kegg_db"]]) + } + } + #Matching des mass exactes mesurees avec les masses des compounds KEGG (pas M+H ou M-H) + reactionM = create.reactionM(DB, ionAnnot, ppm.tol=listArguments[["ppm_tol"]]) + ##PROBABILITY RANKING## + # number of masses with candidates inside the fixed mass window + # and masses with more than one candidate + length(unique(reactionM[reactionM[,"id"]!="unknown",1])) + sum(table(reactionM[reactionM[,"id"]!="unknown",1])>1) + #if (listArguments[["useIso"]]){ + #BUG TODO + # Calculate the ratio between observed and theoretical isotopic patterns. + # If you don't have an assessment of carbon offset to carbon number prediction + # skip this step and use the reactionM as input to weigthM function. + #isoPatt < incorporate.isotopes(comb2plus, reactionM, , samp=12:23, DB=DB) + # calculate the likelihood of each mass to compound assignment using mass accuracy,and isotopic pattern, when present + #wl < weightM(isoPatt,intervals=seq(0,1000,by=500), offset=c(3.115712, 3.434146, 2.350798)) + + #isoPatt=incorporate.isotopes(ionAnnot, reactionM,comb=comb,var=listArguments[["var"]],DB=DB) + + #wl = weightM(reactionM, useIso=true) + #} + #else { + #wl = weightM(reactionM, useIso=FALSE) + #} + wl =weightM(reactionM, useIso=FALSE) + w = design.connection(reactionM) + # Probability calculations + x = 1:ncol(wl$wm) + y = 1:nrow(wl$wm) + conn = gibbs.samp(x, y, 5000, w, wl$wm) + ansConn = export.class.table(conn, reactionM, ionAnnot, DB=DB,html=listArguments[["html"]],filename="AnalysisExample",prob=listArguments[["prob"]]) + if(listArguments[["html"]]){ + #Zip the EICS plot + system(paste('zip -r "Analysis_Report.zip" "AnalysisExample_fig"')) + } + + # calculate the correlations and partial correlations and cross reference then with reactions + mw=which(w==1,arr.ind=TRUE) + #reac2cor function : Use the intensity of putative molecules in repeated samples to calculate correlations and partial + #correlation in a user defined threshold of false discovery rate for significance testing. After the + #correlation test the function also overlay significant correlations with all putative reactions between + #two masses. + #It generates a list of estimated correlations and reactions. + corList=reac2cor(mw,ansConn$classTable,listArguments[["opt"]],listArguments[["corths"]],listArguments[["corprob"]],listArguments[["pcorprob"]]) + ans=list("ansConn"=ansConn,"corList"=corList) + #Generate the siff table for CytoScape + cytoscape_output(corList,ansConn) - #Execute the merge_probmetab function to merge the variableMetadata table and annotations from ProbMetab results - - if(listArguments[["mode_acquisition"]]=="one") { - #Retrocompatibility with previous annotateDiffreport variableMetadata dataframe (must replace mzmed column by mz, and rtmed by rt) - names(variableMetadata)[names(variableMetadata)=="mzmed"] <- "mz" - names(variableMetadata)[names(variableMetadata)=="rtmed"] <- "rt" - variableM=merge_probmetab(variableMetadata, ansConn) - write.table(variableM, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata.tsv") - } else if (listArguments[["mode_acquisition"]]=="two") { - #Retrocompatibility with previous annotateDiffreport variableMetadata dataframe (must replace mzmed column by mz, and rtmed by rt) - names(variableMetadataP)[names(variableMetadata)=="mzmed"] <- "mz" - names(variableMetadataP)[names(variableMetadata)=="rtmed"] <- "rt" - names(variableMetadataN)[names(variableMetadata)=="mzmed"] <- "mz" - names(variableMetadataN)[names(variableMetadata)=="rtmed"] <- "rt" - variableMP=merge_probmetab(variableMetadataP, ansConn) - write.table(variableMP, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata_Positive.tsv") - variableMN=merge_probmetab(variableMetadataN, ansConn) - write.table(variableMN, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata_Negative.tsv") - } - - return(ans) - -} - -##Function that generates a siff table for CytoScape -cytoscape_output=function(corList,ansConn){ - signif_cor=as.data.frame(corList$signif.cor) - classTable=as.data.frame(ansConn$classTable) - #Siff table - siff_table=cbind(signif_cor["node1"],signif_cor["cor"],signif_cor["node2"]) - #attribute table output for Cytoscape - - ## START Code part from the export2cytoscape function of ProbMetab written by Ricardo R. Silva - for (i in 1:nrow(classTable)) if (classTable[i, 1] == ""){ - classTable[i, c(1, 4, 6, 7)] <- classTable[i - 1, c(1, 4, 6, 7)] - } - msel <- as.matrix(classTable[, 1:7]) - msel <- cbind(msel[, 6], msel[,-6]) - colnames(msel)[1] <- "Id" - msel[, 1] <- sub("^\\s+", "", msel[, 1]) - colnames(msel)[1] <- "Id" - ids <- unique(msel[, 1]) - attrMatrix <- matrix("", nrow = length(ids), ncol = ncol(msel)-1) - for (i in 1:length(ids)) { - attrMatrix[i, 1] <- unique(msel[msel[, 1] == ids[i], - 2]) - attrMatrix[i, 2] <- paste("[", paste(msel[msel[, - 1] == ids[i], 3], collapse = ", "), "]", sep = "") - attrMatrix[i, 3] <- paste("[", paste(msel[msel[, - 1] == ids[i], 4], collapse = ", "), "]", sep = "") - attrMatrix[i, 4] <- unique(msel[msel[, 1] == ids[i], - 5]) - attrMatrix[i, 5] <- paste("[", paste(msel[msel[, - 1] == ids[i], 6], collapse = ", "), "]", sep = "") - attrMatrix[i, 6] <- unique(msel[msel[, 1] == ids[i], - 7]) - } - ids <- as.numeric(unique(msel[, 1])) - attrMatrix <- cbind(ids, attrMatrix) - colnames(attrMatrix) <- colnames(msel) - ## END Code part from the export2cytoscape function of ProbMetab writieen by Ricardo R. Silva - write.table(attrMatrix, sep="\t", quote=FALSE, row.names=FALSE, file="Analysis_Report.tsv") - write.table(siff_table, sep="\t", quote=FALSE, row.names=FALSE, file="sif.tsv") - - return(attrMatrix) -} - -##Functions written by Jean-Francois Martin - -deter_ioni <- function (aninfo, pm) -{ - # determine ionisation in ProbMetab result file, used in function merge_probmetab - # input : for 1 ion, aninfo = string with m/z rt and CAMERA annotation from ProbMetab result file - # if the difference between m/z and the probmetab proposed mass is ~1 we use the sign (positive or negative) of this diference - # to define the type of ionisation - # If adduct or fragment was detected, therefore diff >>1 and so, we search for substring "+" ou "2+" ou "3+" ou "-"... - # to define the type of ionisation - # aninfo : vecteur of character resulting of the parsing(sep="#") of the probmetab annotation - if (round(abs(as.numeric(aninfo[1]) - pm),0) ==1) { - if (as.numeric(aninfo[1]) - pm <0) {esi <- "n"} else {esi <- "p"} - } else - if (!is.na(aninfo[4])) { - anstr <- aninfo[4] - # cat(anstr) - if ((grepl("]+",anstr,fixed=T)==T) || (grepl("]2+",anstr,fixed=T)==T) || (grepl("]3+",anstr,fixed=T)==T)) { esi <- "p"} - else - if ((grepl("]-",anstr,fixed=T)==T) || (grepl("]2-",anstr,fixed=T)==T) || (grepl("]3-",anstr,fixed=T)==T)) { esi <- "n"} - # cat(" ioni ",esi,"\n") - } else - { esi <- "u"} - - return(esi) -} - - -merge_probmetab <- function(metaVar,ansConn) { - ## Parse ProbMetab information result file and merge in variable_metaData initial file - ## inputs : - ## metaVar : data.frame of metadataVariable input of probmetab function - ## ansConn : data.frame of ProbMetab result - ## output : dataframe with Probmetab results merge with variableMetadata - ## Constante - ## iannot : indice de la colonne annotation dans le resultat de probMetab - iannot <- 4 - - ## definition of an unique identification of ions mz with 3 decimals and rt(sec) with 1 decimal to avoid - ## duplicate ions name in the diffreport result file - ions <- paste ("M",round(metaVar$mz,3),"T",round(metaVar$rt,1),sep="") - metaVar <- data.frame(ions,metaVar) - - ###### Result data.frame from ProbMetab result list - an_ini <- ansConn$classTable - - ## Suppression of rows without mz and rt or unknown and columns of intensities - ## COLUMNS SUBSCRIPTS HAVE TO BE CHECKED WITh DIFFERENT RESULTS FILES - an <- an_ini[(an_ini[,2]!="unknown"),c(1,2,3,7)] - ## initialisation of vectors receiving the result of the parse of the column annotation (subscrip iannot) - mz <- rep(0,dim(an)[1]) - rt <- rep(0,dim(an)[1]) - propmz <- rep(0,dim(an)[1]) - ioni <- rep("u",dim(an)[1]) - - ## parse the column annotation and define ionisation mode - for (i in 1:dim(an)[1]) { - if (an[i,1] != "") { - info_mzrt <- unlist(strsplit(an[i,iannot],"#")) - propmz[i] <- as.numeric(an[i,1]) - mz[i] <- as.numeric(info_mzrt[1]) - rt[i] <- as.numeric(info_mzrt[2]) - ioni[i] <- deter_ioni(info_mzrt,as.numeric(an[i,1])) - } - else { - propmz[i] <- as.numeric(propmz[i-1]) - mz[i] <- as.numeric(mz[i-1]) - rt[i] <- as.numeric(rt[i-1]) - ioni[i] <- ioni[i-1] - } - } - - ## definition of an unique identification of ions : mz with 3 decimals and rt(sec) with 1 decimal - ## The same as for the metadataVariable data.frame to match with. - ions <- paste ("M",round(mz,3),"T",round(rt,1),sep="") - an <- data.frame(ions,ioni,propmz,mz,rt,an) - - ## transposition of the different probmetab annotations which are in different rows in the initial result data.frame - ## on only 1 row separated with a ";" - li <- as.matrix(table(an$propmz)) - li <- data.frame(dimnames(li)[1],li) - dimnames(li)[[2]][1] <- "propmz" - ions <- rep("u",dim(li)[1]) - propmz <- rep(0,dim(li)[1]) - mpc <- rep("c",dim(li)[1]) - proba <- rep("p",dim(li)[1]) - c <- 0 - while (c < dim(li)[1]) { - c <- c + 1 - suban <- an[an$propmz==li[c,1],] - ions[c] <- as.character(suban[1,1]) - propmz[c] <- as.numeric(suban[1,3]) - mpc[c] <- paste(suban[,7],collapse=";") - proba[c] <- paste(as.character(suban[,8]),collapse=";") - } - - ## Creation of the data.frame with 1 row per ions - anc <- data.frame(ions,propmz,mpc,proba) - anc <- anc[order(anc[,1]),] - - metaVarFinal <- merge(metaVar, anc, by.x=1, by.y=1, all.x=T, all.y=T) - metaVarFinal <- metaVarFinal[,-1] - #write.table(metaVarFinal,file="res.txt", sep="\t", row.names=F, quote=F) - - return (metaVarFinal) -} + #Execute the merge_probmetab function to merge the variableMetadata table and annotations from ProbMetab results + + if(listArguments[["mode_acquisition"]]=="one") { + #Retrocompatibility with previous annotateDiffreport variableMetadata dataframe (must replace mzmed column by mz, and rtmed by rt) + names(variableMetadata)[names(variableMetadata)=="mzmed"] <- "mz" + names(variableMetadata)[names(variableMetadata)=="rtmed"] <- "rt" + variableM=merge_probmetab(variableMetadata, ansConn) + write.table(variableM, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata.tsv") + } else if (listArguments[["mode_acquisition"]]=="two") { + #Retrocompatibility with previous annotateDiffreport variableMetadata dataframe (must replace mzmed column by mz, and rtmed by rt) + names(variableMetadataP)[names(variableMetadata)=="mzmed"] <- "mz" + names(variableMetadataP)[names(variableMetadata)=="rtmed"] <- "rt" + names(variableMetadataN)[names(variableMetadata)=="mzmed"] <- "mz" + names(variableMetadataN)[names(variableMetadata)=="rtmed"] <- "rt" + variableMP=merge_probmetab(variableMetadataP, ansConn) + write.table(variableMP, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata_Positive.tsv") + variableMN=merge_probmetab(variableMetadataN, ansConn) + write.table(variableMN, sep="\t", quote=FALSE, row.names=FALSE, file="variableMetadata_Negative.tsv") + } + + return(ans) + +} + +##Function that generates a siff table for CytoScape +cytoscape_output=function(corList,ansConn){ + signif_cor=as.data.frame(corList$signif.cor) + classTable=as.data.frame(ansConn$classTable) + #Siff table + siff_table=cbind(signif_cor["node1"],signif_cor["cor"],signif_cor["node2"]) + #attribute table output for Cytoscape + + ## START Code part from the export2cytoscape function of ProbMetab written by Ricardo R. Silva + for (i in 1:nrow(classTable)) if (classTable[i, 1] == ""){ + classTable[i, c(1, 4, 6, 7)] <- classTable[i - 1, c(1, 4, 6, 7)] + } + msel <- as.matrix(classTable[, 1:7]) + msel <- cbind(msel[, 6], msel[,-6]) + colnames(msel)[1] <- "Id" + msel[, 1] <- sub("^\\s+", "", msel[, 1]) + colnames(msel)[1] <- "Id" + ids <- unique(msel[, 1]) + attrMatrix <- matrix("", nrow = length(ids), ncol = ncol(msel)-1) + for (i in 1:length(ids)) { + attrMatrix[i, 1] <- unique(msel[msel[, 1] == ids[i], + 2]) + attrMatrix[i, 2] <- paste("[", paste(msel[msel[, + 1] == ids[i], 3], collapse = ", "), "]", sep = "") + attrMatrix[i, 3] <- paste("[", paste(msel[msel[, + 1] == ids[i], 4], collapse = ", "), "]", sep = "") + attrMatrix[i, 4] <- unique(msel[msel[, 1] == ids[i], + 5]) + attrMatrix[i, 5] <- paste("[", paste(msel[msel[, + 1] == ids[i], 6], collapse = ", "), "]", sep = "") + attrMatrix[i, 6] <- unique(msel[msel[, 1] == ids[i], + 7]) + } + ids <- as.numeric(unique(msel[, 1])) + attrMatrix <- cbind(ids, attrMatrix) + colnames(attrMatrix) <- colnames(msel) + ## END Code part from the export2cytoscape function of ProbMetab writieen by Ricardo R. Silva + write.table(attrMatrix, sep="\t", quote=FALSE, row.names=FALSE, file="Analysis_Report.tsv") + write.table(siff_table, sep="\t", quote=FALSE, row.names=FALSE, file="sif.tsv") + + return(attrMatrix) +} + +##Functions written by Jean-Francois Martin + +deter_ioni <- function (aninfo, pm) +{ + # determine ionisation in ProbMetab result file, used in function merge_probmetab + # input : for 1 ion, aninfo = string with m/z rt and CAMERA annotation from ProbMetab result file + # if the difference between m/z and the probmetab proposed mass is ~1 we use the sign (positive or negative) of this diference + # to define the type of ionisation + # If adduct or fragment was detected, therefore diff >>1 and so, we search for substring "+" ou "2+" ou "3+" ou "-"... + # to define the type of ionisation + # aninfo : vecteur of character resulting of the parsing(sep="#") of the probmetab annotation + if (round(abs(as.numeric(aninfo[1]) - pm),0) ==1) { + if (as.numeric(aninfo[1]) - pm <0) {esi <- "n"} else {esi <- "p"} + } else + if (!is.na(aninfo[4])) { + anstr <- aninfo[4] + # cat(anstr) + if ((grepl("]+",anstr,fixed=T)==T) || (grepl("]2+",anstr,fixed=T)==T) || (grepl("]3+",anstr,fixed=T)==T)) { esi <- "p"} + else + if ((grepl("]-",anstr,fixed=T)==T) || (grepl("]2-",anstr,fixed=T)==T) || (grepl("]3-",anstr,fixed=T)==T)) { esi <- "n"} + # cat(" ioni ",esi,"\n") + } else + { esi <- "u"} + + return(esi) +} + + +merge_probmetab <- function(metaVar,ansConn) { + ## Parse ProbMetab information result file and merge in variable_metaData initial file + ## inputs : + ## metaVar : data.frame of metadataVariable input of probmetab function + ## ansConn : data.frame of ProbMetab result + ## output : dataframe with Probmetab results merge with variableMetadata + ## Constante + ## iannot : indice de la colonne annotation dans le resultat de probMetab + iannot <- 4 + + ## definition of an unique identification of ions mz with 3 decimals and rt(sec) with 1 decimal to avoid + ## duplicate ions name in the diffreport result file + ions <- paste ("M",round(metaVar$mz,3),"T",round(metaVar$rt,1),sep="") + metaVar <- data.frame(ions,metaVar) + + ###### Result data.frame from ProbMetab result list + an_ini <- ansConn$classTable + + ## Suppression of rows without mz and rt or unknown and columns of intensities + ## COLUMNS SUBSCRIPTS HAVE TO BE CHECKED WITh DIFFERENT RESULTS FILES + an <- an_ini[(an_ini[,2]!="unknown"),c(1,2,3,7)] + ## initialisation of vectors receiving the result of the parse of the column annotation (subscrip iannot) + mz <- rep(0,dim(an)[1]) + rt <- rep(0,dim(an)[1]) + propmz <- rep(0,dim(an)[1]) + ioni <- rep("u",dim(an)[1]) + + ## parse the column annotation and define ionisation mode + for (i in 1:dim(an)[1]) { + if (an[i,1] != "") { + info_mzrt <- unlist(strsplit(an[i,iannot],"#")) + propmz[i] <- as.numeric(an[i,1]) + mz[i] <- as.numeric(info_mzrt[1]) + rt[i] <- as.numeric(info_mzrt[2]) + ioni[i] <- deter_ioni(info_mzrt,as.numeric(an[i,1])) + } + else { + propmz[i] <- as.numeric(propmz[i-1]) + mz[i] <- as.numeric(mz[i-1]) + rt[i] <- as.numeric(rt[i-1]) + ioni[i] <- ioni[i-1] + } + } + + ## definition of an unique identification of ions : mz with 3 decimals and rt(sec) with 1 decimal + ## The same as for the metadataVariable data.frame to match with. + ions <- paste ("M",round(mz,3),"T",round(rt,1),sep="") + an <- data.frame(ions,ioni,propmz,mz,rt,an) + + ## transposition of the different probmetab annotations which are in different rows in the initial result data.frame + ## on only 1 row separated with a ";" + li <- as.matrix(table(an$propmz)) + li <- data.frame(dimnames(li)[1],li) + dimnames(li)[[2]][1] <- "propmz" + ions <- rep("u",dim(li)[1]) + propmz <- rep(0,dim(li)[1]) + mpc <- rep("c",dim(li)[1]) + proba <- rep("p",dim(li)[1]) + c <- 0 + while (c < dim(li)[1]) { + c <- c + 1 + suban <- an[an$propmz==li[c,1],] + ions[c] <- as.character(suban[1,1]) + propmz[c] <- as.numeric(suban[1,3]) + mpc[c] <- paste(suban[,7],collapse=";") + proba[c] <- paste(as.character(suban[,8]),collapse=";") + } + + ## Creation of the data.frame with 1 row per ions + anc <- data.frame(ions,propmz,mpc,proba) + anc <- anc[order(anc[,1]),] + + metaVarFinal <- merge(metaVar, anc, by.x=1, by.y=1, all.x=T, all.y=T) + metaVarFinal <- metaVarFinal[,-1] + #write.table(metaVarFinal,file="res.txt", sep="\t", row.names=F, quote=F) + + return (metaVarFinal) +} # RETROCOMPATIBILITE avec ancienne version de annotate getVariableMetadata = function(xa) { - # --- variableMetadata --- - peakList=getPeaklist(xa) - peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name"); - variableMetadata=peakList[,!(colnames(peakList) %in% c(sampnames(xa@xcmsSet)))] - variableMetadata$name= paste("M",round(variableMetadata$mz),"T",round(variableMetadata$rt),sep="") - return (variableMetadata) -} + # --- variableMetadata --- + peakList=getPeaklist(xa) + peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name"); + variableMetadata=peakList[,!(colnames(peakList) %in% c(sampnames(xa@xcmsSet)))] + variableMetadata$name= groupnames(xa@xcmsSet) + return (variableMetadata) +} + + +# This function get the raw file path from the arguments +getRawfilePathFromArguments <- function(singlefile, zipfile, listArguments) { + if (!is.null(listArguments[["zipfile"]])) zipfile = listArguments[["zipfile"]] + if (!is.null(listArguments[["zipfilePositive"]])) zipfile = listArguments[["zipfilePositive"]] + if (!is.null(listArguments[["zipfileNegative"]])) zipfile = listArguments[["zipfileNegative"]] + + if (!is.null(listArguments[["singlefile_galaxyPath"]])) { + singlefile_galaxyPaths = listArguments[["singlefile_galaxyPath"]]; + singlefile_sampleNames = listArguments[["singlefile_sampleName"]] + } + if (!is.null(listArguments[["singlefile_galaxyPathPositive"]])) { + singlefile_galaxyPaths = listArguments[["singlefile_galaxyPathPositive"]]; + singlefile_sampleNames = listArguments[["singlefile_sampleNamePositive"]] + } + if (!is.null(listArguments[["singlefile_galaxyPathNegative"]])) { + singlefile_galaxyPaths = listArguments[["singlefile_galaxyPathNegative"]]; + singlefile_sampleNames = listArguments[["singlefile_sampleNameNegative"]] + } + if (exists("singlefile_galaxyPaths")){ + singlefile_galaxyPaths = unlist(strsplit(singlefile_galaxyPaths,",")) + singlefile_sampleNames = unlist(strsplit(singlefile_sampleNames,",")) + + singlefile=NULL + for (singlefile_galaxyPath_i in seq(1:length(singlefile_galaxyPaths))) { + singlefile_galaxyPath=singlefile_galaxyPaths[singlefile_galaxyPath_i] + singlefile_sampleName=singlefile_sampleNames[singlefile_galaxyPath_i] + singlefile[[singlefile_sampleName]] = singlefile_galaxyPath + } + } + return(list(zipfile=zipfile, singlefile=singlefile)) +} + + +# This function retrieve the raw file in the working directory +# - if zipfile: unzip the file with its directory tree +# - if singlefiles: set symlink with the good filename +retrieveRawfileInTheWorkingDirectory <- function(singlefile, zipfile) { + + if(!is.null(singlefile) && (length("singlefile")>0)) { + for (singlefile_sampleName in names(singlefile)) { + singlefile_galaxyPath = singlefile[[singlefile_sampleName]] + if(!file.exists(singlefile_galaxyPath)){ + error_message=paste("Cannot access the sample:",singlefile_sampleName,"located:",singlefile_galaxyPath,". Please, contact your administrator ... if you have one!") + print(error_message); stop(error_message) + } + + file.symlink(singlefile_galaxyPath,singlefile_sampleName) + } + directory = "." + + } + if(!is.null(zipfile) && (zipfile!="")) { + if(!file.exists(zipfile)){ + error_message=paste("Cannot access the Zip file:",zipfile,". Please, contact your administrator ... if you have one!") + print(error_message) + stop(error_message) + } + + #list all file in the zip file + #zip_files=unzip(zipfile,list=T)[,"Name"] + + #unzip + suppressWarnings(unzip(zipfile, unzip="unzip")) + + #get the directory name + filesInZip=unzip(zipfile, list=T); + directories=unique(unlist(lapply(strsplit(filesInZip$Name,"/"), function(x) x[1]))); + directories=directories[!(directories %in% c("__MACOSX")) & file.info(directories)$isdir] + directory = "." + if (length(directories) == 1) directory = directories + + cat("files_root_directory\t",directory,"\n") + + } +}
--- a/macros.xml Mon Jul 04 11:58:10 2016 -0400 +++ b/macros.xml Fri Apr 07 07:14:12 2017 -0400 @@ -16,29 +16,102 @@ <token name="@COMMAND_CAMERA_SCRIPT@"> LANG=C Rscript $__tool_directory__/probmetab.r </token> - - <!-- zipfile load for planemo test --> - <token name="@COMMAND_ZIPFILE_LOAD@"> - #if $zipfile_load_conditional.zipfile_load_select == "yes": - #if $zipfile_load_conditional.zip_file: - zipfile $zipfile_load_conditional.zip_file + + <token name="@COMMAND_LOG_EXIT@"> + ; + return=\$?; + cat 'log.txt'; + sh -c "exit \$return" + </token> + + <!-- raw file load for planemo test --> + <token name="@COMMAND_FILE_LOAD_NEUTRAL@"> + #if $file_load_section_selected.file_load_conditional.file_load_select == "yes": + #if $file_load_section_selected.file_load_conditional.input[0].is_of_type("mzxml") or $file_load_section_selected.file_load_conditional.input[0].is_of_type("mzml") or $file_load_section_selected.file_load_conditional.input[0].is_of_type("mzdata") or $file_load_section_selected.file_load_conditional.input[0].is_of_type("netcdf"): + #set singlefile_galaxyPath = ','.join( [ str( $single_file ) for $single_file in $file_load_section_selected.file_load_conditional.input ] ) + #set singlefile_sampleName = ','.join( [ str( $single_file.name ) for $single_file in $file_load_section_selected.file_load_conditional.input ] ) + singlefile_galaxyPath$polarity '$singlefile_galaxyPath' singlefile_sampleName$polarity '$singlefile_sampleName' + #else + zipfile$polarity '$file_load_section_selected.file_load_conditional.input' #end if - #end if + #end if + </token> + + <token name="@COMMAND_FILE_LOAD_ONE@"> + #set file_load_section_selected = $acquisition_options.file_load_section + #set polarity="" + @COMMAND_FILE_LOAD_NEUTRAL@ + </token> + + <token name="@COMMAND_FILE_LOAD_POSITIVE@"> + #set file_load_section_selected = $acquisition_options.file_load_sectionPositive + #set polarity="Positive" + @COMMAND_FILE_LOAD_NEUTRAL@ + </token> + + <token name="@COMMAND_FILE_LOAD_NEGATIVE@"> + #set file_load_section_selected = $acquisition_options.file_load_sectionNegative + #set polarity="Negative" + @COMMAND_FILE_LOAD_NEUTRAL@ </token> - <xml name="zipfile_load"> - <conditional name="zipfile_load_conditional"> - <param name="zipfile_load_select" type="select" label="Resubmit your zip file" help="Use only if you get a message which say that your original zip file have been deleted on the server." > - <option value="no" >no need</option> - <option value="yes" selected="peakgroups">yes</option> - </param> - <when value="no"> - </when> - <when value="yes"> - <param name="zip_file" type="data" format="no_unzip.zip,zip" label="Zip file" /> - </when> - </conditional> + + <xml name="input_file_load" token_polarity=""> + <section name="file_load_section@POLARITY@" title="@POLARITY@ Resubmit your raw dataset or your zip file"> + <conditional name="file_load_conditional"> + <param name="file_load_select" type="select" label="@POLARITY@ Resubmit your dataset or your zip file" help="Use only if you get a message which say that your original dataset or zip file have been deleted on the server." > + <option value="no" >no need</option> + <option value="yes" >yes</option> + </param> + <when value="no"> + </when> + <when value="yes"> + <param name="input" type="data" format="mzxml,mzml,mzdata,netcdf,no_unzip.zip,zip" multiple="true" label="File(s) from your history containing your chromatograms" help="Single file mode for the format: mzxml, mzml, mzdata and netcdf. Zip file mode for the format: no_unzip.zip, zip. See the help section below." /> + </when> + </conditional> + </section> </xml> - + + <xml name="test_commun"> + <section name="getannot"> + <param name="allowMiss" value="TRUE" /> + <conditional name="option_toexclude"> + <param name="option" value="hide" /> + </conditional> + </section> + <section name="db"> + <param name="kegg_db" value="KEGG" /> + <param name="ppm_tol" value="8" /> + </section> + <section name="export"> + <param name="prob" value="count" /> + <param name="html" value="FALSE" /> + </section> + <section name="reac2cor"> + <param name="opt" value="cor" /> + <param name="corprob" value="0.8" /> + <param name="pcorprob" value="0.8" /> + <param name="corths" value="0.75" /> + </section> + </xml> + + <xml name="test_file_load_zip"> + <section name="file_load_section"> + <conditional name="file_load_conditional"> + <param name="file_load_select" value="yes" /> + <param name="input" value="faahKO_reduce.zip" ftype="zip" /> + </conditional> + </section> + </xml> + + <xml name="test_file_load_single" token_polarity="">> + <section name="file_load_section@POLARITY@"> + <conditional name="file_load_conditional"> + <param name="file_load_select" value="yes" /> + <param name="input" value="wt15.CDF,ko16.CDF,ko15.CDF,wt16.CDF" ftype="netcdf" /> + </conditional> + </section> + </xml> + <token name="@HELP_AUTHORS@"> .. class:: infomark @@ -49,7 +122,7 @@ .. class:: infomark -**Galaxy integration** Misharl Monsoor misharl.monsoor@sb-roscoff.fr ABIMS TEAM, Station biologique de Roscoff. +**Galaxy integration** Misharl Monsoor misharl.monsoor@sb-roscoff.fr (and Gildas Le Corguillé) from ABIMS TEAM, Station biologique de Roscoff. | Contact support@workflow4metabolomics.org for any questions or concerns about the Galaxy implementation of this tool.
--- a/probmetab.r Mon Jul 04 11:58:10 2016 -0400 +++ b/probmetab.r Fri Apr 07 07:14:12 2017 -0400 @@ -4,7 +4,7 @@ # ----- LOG ----- -log_file=file("probmetab.log", open = "wt") +log_file=file("log.txt", open = "wt") sink(log_file) sink(log_file, type = "out") @@ -12,8 +12,8 @@ cat("\tPACKAGE INFO\n") pkgs=c("parallel","BiocGenerics", "Biobase", "Rcpp", "mzR", "igraph", "xcms","snow","CAMERA","batch","ProbMetab") for(p in pkgs) { - suppressWarnings( suppressPackageStartupMessages( stopifnot( library(p, quietly=TRUE, logical.return=TRUE, character.only=TRUE)))) - cat(p,"\t",as.character(packageVersion(p)),"\n",sep="") + suppressWarnings( suppressPackageStartupMessages( stopifnot( library(p, quietly=TRUE, logical.return=TRUE, character.only=TRUE)))) + cat(p,"\t",as.character(packageVersion(p)),"\n",sep="") } source_local <- function(fname){ @@ -21,9 +21,14 @@ base_dir <- dirname(substring(argv[grep("--file=", argv)], 8)) source(paste(base_dir, fname, sep="/")) } + +options(bitmapType='cairo') + cat("\n\n") + + # ----- ARGUMENTS ----- -cat("\tARGUMENTS INFO\n") +cat("\tARGUMENTS INFO\n") listArguments = parseCommandArgs(evaluate=FALSE) #interpretation of arguments given in command line as an R list of objects write.table(as.matrix(listArguments), col.names=F, quote=F, sep='\t') @@ -35,68 +40,59 @@ # ----- INFILE PROCESSING ----- if(listArguments[["mode_acquisition"]]=="one") { - load(listArguments[["xa"]]) - - if (!is.null(listArguments[["zipfile"]])){ - zipfile= listArguments[["zipfile"]]; listArguments[["zipfile"]]=NULL - } - - #Unzip the chromatograms file for plotting EIC pour the HTML file - if(exists("zipfile")) - { - if (zipfile!=""){ - directory=unzip(zipfile) - } - } - if (!exists("xa")) { - xa=xsAnnotate_object - } - source_local("lib.r") - if (!exists("variableMetadata")) variableMetadata= getVariableMetadata(xa); - + load(listArguments[["image"]]) + cat("\t\tXA OBJECT INFO\n") + print(xa) + + source_local("lib.r") + if (!exists("zipfile")) zipfile=NULL + if (!exists("singlefile")) singlefile=NULL + rawFilePath = getRawfilePathFromArguments(singlefile, zipfile, listArguments) + zipfile = rawFilePath$zipfile + singlefile = rawFilePath$singlefile + retrieveRawfileInTheWorkingDirectory(singlefile, zipfile) + + if (!exists("variableMetadata")) variableMetadata= getVariableMetadata(xa); + } else if(listArguments[["inputs_mode"]]=="two"){ - load(listArguments[["image_pos"]]) - - if (!is.null(listArguments[["zipfile"]])){ - zipfile= listArguments[["zipfile"]]; listArguments[["zipfile"]]=NULL - } - - #Unzip the chromatograms file for plotting EIC pour the HTML file - if(exists("zipfile")) { - if (zipfile!=""){ - directory=unzip(zipfile) - } - } - if (!exists("xa")) { - xa=xsAnnotate_object - } - xaP=xa - source_local("lib.r") - if (!exists("variableMetadata")) variableMetadataP= getVariableMetadata(xa) - else variableMetadataP=variableMetadata + + # POSITIVE + load(listArguments[["image_pos"]]) + if (!exists("xa")) xaP=xsAnnotate_object + else xaP=xa + cat("\t\tXA-POSITIVE OBJECT INFO\n") + print(xaP) + + if (!exists("variableMetadata")) variableMetadataP = getVariableMetadata(xa) + else variableMetadataP = variableMetadata - load(listArguments[["image_neg"]]) - - if (!is.null(listArguments[["zipfile"]])){ - zipfile= listArguments[["zipfile"]]; listArguments[["zipfile"]]=NULL - } - - #Unzip the chromatograms file for plotting EIC pour the HTML file - if(exists("zipfile")) { - - if (zipfile!=""){ - directory=unzip(zipfile) - } - } - if (!exists("xa")) { - xa=xsAnnotate_object - } - xaN=xa - source_local("lib.r") - - if (!exists("variableMetadata")) variableMetadataN= getVariableMetadata(xa) - else variableMetadataN=variableMetadata + source_local("lib.r") + if (!exists("zipfile")) zipfile=NULL + if (!exists("singlefile")) singlefile=NULL + rawFilePath = getRawfilePathFromArguments(singlefile, zipfile, listArguments) + zipfilePos = rawFilePath$zipfile + singlefilePos = rawFilePath$singlefile + retrieveRawfileInTheWorkingDirectory(singlefilePos, zipfilePos) + + # NEGATIVE + load(listArguments[["image_neg"]]) + + if (!exists("xa")) xaN=xsAnnotate_object + else xaN=xa + cat("\t\tXA-NEGATIVE OBJECT INFO\n") + print(xaP) + + if (!exists("variableMetadata")) variableMetadataN = getVariableMetadata(xa) + else variableMetadataN = variableMetadata + + source_local("lib.r") + if (!exists("zipfile")) zipfile=NULL + if (!exists("singlefile")) singlefile=NULL + rawFilePath = getRawfilePathFromArguments(singlefile, zipfile, listArguments) + zipfileNeg = rawFilePath$zipfile + singlefileNeg = rawFilePath$singlefile + retrieveRawfileInTheWorkingDirectory(singlefileNeg, zipfileNeg) } #Import the different functions @@ -107,13 +103,12 @@ cat("\tMAIN PROCESSING INFO\n") if(listArguments[["mode_acquisition"]]=="one") { - results=probmetab(xa=xa, variableMetadata=variableMetadata,listArguments=listArguments) + results=probmetab(xa=xa, variableMetadata=variableMetadata,listArguments=listArguments) } else if(listArguments[["inputs_mode"]]=="two"){ - results=probmetab(xaP=xaP, xaN=xaN,variableMetadataP=variableMetadataP, variableMetadataN=variableMetadataN, listArguments=listArguments) + results=probmetab(xaP=xaP, xaN=xaN,variableMetadataP=variableMetadataP, variableMetadataN=variableMetadataN, listArguments=listArguments) } #delete the parameters to avoid the passage to the next tool in .RData image #rm(listArguments) cat("\tDONE\n") #saving R data in .Rdata file to save the variables used in the present tool #save.image(paste("probmetab","RData",sep=".")) -