From 880cd6b1e82f0ff88942e5cdcc46f3992984c06b Mon Sep 17 00:00:00 2001 From: VinzentRisch Date: Wed, 15 Jan 2025 18:56:33 +0100 Subject: [PATCH] added utils find peaks --- conda-recipe/meta.yaml | 2 ++ q2_ms/assets/find_peaks_xcms.R | 55 ++++++++++++++++++++++++++++++++++ q2_ms/tests/__init__.py | 7 +++++ q2_ms/tests/tests_utils.py | 47 +++++++++++++++++++++++++++++ q2_ms/utils.py | 17 +++++++++++ q2_ms/xcms/__init__.py | 7 +++++ q2_ms/xcms/find_peaks_xcms.py | 49 ++++++++++++++++++++++++++++++ 7 files changed, 184 insertions(+) create mode 100644 q2_ms/assets/find_peaks_xcms.R create mode 100644 q2_ms/tests/__init__.py create mode 100644 q2_ms/tests/tests_utils.py create mode 100644 q2_ms/utils.py create mode 100644 q2_ms/xcms/__init__.py create mode 100644 q2_ms/xcms/find_peaks_xcms.py diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index cd91f32..cf79226 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -20,10 +20,12 @@ requirements: - versioningit run: + - bioconductor-xcms - pymzml - qiime2 {{ qiime2_epoch }}.* - q2-types {{ qiime2_epoch }}.* - q2templates {{ qiime2_epoch }}.* + - r-base test: requires: diff --git a/q2_ms/assets/find_peaks_xcms.R b/q2_ms/assets/find_peaks_xcms.R new file mode 100644 index 0000000..4a620a1 --- /dev/null +++ b/q2_ms/assets/find_peaks_xcms.R @@ -0,0 +1,55 @@ +#!/usr/bin/env Rscript + +library(xcms) +library(MsExperiment) +library(MsIO) + + +update_params <- function(params, opt) { + for (key in names(params)) { + if (!is.null(opt[[key]])) { + params[[key]] <- opt[[key]] + } + } + return(params) +} + +# Define command-line options +option_list <- list( + make_option(opt_str = "--mzml_files", type = "character"), + make_option(opt_str = "--ppm", type = "numeric"), + make_option(opt_str = "--min_peakwidth", type = "numeric"), + make_option(opt_str = "--max_peakwidth", type = "numeric"), + make_option(opt_str = "--snthresh", type = "numeric"), + make_option(opt_str = "--prefilter_k", type = "numeric"), + make_option(opt_str = "--prefilter_i", type = "numeric"), + make_option(opt_str = "--mz_center_fun", type = "numeric"), + make_option(opt_str = "--integrate", type = "integer"), + make_option(opt_str = "--mzdiff", type = "numeric"), + make_option(opt_str = "--fitgauss", type = "logical"), + make_option(opt_str = "--noise", type = "numeric"), + make_option(opt_str = "--first_baseline_check", type = "logical"), + make_option(opt_str = "--ms_level", type = "integer"), + make_option(opt_str = "--threads", type = "integer"), + make_option(opt_str = "--output_path", type = "character") +) + +# Parse arguments +opt_parser <- OptionParser(option_list = option_list) +opt <- parse_args(opt_parser) + +# Get full paths to mzML files and read them into an MsExperiment object +mzml_files <- list.files(opt$mzml_files, pattern = "\\.mzML$", full.names = TRUE) +msexperiment <- readMsExperiment(spectraFiles = mzml_files, sampleData = pd, BPPARAM = SerialParam()) + +# Load default parameters for CentWave +CentWaveParams <- CentWaveParam(method="centWave") + +# Update params with command-line arguments +CentWaveParams <- update_params(CentWaveParams, opt) + +# Find peaks using the CentWave algorithm +XCMSExperiment <- findChromPeaks(chr_raw, param = CentWaveParams) + +# Export the XCMSExperiment object to the directory format +saveMsObject(XCMSExperiment, param = PlainTextParam(path = opt$output_path)) diff --git a/q2_ms/tests/__init__.py b/q2_ms/tests/__init__.py new file mode 100644 index 0000000..688e933 --- /dev/null +++ b/q2_ms/tests/__init__.py @@ -0,0 +1,7 @@ +# ---------------------------------------------------------------------------- +# Copyright (c) 2024, QIIME 2 development team. +# +# Distributed under the terms of the Modified BSD License. +# +# The full license is in the file LICENSE, distributed with this software. +# ---------------------------------------------------------------------------- diff --git a/q2_ms/tests/tests_utils.py b/q2_ms/tests/tests_utils.py new file mode 100644 index 0000000..7ab96db --- /dev/null +++ b/q2_ms/tests/tests_utils.py @@ -0,0 +1,47 @@ +from unittest.mock import call, patch + +from qiime2.plugin.testing import TestPluginBase + +from q2_ms.utils import EXTERNAL_CMD_WARNING, run_command + + +class TestRunCommand(TestPluginBase): + package = "q2_ms.tests" + + @patch("subprocess.run") + @patch("builtins.print") + def test_run_command_verbose(self, mock_print, mock_subprocess_run): + # Mock command and working directory + cmd = ["echo", "Hello"] + cwd = "/test/directory" + + # Run the function with verbose=True + run_command(cmd, cwd=cwd, verbose=True) + + # Check if subprocess.run was called with the correct arguments + mock_subprocess_run.assert_called_once_with(cmd, check=True, cwd=cwd) + + # Check if the correct print statements were called + mock_print.assert_has_calls( + [ + call(EXTERNAL_CMD_WARNING), + call("\nCommand:", end=" "), + call("echo Hello", end="\n\n"), + ] + ) + + @patch("subprocess.run") + @patch("builtins.print") + def test_run_command_non_verbose(self, mock_print, mock_subprocess_run): + # Mock command and working directory + cmd = ["echo", "Hello"] + cwd = "/test/directory" + + # Run the function with verbose=False + run_command(cmd, cwd=cwd, verbose=False) + + # Check if subprocess.run was called with the correct arguments + mock_subprocess_run.assert_called_once_with(cmd, check=True, cwd=cwd) + + # Ensure no print statements were made + mock_print.assert_not_called() diff --git a/q2_ms/utils.py b/q2_ms/utils.py new file mode 100644 index 0000000..35c9fca --- /dev/null +++ b/q2_ms/utils.py @@ -0,0 +1,17 @@ +import subprocess + +EXTERNAL_CMD_WARNING = ( + "Running external command line application(s). " + "This may print messages to stdout and/or stderr.\n" + "The command(s) being run are below. These commands " + "cannot be manually re-run as they will depend on " + "temporary files that no longer exist." +) + + +def run_command(cmd, cwd, verbose=True): + if verbose: + print(EXTERNAL_CMD_WARNING) + print("\nCommand:", end=" ") + print(" ".join(cmd), end="\n\n") + subprocess.run(cmd, check=True, cwd=cwd) diff --git a/q2_ms/xcms/__init__.py b/q2_ms/xcms/__init__.py new file mode 100644 index 0000000..688e933 --- /dev/null +++ b/q2_ms/xcms/__init__.py @@ -0,0 +1,7 @@ +# ---------------------------------------------------------------------------- +# Copyright (c) 2024, QIIME 2 development team. +# +# Distributed under the terms of the Modified BSD License. +# +# The full license is in the file LICENSE, distributed with this software. +# ---------------------------------------------------------------------------- diff --git a/q2_ms/xcms/find_peaks_xcms.py b/q2_ms/xcms/find_peaks_xcms.py new file mode 100644 index 0000000..2a7d436 --- /dev/null +++ b/q2_ms/xcms/find_peaks_xcms.py @@ -0,0 +1,49 @@ +import copy +import subprocess + +from q2_ms.types import XCMSExperimentDirFmt, mzMLDirFmt +from q2_ms.utils import run_command + + +def find_peaks_xcms( + mzML: mzMLDirFmt, + ppm: float, + min_peakwidth: str, + max_peakwidth: str, + snthresh: float, + prefilter_k: str, + prefilter_i: int, + mz_center_fun: float, + integrate: int, + mzdiff: float, + fitgauss: bool, + noise: float, + first_baseline_check: float, + ms_level: int, + threads: int, +) -> XCMSExperimentDirFmt: + # Create parameters dict + params = copy.deepcopy(locals()) + + # Innit XCMSExperimentDirFmt + xcms_experiment = XCMSExperimentDirFmt() + + # Run R script + run_find_chrom_peaks(params) + + return xcms_experiment + + +def run_find_chrom_peaks(params): + cmd = ["find_peaks_xcms.R"] + for key, value in params.items(): + cmd.extend([f"--{key}", str(value)]) + + try: + run_command(cmd, verbose=True, cwd=None) + except subprocess.CalledProcessError as e: + raise Exception( + "An error was encountered while running XCMS, " + f"(return code {e.returncode}), please inspect " + "stdout and stderr to learn more." + )