Merge pull request #111 from adrientaudiere/dev

v 0.12.1

Author: adrientaudiere

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: MiscMetabar
 Type: Package
 Title: Miscellaneous Functions for Metabarcoding Analysis
-Version: 0.12.0
+Version: 0.12.1
 Authors@R: person("Adrien", "Taudière", email = "adrien.taudiere@zaclys.net",
   role = c("aut", "cre", "cph"), comment = c(ORCID = "0000-0003-1088-1182"))
 Description: Facilitate the description, transformation, exploration, and reproducibility of metabarcoding analyses. 'MiscMetabar' is mainly built on top of the 'phyloseq', 'dada2' and 'targets' R packages. It helps to build reproducible and robust bioinformatics pipelines in R. 'MiscMetabar' makes ecological analysis of alpha and beta-diversity easier, more reproducible and more powerful by integrating a large number of tools. Important features are described in Taudière A. (2023) <doi:10.21105/joss.06038>.

NEWS.md

@@ -1,12 +1,18 @@
-# MiscMetabar 0.2.0 (in development) 
+# MiscMetabar 0.12.1 (in development) 
+
+- Add parameter name `min_bootstrap` in [add_new_taxonomy_pq()]
+- Bug fix in [assign_idtaxa()]
+- Add parameters `pattern_to_remove` and `remove_NA` to [simplify_taxo()]
+
+# MiscMetabar 0.12.0 
 
 - Add function [assign_idtaxa()] and [learn_idtaxa()] to facilitate the taxonomic assignation using the idtaxa algorithm from the  DECIPHER R package.
 - Add option `idtaxa` to method in [add_new_taxonomy_pq()]
 - Add function [tbl_sum_taxtable()] to summarize tax_table from a phyloseq object 
 - In function [assign_sintax()], add params `too_few` (default value "align_start") and `too_many` (default "merge") to authorize db with variable numbers of rank and parenthesis in taxonomic name, 
 
 
-# MiscMetabar 0.11.1 (in development) 
+# MiscMetabar 0.11.1 
 
 - Add param `suffix` to `add_blast_info()` allowing multiple use of the function on the same phyloseq object (e.g. in order to used different database)
 - Add param `return_DNAStringSet` to `write_temp_fasta()` function to return a DNAStringSet object in place of a temporary file. 

R/dada_phyloseq.R

@@ -1767,55 +1767,63 @@ select_one_sample <- function(physeq, sam_name, silent = FALSE) {
 #'   Note that if trainingSet is not NULL, the ref_fasta is overwrite by the
 #'   trainingSet parameter. To customize learning parameters of the idtaxa
 #'   algorithm you must use trainingSet computed by the function [learn_idtaxa()].
-#' @param min_boostrap  (Int. \[0:1\])
-#' 
+#' @param min_bootstrap  (Int. \[0:1\])
+#'
 #'   Minimum bootstrap value to inform taxonomy. For each bootstrap
-#'   below the min_boostrap value, the taxonomy information is set to NA.
-#' 
+#'   below the min_bootstrap value, the taxonomy information is set to NA.
+#'
 #'   Correspond to parameters :
-#' 
+#'
 #'  - dada2: `minBoot`, default value = 0.5
-#'  
-#'  - sintax: `min_boostrap`, default value = 60
-#' 
-#'  - lca: `threshold`, default value = 60
-#' 
-#'  - idtaxa: `threshold`, default value = 60
-#' 
+#'
+#'  - sintax: `min_bootstrap`, default value = 0.5
+#'
+#'  - lca: `id`, default value = 0.5. Note in that case, the bootstrap value is different.
+#'    See the id parameter in [assign_vsearch_lca()]
+#'
+#'  - idtaxa: `threshold`, default value = 0.6
+#'
 #' @param ... Other arguments passed on to the taxonomic assignation method.
 #' @return A new \code{\link[phyloseq]{phyloseq-class}} object with a larger slot tax_table"
 #' @seealso [dada2::assignTaxonomy()], [assign_sintax()], [assign_vsearch_lca()], [assign_sintax()]
 #' @export
 #'
 #' @author Adrien Taudière
 #'
-add_new_taxonomy_pq <- function(physeq, ref_fasta, suffix = NULL, method = c("dada2", "sintax", "lca", "idtaxa"), trainingSet = NULL, min_boostrap = NULL, ...) {
+add_new_taxonomy_pq <- function(
+    physeq,
+    ref_fasta,
+    suffix = NULL,
+    method = c("dada2", "sintax", "lca", "idtaxa"),
+    trainingSet = NULL,
+    min_bootstrap = NULL,
+    ...) {
   method <- match.arg(method)
 
-  if(is.null(min_boostrap)){
-    min_boostrap <- ifelse(method=="idtaxa", 0.6, 0.5)
+  if (is.null(min_bootstrap)) {
+    min_bootstrap <- ifelse(method == "idtaxa", 0.6, 0.5)
   }
 
   if (is.null(suffix)) {
     suffix <- paste0("_", basename(ref_fasta), "_", method)
   }
   if (method == "dada2") {
     tax_tab <-
-      dada2::assignTaxonomy(physeq@refseq, refFasta = ref_fasta, minBoot = 100 * min_boostrap, ...)
+      dada2::assignTaxonomy(physeq@refseq, refFasta = ref_fasta, minBoot = 100 * min_bootstrap, ...)
     colnames(tax_tab) <-
       make.unique(paste0(colnames(tax_tab), suffix))
     new_tax_tab <- tax_table(cbind(physeq@tax_table, tax_tab))
     new_physeq <- physeq
     tax_table(new_physeq) <- new_tax_tab
   } else if (method == "sintax") {
-    new_physeq <- assign_sintax(physeq, ref_fasta = ref_fasta, suffix = suffix, behavior = "add_to_phyloseq", min_boostrap = min_boostrap, ...)
+    new_physeq <- assign_sintax(physeq, ref_fasta = ref_fasta, suffix = suffix, behavior = "add_to_phyloseq", min_bootstrap = min_bootstrap, ...)
   } else if (method == "lca") {
     new_physeq <- assign_vsearch_lca(physeq, ref_fasta = ref_fasta, suffix = suffix, behavior = "add_to_phyloseq", ...)
   } else if (method == "idtaxa") {
     if (is.null(trainingSet)) {
-      new_physeq <- assign_idtaxa(physeq, seq2search = ref_fasta, suffix = suffix, threshold = 100 * min_boostrap, ...)
+      new_physeq <- assign_idtaxa(physeq, fasta_for_training = ref_fasta, behavior = "add_to_phyloseq", suffix = suffix, threshold = 100 * min_bootstrap, ...)
     } else {
-      new_physeq <- assign_idtaxa(physeq, trainingSet = trainingSet, suffix = suffix, ...)
+      new_physeq <- assign_idtaxa(physeq, trainingSet = trainingSet, behavior = "add_to_phyloseq", suffix = suffix, threshold = 100 * min_bootstrap, ...)
     }
   }
 
@@ -1927,7 +1935,8 @@ tbl_sum_taxtable <- function(physeq, taxonomic_ranks = NULL, ...) {
 #'
 #' @inheritParams clean_pq
 #' @param taxLevels Name of the 7 columns in tax_table required by funguild
-#'
+#' @param db_url a length 1 character string giving the URL to retrieve the database
+#'     from
 #' @return A new object of class `physeq` with Guild information added to
 #'   `tax_table` slot
 #' @export
@@ -1952,7 +1961,7 @@ tbl_sum_taxtable <- function(physeq, taxonomic_ranks = NULL, ...) {
 #' @details
 #' This function is mainly a wrapper of the work of others.
 #'   Please make a reference to `FUNGuildR` package and the associate
-#'   publication (\doi{10.1016/j.funeco.2015.06.006}) if you
+#'   publication (\doi{doi:10.1016/j.funeco.2015.06.006}) if you
 #'   use this function.
 #' @seealso [plot_guild_pq()]
 
@@ -1965,13 +1974,18 @@ add_funguild_info <- function(physeq,
                                 "Family",
                                 "Genus",
                                 "Species"
-                              )) {
+                              ),
+                              db_url = "http://www.stbates.org/funguild_db_2.php") {
   tax_tab <- physeq@tax_table
   FUNGuild_assign <-
     funguild_assign(data.frame(
       "Taxonomy" =
         apply(tax_tab[, taxLevels], 1, paste, collapse = ";")
-    ))
+    ), db_url = db_url)
+  if (is.null(FUNGuild_assign)) {
+    message("No http access to the funguild database. No information were added.")
+    return(physeq)
+  }
   tax_tab <-
     as.matrix(cbind(tax_tab, FUNGuild_assign))
   physeq@tax_table <- tax_table(tax_tab)
@@ -3156,6 +3170,11 @@ rarefy_sample_count_by_modality <-
 #'
 #'  - "add_to_phyloseq" return a phyloseq object with amended slot `@taxtable`.
 #'    Only available if using physeq input and not seq2search input.
+#' @param threshold (Int, default 60) Numeric specifying the confidence at which
+#'    to truncate the output taxonomic classifications.
+#'   Lower values of threshold will classify deeper into the taxonomic tree at
+#'   the expense of accuracy, and vise-versa for higher values of threshold. See
+#'   [DECIPHER::IdTaxa()] man page.
 #' @param column_names (vector of character) names for the column of the
 #'   taxonomy
 #' @param suffix (character) The suffix to name the new columns.
@@ -3199,11 +3218,12 @@ rarefy_sample_count_by_modality <-
 assign_idtaxa <- function(physeq,
                           seq2search = NULL,
                           trainingSet = NULL,
-                          fasta_for_training,
+                          fasta_for_training = NULL,
                           behavior = "return_matrix",
+                          threshold = 60,
                           column_names = c(
                             "Kingdom",
-                            "Phyla",
+                            "Phylum",
                             "Class",
                             "Order",
                             "Family",
@@ -3215,7 +3235,11 @@ assign_idtaxa <- function(physeq,
                           unite = FALSE,
                           verbose = TRUE,
                           ...) {
-  if (is.null(trainingSet) && !is.null(fasta_for_training)) {
+  if (!is.null(trainingSet) && !is.null(fasta_for_training)) {
+    stop("Please provide either trainingSet or fasta_for_training parameters, not both.")
+  } else if (is.null(trainingSet) && is.null(fasta_for_training)) {
+    stop("Please provide either trainingSet or fasta_for_training parameters.")
+  } else if (is.null(trainingSet) && !is.null(fasta_for_training)) {
     if (verbose) {
       message("Training using fasta_for_training file.")
     }
@@ -3230,7 +3254,7 @@ assign_idtaxa <- function(physeq,
     return_DNAStringSet = TRUE
   )
 
-  fasta2search <- OrientNucleotides(RemoveGaps(fasta2search))
+  fasta2search <- DECIPHER::OrientNucleotides(DECIPHER::RemoveGaps(fasta2search))
 
   if (verbose) {
     message("Classifing using training Set with IdTaxa.")
@@ -3240,6 +3264,7 @@ assign_idtaxa <- function(physeq,
     test = fasta2search,
     trainingSet = trainingSet_idtaxa,
     processors = nproc,
+    threshold = threshold,
     ...
   )
 
@@ -3352,9 +3377,9 @@ assign_idtaxa <- function(physeq,
 #'   Please make a reference to [DECIPHER::LearnTaxa()] if you
 #'   use this function.
 learn_idtaxa <- function(fasta_for_training, output_Rdata = NULL, output_path_only = FALSE, unite = FALSE, ...) {
-  seqs <- readDNAStringSet(fasta_for_training)
-  seqs <- RemoveGaps(seqs)
-  seqs <- OrientNucleotides(seqs)
+  seqs <- Biostrings::readDNAStringSet(fasta_for_training)
+  seqs <- DECIPHER::RemoveGaps(seqs)
+  seqs <- DECIPHER::OrientNucleotides(seqs)
 
   taxo_for_learning <- names(seqs)
   if (unite) {
@@ -3374,7 +3399,7 @@ learn_idtaxa <- function(fasta_for_training, output_Rdata = NULL, output_path_on
   taxo_for_learning[!grepl("^Root;", taxo_for_learning)] <-
     paste0("Root;", taxo_for_learning[!grepl("^Root;", taxo_for_learning)])
 
-  train_idtaxa <- LearnTaxa(seqs, taxonomy = taxo_for_learning, ...)
+  train_idtaxa <- DECIPHER::LearnTaxa(seqs, taxonomy = taxo_for_learning, ...)
   if (!is.null(output_Rdata)) {
     save(train_idtaxa, output_Rdata)
   }

R/funguild.R

@@ -24,6 +24,10 @@
 #' modified by Adrien Taudière
 #'
 get_funguild_db <- function(db_url = "http://www.stbates.org/funguild_db_2.php") {
+  if (httr::http_error(db_url)) {
+    return(NULL)
+  }
+
   httr::GET(url = db_url) %>%
     httr::content(as = "text") %>%
     stringr::str_split("\n") %>%
@@ -81,13 +85,17 @@ get_funguild_db <- function(db_url = "http://www.stbates.org/funguild_db_2.php")
 #' "`p__`", ...) are also allowed.
 #' A `character` vector, representing only the taxonomic classification,
 #' is also accepted.
-#' @param tax_col A `character` string, optionally giving an alternate
-#' column name in `otu_table` to use instead of `otu_table$Taxonomy`.
+#'
+#' @param db_url a length 1 character string giving the URL to retrieve the database
+#'     from
 #'
 #' @param db_funguild A `data.frame` representing the FUNGuild as returned by
 #' [get_funguild_db()]
 #' If not supplied, the default database will be downloaded.
 #'
+#' @param tax_col A `character` string, optionally giving an alternate
+#' column name in `otu_table` to use instead of `otu_table$Taxonomy`.
+#'
 #' @return A [`tibble::tibble`] containing all columns of
 #' `otu_table`, plus relevant columns of information from the FUNGuild
 #' @export
@@ -98,8 +106,18 @@ get_funguild_db <- function(db_url = "http://www.stbates.org/funguild_db_2.php")
 #' 20:241-248.
 #' @author Brendan Furneaux (orcid: [0000-0003-3522-7363](https://orcid.org/0000-0003-3522-7363)),
 #' modified by Adrien Taudière
-funguild_assign <- function(otu_table, db_funguild = get_funguild_db(),
-                            tax_col = "Taxonomy") {
+funguild_assign <- function(
+    otu_table,
+    db_url = NULL,
+    db_funguild = NULL,
+    tax_col = "Taxonomy") {
+  if (is.null(db_funguild)) {
+    db_funguild <- get_funguild_db(db_url = db_url)
+  }
+
+  if (is.null(db_funguild)) {
+    return(NULL)
+  }
   if (is.character(otu_table)) {
     otu_table <- tibble::tibble(otu_table)
     names(otu_table) <- tax_col
@@ -120,7 +138,8 @@ funguild_assign <- function(otu_table, db_funguild = get_funguild_db(),
 
   otu_table$taxkey <- make_taxkey(otu_table[[tax_col]])
   all_taxkey <- unique(otu_table$taxkey) %>% na.omit()
-  `.` <- taxon <- taxkey <- searchkey <- taxonomicLevel <- NULL # to pass R CMD check
+  `.` <- taxon <- taxkey <- searchkey <- taxonomicLevel <- NULL
+
   db_funguild <- dplyr::mutate(
     db_funguild,
     searchkey = paste0("@", stringr::str_replace(taxon, "[ _]", "@"), "@")

R/miscellanous.R

@@ -126,18 +126,49 @@ all_object_size <- function() {
 #'
 #' Internally used in [clean_pq()]
 #' @inheritParams clean_pq
+#' @param pattern_to_remove (a vector of character) the pattern to remove using [base::gsub()] function.
 #' @param remove_space (logical; default TRUE): do we remove space?
+#' @param remove_NA (logical; default FALSE): do we remove NA (in majuscule)?
 #' @author Adrien Taudière
 #'
 #' @return A  \code{\link[phyloseq]{phyloseq-class}} object with simplified taxonomy
 #' @export
-simplify_taxo <- function(physeq, remove_space = TRUE) {
+#' @examples
+#' d_fm <- data_fungi_mini
+#' d_fm@tax_table[, "Species"] <- paste0(rep(
+#'   c("s__", "s:"),
+#'   ntaxa(d_fm) / 2
+#' ), d_fm@tax_table[, "Species"])
+#'
+#' # First column is the new vector of Species,
+#' # second column is the column before simplification
+#' cbind(
+#'   simplify_taxo(d_fm)@tax_table[, "Species"],
+#'   d_fm@tax_table[, "Species"]
+#' )
+#' cbind(
+#'   simplify_taxo(d_fm, remove_NA = TRUE)@tax_table[, "Species"],
+#'   d_fm@tax_table[, "Species"]
+#' )
+simplify_taxo <- function(
+    physeq,
+    pattern_to_remove = c(".__", ".*:"),
+    remove_space = TRUE,
+    remove_NA = FALSE) {
   taxo <- physeq@tax_table
-  taxo <- gsub(".__", "", taxo, perl = TRUE)
+  for (p in pattern_to_remove) {
+    taxo <- gsub(p, "", taxo)
+  }
+
   if (remove_space) {
     taxo <- gsub(" ", "", taxo)
     taxo <- gsub("\u00a0", "", taxo)
   }
+
+  if (remove_NA) {
+    taxo <- gsub("NA", "", taxo, ignore.case = FALSE)
+  }
+
   physeq@tax_table <- tax_table(taxo)
   return(physeq)
 }

R/plot_functions.R

@@ -2502,7 +2502,7 @@ biplot_pq <- function(physeq,
 #' @inheritParams clean_pq
 #' @param split_by (required if pairs is NULL) the name of the factor to make all combination
 #'   of couples of values
-#' @param pairs (required if pairs is NULL) the name of the factor in physeq@sam_data` slot
+#' @param pairs (required if split_by is NULL) the name of the factor in physeq@sam_data` slot
 #'   to make plot by pairs of samples. Each level must be present only two times.
 #'   Note that if you set pairs, you also must set fact arguments to passed on to [biplot_pq()].
 #' @param na_remove (logical, default TRUE) if TRUE remove all the samples
@@ -4856,14 +4856,14 @@ hill_curves_pq <- function(physeq,
 #' #      ggtitle("NMDS"))) +
 #' #  patchwork::plot_layout(guides = "collect")
 #'
-#' df_uwot <- umap_pq(data_fungi_mini, pkg = "uwot")
+#' # df_uwot <- umap_pq(data_fungi_mini, pkg = "uwot")
 #'
-#' (ggplot(df_umap, aes(x = x_umap, y = y_umap, col = Height)) +
-#'   geom_point(size = 2) +
-#'   ggtitle("umap::umap")) /
-#'   (ggplot(df_uwot, aes(x = x_umap, y = y_umap, col = Height)) +
-#'     geom_point(size = 2) +
-#'     ggtitle("uwot::umap2"))
+#' #   (ggplot(df_umap, aes(x = x_umap, y = y_umap, col = Height)) +
+#' #     geom_point(size = 2) +
+#' #    ggtitle("umap::umap")) /
+#' #     (ggplot(df_uwot, aes(x = x_umap, y = y_umap, col = Height)) +
+#' #       geom_point(size = 2) +
+#' #       ggtitle("uwot::umap2"))
 #'
 #' @details
 #' This function is mainly a wrapper of the work of others.