Mercurial > repos > iuc > ancombc

comparison ancombc.R @ 0:939c59ab61cf draft

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/ancombc commit 045979180e44c683b5e0760f802af66c05abcae8
author iuc
date Sat, 16 Jul 2022 21:40:04 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:939c59ab61cf
1 #!/usr/bin/env Rscript
2
3 suppressPackageStartupMessages(library("ANCOMBC"))
4 suppressPackageStartupMessages(library("data.table"))
5 suppressPackageStartupMessages(library("optparse"))
6
7 option_list <- list(
8 make_option(c("--phyloseq"), action = "store", dest = "phyloseq", help = "File containing a phyloseq object"),
9 make_option(c("--formula"), action = "store", dest = "formula", help = "Formula"),
10 make_option(c("--p_adj_method"), action = "store", dest = "p_adj_method", help = "Method to adjust p-values"),
11 make_option(c("--zero_cut"), action = "store", dest = "zero_cut", type = "double", help = "Minimum taxa prevalence"),
12 make_option(c("--lib_cut"), action = "store", dest = "lib_cut", type = "integer", help = "Thrshold for filtering samples based on library sizes"),
13 make_option(c("--group"), action = "store", dest = "group", help = "Name of the group variable in the metadata"),
14 make_option(c("--struc_zero"), action = "store", dest = "struc_zero", help = "Detect structural zeros based on group"),
15 make_option(c("--neg_lb"), action = "store", dest = "neg_lb", help = "Classify a taxon as a structural zero using its asymptotic lower bound"),
16 make_option(c("--tol"), action = "store", dest = "tol", type = "double", help = "Iteration convergence tolerance for the E-M algorithm"),
17 make_option(c("--max_iter"), action = "store", dest = "max_iter", help = "Maximum number of iterations for the E-M algorithm"),
18 make_option(c("--conserve"), action = "store", dest = "conserve", help = "Use a conservative variance estimator for the test statistic"),
19 make_option(c("--alpha"), action = "store", dest = "alpha", help = "Level of significance"),
20 make_option(c("--global"), action = "store", dest = "global", help = "Perform global test"),
21 make_option(c("--output_dir"), action = "store", dest = "output_dir", help = "Output directory")
22 )
23
24 parser <- OptionParser(usage = "%prog [options] file", option_list = option_list)
25 args <- parse_args(parser, positional_arguments = TRUE)
26 opt <- args$options
27
28 get_boolean_value <- function(val) {
29 if (val == "true") {
30 return(TRUE)
31 } else {
32 return(FALSE)
33 }
34 }
35
36 get_file_path <- function(dir, file_name) {
37 file_path <- paste(dir, file_name, sep = "/")
38 return(file_path)
39 }
40
41 write_data_frame <- function(dir, file_name, data_frame) {
42 file_path <- get_file_path(dir, file_name)
43 write.table(data_frame, file = file_path, quote = FALSE, row.names = TRUE, col.names = TRUE, sep = "\t")
44 }
45
46 # Convert boolean values to boolean.
47 struc_zero <- get_boolean_value(opt$struc_zero)
48 neg_lb <- get_boolean_value(opt$neg_lb)
49 conserve <- get_boolean_value(opt$conserve)
50 global <- get_boolean_value(opt$global)
51
52 # Construct a phyloseq object.
53 phyloseq_obj <- readRDS(opt$phyloseq)
54
55 # Construct an ANCOM-BC object.
56 ancombc_obj <- ancombc(phyloseq = phyloseq_obj,
57 formula = opt$formula,
58 p_adj_method = opt$p_adj_method,
59 zero_cut = opt$zero_cut,
60 lib_cut = opt$lib_cut,
61 group = opt$group,
62 struc_zero = struc_zero,
63 neg_lb = neg_lb,
64 tol = opt$tol,
65 max_iter = opt$max_iter,
66 conserve = conserve,
67 alpha = opt$alpha,
68 global = global)
69
70 res <- ancombc_obj$res
71
72 # Write the outputs.
73 write_data_frame(opt$output_dir, "feature_table.tabular", ancombc_obj$feature_table)
74 write_data_frame(opt$output_dir, "zero_ind.tabular", ancombc_obj$zero_ind)
75 write.csv2(ancombc_obj$samp_frac, file = get_file_path(opt$output_dir, "samp_frac.tabular"), row.names = FALSE, col.names = FALSE, sep = "\t")
76 write_data_frame(opt$output_dir, "resid.tabular", ancombc_obj$resid)
77 write(ancombc_obj$delta_em, file = get_file_path(opt$output_dir, "delta_em.tabular"))
78 write(ancombc_obj$delta_wls, file = get_file_path(opt$output_dir, "delta_wls.tabular"))
79 write_data_frame(opt$output_dir, "res_beta.tabular", res$beta)
80 write_data_frame(opt$output_dir, "res_se.tabular", res$se)
81 write_data_frame(opt$output_dir, "res_W.tabular", res$W)
82 write_data_frame(opt$output_dir, "res_p_val.tabular", res$p_val)
83 write_data_frame(opt$output_dir, "res_q_val.tabular", res$q_val)
84 write_data_frame(opt$output_dir, "res_diff_abn.tabular", res$diff_abn)
85 write_data_frame(opt$output_dir, "res_global.tabular", ancombc_obj$res_global)