Mercurial > repos > iuc > phyloseq_plot_ordination

comparison phyloseq_tax_glom.R @ 10:cee4982a717b draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/phyloseq commit fa72d860839082a926004d8a97a03a3e27701333
author iuc
date Fri, 04 Apr 2025 10:16:43 +0000
parents
children
comparison
equal deleted inserted replaced
9:32f892954ef6 10:cee4982a717b
1 suppressPackageStartupMessages(library("phyloseq"))
2 suppressPackageStartupMessages(library("dplyr"))
3 suppressPackageStartupMessages(library("optparse"))
4
5 # Define command-line options
6 option_list <- list(
7 make_option(c("-i", "--input"), type = "character", help = "Path to the phyloseq RDS file", metavar = "FILE"),
8 make_option(c("-r", "--rank"), type = "character", help = "Taxonomic rank for aggregation"),
9 make_option("--exclude_otu_ids", action = "store_true", default = FALSE, help = "Exclude OTU IDs from output"),
10 make_option("--single_rank", action = "store_true", default = FALSE, help = "Only output the specified rank column"),
11 make_option("--exclude_na_values", action = "store_true", default = FALSE, help = "Exclude NA values during tax_glom")
12 )
13
14 # Parse arguments
15 opt <- parse_args(OptionParser(option_list = option_list))
16
17 # Validate arguments
18 if (is.null(opt$input) || is.null(opt$rank)) {
19 stop("Error: --input and --rank are required arguments.")
20 }
21
22 if (opt$single_rank && !opt$exclude_otu_ids) {
23 stop("Error: --single_rank can only be used if --exclude_otu_ids is also specified.")
24 }
25
26 # Load the phyloseq object
27 physeq <- readRDS(opt$input)
28
29 # Print available taxonomic ranks
30 cat("Available taxonomic ranks:\n")
31 print(rank_names(physeq))
32
33 # Print original number of OTUs
34 cat("Original number of OTUs:", ntaxa(physeq), "\n")
35
36 # Perform tax_glom
37 physeq_agg <- tax_glom(physeq, taxrank = opt$rank, NArm = opt$exclude_na_values)
38
39 # Print new number of taxa after agglomeration
40 cat("Number of taxa after agglomeration at", opt$rank, "level:", ntaxa(physeq_agg), "\n")
41
42 # Extract the taxonomy table after agglomeration
43 tax_table_agg <- as.data.frame(tax_table(physeq_agg))
44
45 # Convert taxonomic columns to character to preserve NA values
46 tax_table_agg[] <- lapply(tax_table_agg, as.character)
47
48 # Add OTU ID column unless excluded
49 if (!opt$exclude_otu_ids) {
50 tax_table_agg <- cbind("OTU ID" = rownames(tax_table_agg), tax_table_agg)
51 }
52
53 # Extract OTU abundance table and convert to data frame
54 otu_table_agg <- as.data.frame(otu_table(physeq_agg))
55
56 # Append taxonomic information to output
57 otu_table_agg <- cbind(tax_table_agg, otu_table_agg)
58
59 tax_table_agg <- otu_table_agg
60
61 if (opt$single_rank) {
62 # Keep only the specified taxonomic rank column and numeric count columns
63 tax_table_agg <- tax_table_agg %>% select(all_of(opt$rank), where(is.numeric))
64
65 # Group by taxonomic rank and sum the counts
66 tax_table_agg <- tax_table_agg %>%
67 group_by(across(all_of(opt$rank))) %>%
68 summarise(across(where(is.numeric), sum), .groups = "drop")
69 }
70
71 # Save the output as a TSV file
72 output_file <- paste0("physeq_", opt$rank, "_table.tsv")
73 write.table(tax_table_agg, file = output_file, sep = "\t", quote = FALSE, row.names = FALSE, col.names = TRUE)