Skip to content

Commit

Permalink
added utils find peaks
Browse files Browse the repository at this point in the history
  • Loading branch information
VinzentRisch committed Jan 15, 2025
1 parent 6915cd2 commit 880cd6b
Show file tree
Hide file tree
Showing 7 changed files with 184 additions and 0 deletions.
2 changes: 2 additions & 0 deletions conda-recipe/meta.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -20,10 +20,12 @@ requirements:
- versioningit

run:
- bioconductor-xcms
- pymzml
- qiime2 {{ qiime2_epoch }}.*
- q2-types {{ qiime2_epoch }}.*
- q2templates {{ qiime2_epoch }}.*
- r-base

test:
requires:
Expand Down
55 changes: 55 additions & 0 deletions q2_ms/assets/find_peaks_xcms.R
Original file line number Diff line number Diff line change
@@ -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))
7 changes: 7 additions & 0 deletions q2_ms/tests/__init__.py
Original file line number Diff line number Diff line change
@@ -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.
# ----------------------------------------------------------------------------
47 changes: 47 additions & 0 deletions q2_ms/tests/tests_utils.py
Original file line number Diff line number Diff line change
@@ -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()
17 changes: 17 additions & 0 deletions q2_ms/utils.py
Original file line number Diff line number Diff line change
@@ -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)
7 changes: 7 additions & 0 deletions q2_ms/xcms/__init__.py
Original file line number Diff line number Diff line change
@@ -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.
# ----------------------------------------------------------------------------
49 changes: 49 additions & 0 deletions q2_ms/xcms/find_peaks_xcms.py
Original file line number Diff line number Diff line change
@@ -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."
)

0 comments on commit 880cd6b

Please sign in to comment.