Mercurial > repos > lecorguille > xcms_retcor
changeset 10:8828cba9aedd draft
planemo upload for repository https://github.com/workflow4metabolomics/xcms commit 9f72e947d9c241d11221cad561f3525d27231857
| author | lecorguille |
|---|---|
| date | Tue, 18 Sep 2018 16:12:29 -0400 |
| parents | e4e0254a3c0a |
| children | 063b2f99c266 |
| files | README.rst abims_xcms_retcor.xml lib-xcms3.x.x.r lib.r macros.xml macros_xcms.xml static/images/xcms_retcor_workflow.png test-data/faahKO-single-class.xset.group.RData test-data/faahKO.xset.group.RData xcms.r xcms_retcor.r |
| diffstat | 11 files changed, 1066 insertions(+), 913 deletions(-) [+] |
line wrap: on
line diff
--- a/README.rst Wed Nov 29 09:46:41 2017 -0500 +++ b/README.rst Tue Sep 18 16:12:29 2018 -0400 @@ -2,26 +2,42 @@ Changelog/News -------------- +**Version 3.0.0.0 - 08/03/2018** + +- UPGRADE: upgrade the xcms version from 1.46.0 to 3.0.0. So refactoring of a lot of underlying codes and methods. Some parameters may have been renamed. + +- NEW: a bunch of new options: Obiwarp.(centerSample, response, distFun, gapInit, gapExtend, factorDiag, factorGap, localAlignment, initPenalty) + +- IMPROVEMENT: the advanced options are now in sections. It will allow you to access to all the parameters and to know their default values. + +- CHANGE: removing of the TIC and BPC plots. You can new use the dedicated tool "xcms plot chromatogram" + + **Version 2.1.1 - 29/11/2017** - BUGFIX: To avoid issues with accented letter in the parentFile tag of the mzXML files, we changed a hidden mechanim to LC_ALL=C + **Version 2.1.0 - 03/02/2017** - IMPROVEMENT: xcms.retcor can deal with merged individual data + **Version 2.0.8 - 22/12/2016** - BUGFIX: when having only one group (i.e. one folder of raw data) the BPC and TIC pdf files do not contain any graph + **Version 2.0.7 - 06/07/2016** - UPGRADE: upgrate the xcms version from 1.44.0 to 1.46.0 + **Version 2.0.6 - 04/04/2016** - TEST: refactoring to pass planemo test using conda dependencies + **Version 2.0.5 - 10/02/2016** - BUGFIX: better management of errors. Datasets remained green although the process failed
--- a/abims_xcms_retcor.xml Wed Nov 29 09:46:41 2017 -0500 +++ b/abims_xcms_retcor.xml Tue Sep 18 16:12:29 2018 -0400 @@ -1,36 +1,41 @@ -<tool id="abims_xcms_retcor" name="xcms.retcor" version="2.1.1"> +<tool id="abims_xcms_retcor" name="xcms adjustRtime (retcor)" version="@WRAPPER_VERSION@.0"> - <description>Retention Time Correction using retcor function from xcms R package </description> + <description>Retention Time Correction</description> <macros> <import>macros.xml</import> + <import>macros_xcms.xml</import> </macros> <expand macro="requirements"/> <expand macro="stdio"/> <command><![CDATA[ - @COMMAND_XCMS_SCRIPT@ + @COMMAND_RSCRIPT@/xcms_retcor.r image '$image' - xfunction retcor - - xsetRdataOutput '$xsetRData' - ticspdf '$ticsCorPdf' - bicspdf '$bpcsCorPdf' - rplotspdf '$rplotsPdf' method $methods.method - #if $methods.method == "obiwarp": - profStep $methods.profStep + #if $methods.method == "PeakGroups": + minFraction $methods.minFraction + extraPeaks $methods.extraPeaks + smooth $methods.smooth_cond.smooth + ## PeakGroupsSmoothLoess Advanced + span $methods.smooth_cond.PeakGroupsSmoothLoessAdv.span + family $methods.smooth_cond.PeakGroupsSmoothLoessAdv.family #else - smooth $methods.smooth - extra $methods.extra - missing $methods.missing - #if $methods.options.option == "show": - span $methods.options.span - family $methods.options.family - plottype $methods.options.plottype + binSize $methods.binSize + ## Advanced + #if $methods.ObiwarpAdv.centerSample != "": + centerSample $methods.ObiwarpAdv.centerSample #end if + response $methods.ObiwarpAdv.response + distFun $methods.ObiwarpAdv.distFunCond.distFun + gapInit $methods.ObiwarpAdv.distFunCond.gapInit + gapExtend $methods.ObiwarpAdv.distFunCond.gapExtend + factorDiag $methods.ObiwarpAdv.factorDiag + factorGap $methods.ObiwarpAdv.factorGap + localAlignment $methods.ObiwarpAdv.localAlignmentCond.localAlignment + initPenalty $methods.ObiwarpAdv.localAlignmentCond.initPenalty #end if @COMMAND_FILE_LOAD@ @@ -39,46 +44,81 @@ ]]></command> <inputs> - <param name="image" type="data" format="rdata.xcms.raw,rdata.xcms.group,rdata" label="xset RData file" help="output file from another function xcms (xcmsSet, retcor etc.)" /> + <param name="image" type="data" format="rdata.xcms.findchrompeaks,rdata.xcms.group,rdata" label="@INPUT_IMAGE_LABEL@" help="@INPUT_IMAGE_HELP@ from: findChromPeaks, groupChromPeaks" /> <conditional name="methods"> - <param name="method" type="select" label="Method to use for retention time correction" help="[method] See the help section below" > - <option value="obiwarp" >obiwarp</option> - <option value="peakgroups" selected="true">peakgroups</option> + <param name="method" type="select" label="Method to use for retention time correction" help="See the help section below" > + <option value="PeakGroups" selected="true">PeakGroups - retention time correction based on aligment of features (peak groups) present in most/all samples.</option> + <option value="Obiwarp">Obiwarp - alignment based on the complete mz-rt data.</option> </param> - <when value="obiwarp"> - <param name="profStep" type="float" value="1" label="Step size (in m/z)" help="[profStep] to use for profile generation from the raw data files" /> + <when value="PeakGroups"> + <param argument="minFraction" type="float" value="0.9" min="0" max="1" label="Minimum required fraction of samples in which peaks for the peak group were identified" help="(previously missing)"/> + <param argument="extraPeaks" type="integer" value="1" label="Maximal number of additional peaks for all samples to be assigned to a peak group for retention time correction" help="For a data set with 6 samples, ‘extraPeaks = 1’ uses all peak groups with a total peak count lower or equal to ‘6 + 1’. The total peak count is the total number of peaks being assigned to a peak group and considers also multiple peaks within a sample being assigned to the group. (previously extra)" /> + <conditional name="smooth_cond"> + <param argument="smooth" type="select" label="Smooth method" > + <option value="loess" selected="true">loess - non-linear alignment</option> + <option value="linear">linear - linear alignment</option> + </param> + <when value="loess"> + <section name="PeakGroupsSmoothLoessAdv" title="Advanced Options" expanded="False"> + <param argument="span" type="float" value="0.2" label="Degree of smoothing for the loess fitting" /> + <param argument="family" type="select" label="Family" help="if gaussian fitting is by least-squares with no outlier removal, and if symmetric a re descending M estimator is used with Tukey's biweight function, allowing outlier removal"> + <option value="gaussian" selected="true">gaussian</option> + <option value="symmetric">symmetric</option> + </param> + </section> + </when> + <when value="linear" /> + </conditional> </when> - <when value="peakgroups"> - <param name="smooth" type="select" label="Smooth method" help="[smooth] either 'loess’ for non-linear alignment or ‘linear’ for linear alignment" > - <option value="loess">loess</option> - <option value="linear">linear</option> - </param> - <param name="extra" type="integer" value="1" label="Number of extra peaks to allow in retention time correction correction groups" help="[extra]" /> - <param name="missing" type="integer" value="1" label="Number of missing samples to allow in retention time correction groups" help="[missing]" /> - - <conditional name="options"> - <param name="option" type="select" label="Advanced options"> - <option value="show">show</option> - <option value="hide" selected="true">hide</option> - </param> - <when value="show"> - <param name="span" type="float" value="0.2" label="Degree of smoothing for local polynomial regression fitting" help="[span]"/> - - <param name="family" type="select" label="Family" help="[family] if gaussian fitting is by least-squares with no outlier removal, and if symmetric a re descending M estimator is used with Tukey's biweight function, allowing outlier removal"> - <option value="gaussian" selected="true">gaussian</option> - <option value="symmetric">symmetric</option> + <when value="Obiwarp"> + <param argument="binSize" type="float" value="1" label="Bin size (in mz dimension) to be used for the profile matrix generation" help="See ‘step’ parameter in ‘profile-matrix’ documentation for more details. (previously profStep)" /> + <section name="ObiwarpAdv" title="Advanced Options" expanded="False"> + <param argument="centerSample" type="integer" value="" optional="true" label="Index of the center sample in the experiment" help="It defaults to ‘floor(median(1:length(fileNames(object))))’" /> + <param argument="response" type="integer" value="1" label="Defining the responsiveness of warping" help="with ‘response = 0’ giving linear warping on start and end points and ‘response = 100’ warping using all bijective anchors." /> + <conditional name="distFunCond"> + <param argument="distFun" type="select" label="Distance function to be used"> + <option value="cor_opt" selected="true">cor_opt - calculate only 10% diagonal band of distance matrix; better runtime</option> + <option value="cor">cor - Pearson's correlation</option> + <option value="cov">cov - covariance</option> + <option value="prd">prd - product</option> + <option value="euc">euc - Euclidian distance</option> </param> - - <param name="plottype" type="select" help="[plottype] if deviation plot retention time deviation points and regression fit, and if mdevden also plot peak overall peak density and retention time correction peak density"> - <option value="none" selected="true">none</option> - <option value="deviation">deviation</option> - <option value="mdevden">mdevden</option> + <when value="cor_opt"> + <param argument="gapInit" type="float" value="0.3" label="Penalty for gap opening" /> + <param argument="gapExtend" type="float" value="2.4" label="Penalty for gap enlargement" /> + </when> + <when value="cor"> + <param argument="gapInit" type="float" value="0.3" label="Penalty for gap opening" /> + <param argument="gapExtend" type="float" value="2.4" label="Penalty for gap enlargement" /> + </when> + <when value="cov"> + <param argument="gapInit" type="float" value="0.0" label="Penalty for gap opening" /> + <param argument="gapExtend" type="float" value="11.7" label="Penalty for gap enlargement" /> + </when> + <when value="prd"> + <param argument="gapInit" type="float" value="0.0" label="Penalty for gap opening" /> + <param argument="gapExtend" type="float" value="7.8" label="Penalty for gap enlargement" /> + </when> + <when value="euc"> + <param argument="gapInit" type="float" value="0.9" label="Penalty for gap opening" /> + <param argument="gapExtend" type="float" value="1.8" label="Penalty for gap enlargement" /> + </when> + </conditional> + <param argument="factorDiag" type="float" value="2" label="Local weight applied to diagonal moves in the alignment" /> + <param argument="factorGap" type="float" value="1" label="local weight for gap moves in the alignment" /> + <conditional name="localAlignmentCond"> + <param argument="localAlignment" type="select" label="Whether a local alignment should be performed instead of the default global alignment"> + <option value="FALSE" selected="true">FALSE</option> + <option value="TRUE">TRUE</option> </param> - - </when> - <when value="hide"> - </when> - </conditional> + <when value="FALSE"> + <param argument="initPenalty" type="hidden" value="0" label="Penalty for initiating an alignment" /> + </when> + <when value="TRUE"> + <param argument="initPenalty" type="float" value="0" label="Penalty for initiating an alignment" /> + </when> + </conditional> + </section> </when> </conditional> @@ -87,95 +127,98 @@ </inputs> <outputs> - <data name="xsetRData" format="rdata.xcms.retcor" label="${image.name[:-6]}.retcor.RData" /> - <data name="rplotsPdf" format="pdf" label="${image.name[:-6]}.retcor.Rplots.pdf"> - <filter>(methods['method'] == 'peakgroups')</filter> - <filter>(options['option'] == 'show')</filter> - <filter>(family == 'symmetric')</filter> - <filter>(plottype != 'none')</filter> - </data> - <data name="ticsCorPdf" format="pdf" label="${image.name[:-6]}.retcor.TICs_corrected.pdf" /> - <data name="bpcsCorPdf" format="pdf" label="${image.name[:-6]}.retcor.BPCs_corrected.pdf" /> - <data name="log" format="txt" label="xset.log.txt" hidden="true" /> + <data name="xsetRData" format="rdata.xcms.retcor" label="${image.name[:-6]}.adjustRtime.RData" from_work_dir="retcor.RData" /> + <data name="rawVSadjustedPdf" format="pdf" label="${image.name[:-6]}_rawVSadjusted.adjustRtime.Rplots.pdf" from_work_dir="raw_vs_adjusted_rt.pdf" /> </outputs> <tests> - <!--<test> - <param name="image" value="xset.group.RData"/> - <param name="methods|method" value="peakgroups"/> - <param name="methods|smooth" value="loess"/> - <param name="methods|extra" value="1"/> - <param name="methods|missing" value="1"/> - <param name="methods|options|option" value="show"/> - <param name="methods|options|span" value="0.2"/> - <param name="methods|options|family" value="gaussian"/> - <param name="methods|options|plottype" value="deviation"/> - <param name="zipfile_load_conditional|zipfile_load_select" value="yes" /> - <param name="zipfile_load_conditional|zip_file" value="sacuri_dir_root.zip" ftype="zip" /> - <output name="log"> - <assert_contents> - <has_text text="object with 4 samples" /> - <has_text text="Time range: 0.2-1140.1 seconds (0-19 minutes)" /> - <has_text text="Mass range: 50.0021-999.9863 m/z" /> - <has_text text="Peaks: 59359 (about 14840 per sample)" /> - <has_text text="Peak Groups: 0" /> - <has_text text="Sample classes: bio, blank" /> - </assert_contents> - </output> - </test>--> <test> - <param name="image" value="faahKO.xset.group.RData"/> - <param name="methods|method" value="peakgroups"/> - <param name="methods|smooth" value="loess"/> - <param name="methods|extra" value="1"/> - <param name="methods|missing" value="1"/> - <param name="methods|options|option" value="show"/> - <param name="methods|options|span" value="0.2"/> - <param name="methods|options|family" value="gaussian"/> - <param name="methods|options|plottype" value="deviation"/> + <param name="image" value="faahKO-single-class.xset.group.RData" ftype="rdata"/> + <conditional name="methods"> + <param name="method" value="PeakGroups"/> + <param name="extraPeaks" value="1"/> + <param name="minFraction" value="1"/> + <conditional name="smooth_cond"> + <param name="smooth" value="loess"/> + <section name="PeakGroupsSmoothLoessAdv"> + <param name="span" value="0.2"/> + <param name="family" value="gaussian"/> + </section> + </conditional> + </conditional> + <expand macro="test_file_load_single"/> + <assert_stdout> + <has_text text="extraPeaks: 1" /> + <has_text text="minFraction: 1" /> + <has_text text="span: 0.2" /> + <has_text text="object with 4 samples" /> + <has_text text="Time range: 2509.2-4480.3 seconds (41.8-74.7 minutes)" /> + <has_text text="Mass range: 200.1-600 m/z" /> + <has_text text="Peaks: 9251 (about 2313 per sample)" /> + <has_text text="Peak Groups: 0" /> + <has_text text="Sample classes: KO, WT" /> + </assert_stdout> + </test> + <!-- DISABLE FOR TRAVIS : Zip + <test> + <param name="image" value="faahKO.xset.group.RData" ftype="rdata"/> + <conditional name="methods"> + <param name="method" value="PeakGroups"/> + <param name="extraPeaks" value="1"/> + <param name="minFraction" value="1"/> + <conditional name="smooth_cond"> + <param name="smooth" value="loess"/> + <section name="PeakGroupsSmoothLoessAdv"> + <param name="span" value="0.2"/> + <param name="family" value="gaussian"/> + </section> + </conditional> + </conditional> <expand macro="test_file_load_zip"/> - <output name="log"> - <assert_contents> - <has_text text="object with 4 samples" /> - <has_text text="Time range: 2507.7-4481.7 seconds (41.8-74.7 minutes)" /> - <has_text text="Mass range: 200.1-600 m/z" /> - <has_text text="Peaks: 9251 (about 2313 per sample)" /> - <has_text text="Peak Groups: 0" /> - <has_text text="Sample classes: KO, WT" /> - </assert_contents> - </output> + <assert_stdout> + <has_text text="object with 4 samples" /> + <has_text text="Time range: 2509.2-4480.3 seconds (41.8-74.7 minutes)" /> + <has_text text="Mass range: 200.1-600 m/z" /> + <has_text text="Peaks: 9251 (about 2313 per sample)" /> + <has_text text="Peak Groups: 0" /> + <has_text text="Sample classes: KO, WT" /> + </assert_stdout> </test> + --> + <!-- DISABLE FOR TRAVIS + Test to test the different methods parameters <test> - <param name="image" value="faahKO-single-class.xset.group.RData"/> - <param name="methods|method" value="peakgroups"/> - <param name="methods|smooth" value="loess"/> - <param name="methods|extra" value="1"/> - <param name="methods|missing" value="1"/> - <param name="methods|options|option" value="show"/> - <param name="methods|options|span" value="0.2"/> - <param name="methods|options|family" value="gaussian"/> - <param name="methods|options|plottype" value="deviation"/> + <param name="image" value="faahKO-single-class.xset.group.RData" ftype="rdata"/> + <conditional name="methods"> + <param name="method" value="Obiwarp"/> + <section name="ObiwarpAdv"> + <param name="centerSample" value="1"/> + <param name="response" value="0"/> + <conditional name="distFunCond"> + <param name="distFun" value="cov"/> + <param name="gapInit" value="0.1" /> + </conditional> + </section> + </conditional> <expand macro="test_file_load_single"/> - <output name="log"> - <assert_contents> - <has_text text="object with 4 samples" /> - <has_text text="Time range: 2507.7-4481.7 seconds (41.8-74.7 minutes)" /> - <has_text text="Mass range: 200.1-600 m/z" /> - <has_text text="Peaks: 9251 (about 2313 per sample)" /> - <has_text text="Peak Groups: 0" /> - <has_text text="Sample classes: KO, WT" /> - </assert_contents> - </output> + <assert_stdout> + <has_text text="centerSample: 1" /> + <has_text text="response: 0" /> + <has_text text="distFun: cov" /> + <has_text text="gapInit: 0.1" /> + <has_text text="gapExtend: 11.7" /> + </assert_stdout> </test> + --> </tests> <help><![CDATA[ @HELP_AUTHORS@ -=========== -Xcms.retcor -=========== +================ +xcms adjustRtime +================ ----------- Description @@ -191,7 +234,7 @@ .. class:: warningmark -**After an retcor step, it is mandatory to do a group step, otherwise the rest of the workflow will not work with the RData file. (the initial peak grouping becomes invalid and is +**After an adjustRtime step, it is mandatory to do a groupChromPeaks step, otherwise the rest of the workflow will not work with the RData file. (the initial peak grouping becomes invalid and is discarded)** @@ -203,39 +246,29 @@ **Upstream tools** -========================= ================= ======= ========== -Name output file format parameter -========================= ================= ======= ========== -xcms.group xset.group.RData RData RData file -========================= ================= ======= ========== +========================= ============================ ============================== +Name Output file Format +========================= ============================ ============================== +xcms.findChromPeaks raw.xset.RData rdata.xcms.findchrompeaks +------------------------- ---------------------------- ------------------------------ +xcms.groupChromPeaks ``*``.groupChromPeaks.RData rdata.xcms.group +========================= ============================ ============================== **Downstream tools** -+---------------------------+------------------+--------+ -| Name | Output file | Format | -+===========================+==================+========+ -|xcms.group | xset.retcor.RData| RData | -+---------------------------+------------------+--------+ - -The output file **xset.retcor.RData** is an RData file. You can continue your analysis using it in **xcms.group** tool as an next step. +=========================== ============================ ================ +Name Output file Format +=========================== ============================ ================ +xcms.groupChromPeaks ``*``.groupChromPeaks.RData rdata.xcms.group +=========================== ============================ ================ **General schema of the metabolomic workflow** .. image:: xcms_retcor_workflow.png - ------------ -Input files ------------ - -+---------------------------+----------------------+ -| Parameter : num + label | Format | -+===========================+======================+ -| 1 : RData file | rdata.xcms.group | -+---------------------------+----------------------+ - +--------------------------------------------------- ---------- Parameters @@ -244,74 +277,28 @@ Method ------ -**peakgroups** +**PeakGroups** - | xcms ignores those groups by only considering well-behaved peak groups which are missing at most one sample and have at most one extra peak. (Those values can be changed with the **missing** and **extra** arguments.) - | For each of those well-behaved groups, the algorithm calculates a median retention time and, for every sample, a deviation from that median. Within a sample, the observed deviation generally changes over time in a nonlinear fashion. Those changes are approximated using a local polynomial regression technique implemented in the **loess** function. By default, the curve fitting is done using least-squares on all data points. - | However, it is possible to enable outlier detection and removal by setting the **family** argument to **symmetric**. + | This method performs retention time adjustment based on the alignment of chromatographic peak groups present in all/most samples (hence corresponding to house keeping compounds). First the retention time deviation of these peak groups is described by fitting either a polynomial (‘smooth = "loess"’) or a linear ( ‘smooth = "linear"’) model to the data points. These models are subsequently used to adjust the retention time of each spectrum in each sample. + | See the PeakGroups_manual_ + +**Obiwarp** -**obiwarp** + | This method performs retention time adjustment using the Obiwarp method [Prince 2006]. It is based on the code at http://obi-warp.sourceforge.net but supports alignment of multiple samples by aligning each against a _center_ sample. The alignment is performed directly on the ‘profile-matrix’ and can hence be performed independently of the peak detection or peak grouping. + | See the Obiwarp_manual_ - | Calculate retention time deviations for each sample using the obiwarp code at "http://obi-warp.sourceforge.net/". This function is able to align multiple samples by a center-star strategy. Ordered Bijective Interpolated Warping (OBI-Warp) aligns matrices along a single axis using Dynamic Time Warping (DTW) and a one-to-one (bijective) interpolated warp function. OBI-Warp harnesses the non-linear, comprehensive alignment power of DTW and builds on the discrete, non-bijective output of DTW to give natural interpolants that can be used across multiple datasets. - | For the original publication see :**Chromatographic Alignment of ESI-LC-MS Proteomics Data Sets by Ordered Bijective Interpo-lated Warping John T. Prince and, Edward M. Marcotte Analytical Chemistry 2006 78 (17), 6140-6152.** +.. _PeakGroups_manual: https://rdrr.io/bioc/xcms/man/adjustRtime-peakGroups.html#heading-2 +.. _Obiwarp_manual: https://rdrr.io/bioc/xcms/man/adjustRtime-obiwarp.html +@HELP_XCMS_MANUAL@ ------------ Output files ------------ -xset.group.retcor.TICs_corrected.pdf - - | "Total Ion Chromatograms" graph in pdf format,corrected after a retcor step. - -xset.group.retcor.BPCs_corrected.pdf - - | "Total Io"Base Peak Chromatograms" graph in pdf format,corrected after a retcor step - -xset.group.retcor.RData: rdata.xcms.retcor format - - | Rdata file that will be necessary in the **xcms.group** step of the workflow. - - ------- - -.. class:: infomark - -The output file is an xset.retcor.RData file. You can continue your analysis using it in **xcms.fillPeaks** tool. - - ---------------------------------------------------- - ---------------- -Working example ---------------- +xset.groupChromPeaks.adjustRtime.RData: rdata.xcms.retcor format -Input files ------------ - - | RData file -> **xset.group.RData** - -Parameters ----------- - - | Method: -> **peakgroups** - | smooth: -> **loess** - | extra: -> **1** - | missing -> **1** - | Advanced options: -> **show** - | span -> **0.2** - | family -> **gaussian** - | plottype -> **deviation** - - -Output files ------------- - - | **1) xset.group.retcor.RData: RData file** - - | **2) Example of an xset.group.retcor.TICs_corrected pdf file** - -.. image:: xcms_retcor.png + | Rdata file that will be necessary in the **xcms.groupChromPeaks** step of the workflow. --------------------------------------------------- @@ -319,22 +306,37 @@ Changelog/News -------------- +**Version 3.0.0.0 - 08/03/2018** + +- UPGRADE: upgrade the xcms version from 1.46.0 to 3.0.0. So refactoring of a lot of underlying codes and methods. Some parameters may have been renamed. + +- NEW: a bunch of new options: Obiwarp.(centerSample, response, distFun, gapInit, gapExtend, factorDiag, factorGap, localAlignment, initPenalty) + +- IMPROVEMENT: the advanced options are now in sections. It will allow you to access to all the parameters and to know their default values. + +- CHANGE: removing of the TIC and BPC plots. You can now use the dedicated tool "xcms plot chromatogram" + + **Version 2.1.1 - 29/11/2017** - BUGFIX: To avoid issues with accented letter in the parentFile tag of the mzXML files, we changed a hidden mechanim to LC_ALL=C + **Version 2.1.0 - 03/02/2017** - IMPROVEMENT: xcms.retcor can deal with merged individual data + **Version 2.0.8 - 22/12/2016** - BUGFIX: when having only one group (i.e. one folder of raw data) the BPC and TIC pdf files do not contain any graph + **Version 2.0.7 - 06/07/2016** - UPGRADE: upgrate the xcms version from 1.44.0 to 1.46.0 + **Version 2.0.6 - 04/04/2016** - TEST: refactoring to pass planemo test using conda dependencies
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lib-xcms3.x.x.r Tue Sep 18 16:12:29 2018 -0400 @@ -0,0 +1,152 @@ + + +#@TODO: remove this function as soon as we can use xcms 3.x.x from Bioconductor 3.7 +# https://github.com/sneumann/xcms/issues/250 +groupnamesW4M <- function(xdata, mzdec = 0, rtdec = 0) { + mzfmt <- paste("%.", mzdec, "f", sep = "") + rtfmt <- paste("%.", rtdec, "f", sep = "") + + gnames <- paste("M", sprintf(mzfmt, featureDefinitions(xdata)[,"mzmed"]), "T", + sprintf(rtfmt, featureDefinitions(xdata)[,"rtmed"]), sep = "") + + if (any(dup <- duplicated(gnames))) + for (dupname in unique(gnames[dup])) { + dupidx <- which(gnames == dupname) + gnames[dupidx] <- paste(gnames[dupidx], seq(along = dupidx), sep = "_") + } + + return (gnames) +} + +#@TODO: remove this function as soon as we can use xcms 3.x.x from Bioconductor 3.7 +# https://github.com/sneumann/xcms/issues/247 +.concatenate_XCMSnExp <- function(...) { + x <- list(...) + if (length(x) == 0) + return(NULL) + if (length(x) == 1) + return(x[[1]]) + ## Check that all are XCMSnExp objects. + if (!all(unlist(lapply(x, function(z) is(z, "XCMSnExp"))))) + stop("All passed objects should be 'XCMSnExp' objects") + new_x <- as(.concatenate_OnDiskMSnExp(...), "XCMSnExp") + ## If any of the XCMSnExp has alignment results or detected features drop + ## them! + x <- lapply(x, function(z) { + if (hasAdjustedRtime(z)) { + z <- dropAdjustedRtime(z) + warning("Adjusted retention times found, had to drop them.") + } + if (hasFeatures(z)) { + z <- dropFeatureDefinitions(z) + warning("Feature definitions found, had to drop them.") + } + z + }) + ## Combine peaks + fls <- lapply(x, fileNames) + startidx <- cumsum(lengths(fls)) + pks <- lapply(x, chromPeaks) + procH <- lapply(x, processHistory) + for (i in 2:length(fls)) { + pks[[i]][, "sample"] <- pks[[i]][, "sample"] + startidx[i - 1] + procH[[i]] <- lapply(procH[[i]], function(z) { + z@fileIndex <- as.integer(z@fileIndex + startidx[i - 1]) + z + }) + } + pks <- do.call(rbind, pks) + new_x@.processHistory <- unlist(procH) + chromPeaks(new_x) <- pks + if (validObject(new_x)) + new_x +} + +#@TODO: remove this function as soon as we can use xcms 3.x.x from Bioconductor 3.7 +# https://github.com/sneumann/xcms/issues/247 +.concatenate_OnDiskMSnExp <- function(...) { + x <- list(...) + if (length(x) == 0) + return(NULL) + if (length(x) == 1) + return(x[[1]]) + ## Check that all are XCMSnExp objects. + if (!all(unlist(lapply(x, function(z) is(z, "OnDiskMSnExp"))))) + stop("All passed objects should be 'OnDiskMSnExp' objects") + ## Check processingQueue + procQ <- lapply(x, function(z) z@spectraProcessingQueue) + new_procQ <- procQ[[1]] + is_ok <- unlist(lapply(procQ, function(z) + !is.character(all.equal(new_procQ, z)) + )) + if (any(!is_ok)) { + warning("Processing queues from the submitted objects differ! ", + "Dropping the processing queue.") + new_procQ <- list() + } + ## processingData + fls <- lapply(x, function(z) z@processingData@files) + startidx <- cumsum(lengths(fls)) + ## featureData + featd <- lapply(x, fData) + ## Have to update the file index and the spectrum names. + for (i in 2:length(featd)) { + featd[[i]]$fileIdx <- featd[[i]]$fileIdx + startidx[i - 1] + rownames(featd[[i]]) <- MSnbase:::formatFileSpectrumNames( + fileIds = featd[[i]]$fileIdx, + spectrumIds = featd[[i]]$spIdx, + nSpectra = nrow(featd[[i]]), + nFiles = length(unlist(fls)) + ) + } + featd <- do.call(rbind, featd) + featd$spectrum <- 1:nrow(featd) + ## experimentData + expdata <- lapply(x, function(z) { + ed <- z@experimentData + data.frame(instrumentManufacturer = ed@instrumentManufacturer, + instrumentModel = ed@instrumentModel, + ionSource = ed@ionSource, + analyser = ed@analyser, + detectorType = ed@detectorType, + stringsAsFactors = FALSE) + }) + expdata <- do.call(rbind, expdata) + expdata <- new("MIAPE", + instrumentManufacturer = expdata$instrumentManufacturer, + instrumentModel = expdata$instrumentModel, + ionSource = expdata$ionSource, + analyser = expdata$analyser, + detectorType = expdata$detectorType) + + ## protocolData + protodata <- lapply(x, function(z) z@protocolData) + if (any(unlist(lapply(protodata, nrow)) > 0)) + warning("Found non-empty protocol data, but merging protocol data is", + " currently not supported. Skipped.") + ## phenoData + pdata <- do.call(rbind, lapply(x, pData)) + res <- new( + "OnDiskMSnExp", + phenoData = new("NAnnotatedDataFrame", data = pdata), + featureData = new("AnnotatedDataFrame", featd), + processingData = new("MSnProcess", + processing = paste0("Concatenated [", date(), "]"), + files = unlist(fls), smoothed = NA), + experimentData = expdata, + spectraProcessingQueue = new_procQ) + if (validObject(res)) + res +} + +#@TODO: remove this function as soon as we can use xcms 3.x.x from Bioconductor 3.7 +# https://github.com/sneumann/xcms/issues/247 +c.XCMSnExp <- function(...) { + .concatenate_XCMSnExp(...) +} + +#@TODO: remove this function as soon as we can use xcms 3.x.x from Bioconductor 3.7 +# https://github.com/sneumann/xcms/issues/247 +c.MSnbase <- function(...) { + .concatenate_OnDiskMSnExp(...) +}
--- a/lib.r Wed Nov 29 09:46:41 2017 -0500 +++ b/lib.r Tue Sep 18 16:12:29 2018 -0400 @@ -1,319 +1,269 @@ -#Authors ABiMS TEAM -#Lib.r for Galaxy Workflow4Metabolomics xcms tools -# -#version 2.4: lecorguille -# add getPeaklistW4M -#version 2.3: yguitton -# correction for empty PDF when only 1 class -#version 2.2 -# correct bug in Base Peak Chromatogram (BPC) option, not only TIC when scanrange used in xcmsSet -# Note if scanrange is used a warning is prompted in R console but do not stop PDF generation -#version 2.1: yguitton -# Modifications made by Guitton Yann +#@authors ABiMS TEAM, Y. Guitton +# lib.r for Galaxy Workflow4Metabolomics xcms tools + +#@author G. Le Corguille +# solve an issue with batch if arguments are logical TRUE/FALSE +parseCommandArgs <- function(...) { + args <- batch::parseCommandArgs(...) + for (key in names(args)) { + if (args[key] %in% c("TRUE","FALSE")) + args[key] = as.logical(args[key]) + } + return(args) +} +#@author G. Le Corguille +# This function will +# - load the packages +# - display the sessionInfo +loadAndDisplayPackages <- function(pkgs) { + for(pkg in pkgs) suppressPackageStartupMessages( stopifnot( library(pkg, quietly=TRUE, logical.return=TRUE, character.only=TRUE))) + + sessioninfo = sessionInfo() + cat(sessioninfo$R.version$version.string,"\n") + cat("Main packages:\n") + for (pkg in names(sessioninfo$otherPkgs)) { cat(paste(pkg,packageVersion(pkg)),"\t") }; cat("\n") + cat("Other loaded packages:\n") + for (pkg in names(sessioninfo$loadedOnly)) { cat(paste(pkg,packageVersion(pkg)),"\t") }; cat("\n") +} + +#@author G. Le Corguille +# This function merge several chromBPI or chromTIC into one. +mergeChrom <- function(chrom_merged, chrom) { + if (is.null(chrom_merged)) return(NULL) + chrom_merged@.Data <- cbind(chrom_merged@.Data, chrom@.Data) + return(chrom_merged) +} #@author G. Le Corguille -#This function convert if it is required the Retention Time in minutes +# This function merge several xdata into one. +mergeXData <- function(args) { + chromTIC <- NULL + chromBPI <- NULL + chromTIC_adjusted <- NULL + chromBPI_adjusted <- NULL + for(image in args$images) { + + load(image) + # Handle infiles + if (!exists("singlefile")) singlefile <- NULL + if (!exists("zipfile")) zipfile <- NULL + rawFilePath <- getRawfilePathFromArguments(singlefile, zipfile, args) + zipfile <- rawFilePath$zipfile + singlefile <- rawFilePath$singlefile + retrieveRawfileInTheWorkingDirectory(singlefile, zipfile) + + if (exists("raw_data")) xdata <- raw_data + if (!exists("xdata")) stop("\n\nERROR: The RData doesn't contain any object called 'xdata'. This RData should have been created by an old version of XMCS 2.*") + + cat(sampleNamesList$sampleNamesOrigin,"\n") + + if (!exists("xdata_merged")) { + xdata_merged <- xdata + singlefile_merged <- singlefile + md5sumList_merged <- md5sumList + sampleNamesList_merged <- sampleNamesList + chromTIC_merged <- chromTIC + chromBPI_merged <- chromBPI + chromTIC_adjusted_merged <- chromTIC_adjusted + chromBPI_adjusted_merged <- chromBPI_adjusted + } else { + if (is(xdata, "XCMSnExp")) xdata_merged <- c(xdata_merged,xdata) + else if (is(xdata, "OnDiskMSnExp")) xdata_merged <- .concatenate_OnDiskMSnExp(xdata_merged,xdata) + else stop("\n\nERROR: The RData either a OnDiskMSnExp object called raw_data or a XCMSnExp object called xdata") + + singlefile_merged <- c(singlefile_merged,singlefile) + md5sumList_merged$origin <- rbind(md5sumList_merged$origin,md5sumList$origin) + sampleNamesList_merged$sampleNamesOrigin <- c(sampleNamesList_merged$sampleNamesOrigin,sampleNamesList$sampleNamesOrigin) + sampleNamesList_merged$sampleNamesMakeNames <- c(sampleNamesList_merged$sampleNamesMakeNames,sampleNamesList$sampleNamesMakeNames) + chromTIC_merged <- mergeChrom(chromTIC_merged, chromTIC) + chromBPI_merged <- mergeChrom(chromBPI_merged, chromBPI) + chromTIC_adjusted_merged <- mergeChrom(chromTIC_adjusted_merged, chromTIC_adjusted) + chromBPI_adjusted_merged <- mergeChrom(chromBPI_adjusted_merged, chromBPI_adjusted) + } + } + rm(image) + xdata <- xdata_merged; rm(xdata_merged) + singlefile <- singlefile_merged; rm(singlefile_merged) + md5sumList <- md5sumList_merged; rm(md5sumList_merged) + sampleNamesList <- sampleNamesList_merged; rm(sampleNamesList_merged) + + if (!is.null(args$sampleMetadata)) { + cat("\tXSET PHENODATA SETTING...\n") + sampleMetadataFile <- args$sampleMetadata + sampleMetadata <- getDataFrameFromFile(sampleMetadataFile, header=F) + xdata@phenoData@data$sample_group=sampleMetadata$V2[match(xdata@phenoData@data$sample_name,sampleMetadata$V1)] + + if (any(is.na(pData(xdata)$sample_group))) { + sample_missing <- pData(xdata)$sample_name[is.na(pData(xdata)$sample_group)] + error_message <- paste("Those samples are missing in your sampleMetadata:", paste(sample_missing, collapse=" ")) + print(error_message) + stop(error_message) + } + } + + if (!is.null(chromTIC_merged)) { chromTIC <- chromTIC_merged; chromTIC@phenoData <- xdata@phenoData } + if (!is.null(chromBPI_merged)) { chromBPI <- chromBPI_merged; chromBPI@phenoData <- xdata@phenoData } + if (!is.null(chromTIC_adjusted_merged)) { chromTIC_adjusted <- chromTIC_adjusted_merged; chromTIC_adjusted@phenoData <- xdata@phenoData } + if (!is.null(chromBPI_adjusted_merged)) { chromBPI_adjusted <- chromBPI_adjusted_merged; chromBPI_adjusted@phenoData <- xdata@phenoData } + + return(list("xdata"=xdata, "singlefile"=singlefile, "md5sumList"=md5sumList,"sampleNamesList"=sampleNamesList, "chromTIC"=chromTIC, "chromBPI"=chromBPI, "chromTIC_adjusted"=chromTIC_adjusted, "chromBPI_adjusted"=chromBPI_adjusted)) +} + +#@author G. Le Corguille +# This function convert if it is required the Retention Time in minutes RTSecondToMinute <- function(variableMetadata, convertRTMinute) { if (convertRTMinute){ #converting the retention times (seconds) into minutes print("converting the retention times into minutes in the variableMetadata") - variableMetadata[,"rt"]=variableMetadata[,"rt"]/60 - variableMetadata[,"rtmin"]=variableMetadata[,"rtmin"]/60 - variableMetadata[,"rtmax"]=variableMetadata[,"rtmax"]/60 + variableMetadata[,"rt"] <- variableMetadata[,"rt"]/60 + variableMetadata[,"rtmin"] <- variableMetadata[,"rtmin"]/60 + variableMetadata[,"rtmax"] <- variableMetadata[,"rtmax"]/60 } return (variableMetadata) } #@author G. Le Corguille -#This function format ions identifiers +# This function format ions identifiers formatIonIdentifiers <- function(variableMetadata, numDigitsRT=0, numDigitsMZ=0) { - splitDeco = strsplit(as.character(variableMetadata$name),"_") - idsDeco = sapply(splitDeco, function(x) { deco=unlist(x)[2]; if (is.na(deco)) return ("") else return(paste0("_",deco)) }) - namecustom = make.unique(paste0("M",round(variableMetadata[,"mz"],numDigitsMZ),"T",round(variableMetadata[,"rt"],numDigitsRT),idsDeco)) - variableMetadata=cbind(name=variableMetadata$name, namecustom=namecustom, variableMetadata[,!(colnames(variableMetadata) %in% c("name"))]) + splitDeco <- strsplit(as.character(variableMetadata$name),"_") + idsDeco <- sapply(splitDeco, function(x) { deco=unlist(x)[2]; if (is.na(deco)) return ("") else return(paste0("_",deco)) }) + namecustom <- make.unique(paste0("M",round(variableMetadata[,"mz"],numDigitsMZ),"T",round(variableMetadata[,"rt"],numDigitsRT),idsDeco)) + variableMetadata <- cbind(name=variableMetadata$name, namecustom=namecustom, variableMetadata[,!(colnames(variableMetadata) %in% c("name"))]) return(variableMetadata) } #@author G. Le Corguille +# This function convert the remain NA to 0 in the dataMatrix +naTOzeroDataMatrix <- function(dataMatrix, naTOzero) { + if (naTOzero){ + dataMatrix[is.na(dataMatrix)] <- 0 + } + return (dataMatrix) +} + +#@author G. Le Corguille +# Draw the plotChromPeakDensity 3 per page in a pdf file +getPlotChromPeakDensity <- function(xdata, mzdigit=4) { + pdf(file="plotChromPeakDensity.pdf", width=16, height=12) + + par(mfrow = c(3, 1), mar = c(4, 4, 1, 0.5)) + + group_colors <- brewer.pal(3, "Set1")[1:length(unique(xdata$sample_group))] + names(group_colors) <- unique(xdata$sample_group) + + xlim <- c(min(featureDefinitions(xdata)$rtmin), max(featureDefinitions(xdata)$rtmax)) + for (i in 1:nrow(featureDefinitions(xdata))) { + mzmin = featureDefinitions(xdata)[i,]$mzmin + mzmax = featureDefinitions(xdata)[i,]$mzmax + plotChromPeakDensity(xdata, mz=c(mzmin,mzmax), col=group_colors, pch=16, xlim=xlim, main=paste(round(mzmin,mzdigit),round(mzmax,mzdigit))) + legend("topright", legend=names(group_colors), col=group_colors, cex=0.8, lty=1) + } + + dev.off() +} + +#@author G. Le Corguille +# Draw the plotChromPeakDensity 3 per page in a pdf file +getPlotAdjustedRtime <- function(xdata) { + + pdf(file="raw_vs_adjusted_rt.pdf", width=16, height=12) + + # Color by group + group_colors <- brewer.pal(3, "Set1")[1:length(unique(xdata$sample_group))] + if (length(group_colors) > 1) { + names(group_colors) <- unique(xdata$sample_group) + plotAdjustedRtime(xdata, col = group_colors[xdata$sample_group]) + legend("topright", legend=names(group_colors), col=group_colors, cex=0.8, lty=1) + } + + # Color by sample + plotAdjustedRtime(xdata, col = rainbow(length(xdata@phenoData@data$sample_name))) + legend("topright", legend=xdata@phenoData@data$sample_name, col=rainbow(length(xdata@phenoData@data$sample_name)), cex=0.8, lty=1) + + dev.off() +} + +#@author G. Le Corguille # value: intensity values to be used into, maxo or intb -getPeaklistW4M <- function(xset, intval="into",convertRTMinute=F,numDigitsMZ=4,numDigitsRT=0,variableMetadataOutput,dataMatrixOutput) { - variableMetadata_dataMatrix = peakTable(xset, method="medret", value=intval) - variableMetadata_dataMatrix = cbind(name=groupnames(xset),variableMetadata_dataMatrix) +getPeaklistW4M <- function(xdata, intval="into", convertRTMinute=F, numDigitsMZ=4, numDigitsRT=0, naTOzero=T, variableMetadataOutput, dataMatrixOutput) { + dataMatrix <- featureValues(xdata, method="medret", value=intval) + colnames(dataMatrix) <- tools::file_path_sans_ext(colnames(dataMatrix)) + dataMatrix = cbind(name=groupnamesW4M(xdata), dataMatrix) + variableMetadata <- featureDefinitions(xdata) + colnames(variableMetadata)[1] = "mz"; colnames(variableMetadata)[4] = "rt" + variableMetadata = data.frame(name=groupnamesW4M(xdata), variableMetadata) - dataMatrix = variableMetadata_dataMatrix[,(make.names(colnames(variableMetadata_dataMatrix)) %in% c("name", make.names(sampnames(xset))))] - - variableMetadata = variableMetadata_dataMatrix[,!(make.names(colnames(variableMetadata_dataMatrix)) %in% c(make.names(sampnames(xset))))] - variableMetadata = RTSecondToMinute(variableMetadata, convertRTMinute) - variableMetadata = formatIonIdentifiers(variableMetadata, numDigitsRT=numDigitsRT, numDigitsMZ=numDigitsMZ) + variableMetadata <- RTSecondToMinute(variableMetadata, convertRTMinute) + variableMetadata <- formatIonIdentifiers(variableMetadata, numDigitsRT=numDigitsRT, numDigitsMZ=numDigitsMZ) + dataMatrix <- naTOzeroDataMatrix(dataMatrix, naTOzero) write.table(variableMetadata, file=variableMetadataOutput,sep="\t",quote=F,row.names=F) write.table(dataMatrix, file=dataMatrixOutput,sep="\t",quote=F,row.names=F) + } -#@author Y. Guitton -getBPC <- function(file,rtcor=NULL, ...) { - object <- xcmsRaw(file) - sel <- profRange(object, ...) - cbind(if (is.null(rtcor)) object@scantime[sel$scanidx] else rtcor ,xcms:::colMax(object@env$profile[sel$massidx,sel$scanidx,drop=FALSE])) - #plotChrom(xcmsRaw(file), base=T) +#@author G. Le Corguille +# It allow different of field separators +getDataFrameFromFile <- function(filename, header=T) { + myDataFrame <- read.table(filename, header=header, sep=";", stringsAsFactors=F) + if (ncol(myDataFrame) < 2) myDataFrame <- read.table(filename, header=header, sep="\t", stringsAsFactors=F) + if (ncol(myDataFrame) < 2) myDataFrame <- read.table(filename, header=header, sep=",", stringsAsFactors=F) + if (ncol(myDataFrame) < 2) { + error_message="Your tabular file seems not well formatted. The column separators accepted are ; , and tabulation" + print(error_message) + stop(error_message) + } + return(myDataFrame) } -#@author Y. Guitton -getBPCs <- function (xcmsSet=NULL, pdfname="BPCs.pdf",rt=c("raw","corrected"), scanrange=NULL) { - cat("Creating BIC pdf...\n") +#@author G. Le Corguille +# Draw the BPI and TIC graphics +# colored by sample names or class names +getPlotChromatogram <- function(chrom, xdata, pdfname="Chromatogram.pdf", aggregationFun = "max") { - if (is.null(xcmsSet)) { - cat("Enter an xcmsSet \n") - stop() - } else { - files <- filepaths(xcmsSet) - } - - phenoDataClass<-as.vector(levels(xcmsSet@phenoData[,1])) #sometime phenoData have more than 1 column use first as class - - classnames<-vector("list",length(phenoDataClass)) - for (i in 1:length(phenoDataClass)){ - classnames[[i]]<-which( xcmsSet@phenoData[,1]==phenoDataClass[i]) - } + if (aggregationFun == "sum") + type="Total Ion Chromatograms" + else + type="Base Peak Intensity Chromatograms" - N <- dim(phenoData(xcmsSet))[1] + adjusted="Raw" + if (hasAdjustedRtime(xdata)) + adjusted="Adjusted" - TIC <- vector("list",N) - - - for (j in 1:N) { + main <- paste(type,":",adjusted,"data") - TIC[[j]] <- getBPC(files[j]) - #good for raw - # seems strange for corrected - #errors if scanrange used in xcmsSetgeneration - if (!is.null(xcmsSet) && rt == "corrected") - rtcor <- xcmsSet@rt$corrected[[j]] - else - rtcor <- NULL + pdf(pdfname, width=16, height=10) - TIC[[j]] <- getBPC(files[j],rtcor=rtcor) - # TIC[[j]][,1]<-rtcor + # Color by group + group_colors <- brewer.pal(3, "Set1")[1:length(unique(xdata$sample_group))] + if (length(group_colors) > 1) { + names(group_colors) <- unique(xdata$sample_group) + plot(chrom, col = group_colors[chrom$sample_group], main=main) + legend("topright", legend=names(group_colors), col=group_colors, cex=0.8, lty=1) } - - - pdf(pdfname,w=16,h=10) - cols <- rainbow(N) - lty = 1:N - pch = 1:N - #search for max x and max y in BPCs - xlim = range(sapply(TIC, function(x) range(x[,1]))) - ylim = range(sapply(TIC, function(x) range(x[,2]))) - ylim = c(-ylim[2], ylim[2]) - - - ##plot start - - if (length(phenoDataClass)>2){ - for (k in 1:(length(phenoDataClass)-1)){ - for (l in (k+1):length(phenoDataClass)){ - #print(paste(phenoDataClass[k],"vs",phenoDataClass[l],sep=" ")) - plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",phenoDataClass[k]," vs ",phenoDataClass[l], sep=""), xlab = "Retention Time (min)", ylab = "BPC") - colvect<-NULL - for (j in 1:length(classnames[[k]])) { - tic <- TIC[[classnames[[k]][j]]] - # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") - points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") - colvect<-append(colvect,cols[classnames[[k]][j]]) - } - for (j in 1:length(classnames[[l]])) { - # i=class2names[j] - tic <- TIC[[classnames[[l]][j]]] - points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") - colvect<-append(colvect,cols[classnames[[l]][j]]) - } - legend("topright",paste(basename(files[c(classnames[[k]],classnames[[l]])])), col = colvect, lty = lty, pch = pch) - } - } - }#end if length >2 + # Color by sample + plot(chrom, col = rainbow(length(xdata@phenoData@data$sample_name)), main=main) + legend("topright", legend=xdata@phenoData@data$sample_name, col=rainbow(length(xdata@phenoData@data$sample_name)), cex=0.8, lty=1) - if (length(phenoDataClass)==2){ - k=1 - l=2 - colvect<-NULL - plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",phenoDataClass[k],"vs",phenoDataClass[l], sep=""), xlab = "Retention Time (min)", ylab = "BPC") - - for (j in 1:length(classnames[[k]])) { - - tic <- TIC[[classnames[[k]][j]]] - # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") - points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") - colvect<-append(colvect,cols[classnames[[k]][j]]) - } - for (j in 1:length(classnames[[l]])) { - # i=class2names[j] - tic <- TIC[[classnames[[l]][j]]] - points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") - colvect<-append(colvect,cols[classnames[[l]][j]]) - } - legend("topright",paste(basename(files[c(classnames[[k]],classnames[[l]])])), col = colvect, lty = lty, pch = pch) - - }#end length ==2 - - #case where only one class - if (length(phenoDataClass)==1){ - k=1 - ylim = range(sapply(TIC, function(x) range(x[,2]))) - colvect<-NULL - plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",phenoDataClass[k], sep=""), xlab = "Retention Time (min)", ylab = "BPC") - - for (j in 1:length(classnames[[k]])) { - tic <- TIC[[classnames[[k]][j]]] - # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") - points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") - colvect<-append(colvect,cols[classnames[[k]][j]]) - } - - legend("topright",paste(basename(files[c(classnames[[k]])])), col = colvect, lty = lty, pch = pch) - - }#end length ==1 - - dev.off() #pdf(pdfname,w=16,h=10) - - invisible(TIC) + dev.off() } - -#@author Y. Guitton -getTIC <- function(file,rtcor=NULL) { - object <- xcmsRaw(file) - cbind(if (is.null(rtcor)) object@scantime else rtcor, rawEIC(object,mzrange=range(object@env$mz))$intensity) -} - -## -## overlay TIC from all files in current folder or from xcmsSet, create pdf -## -#@author Y. Guitton -getTICs <- function(xcmsSet=NULL,files=NULL, pdfname="TICs.pdf",rt=c("raw","corrected")) { - cat("Creating TIC pdf...\n") - - if (is.null(xcmsSet)) { - filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]", "[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]") - filepattern <- paste(paste("\\.", filepattern, "$", sep = ""), collapse = "|") - if (is.null(files)) - files <- getwd() - info <- file.info(files) - listed <- list.files(files[info$isdir], pattern = filepattern, recursive = TRUE, full.names = TRUE) - files <- c(files[!info$isdir], listed) - } else { - files <- filepaths(xcmsSet) - } - - phenoDataClass<-as.vector(levels(xcmsSet@phenoData[,1])) #sometime phenoData have more than 1 column use first as class - classnames<-vector("list",length(phenoDataClass)) - for (i in 1:length(phenoDataClass)){ - classnames[[i]]<-which( xcmsSet@phenoData[,1]==phenoDataClass[i]) - } - - N <- length(files) - TIC <- vector("list",N) - - for (i in 1:N) { - if (!is.null(xcmsSet) && rt == "corrected") - rtcor <- xcmsSet@rt$corrected[[i]] else - rtcor <- NULL - TIC[[i]] <- getTIC(files[i],rtcor=rtcor) - } - - pdf(pdfname,w=16,h=10) - cols <- rainbow(N) - lty = 1:N - pch = 1:N - #search for max x and max y in TICs - xlim = range(sapply(TIC, function(x) range(x[,1]))) - ylim = range(sapply(TIC, function(x) range(x[,2]))) - ylim = c(-ylim[2], ylim[2]) - - - ##plot start - if (length(phenoDataClass)>2){ - for (k in 1:(length(phenoDataClass)-1)){ - for (l in (k+1):length(phenoDataClass)){ - #print(paste(phenoDataClass[k],"vs",phenoDataClass[l],sep=" ")) - plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",phenoDataClass[k]," vs ",phenoDataClass[l], sep=""), xlab = "Retention Time (min)", ylab = "TIC") - colvect<-NULL - for (j in 1:length(classnames[[k]])) { - tic <- TIC[[classnames[[k]][j]]] - # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") - points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") - colvect<-append(colvect,cols[classnames[[k]][j]]) - } - for (j in 1:length(classnames[[l]])) { - # i=class2names[j] - tic <- TIC[[classnames[[l]][j]]] - points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") - colvect<-append(colvect,cols[classnames[[l]][j]]) - } - legend("topright",paste(basename(files[c(classnames[[k]],classnames[[l]])])), col = colvect, lty = lty, pch = pch) - } - } - }#end if length >2 - if (length(phenoDataClass)==2){ - k=1 - l=2 - - plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",phenoDataClass[k],"vs",phenoDataClass[l], sep=""), xlab = "Retention Time (min)", ylab = "TIC") - colvect<-NULL - for (j in 1:length(classnames[[k]])) { - tic <- TIC[[classnames[[k]][j]]] - # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") - points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") - colvect<-append(colvect,cols[classnames[[k]][j]]) - } - for (j in 1:length(classnames[[l]])) { - # i=class2names[j] - tic <- TIC[[classnames[[l]][j]]] - points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") - colvect<-append(colvect,cols[classnames[[l]][j]]) - } - legend("topright",paste(basename(files[c(classnames[[k]],classnames[[l]])])), col = colvect, lty = lty, pch = pch) - - }#end length ==2 - - #case where only one class - if (length(phenoDataClass)==1){ - k=1 - ylim = range(sapply(TIC, function(x) range(x[,2]))) - - plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",phenoDataClass[k], sep=""), xlab = "Retention Time (min)", ylab = "TIC") - colvect<-NULL - for (j in 1:length(classnames[[k]])) { - tic <- TIC[[classnames[[k]][j]]] - # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") - points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") - colvect<-append(colvect,cols[classnames[[k]][j]]) - } - - legend("topright",paste(basename(files[c(classnames[[k]])])), col = colvect, lty = lty, pch = pch) - - }#end length ==1 - - dev.off() #pdf(pdfname,w=16,h=10) - - invisible(TIC) -} - - - -## -## Get the polarities from all the samples of a condition +# Get the polarities from all the samples of a condition #@author Misharl Monsoor misharl.monsoor@sb-roscoff.fr ABiMS TEAM #@author Gildas Le Corguille lecorguille@sb-roscoff.fr ABiMS TEAM -getSampleMetadata <- function(xcmsSet=NULL, sampleMetadataOutput="sampleMetadata.tsv") { +getSampleMetadata <- function(xdata=NULL, sampleMetadataOutput="sampleMetadata.tsv") { cat("Creating the sampleMetadata file...\n") #Create the sampleMetada dataframe - sampleMetadata=xset@phenoData - sampleNamesOrigin=rownames(sampleMetadata) - sampleNamesMakeNames=make.names(sampleNamesOrigin) + sampleMetadata <- xdata@phenoData@data + rownames(sampleMetadata) <- NULL + colnames(sampleMetadata) <- c("sampleMetadata", "class") + + sampleNamesOrigin <- sampleMetadata$sampleMetadata + sampleNamesMakeNames <- make.names(sampleNamesOrigin) if (any(duplicated(sampleNamesMakeNames))) { write("\n\nERROR: Usually, R has trouble to deal with special characters in its column names, so it rename them using make.names().\nIn your case, at least two columns after the renaming obtain the same name, thus XCMS will collapse those columns per name.", stderr()) @@ -330,63 +280,49 @@ } } - sampleMetadata$sampleMetadata=sampleNamesMakeNames - sampleMetadata=cbind(sampleMetadata["sampleMetadata"],sampleMetadata["class"]) #Reorder columns - rownames(sampleMetadata)=NULL + sampleMetadata$sampleMetadata <- sampleNamesMakeNames + - #Create a list of files name in the current directory - list_files=xset@filepaths #For each sample file, the following actions are done - for (file in list_files){ + for (fileIdx in 1:length(fileNames(xdata))) { #Check if the file is in the CDF format - if (!mzR:::netCDFIsFile(file)){ + if (!mzR:::netCDFIsFile(fileNames(xdata))) { # If the column isn't exist, with add one filled with NA - if (is.null(sampleMetadata$polarity)) sampleMetadata$polarity=NA + if (is.null(sampleMetadata$polarity)) sampleMetadata$polarity <- NA - #Create a simple xcmsRaw object for each sample - xcmsRaw=xcmsRaw(file) #Extract the polarity (a list of polarities) - polarity=xcmsRaw@polarity + polarity <- fData(xdata)[fData(xdata)$fileIdx == fileIdx,"polarity"] #Verify if all the scans have the same polarity - uniq_list=unique(polarity) + uniq_list <- unique(polarity) if (length(uniq_list)>1){ - polarity="mixed" + polarity <- "mixed" } else { - polarity=as.character(uniq_list) + polarity <- as.character(uniq_list) } - #Transforms the character to obtain only the sample name - filename=basename(file) - library(tools) - samplename=file_path_sans_ext(filename) #Set the polarity attribute - sampleMetadata$polarity[sampleMetadata$sampleMetadata==samplename]=polarity - - #Delete xcmsRaw object because it creates a bug for the fillpeaks step - rm(xcmsRaw) + sampleMetadata$polarity[fileIdx] <- polarity } } write.table(sampleMetadata, sep="\t", quote=FALSE, row.names=FALSE, file=sampleMetadataOutput) - return(list("sampleNamesOrigin"=sampleNamesOrigin,"sampleNamesMakeNames"=sampleNamesMakeNames)) + return(list("sampleNamesOrigin"=sampleNamesOrigin, "sampleNamesMakeNames"=sampleNamesMakeNames)) } -## -## This function check if xcms will found all the files -## +# This function check if xcms will found all the files #@author Gildas Le Corguille lecorguille@sb-roscoff.fr ABiMS TEAM checkFilesCompatibilityWithXcms <- function(directory) { cat("Checking files filenames compatibilities with xmcs...\n") # WHAT XCMS WILL FIND filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]","[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]") - filepattern <- paste(paste("\\.", filepattern, "$", sep = ""),collapse = "|") + filepattern <- paste(paste("\\.", filepattern, "$", sep=""),collapse="|") info <- file.info(directory) - listed <- list.files(directory[info$isdir], pattern = filepattern,recursive = TRUE, full.names = TRUE) + listed <- list.files(directory[info$isdir], pattern=filepattern, recursive=TRUE, full.names=TRUE) files <- c(directory[!info$isdir], listed) files_abs <- file.path(getwd(), files) exists <- file.exists(files_abs) @@ -394,8 +330,8 @@ files[exists] <- sub("//","/",files[exists]) # WHAT IS ON THE FILESYSTEM - filesystem_filepaths=system(paste("find $PWD/",directory," -not -name '\\.*' -not -path '*conda-env*' -type f -name \"*\"", sep=""), intern=T) - filesystem_filepaths=filesystem_filepaths[grep(filepattern, filesystem_filepaths, perl=T)] + filesystem_filepaths <- system(paste0("find \"$PWD/",directory,"\" -not -name '\\.*' -not -path '*conda-env*' -type f -name \"*\""), intern=T) + filesystem_filepaths <- filesystem_filepaths[grep(filepattern, filesystem_filepaths, perl=T)] # COMPARISON if (!is.na(table(filesystem_filepaths %in% files)["FALSE"])) { @@ -406,16 +342,26 @@ } +#This function list the compatible files within the directory as xcms did +#@author Gildas Le Corguille lecorguille@sb-roscoff.fr ABiMS TEAM +getMSFiles <- function (directory) { + filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]","[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]") + filepattern <- paste(paste("\\.", filepattern, "$", sep=""),collapse="|") + info <- file.info(directory) + listed <- list.files(directory[info$isdir], pattern=filepattern,recursive=TRUE, full.names=TRUE) + files <- c(directory[!info$isdir], listed) + exists <- file.exists(files) + files <- files[exists] + return(files) +} -## -## This function check if XML contains special caracters. It also checks integrity and completness. -## +# This function check if XML contains special caracters. It also checks integrity and completness. #@author Misharl Monsoor misharl.monsoor@sb-roscoff.fr ABiMS TEAM checkXmlStructure <- function (directory) { cat("Checking XML structure...\n") - cmd=paste("IFS=$'\n'; for xml in $(find",directory,"-not -name '\\.*' -not -path '*conda-env*' -type f -iname '*.*ml*'); do if [ $(xmllint --nonet --noout \"$xml\" 2> /dev/null; echo $?) -gt 0 ]; then echo $xml;fi; done;") - capture=system(cmd,intern=TRUE) + cmd <- paste0("IFS=$'\n'; for xml in $(find '",directory,"' -not -name '\\.*' -not -path '*conda-env*' -type f -iname '*.*ml*'); do if [ $(xmllint --nonet --noout \"$xml\" 2> /dev/null; echo $?) -gt 0 ]; then echo $xml;fi; done;") + capture <- system(cmd, intern=TRUE) if (length(capture)>0){ #message=paste("The following mzXML or mzML file is incorrect, please check these files first:",capture) @@ -427,24 +373,22 @@ } -## -## This function check if XML contain special characters -## +# This function check if XML contain special characters #@author Misharl Monsoor misharl.monsoor@sb-roscoff.fr ABiMS TEAM deleteXmlBadCharacters<- function (directory) { cat("Checking Non ASCII characters in the XML...\n") - processed=F - l=system( paste("find",directory, "-not -name '\\.*' -not -path '*conda-env*' -type f -iname '*.*ml*'"),intern=TRUE) + processed <- F + l <- system( paste0("find '",directory, "' -not -name '\\.*' -not -path '*conda-env*' -type f -iname '*.*ml*'"), intern=TRUE) for (i in l){ - cmd=paste("LC_ALL=C grep '[^ -~]' \"",i,"\"",sep="") - capture=suppressWarnings(system(cmd,intern=TRUE)) + cmd <- paste("LC_ALL=C grep '[^ -~]' \"", i, "\"", sep="") + capture <- suppressWarnings(system(cmd, intern=TRUE)) if (length(capture)>0){ - cmd=paste("perl -i -pe 's/[^[:ascii:]]//g;'",i) + cmd <- paste("perl -i -pe 's/[^[:ascii:]]//g;'",i) print( paste("WARNING: Non ASCII characters have been removed from the ",i,"file") ) - c=system(cmd,intern=TRUE) - capture="" - processed=T + c <- system(cmd, intern=TRUE) + capture <- "" + processed <- T } } if (processed) cat("\n\n") @@ -452,17 +396,15 @@ } -## -## This function will compute MD5 checksum to check the data integrity -## +# This function will compute MD5 checksum to check the data integrity #@author Gildas Le Corguille lecorguille@sb-roscoff.fr getMd5sum <- function (directory) { cat("Compute md5 checksum...\n") # WHAT XCMS WILL FIND filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]","[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]") - filepattern <- paste(paste("\\.", filepattern, "$", sep = ""),collapse = "|") + filepattern <- paste(paste("\\.", filepattern, "$", sep=""),collapse="|") info <- file.info(directory) - listed <- list.files(directory[info$isdir], pattern = filepattern,recursive = TRUE, full.names = TRUE) + listed <- list.files(directory[info$isdir], pattern=filepattern, recursive=TRUE, full.names=TRUE) files <- c(directory[!info$isdir], listed) exists <- file.exists(files) files <- files[exists] @@ -476,80 +418,93 @@ # 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"]] +#@author Gildas Le Corguille lecorguille@sb-roscoff.fr +getRawfilePathFromArguments <- function(singlefile, zipfile, args, prefix="") { + if (!(prefix %in% c("","Positive","Negative","MS1","MS2"))) stop("prefix must be either '', 'Positive', 'Negative', 'MS1' or 'MS2'") + + if (!is.null(args[[paste0("zipfile",prefix)]])) zipfile <- args[[paste0("zipfile",prefix)]] - 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 (!is.null(args[[paste0("singlefile_galaxyPath",prefix)]])) { + singlefile_galaxyPaths <- args[[paste0("singlefile_galaxyPath",prefix)]] + singlefile_sampleNames <- args[[paste0("singlefile_sampleName",prefix)]] + } + 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] + # In case, an url is used to import data within Galaxy + singlefile_sampleName <- tail(unlist(strsplit(singlefile_sampleName,"/")), n=1) + singlefile[[singlefile_sampleName]] <- singlefile_galaxyPath } - 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 - } - } - for (argument in c("zipfile","zipfilePositive","zipfileNegative","singlefile_galaxyPath","singlefile_sampleName","singlefile_galaxyPathPositive","singlefile_sampleNamePositive","singlefile_galaxyPathNegative","singlefile_sampleNameNegative")) { - listArguments[[argument]]=NULL - } - return(list(zipfile=zipfile, singlefile=singlefile, listArguments=listArguments)) + } + 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 +#@author Gildas Le Corguille lecorguille@sb-roscoff.fr retrieveRawfileInTheWorkingDirectory <- function(singlefile, zipfile) { if(!is.null(singlefile) && (length("singlefile")>0)) { for (singlefile_sampleName in names(singlefile)) { - singlefile_galaxyPath = singlefile[[singlefile_sampleName]] + 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!") + 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) + if (!suppressWarnings( try (file.link(singlefile_galaxyPath, singlefile_sampleName), silent=T))) + file.copy(singlefile_galaxyPath, singlefile_sampleName) + } - directory = "." + directory <- "." } - if(!is.null(zipfile) && (zipfile!="")) { + 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!") + 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"] + #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 + suppressWarnings(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") } return (directory) } + + +# This function retrieve a xset like object +#@author Gildas Le Corguille lecorguille@sb-roscoff.fr +getxcmsSetObject <- function(xobject) { + # XCMS 1.x + if (class(xobject) == "xcmsSet") + return (xobject) + # XCMS 3.x + if (class(xobject) == "XCMSnExp") { + # Get the legacy xcmsSet object + suppressWarnings(xset <- as(xobject, 'xcmsSet')) + if (!is.null(xset@phenoData$sample_group)) + sampclass(xset) <- xset@phenoData$sample_group + else + sampclass(xset) <- "." + return (xset) + } +}
--- a/macros.xml Wed Nov 29 09:46:41 2017 -0500 +++ b/macros.xml Tue Sep 18 16:12:29 2018 -0400 @@ -1,141 +1,52 @@ <?xml version="1.0"?> <macros> - <xml name="requirements"> - <requirements> - <requirement type="package" version="0.4_1">r-snow</requirement> - <requirement type="package" version="1.46.0">bioconductor-xcms</requirement> - <requirement type="package" version="1.1_4">r-batch</requirement> - </requirements> - </xml> - <xml name="requirements_light"> - <requirements> - <requirement type="package" version="1.46.0">bioconductor-xcms</requirement> - </requirements> - </xml> <xml name="stdio"> <stdio> <exit_code range="1" level="fatal" /> </stdio> </xml> - <token name="@COMMAND_XCMS_SCRIPT@"> - LC_ALL=C Rscript $__tool_directory__/xcms.r - </token> + <!-- COMMAND --> + <token name="@COMMAND_RSCRIPT@">LC_ALL=C Rscript $__tool_directory__/</token> <token name="@COMMAND_LOG_EXIT@"> ; return=\$?; - mv log.txt '$log'; - cat '$log'; + cat 'log.txt'; sh -c "exit \$return" </token> - <!-- zipfile load for planemo test --> - - <token name="@COMMAND_FILE_LOAD@"> - #if $file_load_section.file_load_conditional.file_load_select == "yes": - #if $file_load_section.file_load_conditional.input[0].is_of_type("mzxml") or $file_load_section.file_load_conditional.input[0].is_of_type("mzml") or $file_load_section.file_load_conditional.input[0].is_of_type("mzdata") or $file_load_section.file_load_conditional.input[0].is_of_type("netcdf"): - #set singlefile_galaxyPath = ','.join( [ str( $single_file ) for $single_file in $file_load_section.file_load_conditional.input ] ) - #set singlefile_sampleName = ','.join( [ str( $single_file.name ) for $single_file in $file_load_section.file_load_conditional.input ] ) - - singlefile_galaxyPath '$singlefile_galaxyPath' singlefile_sampleName '$singlefile_sampleName' - #else - zipfile '$file_load_section.file_load_conditional.input' - #end if - #end if - </token> - - <xml name="input_file_load"> - <section name="file_load_section" title="Resubmit your raw dataset or your zip file"> - <conditional name="file_load_conditional"> - <param name="file_load_select" type="select" label="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> + <!-- INPUT_VALIDATORS --> + <xml name="input_validator_range_integer"> + <validator type="regex" message="The format is 'min,max'" >[0-9]+ *, *[0-9]+</validator> </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 name="input_validator_range_float"> + <validator type="regex" message="The format is 'min,max'" >[0-9]+\.?[0-9]* *, *[0-9]+\.?[0-9]*</validator> </xml> - <xml name="test_file_load_single"> - <section name="file_load_section"> - <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 name="input_validator_list_integer"> + <validator type="regex" message="The format is '1,2,4,6'" >[0-9, ]+</validator> </xml> - <token name="@COMMAND_PEAKLIST@"> - #if $peaklist.peaklistBool - variableMetadataOutput '$variableMetadata' - dataMatrixOutput '$dataMatrix' - convertRTMinute $peaklist.convertRTMinute - numDigitsMZ $peaklist.numDigitsMZ - numDigitsRT $peaklist.numDigitsRT - intval $peaklist.intval - #end if - </token> + + <token name="@INPUT_IMAGE_LABEL@">RData file</token> + <token name="@INPUT_IMAGE_HELP@">It contains a xcms3::XCMSnExp object (named xdata)</token> - <xml name="input_peaklist"> - <conditional name="peaklist"> - <param name="peaklistBool" type="boolean" label="Get a Peak List" /> - <when value="true"> - <param name="convertRTMinute" type="boolean" checked="false" truevalue="TRUE" falsevalue="FALSE" label="Convert retention time (seconds) into minutes" help="Convert the columns rtmed, rtmin and rtmax into minutes"/> - <param name="numDigitsMZ" type="integer" value="4" label="Number of decimal places for mass values reported in ions' identifiers." help="A minimum of 4 decimal places is recommended. Useful to avoid duplicates within identifiers" /> - <param name="numDigitsRT" type="integer" value="0" label="Number of decimal places for retention time values reported in ions' identifiers." help="Useful to avoid duplicates within identifiers" /> - <param name="intval" type="select" label="Reported intensity values" help="[intval] See the help section below"> - <option value="into" selected="true">into</option> - <option value="maxo">maxo</option> - <option value="intb">intb</option> - </param> - </when> - <when value="false" /> - </conditional> - </xml> - <xml name="output_peaklist" token_function=""> - <data name="variableMetadata" format="tabular" label="${image.name[:-6]}.@FUNCTION@.variableMetadata.tsv"> - <filter>(peaklist['peaklistBool'])</filter> - </data> - <data name="dataMatrix" format="tabular" label="${image.name[:-6]}.@FUNCTION@.dataMatrix.tsv" > - <filter>(peaklist['peaklistBool'])</filter> - </data> - </xml> - - <token name="@HELP_AUTHORS@"> -.. class:: infomark - -**Authors** Colin A. Smith csmith@scripps.edu, Ralf Tautenhahn rtautenh@gmail.com, Steffen Neumann sneumann@ipb-halle.de, Paul Benton hpaul.benton08@imperial.ac.uk and Christopher Conley cjconley@ucdavis.edu + <!-- MISC --> + <token name="@HELP_AUTHORS_WRAPPERS@"> .. class:: infomark -**Galaxy integration** ABiMS TEAM - UPMC/CNRS - Station biologique de Roscoff and Yann Guitton yann.guitton@oniris-nantes.fr - part of Workflow4Metabolomics.org [W4M] +**Galaxy integration** ABiMS TEAM - SU/CNRS - Station biologique de Roscoff and Yann Guitton - LABERCA +Part of Workflow4Metabolomics.org [W4M] | Contact support@workflow4metabolomics.org for any questions or concerns about the Galaxy implementation of this tool. ---------------------------------------------------- - </token> - - <xml name="citation"> - <citations> - <citation type="doi">10.1021/ac051437y</citation> + <xml name="citation_w4m"> <citation type="doi">10.1093/bioinformatics/btu813</citation> - </citations> </xml> </macros>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/macros_xcms.xml Tue Sep 18 16:12:29 2018 -0400 @@ -0,0 +1,244 @@ +<?xml version="1.0"?> +<macros> + + <token name="@WRAPPER_VERSION@">3.0.0</token> + <xml name="requirements"> + <requirements> + <requirement type="package" version="@WRAPPER_VERSION@">bioconductor-xcms</requirement> + <requirement type="package" version="1.1_4">r-batch</requirement> + <requirement type="package" version="1.1_2">r-rcolorbrewer</requirement> + <requirement type="package" version="6.0">unzip</requirement> + <yield /> + </requirements> + </xml> + + <!-- FILE_LOAD for planemo test --> + <token name="@COMMAND_FILE_LOAD@"> + #if $file_load_section.file_load_conditional.file_load_select == "yes": + #if $file_load_section.file_load_conditional.input[0].is_of_type("mzxml") or $file_load_section.file_load_conditional.input[0].is_of_type("mzml") or $file_load_section.file_load_conditional.input[0].is_of_type("mzdata") or $file_load_section.file_load_conditional.input[0].is_of_type("netcdf"): + #set singlefile_galaxyPath = '|'.join( [ str( $single_file ) for $single_file in $file_load_section.file_load_conditional.input ] ) + #set singlefile_sampleName = '|'.join( [ str( $single_file.name ) for $single_file in $file_load_section.file_load_conditional.input ] ) + + singlefile_galaxyPath '$singlefile_galaxyPath' singlefile_sampleName '$singlefile_sampleName' + #else + zipfile '$file_load_section.file_load_conditional.input' + #end if + #end if + </token> + + <xml name="input_file_load"> + <section name="file_load_section" title="Resubmit your raw dataset or your zip file"> + <conditional name="file_load_conditional"> + <param name="file_load_select" type="select" label="Resubmit your dataset or your zip file" help="Use only if you get a message saying 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_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_zip_sacuri"> + <section name="file_load_section"> + <conditional name="file_load_conditional"> + <param name="file_load_select" value="yes" /> + <param name="input" value="sacuri_dir_root.zip" ftype="zip" /> + </conditional> + </section> + </xml> + + <xml name="test_file_load_single"> + <section name="file_load_section"> + <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> + + <xml name="test_file_load_single_ko15"> + <section name="file_load_section"> + <conditional name="file_load_conditional"> + <param name="file_load_select" value="yes" /> + <param name="input" value="ko15.CDF" ftype="netcdf" /> + </conditional> + </section> + </xml> + + <!-- PEAKLIST --> + <token name="@COMMAND_PEAKLIST@"> + #if $peaklist.peaklistBool + convertRTMinute $peaklist.convertRTMinute + numDigitsMZ $peaklist.numDigitsMZ + numDigitsRT $peaklist.numDigitsRT + intval $peaklist.intval + naTOzero $peaklist.naTOzero + #end if + </token> + + <xml name="input_peaklist"> + <conditional name="peaklist"> + <param name="peaklistBool" type="boolean" label="Get a Peak List" /> + <when value="true"> + <param name="convertRTMinute" type="boolean" checked="false" truevalue="TRUE" falsevalue="FALSE" label="Convert retention time (seconds) into minutes" help="Convert the columns rtmed, rtmin and rtmax into minutes"/> + <param name="numDigitsMZ" type="integer" value="4" label="Number of decimal places for mass values reported in ions' identifiers." help="A minimum of 4 decimal places is recommended. Useful to avoid duplicates within identifiers" /> + <param name="numDigitsRT" type="integer" value="0" label="Number of decimal places for retention time values reported in ions' identifiers." help="Useful to avoid duplicates within identifiers" /> + <param name="intval" type="select" label="Reported intensity values" help="[intval] See the help section below"> + <option value="into" selected="true">into</option> + <option value="maxo">maxo</option> + <option value="intb">intb</option> + </param> + <param name="naTOzero" type="boolean" checked="true" truevalue="TRUE" falsevalue="FALSE" label="Replace the remain NA by 0 in the dataMatrix" help="Rather mandatory for some downstream statistical steps"/> + </when> + <when value="false" /> + </conditional> + </xml> + + <xml name="output_peaklist" token_function=""> + <data name="variableMetadata" format="tabular" label="${image.name[:-6]}.@FUNCTION@.variableMetadata.tsv" from_work_dir="variableMetadata.tsv" > + <filter>(peaklist['peaklistBool'])</filter> + </data> + <data name="dataMatrix" format="tabular" label="${image.name[:-6]}.@FUNCTION@.dataMatrix.tsv" from_work_dir="dataMatrix.tsv" > + <filter>(peaklist['peaklistBool'])</filter> + </data> + </xml> + + <token name="@HELP_PEAKLIST@"> + +Get a Peak List +--------------- + +If 'true', the module generates two additional files corresponding to the peak list: +- the variable metadata file (corresponding to information about extracted ions such as mass or retention time) +- the data matrix (corresponding to related intensities) + +**decimal places for [mass or retention time] values in identifiers** + + | Ions' identifiers are constructed as MxxxTyyy where 'xxx' is the ion median mass and 'yyy' the ion median retention time. + | Two parameters are used to adjust the number of decimal places wanted in identifiers for mass and retention time respectively. + | Theses parameters do not affect decimal places in columns other than the identifier one. + +**Reported intensity values** + + | This parameter determines which values should be reported as intensities in the dataMatrix table; it correspond to xcms 'intval' parameter: + | - into: integrated area of original (raw) peak + | - maxo: maximum intensity of original (raw) peak + | - intb: baseline corrected integrated peak area (only available if peak detection was done by ‘findPeaks.centWave’) + + </token> + + <token name="@HELP_PEAKLIST_OUTPUT@"> +xset.variableMetadata.tsv : tabular format + + | Table containing information about ions; can be used as one input of **Quality_Metrics** or **Generic_filter** modules. + +xset.dataMatrix.tsv : tabular format + + | Table containing ions' intensities; can be used as one input of **Quality_Metrics** or **Generic_filter** modules. + </token> + + <!-- CENTWAVE --> + <token name="@COMMAND_CENTWAVE@"> + ppm $methods.ppm + peakwidth "c($methods.peakwidth)" + + ## Advanced + snthresh $methods.CentWaveAdv.snthresh + prefilter "c($methods.CentWaveAdv.prefilter)" + mzCenterFun $methods.CentWaveAdv.mzCenterFun + integrate $methods.CentWaveAdv.integrate + mzdiff $methods.CentWaveAdv.mzdiff + fitgauss $methods.CentWaveAdv.fitgauss + noise $methods.CentWaveAdv.noise + verboseColumns $methods.CentWaveAdv.verboseColumns + </token> + + <xml name="input_centwave"> + <param argument="ppm" type="integer" value="25" label="Max tolerated ppm m/z deviation in consecutive scans in ppm" help="for the initial ROI definition." /> + <param argument="peakwidth" type="text" value="20,50" label="Min,Max peak width in seconds" help="with the expected approximate peak width in chromatographic space."> + <expand macro="input_validator_range_float"/> + </param> + </xml> + + <xml name="input_centwaveAdv"> + <param argument="snthresh" type="integer" value="10" label="Signal to Noise ratio cutoff" /> + <param argument="prefilter" type="text" value="3,100" label="Prefilter step for for the first analysis step (ROI detection)" help="Separate by coma k, I. Mass traces are only retained if they contain at least ‘k‘ peaks with intensity ‘>= I‘."> + <expand macro="input_validator_range_integer"/> + </param> + <param argument="mzCenterFun" type="select" label="Name of the function to calculate the m/z center of the chromatographic peak" > + <option value="wMean">intensity weighted mean of the peak's m/z values</option> + <option value="mean">mean of the peak's m/z values</option> + <option value="apex">use the m/z value at the peak apex</option> + <option value="wMeanApex3">ntensity weighted mean of the m/z value at the peak apex and the m/z values left and right of it</option> + <option value="meanApex3">mean of the m/z value of the peak apex and the m/z values left and right of it</option> + </param> + <param argument="integrate" type="select" label="Integration method" > + <option value="1">peak limits are found through descent on the mexican hat filtered data (more robust, but less exact)</option> + <option value="2">peak limits based on real data (more accurate but prone to noise)</option> + </param> + <param argument="mzdiff" type="float" value="-0.001" label="Minimum difference in m/z for peaks with overlapping retention times" help="can be negative to allow overlap" /> + <param argument="fitgauss" type="boolean" checked="false" truevalue="TRUE" falsevalue="FALSE" label="fitgauss" help="whether or not a Gaussian should be fitted to each peak" /> + <param argument="noise" type="integer" value="0" label="Noise filter" help="allowing to set a minimum intensity required for centroids to be considered in the first analysis step (centroids with intensity lower than ‘noise’ are omitted from ROI detection)." /> + <param argument="verboseColumns" type="boolean" checked="false" truevalue="TRUE" falsevalue="FALSE" label="verbose Columns" help="whether additional peak meta data columns should be returned" /> + </xml> + + <token name="@COMMAND_CENTWAVEADVROI@"> + #if $sectionROI.roiList: + roiList '$sectionROI.roiList' + firstBaselineCheck $sectionROI.firstBaselineCheck + #if $sectionROI.roiScales != "": + roiScales "c($sectionROI.roiScales)" + #end if + #end if + </token> + + <xml name="input_centwaveAdvROI" token_optional="true"> + <param argument="roiList" type="data" format="tabular" optional="@OPTIONAL@" label="List of regions-of-interest (ROI) representing detected mass traces" help="If ROIs are submitted the first analysis step is omitted and chromatographic peak detection is performed on the submitted ROIs. Each ROI is expected to have the following elements specified: ‘scmin’ (start scan index), ‘scmax’ (end scan index), ‘mzmin’ (minimum m/z), ‘mzmax’ (maximum m/z), ‘length’ (number of scans), ‘intensity’ (summed intensity)." /> + <param argument="firstBaselineCheck" type="boolean" checked="true" truevalue="TRUE" falsevalue="FALSE" label="Is continuous data within regions of interest is checked to be above the first baseline." /> + <param argument="roiScales" type="text" value="" optional="true" label="Numeric vector defining the scale for each region of interest in ‘roiList’" help="Length equal to ‘roiList’ - Should be used for the centWave-wavelets (format 0.9,1,0.2)"> + <expand macro="input_validator_range_float"/> + </param> + </xml> + + <!-- MISC --> + <token name="@HELP_AUTHORS@"> +.. class:: infomark + +**Authors** Colin A. Smith csmith@scripps.edu, Ralf Tautenhahn rtautenh@gmail.com, Steffen Neumann sneumann@ipb-halle.de, Paul Benton hpaul.benton08@imperial.ac.uk and Christopher Conley cjconley@ucdavis.edu + +@HELP_AUTHORS_WRAPPERS@ + +--------------------------------------------------- + + </token> + + <token name="@HELP_XCMS_MANUAL@"> + +For details and explanations concerning all the parameters and workflow of xcms_ package, see its manual_ and this example_ + +.. _xcms: https://bioconductor.org/packages/release/bioc/html/xcms.html +.. _manual: http://www.bioconductor.org/packages/release/bioc/manuals/xcms/man/xcms.pdf +.. _example: https://bioconductor.org/packages/release/bioc/vignettes/xcms/inst/doc/xcms.html + + </token> + + <xml name="citation"> + <citations> + <citation type="doi">10.1021/ac051437y</citation> + <expand macro="citation_w4m"/> + </citations> + </xml> +</macros>
--- a/xcms.r Wed Nov 29 09:46:41 2017 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,229 +0,0 @@ -#!/usr/bin/env Rscript -# xcms.r version="2.2.0" -#Authors ABIMS TEAM -#BPC Addition from Y.guitton - - -# ----- LOG FILE ----- -log_file=file("log.txt", open = "wt") -sink(log_file) -sink(log_file, type = "output") - - -# ----- PACKAGE ----- -cat("\tPACKAGE INFO\n") -#pkgs=c("xcms","batch") -pkgs=c("parallel","BiocGenerics", "Biobase", "Rcpp", "mzR", "xcms","snow","batch") -for(pkg in pkgs) { - suppressPackageStartupMessages( stopifnot( library(pkg, quietly=TRUE, logical.return=TRUE, character.only=TRUE))) - cat(pkg,"\t",as.character(packageVersion(pkg)),"\n",sep="") -} -source_local <- function(fname){ argv <- commandArgs(trailingOnly = FALSE); base_dir <- dirname(substring(argv[grep("--file=", argv)], 8)); source(paste(base_dir, fname, sep="/")) } -cat("\n\n"); - - - - - -# ----- ARGUMENTS ----- -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') - -cat("\n\n"); - - -# ----- ARGUMENTS PROCESSING ----- -cat("\tINFILE PROCESSING INFO\n") - -#image is an .RData file necessary to use xset variable given by previous tools -if (!is.null(listArguments[["image"]])){ - load(listArguments[["image"]]); listArguments[["image"]]=NULL -} - -#Import the different functions -source_local("lib.r") - -cat("\n\n") - -#Import the different functions - -# ----- PROCESSING INFILE ----- -cat("\tARGUMENTS PROCESSING INFO\n") - -# Save arguments to generate a report -if (!exists("listOFlistArguments")) listOFlistArguments=list() -listOFlistArguments[[paste(format(Sys.time(), "%y%m%d-%H:%M:%S_"),listArguments[["xfunction"]],sep="")]] = listArguments - - -#saving the commun parameters -thefunction = listArguments[["xfunction"]]; listArguments[["xfunction"]]=NULL #delete from the list of arguments - -xsetRdataOutput = paste(thefunction,"RData",sep=".") -if (!is.null(listArguments[["xsetRdataOutput"]])){ - xsetRdataOutput = listArguments[["xsetRdataOutput"]]; listArguments[["xsetRdataOutput"]]=NULL -} - -#saving the specific parameters -rplotspdf = "Rplots.pdf" -if (!is.null(listArguments[["rplotspdf"]])){ - rplotspdf = listArguments[["rplotspdf"]]; listArguments[["rplotspdf"]]=NULL -} -sampleMetadataOutput = "sampleMetadata.tsv" -if (!is.null(listArguments[["sampleMetadataOutput"]])){ - sampleMetadataOutput = listArguments[["sampleMetadataOutput"]]; listArguments[["sampleMetadataOutput"]]=NULL -} -variableMetadataOutput = "variableMetadata.tsv" -if (!is.null(listArguments[["variableMetadataOutput"]])){ - variableMetadataOutput = listArguments[["variableMetadataOutput"]]; listArguments[["variableMetadataOutput"]]=NULL -} -dataMatrixOutput = "dataMatrix.tsv" -if (!is.null(listArguments[["dataMatrixOutput"]])){ - dataMatrixOutput = listArguments[["dataMatrixOutput"]]; listArguments[["dataMatrixOutput"]]=NULL -} -if (!is.null(listArguments[["convertRTMinute"]])){ - convertRTMinute = listArguments[["convertRTMinute"]]; listArguments[["convertRTMinute"]]=NULL -} -if (!is.null(listArguments[["numDigitsMZ"]])){ - numDigitsMZ = listArguments[["numDigitsMZ"]]; listArguments[["numDigitsMZ"]]=NULL -} -if (!is.null(listArguments[["numDigitsRT"]])){ - numDigitsRT = listArguments[["numDigitsRT"]]; listArguments[["numDigitsRT"]]=NULL -} -if (!is.null(listArguments[["intval"]])){ - intval = listArguments[["intval"]]; listArguments[["intval"]]=NULL -} - -if (thefunction %in% c("xcmsSet","retcor")) { - ticspdf = listArguments[["ticspdf"]]; listArguments[["ticspdf"]]=NULL - bicspdf = listArguments[["bicspdf"]]; listArguments[["bicspdf"]]=NULL -} - - -if (thefunction %in% c("xcmsSet","retcor","fillPeaks")) { - if (!exists("singlefile")) singlefile=NULL - if (!exists("zipfile")) zipfile=NULL - rawFilePath = getRawfilePathFromArguments(singlefile, zipfile, listArguments) - zipfile = rawFilePath$zipfile - singlefile = rawFilePath$singlefile - listArguments = rawFilePath$listArguments - directory = retrieveRawfileInTheWorkingDirectory(singlefile, zipfile) - md5sumList=list("origin"=getMd5sum(directory)) -} - -#addition of the directory to the list of arguments in the first position -if (thefunction == "xcmsSet") { - checkXmlStructure(directory) - checkFilesCompatibilityWithXcms(directory) - listArguments=append(directory, listArguments) -} - - -#addition of xset object to the list of arguments in the first position -if (exists("xset")){ - listArguments=append(list(xset), listArguments) -} - -cat("\n\n") - - - - -# ----- MAIN PROCESSING INFO ----- -cat("\tMAIN PROCESSING INFO\n") - - -#Verification of a group step before doing the fillpeaks job. - -if (thefunction == "fillPeaks") { - res=try(is.null(groupnames(xset))) - if (class(res) == "try-error"){ - error<-geterrmessage() - write(error, stderr()) - stop("You must always do a group step after a retcor. Otherwise it won't work for the fillpeaks step") - } - -} - -#change the default display settings -#dev.new(file="Rplots.pdf", width=16, height=12) -pdf(file=rplotspdf, width=16, height=12) -if (thefunction == "group") { - par(mfrow=c(2,2)) -} -#else if (thefunction == "retcor") { -#try to change the legend display -# par(xpd=NA) -# par(xpd=T, mar=par()$mar+c(0,0,0,4)) -#} - - -#execution of the function "thefunction" with the parameters given in "listArguments" - -cat("\t\tCOMPUTE\n") -xset = do.call(thefunction, listArguments) - -# check if there are no peaks -if (nrow(peaks(xset)) == 0) { - stop("No peaks were detected. You should review your settings") -} - - -cat("\n\n") - -dev.off() #dev.new(file="Rplots.pdf", width=16, height=12) - -if (thefunction == "xcmsSet") { - - #transform the files absolute pathways into relative pathways - xset@filepaths<-sub(paste(getwd(),"/",sep="") ,"", xset@filepaths) - if(exists("zipfile") && !is.null(zipfile) && (zipfile!="")) { - - #Modify the samples names (erase the path) - for(i in 1:length(sampnames(xset))){ - - sample_name=unlist(strsplit(sampnames(xset)[i], "/")) - sample_name=sample_name[length(sample_name)] - sample_name= unlist(strsplit(sample_name,"[.]"))[1] - sampnames(xset)[i]=sample_name - - } - - } - -} - -# -- TIC -- -if (thefunction == "xcmsSet") { - cat("\t\tGET TIC GRAPH\n") - sampleNamesList = getSampleMetadata(xcmsSet=xset, sampleMetadataOutput=sampleMetadataOutput) - getTICs(xcmsSet=xset, pdfname=ticspdf,rt="raw") - getBPCs(xcmsSet=xset,rt="raw",pdfname=bicspdf) -} else if (thefunction == "retcor") { - cat("\t\tGET TIC GRAPH\n") - getTICs(xcmsSet=xset, pdfname=ticspdf,rt="corrected") - getBPCs(xcmsSet=xset,rt="corrected",pdfname=bicspdf) -} - -if ((thefunction == "group" || thefunction == "fillPeaks") && exists("intval")) { - getPeaklistW4M(xset,intval,convertRTMinute,numDigitsMZ,numDigitsRT,variableMetadataOutput,dataMatrixOutput) -} - - -cat("\n\n") - -# ----- EXPORT ----- - -cat("\tXSET OBJECT INFO\n") -print(xset) -#delete the parameters to avoid the passage to the next tool in .RData image - - -#saving R data in .Rdata file to save the variables used in the present tool -objects2save = c("xset","zipfile","singlefile","listOFlistArguments","md5sumList","sampleNamesList") -save(list=objects2save[objects2save %in% ls()], file=xsetRdataOutput) - -cat("\n\n") - - -cat("\tDONE\n")
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/xcms_retcor.r Tue Sep 18 16:12:29 2018 -0400 @@ -0,0 +1,102 @@ +#!/usr/bin/env Rscript + +# ----- LOG FILE ----- +log_file=file("log.txt", open = "wt") +sink(log_file) +sink(log_file, type = "output") + + +# ----- PACKAGE ----- +cat("\tSESSION INFO\n") + +#Import the different functions +source_local <- function(fname){ argv <- commandArgs(trailingOnly=FALSE); base_dir <- dirname(substring(argv[grep("--file=", argv)], 8)); source(paste(base_dir, fname, sep="/")) } +source_local("lib.r") +source_local("lib-xcms3.x.x.r") + +pkgs <- c("xcms","batch","RColorBrewer") +loadAndDisplayPackages(pkgs) +cat("\n\n"); + + +# ----- ARGUMENTS ----- +cat("\tARGUMENTS INFO\n") +args = parseCommandArgs(evaluate=FALSE) #interpretation of arguments given in command line as an R list of objects +write.table(as.matrix(args), col.names=F, quote=F, sep='\t') + +cat("\n\n") + +# ----- PROCESSING INFILE ----- +cat("\tARGUMENTS PROCESSING INFO\n") + +#saving the specific parameters +method <- args$method + +cat("\n\n") + + +# ----- ARGUMENTS PROCESSING ----- +cat("\tINFILE PROCESSING INFO\n") + +#image is an .RData file necessary to use xset variable given by previous tools +load(args$image); args$image=NULL +if (!exists("xdata")) stop("\n\nERROR: The RData doesn't contain any object called 'xdata'. This RData should have been created by an old version of XMCS 2.*") + +# Handle infiles +if (!exists("singlefile")) singlefile <- NULL +if (!exists("zipfile")) zipfile <- NULL +rawFilePath <- getRawfilePathFromArguments(singlefile, zipfile, args) +zipfile <- rawFilePath$zipfile +singlefile <- rawFilePath$singlefile +directory <- retrieveRawfileInTheWorkingDirectory(singlefile, zipfile) + +cat("\n\n") + + +# ----- MAIN PROCESSING INFO ----- +cat("\tMAIN PROCESSING INFO\n") + + +cat("\t\tCOMPUTE\n") + +cat("\t\t\tAlignment/Retention Time correction\n") +# clear the arguement list to remove unexpected key/value as singlefile_galaxyPath or method ... +args <- args[names(args) %in% slotNames(do.call(paste0(method,"Param"), list()))] + +adjustRtimeParam <- do.call(paste0(method,"Param"), args) +print(adjustRtimeParam) +xdata <- adjustRtime(xdata, param=adjustRtimeParam) + +cat("\t\t\tCompute and Store TIC and BPI\n") +chromTIC_adjusted = chromatogram(xdata, aggregationFun = "sum") +chromBPI_adjusted = chromatogram(xdata, aggregationFun = "max") + +cat("\n\n") + + +# -- TIC -- +cat("\t\tDRAW GRAPHICS\n") +getPlotAdjustedRtime(xdata) + +cat("\n\n") + +# ----- EXPORT ----- + +cat("\tXCMSnExp OBJECT INFO\n") +print(xdata) +cat("\n\n") + +cat("\txcmsSet OBJECT INFO\n") +# Get the legacy xcmsSet object +xset <- getxcmsSetObject(xdata) +print(xset) +cat("\n\n") + +#saving R data in .Rdata file to save the variables used in the present tool +objects2save = c("xdata","zipfile","singlefile","md5sumList","sampleNamesList", "chromTIC", "chromBPI", "chromTIC_adjusted", "chromBPI_adjusted") +save(list=objects2save[objects2save %in% ls()], file="retcor.RData") + +cat("\n\n") + + +cat("\tDONE\n")