Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add: implement tool to add rank names to phyloseq object #6625

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 4 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
72 changes: 72 additions & 0 deletions tools/phyloseq/add_rank_names_to_phyloseq.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@
#!/usr/bin/env Rscript

suppressPackageStartupMessages(library("optparse"))
suppressPackageStartupMessages(library("phyloseq"))
suppressPackageStartupMessages(library("tidyverse"))

# Option parsing
option_list <- list(
make_option(c("--input"),
action = "store", dest = "input",
help = "Input file containing a phyloseq object"
),
make_option(c("--output"),
action = "store", dest = "output",
help = "Output file for the updated phyloseq object"
),
make_option(c("--ranks"),
action = "store", dest = "ranks",
help = "Comma-separated list of taxonomy ranks (default: Kingdom,Phylum,Class,Order,Family,Genus,Species)"
)
)

parser <- OptionParser(usage = "%prog [options] file", option_list = option_list)
args <- parse_args(parser, positional_arguments = TRUE)
opt <- args$options

cat("Input file: ", opt$input, "\n")
cat("Output file: ", opt$output, "\n")
cat("Ranks provided: ", opt$ranks, "\n")

if (is.null(opt$ranks)) {
opt$ranks <- "Kingdom,Phylum,Class,Order,Family,Genus,Species"
}

# Parse rank names
rank_names <- unlist(strsplit(opt$ranks, ","))

# Load phyloseq object
physeq <- readRDS(opt$input)

# Check if physeq object is loaded successfully
if (is.null(physeq)) {
stop("Error: Failed to load the Phyloseq object. Check the input file.")
}

cat("Phyloseq object successfully loaded.\n")
cat("Class of loaded object: ", class(physeq), "\n")

# Check the current tax_table
cat("Current tax_table:\n")
print(tax_table(physeq))

# Add a column for Species if missing
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this part is needed (the add species part), we cannot really know which rank is missing if it is missing, the user should provide names matching the columns exactly, otherwise the tool should fail, you could change if (ncol(tax_table(physeq)) > length(rank_names)) { to if (ncol(tax_table(physeq)) == length(rank_names)) { instead, and provide a message for further help

if (ncol(tax_table(physeq)) < length(rank_names)) {
cat("Warning: Taxonomy table has fewer columns than the provided ranks. Filling missing columns with NA.\n")
while (ncol(tax_table(physeq)) < length(rank_names)) {
tax_table(physeq) <- cbind(tax_table(physeq), NA)
}
} else if (ncol(tax_table(physeq)) > length(rank_names)) {
stop("Error: Number of columns in tax_table exceeds the number of provided ranks. Please adjust your rank list.")
}

# Set column names to the provided ranks
colnames(tax_table(physeq)) <- rank_names

# Confirm the changes
cat("Updated tax_table:\n")
print(tax_table(physeq))

# Save the updated phyloseq object
saveRDS(physeq, file = opt$output, compress = TRUE)
cat("Updated Phyloseq object saved to: ", opt$output, "\n")
37 changes: 37 additions & 0 deletions tools/phyloseq/add_rank_names_to_phyloseq.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
<tool id="add_rank_names_to_phyloseq" name="Add Rank Names to Phyloseq Object" version="@TOOL_VERSION@+galaxy@VERSION_SUFFIX@" profile="@PROFILE@">
<description>Add taxonomy rank names to a phyloseq object</description>
<macros>
<import>macros.xml</import>
</macros>
<expand macro="bio_tools"/>
<expand macro="requirements"/>
<command detect_errors="exit_code"><![CDATA[
Rscript '${__tool_directory__}/add_rank_names_to_phyloseq.R' --input '$input' --output '$output' --ranks '$ranks'
]]></command>

<inputs>
<expand macro="phyloseq_input"/>
<param name="ranks" type="text" label="Comma-separated list of taxonomy ranks" value="Kingdom,Phylum,Class,Order,Family,Genus,Species"/>
</inputs>

<outputs>
<data name="output" format="phyloseq"/>
</outputs>

<tests>
<test>
<param name="input" value="output.phyloseq" ftype="phyloseq"/>
<param name="ranks" value="Kingdom,Phylum,Class,Order,Family,Genus,Species"/>
<output name="output" ftype="phyloseq">
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you could compare the new phylosq object to expected output, like:

<output name="ampvis" value="AalborgWWTPs-subset_samples.rds" ftype="ampvis2" compare="sim_size"/>

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Better just use an has_size assertion and remove value="AalborgWWTPs-subset_samples.rds".

<assert_contents>
<has_size value="87125" delta="1000"/>
</assert_contents>
</output>
</test>
</tests>

<help>
This tool adds taxonomy rank names to a phyloseq object in the `tax_table` slot.
</help>
<expand macro="citations"/>
</tool>
Loading