|
| 1 | +#!/usr/bin/env Rscript |
| 2 | + |
| 3 | +suppressPackageStartupMessages(library("optparse")) |
| 4 | +suppressPackageStartupMessages(library("phyloseq")) |
| 5 | +suppressPackageStartupMessages(library("tidyverse")) |
| 6 | + |
| 7 | +# Option parsing |
| 8 | +option_list <- list( |
| 9 | + make_option(c("--input"), |
| 10 | + action = "store", dest = "input", |
| 11 | + help = "Input file containing a phyloseq object" |
| 12 | + ), |
| 13 | + make_option(c("--output"), |
| 14 | + action = "store", dest = "output", |
| 15 | + help = "Output file for the updated phyloseq object" |
| 16 | + ), |
| 17 | + make_option(c("--ranks"), |
| 18 | + action = "store", dest = "ranks", |
| 19 | + help = "Comma-separated list of taxonomy ranks (default: Kingdom,Phylum,Class,Order,Family,Genus,Species)" |
| 20 | + ) |
| 21 | +) |
| 22 | + |
| 23 | +parser <- OptionParser(usage = "%prog [options] file", option_list = option_list) |
| 24 | +args <- parse_args(parser, positional_arguments = TRUE) |
| 25 | +opt <- args$options |
| 26 | + |
| 27 | +cat("Input file: ", opt$input, "\n") |
| 28 | +cat("Output file: ", opt$output, "\n") |
| 29 | +cat("Ranks provided: ", opt$ranks, "\n") |
| 30 | + |
| 31 | +if (is.null(opt$ranks)) { |
| 32 | + opt$ranks <- "Kingdom,Phylum,Class,Order,Family,Genus,Species" |
| 33 | +} |
| 34 | + |
| 35 | +# Parse rank names |
| 36 | +rank_names <- unlist(strsplit(opt$ranks, ",")) |
| 37 | + |
| 38 | +# Load phyloseq object |
| 39 | +physeq <- readRDS(opt$input) |
| 40 | + |
| 41 | +# Check if physeq object is loaded successfully |
| 42 | +if (is.null(physeq)) { |
| 43 | + stop("Error: Failed to load the Phyloseq object. Check the input file.") |
| 44 | +} |
| 45 | + |
| 46 | +cat("Phyloseq object successfully loaded.\n") |
| 47 | +cat("Class of loaded object: ", class(physeq), "\n") |
| 48 | + |
| 49 | +# Check the current tax_table |
| 50 | +cat("Current tax_table:\n") |
| 51 | +print(tax_table(physeq)) |
| 52 | + |
| 53 | + |
| 54 | +# Strict check for taxonomy table and provided ranks |
| 55 | +if (ncol(tax_table(physeq)) != length(rank_names)) { |
| 56 | + stop( |
| 57 | + "Error: Number of columns in tax_table does not match the number of provided ranks. ", |
| 58 | + "Please ensure the taxonomy table matches the ranks exactly." |
| 59 | + ) |
| 60 | +} |
| 61 | + |
| 62 | +# Set column names to the provided ranks |
| 63 | +colnames(tax_table(physeq)) <- rank_names |
| 64 | + |
| 65 | +# Confirm the changes |
| 66 | +cat("Updated tax_table:\n") |
| 67 | +print(tax_table(physeq)) |
| 68 | + |
| 69 | +# Save the updated phyloseq object |
| 70 | +saveRDS(physeq, file = opt$output, compress = TRUE) |
| 71 | +cat("Updated Phyloseq object saved to: ", opt$output, "\n") |
0 commit comments