Mercurial > repos > iuc > phyloseq_plot_richness

comparison phyloseq_from_dada2.R @ 2:fb7c4bbe8994 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 5ec9f9e81bb9a42dec5c331dd23215ca0b027b2b
author iuc
date Sat, 16 Mar 2024 07:55:57 +0000
parents 1ff178d1757e
children
comparison
equal deleted inserted replaced
1:1ff178d1757e 2:fb7c4bbe8994
5 suppressPackageStartupMessages(library("tidyverse")) 5 suppressPackageStartupMessages(library("tidyverse"))
6 6
7 option_list <- list( 7 option_list <- list(
8 make_option(c("--sequence_table"), action = "store", dest = "sequence_table", help = "Input sequence table"), 8 make_option(c("--sequence_table"), action = "store", dest = "sequence_table", help = "Input sequence table"),
9 make_option(c("--taxonomy_table"), action = "store", dest = "taxonomy_table", help = "Input taxonomy table"), 9 make_option(c("--taxonomy_table"), action = "store", dest = "taxonomy_table", help = "Input taxonomy table"),
10 make_option(c("--sample_table"), action = "store", default = NULL, dest = "sample_table", help = "Input sample table"),
10 make_option(c("--output"), action = "store", dest = "output", help = "RDS output") 11 make_option(c("--output"), action = "store", dest = "output", help = "RDS output")
11 ) 12 )
12 13
13 parser <- OptionParser(usage = "%prog [options] file", option_list = option_list) 14 parser <- OptionParser(usage = "%prog [options] file", option_list = option_list)
14 args <- parse_args(parser, positional_arguments = TRUE) 15 args <- parse_args(parser, positional_arguments = TRUE)
15 opt <- args$options 16 opt <- args$options
16 # The input sequence_table is an integer matrix 17 # The input sequence_table is an integer matrix
17 # stored as tabular (rows = samples, columns = ASVs). 18 # stored as tabular (rows = samples, columns = ASVs).
18 seq_table_numeric_matrix <- data.matrix(read.table(opt$sequence_table, sep = "\t")) 19 seq_table_numeric_matrix <- data.matrix(read.table(opt$sequence_table, header = T, sep = "\t", row.names = 1, check.names = FALSE))
19 # The input taxonomy_table is a table containing 20 # The input taxonomy_table is a table containing
20 # the assigned taxonomies exceeding the minBoot 21 # the assigned taxonomies exceeding the minBoot
21 # level of bootstrapping confidence. Rows correspond 22 # level of bootstrapping confidence. Rows correspond
22 # to sequences, columns to taxonomic levels. NA 23 # to sequences, columns to taxonomic levels. NA
23 # indicates that the sequence was not consistently 24 # indicates that the sequence was not consistently
24 # classified at that level at the minBoot threshold. 25 # classified at that level at the minBoot threshold.
25 tax_table_matrix <- as.matrix(read.table(opt$taxonomy_table, header = FALSE, sep = "\t")) 26 tax_table_matrix <- as.matrix(read.table(opt$taxonomy_table, header = T, sep = "\t", row.names = 1, check.names = FALSE))
26 # Construct a tax_table object. The rownames of 27 # Construct a tax_table object. The rownames of
27 # tax_tab must match the OTU names (taxa_names) 28 # tax_tab must match the OTU names (taxa_names)
28 # of the otu_table defined below. 29 # of the otu_table defined below.
29 tax_tab <- tax_table(tax_table_matrix) 30 tax_tab <- tax_table(tax_table_matrix)
31
30 # Construct an otu_table object. 32 # Construct an otu_table object.
31 otu_tab <- otu_table(seq_table_numeric_matrix, taxa_are_rows = TRUE) 33 otu_tab <- otu_table(seq_table_numeric_matrix, taxa_are_rows = TRUE)
34
32 # Construct a phyloseq object. 35 # Construct a phyloseq object.
33 phyloseq_obj <- phyloseq(otu_tab, tax_tab) 36 phyloseq_obj <- phyloseq(otu_tab, tax_tab)
37 if (!is.null(opt$sample_table)) {
38 sample_tab <- sample_data(
39 read.table(opt$sample_table, header = T, sep = "\t", row.names = 1, check.names = FALSE)
40 )
41 phyloseq_obj <- merge_phyloseq(phyloseq_obj, sample_tab)
42 }
43
44 # use short names for our ASVs and save the ASV sequences
45 # refseq slot of the phyloseq object as described in
46 # https://benjjneb.github.io/dada2/tutorial.html
47 dna <- Biostrings::DNAStringSet(taxa_names(phyloseq_obj))
48 names(dna) <- taxa_names(phyloseq_obj)
49 phyloseq_obj <- merge_phyloseq(phyloseq_obj, dna)
50 taxa_names(phyloseq_obj) <- paste0("ASV", seq(ntaxa(phyloseq_obj)))
51
52 print(phyloseq_obj)
53
54 # save R object to file
34 saveRDS(phyloseq_obj, file = opt$output, compress = TRUE) 55 saveRDS(phyloseq_obj, file = opt$output, compress = TRUE)