From dfc1a1083eee5dc953bf8293a676f323b3e150ac Mon Sep 17 00:00:00 2001 From: raivo-otus Date: Thu, 2 Apr 2026 14:09:58 +0300 Subject: [PATCH 1/4] refactor to fit mia style Refactored from previous "draft implementation" with assistance from Claude Opus 4.6 to fit mia style. In the process updated to use httr2 for downloads and some optimizations to logic. --- DESCRIPTION | 3 + NAMESPACE | 19 ++ R/fetchMetalogTSE.R | 380 ++++++++++++++++++++++++++ man/fetchMetalogTSE.Rd | 90 ++++++ tests/testthat/test-fetchMetalogTSE.R | 307 +++++++++++++++++++++ 5 files changed, 799 insertions(+) create mode 100644 R/fetchMetalogTSE.R create mode 100644 man/fetchMetalogTSE.Rd create mode 100644 tests/testthat/test-fetchMetalogTSE.R diff --git a/DESCRIPTION b/DESCRIPTION index a6477e256..75e1389f8 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 78b47e87a..95e3ac9e8 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) @@ -361,6 +362,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) @@ -400,6 +403,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) @@ -413,6 +417,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,"%>%") @@ -431,6 +440,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..33fcfa8d5 --- /dev/null +++ b/R/fetchMetalogTSE.R @@ -0,0 +1,380 @@ +#' 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 scalar}. The Metalog collection to +#' download. Must be one of \code{"human"}, \code{"animal"}, \code{"ocean"}, +#' or \code{"environmental"}. +#' +#' @param meta.type \code{Character scalar}. The metadata scope to download. +#' Must be one of \code{"core"}, \code{"extended"}, or \code{"all"}. +#' (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}). +#' +#' @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") +#' } +#' +NULL + +#' @rdname fetchMetalogTSE +#' @importFrom Matrix sparseMatrix rowSums +#' @importFrom TreeSummarizedExperiment TreeSummarizedExperiment +#' @importFrom S4Vectors SimpleList DataFrame metadata metadata<- +#' @importFrom data.table fread setnames setkey dcast tstrsplit +#' @export +fetchMetalogTSE <- function( + collection, + meta.type = "core", + samplelist = NULL, + use.cache = TRUE) { + ################################ Input check ############################### + allowed_collections <- c("human", "animal", "ocean", "environmental") + if (!.is_non_empty_string(collection) || + !collection %in% allowed_collections) { + stop( + "'collection' must be one of: ", + paste(dQuote(allowed_collections), collapse = ", "), + call. = FALSE + ) + } + allowed_meta_types <- c("core", "extended", "all") + 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) + } + ############################## Input check end ############################# + # Construct download URLs, download and cache + data_files <- .resolve_metalog_url(collection, meta.type, use.cache) + # Latest database file for taxonomy mapping + mapping_db <- .download_metalog_file( + "https://metalog.embl.de/static/download/profiles/metaphlan4_clades.tsv.gz", + use.cache = use.cache + ) + # Load assay as sparse matrix + assay_list <- .load_metalog_assay(data_files[["assay"]]) + # Optional filtering to requested samples + if (!is.null(samplelist)) { + message("Filtering to requested samples...") + assay_list <- .filter_metalog_samples(assay_list, samplelist) + } + # Load metadata (pivoted to wide format) + md_df <- .load_metalog_metadata( + data_files[["md"]], assay_list[["samples"]]) + # Map SGBs to full taxonomic lineage + tax <- .construct_metalog_taxmap(mapping_db, assay_list[["taxa"]]) + + tse <- TreeSummarizedExperiment( + assays = SimpleList("relabundance" = assay_list[["assay"]]), + colData = DataFrame(md_df), + rowData = DataFrame(tax) + ) + + # Store provenance information + metadata(tse)$metalog <- list( + source = "https://metalog.embl.de/", + license = "Open Database License (ODbL) v1.0", + collection = collection, + meta.type = meta.type, + date_profile = .parse_metalog_date(data_files[["assay"]]), + date_metadata = .parse_metalog_date(data_files[["md"]]), + date_fetched = Sys.Date() + ) + return(tse) +} + +################################ 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 <- .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) +.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)) + i <- match(dt$clade_name, taxa) + j <- match(dt$sample_alias, samples) + X <- Matrix::sparseMatrix( + i = i, j = j, 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 +} + +# 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 <- Matrix::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" + ) +}) From c795a720d921bab86fc8bfb4b90360405bddadb5 Mon Sep 17 00:00:00 2001 From: raivo-otus Date: Tue, 7 Apr 2026 14:11:06 +0300 Subject: [PATCH 2/4] Add phylogeny tree --- ..Rcheck/00check.log | 13 +++++++++++++ R/fetchMetalogTSE.R | 41 ++++++++++++++++++++++++++++++++++++++++- 2 files changed, 53 insertions(+), 1 deletion(-) create mode 100644 ..Rcheck/00check.log 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/R/fetchMetalogTSE.R b/R/fetchMetalogTSE.R index 33fcfa8d5..cf64c3ea0 100644 --- a/R/fetchMetalogTSE.R +++ b/R/fetchMetalogTSE.R @@ -82,6 +82,7 @@ NULL #' @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, @@ -149,11 +150,15 @@ fetchMetalogTSE <- function( data_files[["md"]], assay_list[["samples"]]) # Map SGBs to full taxonomic lineage tax <- .construct_metalog_taxmap(mapping_db, assay_list[["taxa"]]) + # Download and prune the MetaPhlAn4 SGB phylogeny to taxa in the data + tree_info <- .construct_metalog_tree(assay_list[["taxa"]], use.cache) tse <- TreeSummarizedExperiment( assays = SimpleList("relabundance" = assay_list[["assay"]]), colData = DataFrame(md_df), - rowData = DataFrame(tax) + rowData = DataFrame(tax), + rowTree = tree_info[["tree"]], + rowNodeLab = tree_info[["node_lab"]] ) # Store provenance information @@ -350,6 +355,40 @@ fetchMetalogTSE <- function( 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)) From a672a2089fa1447bf205a1836812376550f67029 Mon Sep 17 00:00:00 2001 From: raivo-otus Date: Tue, 7 Apr 2026 15:07:34 +0300 Subject: [PATCH 3/4] use dense matrix assays for better support --- R/fetchMetalogTSE.R | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/R/fetchMetalogTSE.R b/R/fetchMetalogTSE.R index cf64c3ea0..221500600 100644 --- a/R/fetchMetalogTSE.R +++ b/R/fetchMetalogTSE.R @@ -78,7 +78,6 @@ NULL #' @rdname fetchMetalogTSE -#' @importFrom Matrix sparseMatrix rowSums #' @importFrom TreeSummarizedExperiment TreeSummarizedExperiment #' @importFrom S4Vectors SimpleList DataFrame metadata metadata<- #' @importFrom data.table fread setnames setkey dcast tstrsplit @@ -275,7 +274,7 @@ fetchMetalogTSE <- function( list(assay = assay_file, md = md_file) } -# Load MetaPhlAn4 profiles into a sparse matrix (rows = taxa, cols = samples) +# Load MetaPhlAn4 profiles into a dense matrix (rows = taxa, cols = samples) .load_metalog_assay <- function(path, sep = "\t") { # data.table NSE bindings clade_name <- rel_abund <- sample_alias <- NULL @@ -291,13 +290,15 @@ fetchMetalogTSE <- function( by = .(clade_name, sample_alias)] taxa <- sort(unique(dt$clade_name)) samples <- sort(unique(dt$sample_alias)) - i <- match(dt$clade_name, taxa) - j <- match(dt$sample_alias, samples) - X <- Matrix::sparseMatrix( - i = i, j = j, x = dt$rel_abund, - dims = c(length(taxa), length(samples)), + X <- matrix( + 0, nrow = length(taxa), ncol = length(samples), dimnames = list(taxa, samples) ) + idx <- cbind( + match(dt$clade_name, taxa), + match(dt$sample_alias, samples) + ) + X[idx] <- dt$rel_abund list(assay = X, taxa = taxa, samples = samples) } @@ -411,7 +412,7 @@ fetchMetalogTSE <- function( assay_list$assay <- assay_list$assay[, keep_samples, drop = FALSE] assay_list$samples <- keep_samples # Drop taxa with zero abundance - row_sums <- Matrix::rowSums(assay_list$assay) + 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 From 156e4e1ac424fea6ceee84827627ac7cb8ce2ab2 Mon Sep 17 00:00:00 2001 From: raivo-otus Date: Wed, 8 Apr 2026 15:13:34 +0300 Subject: [PATCH 4/4] feature to combine multiple as sparse and make dense later --- R/fetchMetalogTSE.R | 196 +++++++++++++++++++++++++++++++++++--------- 1 file changed, 157 insertions(+), 39 deletions(-) diff --git a/R/fetchMetalogTSE.R b/R/fetchMetalogTSE.R index 221500600..a6ec20a6f 100644 --- a/R/fetchMetalogTSE.R +++ b/R/fetchMetalogTSE.R @@ -8,12 +8,18 @@ #' provided to retrieve only a subset of samples (exported from the Metalog #' web UI). #' -#' @param collection \code{Character scalar}. The Metalog collection to -#' download. Must be one of \code{"human"}, \code{"animal"}, \code{"ocean"}, -#' or \code{"environmental"}. +#' @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"}, or \code{"all"}. +#' 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 @@ -27,6 +33,13 @@ #' 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 @@ -73,6 +86,9 @@ #' #' # 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 @@ -87,18 +103,21 @@ fetchMetalogTSE <- function( collection, meta.type = "core", samplelist = NULL, - use.cache = TRUE) { + use.cache = TRUE, + make.dense = TRUE) { ################################ Input check ############################### allowed_collections <- c("human", "animal", "ocean", "environmental") - if (!.is_non_empty_string(collection) || - !collection %in% allowed_collections) { + if (!is.character(collection) || length(collection) < 1 || + anyNA(collection) || !all(nzchar(collection)) || + !all(collection %in% allowed_collections)) { stop( - "'collection' must be one of: ", + "'collection' must be a character vector with values from: ", paste(dQuote(allowed_collections), collapse = ", "), call. = FALSE ) } - allowed_meta_types <- c("core", "extended", "all") + 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( @@ -129,31 +148,50 @@ fetchMetalogTSE <- function( 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 ############################# - # Construct download URLs, download and cache - data_files <- .resolve_metalog_url(collection, meta.type, use.cache) - # Latest database file for taxonomy mapping + # 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 ) - # Load assay as sparse matrix - assay_list <- .load_metalog_assay(data_files[["assay"]]) - # Optional filtering to requested samples - if (!is.null(samplelist)) { - message("Filtering to requested samples...") - assay_list <- .filter_metalog_samples(assay_list, samplelist) + # 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"]]) } - # Load metadata (pivoted to wide format) - md_df <- .load_metalog_metadata( - data_files[["md"]], assay_list[["samples"]]) + 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, assay_list[["taxa"]]) + 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(assay_list[["taxa"]], use.cache) + tree_info <- .construct_metalog_tree(merged[["taxa"]], use.cache) tse <- TreeSummarizedExperiment( - assays = SimpleList("relabundance" = assay_list[["assay"]]), + assays = SimpleList("relabundance" = merged[["assay"]]), colData = DataFrame(md_df), rowData = DataFrame(tax), rowTree = tree_info[["tree"]], @@ -161,18 +199,95 @@ fetchMetalogTSE <- function( ) # 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 = .parse_metalog_date(data_files[["assay"]]), - date_metadata = .parse_metalog_date(data_files[["md"]]), + 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. @@ -266,15 +381,20 @@ fetchMetalogTSE <- function( download_dir = cache_dir, use.cache = use.cache ) - md_file <- .download_metalog_file( - target_url = md_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 dense matrix (rows = taxa, cols = samples) +# 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 @@ -290,15 +410,13 @@ fetchMetalogTSE <- function( by = .(clade_name, sample_alias)] taxa <- sort(unique(dt$clade_name)) samples <- sort(unique(dt$sample_alias)) - X <- matrix( - 0, nrow = length(taxa), ncol = length(samples), + 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) ) - idx <- cbind( - match(dt$clade_name, taxa), - match(dt$sample_alias, samples) - ) - X[idx] <- dt$rel_abund list(assay = X, taxa = taxa, samples = samples) }