diff --git a/..Rcheck/00check.log b/..Rcheck/00check.log new file mode 100644 index 000000000..5b9f86eac --- /dev/null +++ b/..Rcheck/00check.log @@ -0,0 +1,13 @@ +* using log directory ‘/home/rpth/Projects/mia/..Rcheck’ +* using R version 4.5.3 (2026-03-11) +* using platform: x86_64-pc-linux-gnu +* R was compiled by + gcc (Ubuntu 11.4.0-1ubuntu1~22.04.3) 11.4.0 + GNU Fortran (Ubuntu 11.4.0-1ubuntu1~22.04.3) 11.4.0 +* running under: Pop!_OS 22.04 LTS +* using session charset: UTF-8 +* checking for file ‘./DESCRIPTION’ ... ERROR +Required fields missing or empty: + ‘Author’ ‘Maintainer’ +* DONE +Status: 1 ERROR diff --git a/DESCRIPTION b/DESCRIPTION index 9bb0d5ec7..09c6edaf1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -67,14 +67,17 @@ Imports: BiocParallel, Biostrings, bluster, + data.table, DECIPHER, decontam, DelayedArray, DelayedMatrixStats, DirichletMultinomial, dplyr, + httr2, IRanges, MASS, + Matrix, MatrixGenerics, methods, ecodive, diff --git a/NAMESPACE b/NAMESPACE index b0f050187..fd3d715e2 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -57,6 +57,7 @@ export(estimateDominance) export(estimateEvenness) export(estimateFaith) export(estimateRichness) +export(fetchMetalogTSE) export(full_join) export(getAbundanceClass) export(getAbundant) @@ -362,6 +363,8 @@ importFrom(IRanges,LogicalList) importFrom(IRanges,NumericList) importFrom(IRanges,relist) importFrom(MASS,isoMDS) +importFrom(Matrix,rowSums) +importFrom(Matrix,sparseMatrix) importFrom(MatrixGenerics,rowSums2) importFrom(MultiAssayExperiment,ExperimentList) importFrom(MultiAssayExperiment,MultiAssayExperiment) @@ -401,6 +404,7 @@ importFrom(SummarizedExperiment,assays) importFrom(SummarizedExperiment,colData) importFrom(SummarizedExperiment,rowData) importFrom(SummarizedExperiment,rowRanges) +importFrom(TreeSummarizedExperiment,TreeSummarizedExperiment) importFrom(TreeSummarizedExperiment,changeTree) importFrom(TreeSummarizedExperiment,subsetByLeaf) importFrom(ape,as.phylo) @@ -414,6 +418,11 @@ importFrom(ape,is.rooted) importFrom(ape,read.tree) importFrom(ape,reorder.phylo) importFrom(bluster,clusterRows) +importFrom(data.table,dcast) +importFrom(data.table,fread) +importFrom(data.table,setkey) +importFrom(data.table,setnames) +importFrom(data.table,tstrsplit) importFrom(decontam,isContaminant) importFrom(decontam,isNotContaminant) importFrom(dplyr,"%>%") @@ -432,6 +441,16 @@ importFrom(dplyr,sym) importFrom(dplyr,tally) importFrom(ecodive,unweighted_unifrac) importFrom(ecodive,weighted_unifrac) +importFrom(httr2,req_error) +importFrom(httr2,req_options) +importFrom(httr2,req_perform) +importFrom(httr2,req_progress) +importFrom(httr2,req_retry) +importFrom(httr2,req_timeout) +importFrom(httr2,req_user_agent) +importFrom(httr2,request) +importFrom(httr2,resp_header) +importFrom(httr2,resp_status) importFrom(rlang,":=") importFrom(rlang,sym) importFrom(scater,calculateMDS) diff --git a/R/fetchMetalogTSE.R b/R/fetchMetalogTSE.R new file mode 100644 index 000000000..a6ec20a6f --- /dev/null +++ b/R/fetchMetalogTSE.R @@ -0,0 +1,538 @@ +#' Fetch data from the Metalog database as +#' \code{TreeSummarizedExperiment} +#' +#' \code{fetchMetalogTSE} downloads MetaPhlAn4 taxonomic profiles and +#' associated sample metadata from the +#' \href{https://metalog.embl.de/}{Metalog} database and returns them as a +#' \code{TreeSummarizedExperiment} object. An optional sample list can be +#' provided to retrieve only a subset of samples (exported from the Metalog +#' web UI). +#' +#' @param collection \code{Character vector}. One or more Metalog collections +#' to download. Each value must be one of \code{"human"}, \code{"animal"}, +#' \code{"ocean"}, or \code{"environmental"}. When multiple collections are +#' given, their assays and metadata are merged before constructing the +#' final object; a \code{collection} column is added to \code{colData} +#' identifying the source collection of each sample. +#' +#' @param meta.type \code{Character scalar}. The metadata scope to download. +#' Must be one of \code{"core"}, \code{"extended"}, \code{"all"}, or +#' \code{"none"}. When \code{"none"}, no sample metadata is downloaded and +#' \code{colData} is empty (in the multi-collection case it retains only +#' the \code{collection} column identifying the source of each sample). +#' (Default: \code{"core"}). +#' +#' @param samplelist \code{Character scalar} or \code{NULL}. File path to a +#' sample list exported from the Metalog web UI. Supported formats are +#' \code{csv}, \code{tsv}, and \code{json}. When provided, the returned +#' object is filtered to include only the listed samples; taxa with zero +#' abundance across the retained samples are also removed. +#' (Default: \code{NULL}). +#' +#' @param use.cache \code{Logical scalar}. Should previously downloaded files +#' be reused? When \code{TRUE}, cached files in the download directory are +#' used if available. (Default: \code{TRUE}). +#' +#' @param make.dense \code{Logical scalar}. Should the assay be returned as a +#' dense base R matrix? Internally the assay is built as a sparse matrix to +#' conserve memory while loading and merging collections. When \code{TRUE} +#' (the default), it is converted to a dense matrix as a final step for +#' broader compatibility with downstream tools. Set to \code{FALSE} to +#' retain the sparse representation. (Default: \code{TRUE}). +#' +#' @details +#' Data is downloaded from the Metalog database +#' (\url{https://metalog.embl.de/}) and returned as a +#' \code{TreeSummarizedExperiment}. The assay is stored as a sparse matrix +#' named \code{"relabundance"} containing MetaPhlAn4 relative abundances at +#' SGB (species-level genome bin) resolution. Full taxonomic lineages are +#' mapped to standard ranks (Kingdom through SGB) and stored in +#' \code{rowData}. Sample metadata (in long format on the server) is pivoted +#' to wide format and stored in \code{colData}. +#' +#' This function requires an internet connection to download data from the +#' Metalog server. +#' +#' Provenance information is stored in \code{metadata(tse)$metalog} as a +#' list containing the source URL, license, collection and metadata type +#' used, the date stamps parsed from the downloaded file names +#' (\code{date_profile}, \code{date_metadata}), and the date the data was +#' fetched (\code{date_fetched}). +#' +#' @return A +#' \code{\link[TreeSummarizedExperiment:TreeSummarizedExperiment-class]{TreeSummarizedExperiment}} +#' object +#' +#' @name fetchMetalogTSE +#' @seealso +#' \code{\link[=importMetaPhlAn]{importMetaPhlAn}} +#' +#' @export +#' +#' @author Rasmus Hindström +#' +#' @references +#' Metalog database: \url{https://metalog.embl.de/} +#' +#' The data is made available under the Open Database License (ODbL) v1.0. +#' +#' @examples +#' \dontrun{ +#' # Fetch the human collection with core metadata +#' tse <- fetchMetalogTSE("human") +#' +#' # Fetch with extended metadata +#' tse <- fetchMetalogTSE("human", meta.type = "extended") +#' +#' # Fetch a subset of samples using a sample list +#' tse <- fetchMetalogTSE("human", samplelist = "my_samples.csv") +#' +#' # Fetch and merge multiple collections +#' tse <- fetchMetalogTSE(c("human", "animal")) +#' } +#' +NULL + +#' @rdname fetchMetalogTSE +#' @importFrom TreeSummarizedExperiment TreeSummarizedExperiment +#' @importFrom S4Vectors SimpleList DataFrame metadata metadata<- +#' @importFrom data.table fread setnames setkey dcast tstrsplit +#' @importFrom ape read.tree keep.tip +#' @export +fetchMetalogTSE <- function( + collection, + meta.type = "core", + samplelist = NULL, + use.cache = TRUE, + make.dense = TRUE) { + ################################ Input check ############################### + allowed_collections <- c("human", "animal", "ocean", "environmental") + if (!is.character(collection) || length(collection) < 1 || + anyNA(collection) || !all(nzchar(collection)) || + !all(collection %in% allowed_collections)) { + stop( + "'collection' must be a character vector with values from: ", + paste(dQuote(allowed_collections), collapse = ", "), + call. = FALSE + ) + } + collection <- unique(collection) + allowed_meta_types <- c("core", "extended", "all", "none") + if (!.is_non_empty_string(meta.type) || + !meta.type %in% allowed_meta_types) { + stop( + "'meta.type' must be one of: ", + paste(dQuote(allowed_meta_types), collapse = ", "), + call. = FALSE + ) + } + if (!is.null(samplelist) && !.is_non_empty_string(samplelist)) { + stop("'samplelist' must be a single character value or NULL.", + call. = FALSE) + } + if (!is.null(samplelist)) { + if (!file.exists(samplelist)) { + stop("'samplelist' file does not exist: ", samplelist, + call. = FALSE) + } + ext <- tolower(tools::file_ext(samplelist)) + allowed_exts <- c("csv", "tsv", "json") + if (!ext %in% allowed_exts) { + stop( + "'samplelist' file type must be one of: ", + paste(allowed_exts, collapse = ", "), + call. = FALSE + ) + } + } + if (!.is_a_bool(use.cache)) { + stop("'use.cache' must be TRUE or FALSE.", call. = FALSE) + } + if (!.is_a_bool(make.dense)) { + stop("'make.dense' must be TRUE or FALSE.", call. = FALSE) + } + ############################## Input check end ############################# + # Latest database file for taxonomy mapping (shared across collections) + mapping_db <- .download_metalog_file( + "https://metalog.embl.de/static/download/profiles/metaphlan4_clades.tsv.gz", + use.cache = use.cache + ) + # Per-collection: download, load assay, optional sample filter, load md + data_files_list <- lapply(collection, .resolve_metalog_url, + meta.type = meta.type, use.cache = use.cache) + names(data_files_list) <- collection + per_coll <- lapply(collection, function(co) { + df <- data_files_list[[co]] + al <- .load_metalog_assay(df[["assay"]]) + if (!is.null(samplelist)) { + message("Filtering to requested samples in '", co, "'...") + al <- .filter_metalog_samples(al, samplelist) + } + md <- if (identical(meta.type, "none")) { + data.frame(row.names = al[["samples"]]) + } else { + .load_metalog_metadata(df[["md"]], al[["samples"]]) + } + list(assay_list = al, md = md) + }) + names(per_coll) <- collection + # Merge assays and metadata across collections + merged <- .merge_metalog_assays(lapply(per_coll, `[[`, "assay_list")) + if (make.dense) { + merged[["assay"]] <- as.matrix(merged[["assay"]]) + } + md_df <- .merge_metalog_metadata( + lapply(per_coll, `[[`, "md"), collection, merged[["samples"]], + add.collection = !(identical(meta.type, "none") && + length(collection) == 1)) + # Map SGBs to full taxonomic lineage + tax <- .construct_metalog_taxmap(mapping_db, merged[["taxa"]]) + # Download and prune the MetaPhlAn4 SGB phylogeny to taxa in the data + tree_info <- .construct_metalog_tree(merged[["taxa"]], use.cache) + + tse <- TreeSummarizedExperiment( + assays = SimpleList("relabundance" = merged[["assay"]]), + colData = DataFrame(md_df), + rowData = DataFrame(tax), + rowTree = tree_info[["tree"]], + rowNodeLab = tree_info[["node_lab"]] + ) + + # Store provenance information + date_profile <- vapply(data_files_list, + function(df) .parse_metalog_date(df[["assay"]]), character(1)) + date_metadata <- vapply(data_files_list, function(df) { + if (is.na(df[["md"]])) NA_character_ + else .parse_metalog_date(df[["md"]]) + }, character(1)) + metadata(tse)$metalog <- list( + source = "https://metalog.embl.de/", + license = "Open Database License (ODbL) v1.0", + collection = collection, + meta.type = meta.type, + date_profile = date_profile, + date_metadata = date_metadata, + date_fetched = Sys.Date() + ) + return(tse) +} + +# Merge per-collection assay matrices into a single matrix. +# Rows (taxa) are the union; columns (samples) are concatenated. Sample +# alias collisions across collections trigger a hard error. +#' @importFrom Matrix sparseMatrix +.merge_metalog_assays <- function(assay_lists) { + if (length(assay_lists) == 1) { + al <- assay_lists[[1]] + return(list( + assay = al[["assay"]], + taxa = al[["taxa"]], + samples = al[["samples"]] + )) + } + taxa <- sort(unique(unlist( + lapply(assay_lists, `[[`, "taxa"), use.names = FALSE))) + samples <- unlist( + lapply(assay_lists, `[[`, "samples"), use.names = FALSE) + dups <- unique(samples[duplicated(samples)]) + if (length(dups) > 0) { + stop( + "Duplicate sample aliases across collections: ", + paste(dQuote(utils::head(dups, 5)), collapse = ", "), + if (length(dups) > 5) ", ..." else "", + call. = FALSE + ) + } + # Collect (i, j, x) triplets from each per-collection sparse matrix and + # remap them into the merged taxa/sample index space. + sample_offset <- 0L + triplets <- lapply(assay_lists, function(al) { + m <- methods::as(al[["assay"]], "TsparseMatrix") + i <- match(rownames(m)[m@i + 1L], taxa) + j <- m@j + 1L + sample_offset + sample_offset <<- sample_offset + ncol(m) + list(i = i, j = j, x = m@x) + }) + X <- Matrix::sparseMatrix( + i = unlist(lapply(triplets, `[[`, "i"), use.names = FALSE), + j = unlist(lapply(triplets, `[[`, "j"), use.names = FALSE), + x = unlist(lapply(triplets, `[[`, "x"), use.names = FALSE), + dims = c(length(taxa), length(samples)), + dimnames = list(taxa, samples) + ) + list(assay = X, taxa = taxa, samples = samples) +} + +# Merge per-collection metadata data.frames. Adds a 'collection' column +# identifying the source collection of each sample. Non-overlapping +# columns are filled with NA. +.merge_metalog_metadata <- function(md_list, collections, samples, + add.collection = TRUE) { + if (add.collection) { + for (i in seq_along(md_list)) { + md_list[[i]][["collection"]] <- collections[[i]] + } + } + combined <- data.table::rbindlist( + lapply(md_list, function(x) { + x[["sample_alias"]] <- rownames(x) + x + }), + fill = TRUE, use.names = TRUE + ) + df <- as.data.frame(combined) + rownames(df) <- df[["sample_alias"]] + df[["sample_alias"]] <- NULL + # Reorder to match merged sample order + df <- df[samples, , drop = FALSE] + df +} + +################################ HELP FUNCTIONS ################################ + +# Extract the YYYY-MM-DD date stamp from a Metalog filename. +# Returns NA_character_ if no date is found. +.parse_metalog_date <- function(path) { + fname <- basename(path) + m <- regmatches(fname, regexpr("[0-9]{4}-[0-9]{2}-[0-9]{2}", fname)) + if (length(m) == 0) NA_character_ else m +} + +# Download a file from Metalog, handling HTTP -> HTTPS redirect issues. +# Adapted from Metalog's example script. +#' @importFrom httr2 request req_user_agent req_options req_error +#' req_timeout req_retry req_progress req_perform resp_status resp_header +.download_metalog_file <- function( + target_url, + download_dir = tools::R_user_dir("mia", "cache"), + use.cache = TRUE) { + base_filename <- basename(target_url) + # Check cache + if (use.cache) { + pattern <- sub( + "latest", "[0-9]{4}-[0-9]{2}-[0-9]{2}", base_filename) + matching_files <- list.files( + download_dir, pattern = pattern, full.names = TRUE) + if (length(matching_files) > 0) { + latest_file <- max(matching_files) + message("Loaded cached file: ", latest_file) + return(latest_file) + } + } + if (!dir.exists(download_dir)) { + dir.create(download_dir, recursive = TRUE, showWarnings = FALSE) + } + message("Fetching file from: ", target_url) + pkg_version <- utils::packageDescription("mia", fields = "Version") + ua <- paste0("mia/", pkg_version, " (R; Bioconductor)") + # Metalog server may downgrade HTTPS to HTTP on redirect. Intercept + # the redirect and force HTTPS. + initial_resp <- httr2::request(target_url) |> + httr2::req_user_agent(ua) |> + httr2::req_options(followlocation = FALSE) |> + httr2::req_error(is_error = ~ FALSE) |> + httr2::req_timeout(300) |> + httr2::req_perform() + status <- httr2::resp_status(initial_resp) + if (status >= 300 && status < 400) { + final_url <- httr2::resp_header(initial_resp, "location") + if (is.null(final_url)) { + stop("Redirect response missing Location header.", + call. = FALSE) + } + final_url <- sub("^http://", "https://", final_url) + message("Intercepted redirect. Forcing HTTPS: ", final_url) + } else if (status == 200) { + final_url <- target_url + } else { + stop("Initial request failed with status: ", status, + call. = FALSE) + } + # Download the file with retry and timeout + filename <- basename(final_url) + destfile <- file.path(download_dir, filename) + message("Downloading to: ", destfile) + tryCatch({ + httr2::request(final_url) |> + httr2::req_user_agent(ua) |> + httr2::req_timeout(300) |> + httr2::req_retry(max_tries = 3) |> + httr2::req_progress() |> + httr2::req_perform(path = destfile) + }, error = function(e) { + if (file.exists(destfile)) file.remove(destfile) + stop("Error downloading the file: ", conditionMessage(e), + call. = FALSE) + }) + return(destfile) +} + +# Construct Metalog download URLs and fetch assay + metadata files +.resolve_metalog_url <- function(collection, meta.type, use.cache) { + cache_dir <- tools::R_user_dir("mia", "cache") + base_url <- "https://metalog.embl.de/static/download" + assay_url <- sprintf( + "%s/profiles/%s_metaphlan4_latest.tsv.gz", base_url, collection) + md_url <- sprintf( + "%s/metadata/%s_%s_long_latest.tsv.gz", + base_url, collection, meta.type) + assay_file <- .download_metalog_file( + target_url = assay_url, + download_dir = cache_dir, + use.cache = use.cache + ) + md_file <- if (identical(meta.type, "none")) { + NA_character_ + } else { + .download_metalog_file( + target_url = md_url, + download_dir = cache_dir, + use.cache = use.cache + ) + } + list(assay = assay_file, md = md_file) +} + +# Load MetaPhlAn4 profiles into a sparse matrix (rows = taxa, cols = samples) +#' @importFrom Matrix sparseMatrix +.load_metalog_assay <- function(path, sep = "\t") { + # data.table NSE bindings + clade_name <- rel_abund <- sample_alias <- NULL + dt <- data.table::fread(path, sep = sep) + data.table::setnames( + dt, seq_len(3), c("sample_alias", "clade_name", "rel_abund")) + dt <- dt[startsWith(clade_name, "t__SGB"), ] + dt[, rel_abund := as.numeric(rel_abund)] + dt <- dt[!is.na(rel_abund) & rel_abund != 0] + # Aggregate duplicates + data.table::setkey(dt, clade_name, sample_alias) + dt <- dt[, .(rel_abund = sum(rel_abund)), + by = .(clade_name, sample_alias)] + taxa <- sort(unique(dt$clade_name)) + samples <- sort(unique(dt$sample_alias)) + X <- Matrix::sparseMatrix( + i = match(dt$clade_name, taxa), + j = match(dt$sample_alias, samples), + x = dt$rel_abund, + dims = c(length(taxa), length(samples)), + dimnames = list(taxa, samples) + ) + list(assay = X, taxa = taxa, samples = samples) +} + +# Load Metalog metadata, pivot to wide, and subset to samples present in assay +.load_metalog_metadata <- function(meta_path, samples, sep = "\t") { + sample_alias <- metadata_item <- NULL + dt <- data.table::fread(meta_path, sep = sep, na.strings = c("", "NA")) + wide <- data.table::dcast( + dt, + sample_alias ~ metadata_item, + value.var = "value", + fill = NA_character_ + ) + missing <- setdiff(samples, wide$sample_alias) + if (length(missing) > 0) { + warning( + length(missing), " sample(s) present in assay data but missing ", + "from metadata.", call. = FALSE + ) + } + # Subset and reorder via data.table key lookup + data.table::setkey(wide, sample_alias) + wide <- wide[.(samples)] + meta_df <- as.data.frame(wide) + rownames(meta_df) <- meta_df$sample_alias + meta_df$sample_alias <- NULL + meta_df +} + +# Map SGB clade names to full taxonomic lineages +.construct_metalog_taxmap <- function(database, taxa) { + clade_name <- lineage <- NULL + taxmap <- data.table::fread(database, sep = "\t", header = TRUE) + taxmap <- taxmap[startsWith(clade_name, "t__SGB")] + taxmap <- taxmap[!duplicated(clade_name)] + # Align to taxa order; unmatched SGBs get NA + idx <- match(taxa, taxmap$clade_name) + n_missing <- sum(is.na(idx)) + if (n_missing > 0) { + warning( + n_missing, " of ", length(taxa), + " taxa could not be mapped to a full lineage.", + call. = FALSE + ) + } + taxmap <- taxmap[idx] + # Parse lineage into standard taxonomy ranks + lineage_cols <- c( + "Kingdom", "Phylum", "Class", "Order", + "Family", "Genus", "Species", "SGB" + ) + taxmap[, (lineage_cols) := data.table::tstrsplit(lineage, "\\|")] + result <- as.data.frame(taxmap[, lineage_cols, with = FALSE]) + rownames(result) <- taxa + result +} + +# Download the MetaPhlAn4 SGB phylogeny and prune to taxa in the dataset. +# Returns a list with the pruned tree and a per-taxon tip label vector +# (NA for taxa absent from the tree) suitable for rowNodeLab. +.construct_metalog_tree <- function(taxa, use.cache) { + tree_url <- paste0( + "http://cmprod1.cibio.unitn.it/biobakery4/metaphlan_databases/", + "mpa_vJun23_CHOCOPhlAnSGB_202307.nwk" + ) + tree_file <- .download_metalog_file(tree_url, use.cache = use.cache) + tree <- ape::read.tree(tree_file) + # Tree tips are bare numeric SGB ids (e.g. "122987"); taxa are + # MetaPhlAn clade names like "t__SGB1234" or "t__SGB1234_group". + # Extract the numeric id from each taxon to match tips. + taxa_id <- vapply(taxa, function(x) { + m <- regmatches(x, regexpr("SGB[0-9]+", x)) + if (length(m) == 0) NA_character_ else sub("^SGB", "", m) + }, character(1)) + node_lab <- ifelse(taxa_id %in% tree$tip.label, taxa_id, NA_character_) + keep <- node_lab[!is.na(node_lab)] + n_missing <- sum(is.na(node_lab)) + if (n_missing > 0) { + warning( + n_missing, " of ", length(taxa), + " taxa could not be matched to a tip in the SGB tree.", + call. = FALSE + ) + } + if (length(keep) == 0) { + stop("No taxa matched any tip in the SGB tree.", call. = FALSE) + } + tree <- ape::keep.tip(tree, keep) + list(tree = tree, node_lab = node_lab) +} + +# Filter assay data to samples listed in a sample list file +.filter_metalog_samples <- function(assay_list, samplelist) { + ext <- tolower(tools::file_ext(samplelist)) + if (ext %in% c("csv", "tsv")) { + sl_df <- data.table::fread(samplelist) + } else if (ext == "json") { + .require_package("jsonlite") + sl_df <- as.data.frame(jsonlite::fromJSON(samplelist)) + } + target_samples <- sl_df[["sample_alias"]] + available_samples <- assay_list[["samples"]] + keep_samples <- intersect(target_samples, available_samples) + if (length(keep_samples) == 0) { + stop( + "None of the samples in 'samplelist' were found in the dataset.", + call. = FALSE + ) + } + # Subset columns (samples) + assay_list$assay <- assay_list$assay[, keep_samples, drop = FALSE] + assay_list$samples <- keep_samples + # Drop taxa with zero abundance + row_sums <- rowSums(assay_list$assay) + keep_taxa <- names(row_sums[row_sums > 0]) + assay_list$assay <- assay_list$assay[keep_taxa, , drop = FALSE] + assay_list$taxa <- keep_taxa + assay_list +} diff --git a/man/fetchMetalogTSE.Rd b/man/fetchMetalogTSE.Rd new file mode 100644 index 000000000..aeffc76b9 --- /dev/null +++ b/man/fetchMetalogTSE.Rd @@ -0,0 +1,90 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/fetchMetalogTSE.R +\name{fetchMetalogTSE} +\alias{fetchMetalogTSE} +\title{Fetch data from the Metalog database as +\code{TreeSummarizedExperiment}} +\usage{ +fetchMetalogTSE( + collection, + meta.type = "core", + samplelist = NULL, + use.cache = TRUE +) +} +\arguments{ +\item{collection}{\code{Character scalar}. The Metalog collection to +download. Must be one of \code{"human"}, \code{"animal"}, \code{"ocean"}, +or \code{"environmental"}.} + +\item{meta.type}{\code{Character scalar}. The metadata scope to download. +Must be one of \code{"core"}, \code{"extended"}, or \code{"all"}. +(Default: \code{"core"}).} + +\item{samplelist}{\code{Character scalar} or \code{NULL}. File path to a +sample list exported from the Metalog web UI. Supported formats are +\code{csv}, \code{tsv}, and \code{json}. When provided, the returned +object is filtered to include only the listed samples; taxa with zero +abundance across the retained samples are also removed. +(Default: \code{NULL}).} + +\item{use.cache}{\code{Logical scalar}. Should previously downloaded files +be reused? When \code{TRUE}, cached files in the download directory are +used if available. (Default: \code{TRUE}).} +} +\value{ +A +\code{\link[TreeSummarizedExperiment:TreeSummarizedExperiment-class]{TreeSummarizedExperiment}} +object +} +\description{ +\code{fetchMetalogTSE} downloads MetaPhlAn4 taxonomic profiles and +associated sample metadata from the +\href{https://metalog.embl.de/}{Metalog} database and returns them as a +\code{TreeSummarizedExperiment} object. An optional sample list can be +provided to retrieve only a subset of samples (exported from the Metalog +web UI). +} +\details{ +Data is downloaded from the Metalog database +(\url{https://metalog.embl.de/}) and returned as a +\code{TreeSummarizedExperiment}. The assay is stored as a sparse matrix +named \code{"relabundance"} containing MetaPhlAn4 relative abundances at +SGB (species-level genome bin) resolution. Full taxonomic lineages are +mapped to standard ranks (Kingdom through SGB) and stored in +\code{rowData}. Sample metadata (in long format on the server) is pivoted +to wide format and stored in \code{colData}. + +This function requires an internet connection to download data from the +Metalog server. + +Provenance information is stored in \code{metadata(tse)$metalog} as a +list containing the source URL, license, collection and metadata type +used, the date stamps parsed from the downloaded file names +(\code{date_profile}, \code{date_metadata}), and the date the data was +fetched (\code{date_fetched}). +} +\examples{ +\dontrun{ +# Fetch the human collection with core metadata +tse <- fetchMetalogTSE("human") + +# Fetch with extended metadata +tse <- fetchMetalogTSE("human", meta.type = "extended") + +# Fetch a subset of samples using a sample list +tse <- fetchMetalogTSE("human", samplelist = "my_samples.csv") +} + +} +\references{ +Metalog database: \url{https://metalog.embl.de/} + +The data is made available under the Open Database License (ODbL) v1.0. +} +\seealso{ +\code{\link[=importMetaPhlAn]{importMetaPhlAn}} +} +\author{ +Rasmus Hindström +} diff --git a/tests/testthat/test-fetchMetalogTSE.R b/tests/testthat/test-fetchMetalogTSE.R new file mode 100644 index 000000000..369be9e91 --- /dev/null +++ b/tests/testthat/test-fetchMetalogTSE.R @@ -0,0 +1,307 @@ +################################################################################ +# Helper: create small fixture files for the data-processing helpers +################################################################################ + +# Minimal long-format MetaPhlAn4 profile (3 cols: sample, clade, abundance) +.make_assay_fixture <- function(dir) { + path <- file.path(dir, "assay.tsv") + lines <- c( + "sample_alias\tclade_name\trel_abund", + "S1\tk__Bacteria\t100.0", + "S1\tt__SGB1234\t60.5", + "S1\tt__SGB5678\t39.5", + "S2\tt__SGB1234\t80.0", + "S2\tt__SGB5678\t20.0", + "S3\tt__SGB1234\t50.0", + "S3\tt__SGB9999\t50.0" + ) + writeLines(lines, path) + path +} + +# Minimal long-format metadata +.make_metadata_fixture <- function(dir) { + path <- file.path(dir, "metadata.tsv") + lines <- c( + "sample_alias\tmetadata_item\tvalue", + "S1\tage\t30", + "S1\tcountry\tFI", + "S2\tage\t45", + "S2\tcountry\tDE", + "S3\tage\t25", + "S3\tcountry\tSE" + ) + writeLines(lines, path) + path +} + +# Minimal taxonomy mapping database +.make_taxdb_fixture <- function(dir) { + path <- file.path(dir, "taxdb.tsv") + lines <- c( + paste("clade_name", "NCBI_taxids", "lineage", sep = "\t"), + paste("t__SGB1234", "12345", + "k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Lactobacillaceae|g__Lactobacillus|s__L_acidophilus|t__SGB1234", + sep = "\t"), + paste("t__SGB5678", "56789", + "k__Bacteria|p__Firmicutes|c__Bacilli|o__Lactobacillales|f__Streptococcaceae|g__Streptococcus|s__S_thermophilus|t__SGB5678", + sep = "\t"), + paste("t__SGB9999", "99999", + "k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Enterobacterales|f__Enterobacteriaceae|g__Escherichia|s__E_coli|t__SGB9999", + sep = "\t"), + paste("k__Bacteria", "2", "k__Bacteria", sep = "\t") + ) + writeLines(lines, path) + path +} + +# Minimal sample list (csv) +.make_samplelist_fixture <- function(dir, samples = c("S1", "S2")) { + path <- file.path(dir, "samplelist.csv") + df <- data.frame(sample_alias = samples) + write.csv(df, path, row.names = FALSE) + path +} + +################################################################################ +# Input validation tests +################################################################################ + +test_that("fetchMetalogTSE rejects invalid collection", { + expect_error(fetchMetalogTSE("invalid_collection"), + "'collection' must be one of") + expect_error(fetchMetalogTSE(123), + "'collection' must be one of") + expect_error(fetchMetalogTSE(c("human", "animal")), + "'collection' must be one of") +}) + +test_that("fetchMetalogTSE rejects invalid meta.type", { + expect_error(fetchMetalogTSE("human", meta.type = "bad"), + "'meta.type' must be one of") + expect_error(fetchMetalogTSE("human", meta.type = 42), + "'meta.type' must be one of") +}) + +test_that("fetchMetalogTSE rejects invalid use.cache", { + expect_error(fetchMetalogTSE("human", use.cache = "yes"), + "'use.cache' must be TRUE or FALSE") + expect_error(fetchMetalogTSE("human", use.cache = NA), + "'use.cache' must be TRUE or FALSE") +}) + +test_that("fetchMetalogTSE rejects non-string samplelist", { + expect_error(fetchMetalogTSE("human", samplelist = 123), + "'samplelist' must be a single character value or NULL") +}) + +test_that("fetchMetalogTSE rejects non-existent samplelist", { + expect_error( + fetchMetalogTSE("human", samplelist = "no_such_file.csv"), + "'samplelist' file does not exist") +}) + +test_that("fetchMetalogTSE rejects unsupported samplelist extension", { + tmp <- tempfile(fileext = ".xlsx") + writeLines("placeholder", tmp) + on.exit(unlink(tmp)) + expect_error(fetchMetalogTSE("human", samplelist = tmp), + "'samplelist' file type must be one of") +}) + +################################################################################ +# .load_metalog_assay +################################################################################ + +test_that(".load_metalog_assay returns correct structure", { + dir <- tempdir() + path <- .make_assay_fixture(dir) + on.exit(unlink(path)) + result <- mia:::.load_metalog_assay(path) + expect_type(result, "list") + expect_named(result, c("assay", "taxa", "samples")) + expect_s4_class(result$assay, "dgCMatrix") +}) + +test_that(".load_metalog_assay filters to SGB rows only", { + dir <- tempdir() + path <- .make_assay_fixture(dir) + on.exit(unlink(path)) + result <- mia:::.load_metalog_assay(path) + # k__Bacteria row should be excluded + expect_true(all(startsWith(result$taxa, "t__SGB"))) + expect_equal(length(result$taxa), 3L) +}) + +test_that(".load_metalog_assay has correct dimensions", { + dir <- tempdir() + path <- .make_assay_fixture(dir) + on.exit(unlink(path)) + result <- mia:::.load_metalog_assay(path) + # 3 taxa (SGB1234, SGB5678, SGB9999), 3 samples (S1, S2, S3) + expect_equal(nrow(result$assay), 3L) + expect_equal(ncol(result$assay), 3L) + expect_equal(length(result$samples), 3L) +}) + +test_that(".load_metalog_assay aggregates duplicate entries", { + dir <- tempdir() + path <- file.path(dir, "assay_dup.tsv") + lines <- c( + "sample_alias\tclade_name\trel_abund", + "S1\tt__SGB1234\t30.0", + "S1\tt__SGB1234\t20.0" + ) + writeLines(lines, path) + on.exit(unlink(path)) + result <- mia:::.load_metalog_assay(path) + # Should sum to 50.0 + expect_equal(as.numeric(result$assay["t__SGB1234", "S1"]), 50.0) +}) + +################################################################################ +# .load_metalog_metadata +################################################################################ + +test_that(".load_metalog_metadata returns wide data.frame", { + dir <- tempdir() + path <- .make_metadata_fixture(dir) + on.exit(unlink(path)) + result <- mia:::.load_metalog_metadata(path, c("S1", "S2")) + expect_s3_class(result, "data.frame") + expect_true("age" %in% colnames(result)) + expect_true("country" %in% colnames(result)) + expect_equal(nrow(result), 2L) +}) + +test_that(".load_metalog_metadata subsets to requested samples", { + dir <- tempdir() + path <- .make_metadata_fixture(dir) + on.exit(unlink(path)) + result <- mia:::.load_metalog_metadata(path, c("S2")) + expect_equal(nrow(result), 1L) + expect_equal(rownames(result), "S2") +}) + +test_that(".load_metalog_metadata warns on missing samples", { + dir <- tempdir() + path <- .make_metadata_fixture(dir) + on.exit(unlink(path)) + expect_warning( + mia:::.load_metalog_metadata(path, c("S1", "MISSING")), + "sample\\(s\\) present in assay data but missing" + ) +}) + +################################################################################ +# .construct_metalog_taxmap +################################################################################ + +test_that(".construct_metalog_taxmap returns taxonomy data.frame", { + dir <- tempdir() + db_path <- .make_taxdb_fixture(dir) + on.exit(unlink(db_path)) + taxa <- c("t__SGB1234", "t__SGB5678") + result <- mia:::.construct_metalog_taxmap(db_path, taxa) + expect_s3_class(result, "data.frame") + expect_equal(nrow(result), 2L) + expect_equal(rownames(result), taxa) + expect_equal( + colnames(result), + c("Kingdom", "Phylum", "Class", "Order", + "Family", "Genus", "Species", "SGB") + ) +}) + +test_that(".construct_metalog_taxmap parses lineage correctly", { + dir <- tempdir() + db_path <- .make_taxdb_fixture(dir) + on.exit(unlink(db_path)) + result <- mia:::.construct_metalog_taxmap(db_path, "t__SGB1234") + expect_equal(result["t__SGB1234", "Kingdom"], "k__Bacteria") + expect_equal(result["t__SGB1234", "Phylum"], "p__Firmicutes") + expect_equal(result["t__SGB1234", "Genus"], "g__Lactobacillus") +}) + +test_that(".construct_metalog_taxmap preserves taxa order", { + dir <- tempdir() + db_path <- .make_taxdb_fixture(dir) + on.exit(unlink(db_path)) + taxa <- c("t__SGB5678", "t__SGB1234") + result <- mia:::.construct_metalog_taxmap(db_path, taxa) + expect_equal(rownames(result), taxa) +}) + +test_that(".construct_metalog_taxmap warns on unmatched taxa", { + dir <- tempdir() + db_path <- .make_taxdb_fixture(dir) + on.exit(unlink(db_path)) + taxa <- c("t__SGB1234", "t__SGB0000") + expect_warning( + mia:::.construct_metalog_taxmap(db_path, taxa), + "1 of 2 taxa could not be mapped" + ) +}) + +################################################################################ +# .parse_metalog_date +################################################################################ + +test_that(".parse_metalog_date extracts date from filename", { + expect_equal( + mia:::.parse_metalog_date( + "/cache/human_metaphlan4_2025-03-15.tsv.gz"), + "2025-03-15" + ) + expect_equal( + mia:::.parse_metalog_date( + "/cache/human_core_long_2024-12-01.tsv.gz"), + "2024-12-01" + ) +}) + +test_that(".parse_metalog_date returns NA when no date present", { + expect_true(is.na( + mia:::.parse_metalog_date("some_file_without_date.tsv.gz") + )) +}) + +################################################################################ +# .filter_metalog_samples +################################################################################ + +test_that(".filter_metalog_samples subsets to requested samples", { + dir <- tempdir() + assay_path <- .make_assay_fixture(dir) + sl_path <- .make_samplelist_fixture(dir, samples = c("S1", "S2")) + on.exit(unlink(c(assay_path, sl_path))) + assay_list <- mia:::.load_metalog_assay(assay_path) + result <- mia:::.filter_metalog_samples(assay_list, sl_path) + expect_equal(sort(result$samples), c("S1", "S2")) + expect_equal(ncol(result$assay), 2L) +}) + +test_that(".filter_metalog_samples drops zero-abundance taxa", { + dir <- tempdir() + assay_path <- .make_assay_fixture(dir) + # S1 and S2 have SGB1234 and SGB5678 but NOT SGB9999 + sl_path <- .make_samplelist_fixture(dir, samples = c("S1", "S2")) + on.exit(unlink(c(assay_path, sl_path))) + assay_list <- mia:::.load_metalog_assay(assay_path) + result <- mia:::.filter_metalog_samples(assay_list, sl_path) + # SGB9999 only in S3, so it should be dropped + expect_false("t__SGB9999" %in% result$taxa) + expect_equal(length(result$taxa), 2L) +}) + +test_that(".filter_metalog_samples errors when no samples match", { + dir <- tempdir() + assay_path <- .make_assay_fixture(dir) + sl_path <- .make_samplelist_fixture(dir, samples = c("NONEXISTENT")) + on.exit(unlink(c(assay_path, sl_path))) + assay_list <- mia:::.load_metalog_assay(assay_path) + expect_error( + mia:::.filter_metalog_samples(assay_list, sl_path), + "None of the samples" + ) +})