diff --git a/DESCRIPTION b/DESCRIPTION index 2fd252a3d..729179fc5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -78,7 +78,6 @@ Imports: MASS, MatrixGenerics, methods, - ecodive, rlang, S4Vectors, scater, @@ -94,6 +93,7 @@ Suggests: BiocStyle, biomformat, dada2, + ecodive, knitr, mediation, miaTime, diff --git a/NAMESPACE b/NAMESPACE index b0f050187..6e3059057 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -57,6 +57,9 @@ export(estimateDominance) export(estimateEvenness) export(estimateFaith) export(estimateRichness) +export(exportMothur) +export(exportQIIME2) +export(exportRaw) export(full_join) export(getAbundanceClass) export(getAbundant) @@ -413,6 +416,8 @@ importFrom(ape,is.binary) importFrom(ape,is.rooted) importFrom(ape,read.tree) importFrom(ape,reorder.phylo) +importFrom(ape,write.FASTA) +importFrom(ape,write.tree) importFrom(bluster,clusterRows) importFrom(decontam,isContaminant) importFrom(decontam,isNotContaminant) @@ -430,8 +435,6 @@ importFrom(dplyr,select) importFrom(dplyr,summarise) importFrom(dplyr,sym) importFrom(dplyr,tally) -importFrom(ecodive,unweighted_unifrac) -importFrom(ecodive,weighted_unifrac) importFrom(rlang,":=") importFrom(rlang,sym) importFrom(scater,calculateMDS) diff --git a/NEWS b/NEWS index f5f8ec97c..434e3bae4 100644 --- a/NEWS +++ b/NEWS @@ -190,5 +190,6 @@ Changes in version 1.19.x Changes in version 1.21.x + Fixed explained percentages on PCA axes for jointRPCA and hardened rank-deficient SVD handling in joint-RPCA (1.21.1, 2026-05-12) -* Generalise agglomerateByModule to non-binary numerical modules (1.21.2, 2026-05-31) -* Modify joint-RPCA to match with Gemelli implementation (1.21.3, 2026-06-16) ++ Generalise agglomerateByModule to non-binary numerical modules (1.21.2, 2026-05-31) ++ Modify joint-RPCA to match with Gemelli implementation (1.21.3, 2026-06-16) ++ Implement Striped Unifrac algorithm (1.21.3, 2026-06-12) diff --git a/R/AllGenerics.R b/R/AllGenerics.R index 241789a42..5ea880f5d 100644 --- a/R/AllGenerics.R +++ b/R/AllGenerics.R @@ -105,6 +105,21 @@ setGeneric("convertToBIOM", signature = c("x"), setGeneric("convertToPhyloseq", signature = c("x"), function(x, ...) standardGeneric("convertToPhyloseq")) +#' @rdname export-methods +#' @export +setGeneric("exportRaw", signature = c("x"), function(x, ...) + standardGeneric("exportRaw")) + +#' @rdname export-methods +#' @export +setGeneric("exportQIIME2", signature = c("x"), function(x, ...) + standardGeneric("exportQIIME2")) + +#' @rdname export-methods +#' @export +setGeneric("exportMothur", signature = c("x"), function(x, ...) + standardGeneric("exportMothur")) + #' @rdname isContaminant #' @export setGeneric("addContaminantQC", signature = c("x"), diff --git a/R/RcppExports.R b/R/RcppExports.R index 70c294092..1f301c944 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -9,3 +9,7 @@ .Call(`_mia_apply_transformation_difference_or_division`, mat, method) } +.unifrac_cpp <- function(assay, rowTree, weighted) { + .Call(`_mia_unifrac_cpp`, assay, rowTree, weighted) +} + diff --git a/R/addAlpha.R b/R/addAlpha.R index 9928c3d69..5b2e6d3af 100644 --- a/R/addAlpha.R +++ b/R/addAlpha.R @@ -244,9 +244,9 @@ #' evenly the abundances of different species are distributed. The following #' evenness indices are provided: #' -#' By default, four indices are returned, each taking into account different -#' aspects: richness (the number of observed unique features), -#' dominance (Berger-Parker), information (Shannon), and phylogenetics (Faith) +#' By default, four indices are returned, each taking into account different +#' aspects: richness (the number of observed unique features), +#' dominance (Berger-Parker), information (Shannon), and phylogenetics (Faith) #' (Cassol et al., 2025). #' #' The available evenness indices include the following (all in lowercase): @@ -360,6 +360,11 @@ #' in characterizing microbiomes. #' _Genome Res._ 31(11):2131-2137. doi: 10.1101/gr.275777.121 #' +#' Armstrong G. et al. (2021) +#' Efficient computation of Faith's phylogenetic diversity with applications +#' in characterizing microbiomes. +#' _Genome Res._ 31(11):2131-2137. doi: 10.1101/gr.275777.121 +#' #' Beisel J-N. et al. (2003) #' A Comparative Analysis of Diversity Index Sensitivity. #' _Internal Rev. Hydrobiol._ 88(1):3-15. @@ -447,9 +452,9 @@ #' A tribute to Claude Shannon (1916 –2001) and a plea for more rigorous use of #' species richness, species diversity and the ‘Shannon–Wiener’ Index. #' _Alpha Ecology & Biogeography_ 12, 177–197. -#' -#' Cassol, I., Ibañez, M. & Bustamante, J.P. (2025) -#' Key features and guidelines for the application of microbial alpha diversity +#' +#' Cassol, I., Ibañez, M. & Bustamante, J.P. (2025) +#' Key features and guidelines for the application of microbial alpha diversity #' metrics. _Sci Rep_ 15, 622. doi:10.1038/s41598-024-77864-y #' #' @seealso diff --git a/R/addDissimilarity.R b/R/addDissimilarity.R index e4d8e7661..844d80478 100644 --- a/R/addDissimilarity.R +++ b/R/addDissimilarity.R @@ -91,9 +91,8 @@ #' all the abundances of features are equal between two samples, and 0 means #' that samples have completely different relative abundances. #' -#' Unifrac is calculated with \code{ -#' \link[ecodive:unweighted_unifrac]{ecodive:unweighted_unifrac()}} -#' or \code{\link[ecodive:weighted_unifrac]{ecodive:weighted_unifrac()}}. +#' Unifrac is calculated with a C++ implementation of the Striped Unifrac +#' algorithm (McDonald et al. 2021). #' #' If rarefaction is enabled, \code{\link[vegan:avgdist]{vegan:avgdist()}} is #' utilized. @@ -131,6 +130,11 @@ #' #' Lozupone C, Knight R. ``Unifrac: a new phylogenetic method for comparing #' microbial communities.'' Appl Environ Microbiol. 2005 71 (12):8228-35. +#' +#' McDonald D, Vázquez-Baeza Y, Koslicki D, McClelland J, Reeve N, Xu Z, +#' Gonzalez A, Knight R. ``Striped UniFrac: enabling microbiome analysis at +#' unprecedented scale.'' Nat Methods. 2018 15 (11):847-848. +#' doi: 10.1038/s41592-018-0187-8. #' #' For JSD dissimilarity: #' Jensen-Shannon Divergence and Hilbert space embedding. diff --git a/R/calculateUnifrac.R b/R/calculateUnifrac.R index 903daae44..bb9320cd5 100644 --- a/R/calculateUnifrac.R +++ b/R/calculateUnifrac.R @@ -1,8 +1,6 @@ #' @importFrom ape drop.tip -#' @importFrom ecodive weighted_unifrac -#' @importFrom ecodive unweighted_unifrac .get_unifrac <- function( - x, tree, weighted = FALSE, node.label = nodeLab, nodeLab = NULL, ...){ + x, tree, new = FALSE, weighted = FALSE, node.label = nodeLab, nodeLab = NULL, ...){ # Transpose the matrix so that the orientation is the same as in other # dissimilatity methods x <- t(x) @@ -72,9 +70,8 @@ # multiple rows are linked to single tip. x <- .merge_assay_by_rows(x, node.label, ...) - # Calculate unifrac. Use implementation from ecodive package - FUN <- if( weighted ) weighted_unifrac else unweighted_unifrac - res <- FUN(t(x), tree = tree) + # Calculate unifrac with C++ algorithm + res <- .unifrac_cpp(x, tree, weighted) return(res) } diff --git a/R/exporters.R b/R/exporters.R new file mode 100644 index 000000000..4ff078b6b --- /dev/null +++ b/R/exporters.R @@ -0,0 +1,296 @@ +#' Exporters to common formats for microbiome data outside of R +#' +#' @description +#' There are a few very popular external tools for microbiome analysis, +#' including QIIME2 and mothur. However, R does not currently provide any class +#' to accommodate those data formats. When exporting data from mia to external +#' tools, the best approach is therefore to break a data container into its +#' building blocks (assays, side information, trees, etc.). +#' +#' Thanks to \code{exportRaw}, \code{exportQIIME2} and \code{exportMothur}, +#' it is now possible to export a +#' \code{\link[TreeSummarizedExperiment]{TreeSummarizedExperiment}} object as +#' raw elements or near-ready QIIME2 and mothur formats, respectively. This way, +#' migrating from mia to an external system is still a bad idea, but at least it +#' is fairly straightforward. +#' +#' @param x a \code{\link[TreeSummarizedExperiment]{TreeSummarizedExperiment}} +#' object. +#' +#' @param dpath \code{Character scalar}. +#' +#' @param rowdata.file \code{Character scalar}. (Default: \code{"rowdata"}) +#' +#' @param coldata.file \code{Character scalar}. (Default: \code{"coldata"}) +#' +#' @param assay.dir \code{Character scalar}. (Default: \code{"assays"}) +#' +#' @param rowtree.dir \code{Character scalar}. (Default: \code{"row_trees"}) +#' +#' @param coltree.dir \code{Character scalar}. (Default: \code{"col_trees"}) +#' +#' @param dimred.dir \code{Character scalar}. (Default: \code{"dim_reds"}) +#' +#' @param altexp.dir \code{Character scalar}. (Default: \code{"alt_exps"}) +#' +#' @param assay.type \code{Character scalar}. (Default: \code{"counts"}) +#' +#' @param tree.name \code{Character scalar}. (Default: \code{"phylo"}) +#' +#' @param group.var \code{Character scalar}. (Default: \code{NULL}) +#' +#' @param ... Unused. +#' +#' @returns Directory at \code{dpath} with components of \code{x} each stored +#' as a file in the proper format. +#' +#' @details The output directory contains the elements of \code{x}. For +#' \code{exportQIIME2} and \code{exportMothur}, data will need some more +#' processing using the target tool. For some tips, check the rbiom package +#' vignettes on converting data: +#' \url{https://cmmr.github.io/rbiom/articles/convert.html} +#' +#' @examples +#' library(TreeSummarizedExperiment) +#' +#' tse <- makeTSE() +#' assayNames(tse) <- "counts" +#' names(rowData(tse))[1] <- "Genus" +#' +#' # Export raw TreeSE components in custom directory +#' exportRaw(tse, "out") +#' +#' # Export TreeSE components in near-ready QIIME2 format +#' exportQIIME2(tse, "qiime2_dir") +#' +#' # Export TreeSE components in near-ready mothur format +#' exportMothur(tse, "mothur_dir") +#' +#' @name export-methods +#' @aliases exportRaw exportQIIME2 exportMothur +NULL + + +#' @rdname export-methods +#' @importFrom ape write.tree +setMethod("exportRaw", signature = c(x = "TreeSummarizedExperiment"), + function(x, dpath, rowdata.file = "rowdata", coldata.file = "coldata", + assay.dir = "assays", rowtree.dir = "row_trees", coltree.dir = "col_trees", + dimred.dir = "dim_reds", altexp.dir = "alt_exps"){ + + if( !endsWith(dpath, "/") ) dpath <- paste0(dpath, "/") + if( !dir.exists(dpath) ) dir.create(dpath) + + write.table(rowData(x), paste0(dpath, rowdata.file, ".tsv"), sep = "\t") + + write.table(colData(x), paste0(dpath, coldata.file, ".tsv"), sep = "\t") + + assay.dir <- .create_slot_dir(x, assays, assay.dir, dpath) + + for( assay_name in assayNames(x) ){ + write.table( + assay(x, assay_name), + paste0(assay.dir, assay_name, ".tsv"), + sep = "\t" + ) + } + + rowtree.dir <- .create_slot_dir(x, rowTreeNames, rowtree.dir, dpath) + + for( tree_name in rowTreeNames(x) ){ + write.tree( + rowTree(x, tree_name), paste0(rowtree.dir, tree_name, ".nwk") + ) + } + + coltree.dir <- .create_slot_dir(x, colTreeNames, coltree.dir, dpath) + + for( tree_name in colTreeNames(x) ){ + write.tree( + colTree(x, tree_name), paste0(coltree.dir, tree_name, ".nwk") + ) + } + + dimred.dir <- .create_slot_dir(x, reducedDims, dimred.dir, dpath) + + for( dimred_name in reducedDimNames(x) ){ + write.table( + reducedDim(x, dimred_name), + paste0(dimred.dir, dimred_name, ".tsv"), + sep = "\t" + ) + } + + altexp.dir <- .create_slot_dir(x, altExps, altexp.dir, dpath) + + for( altexp_name in altExpNames(x) ){ + + altexp_path <- paste0(altexp.dir, altexp_name) + + if( is(altExp(x, altexp_name), "SummarizedExperiment") ){ + + exportToRaw(altExp(x, altexp_name), altexp_path) + + }else{ + + write.table(altExp(x, altexp_name), altexp_path, sep = "\t") + + } + } + invisible(NULL) +}) + + +.create_slot_dir <- function(x, FUN, slot.path, main.path = ""){ + + if( length(FUN(x)) != 0L ){ + slot.path <- paste0(main.path, slot.path) + if( !endsWith(slot.path, "/") ) slot.path <- paste0(slot.path, "/") + dir.create(slot.path) + } + return(slot.path) +} + + +#' @rdname export-methods +#' @importFrom ape write.tree write.FASTA +setMethod("exportQIIME2", signature = c(x = "TreeSummarizedExperiment"), + function(x, dpath, assay.type = "counts", tree.name = "phylo", + group.var = NULL){ + + if( !endsWith(dpath, "/") ) dpath <- paste0(dpath, "/") + if( !dir.exists(dpath) ) dir.create(dpath) + + if( !is.null(group.var) ){ + + group <- rowData(x)[[group.var]] + group <- gsub("-", "_", group) + group <- cbind(rownames(x), group) + colnames(group) <- c("Feature ID", group.var) + + write.table( + group, paste0(dpath, group.var, ".tsv"), + sep = "\t", quote = FALSE, row.names = FALSE + ) + } + + row_data <- apply(rowData(x)[taxonomyRanks(x)], 1L, paste, collapse = ";_") + row_data <- gsub("(;_|;_NA)+$", "", row_data) + + row_data <- data.frame(rownames(x), row_data, 1L, row.names = NULL) + colnames(row_data) <- c("Feature ID", "Taxon", "Confidence") + + write.table( + row_data, paste0(dpath, "taxonomy.tsv"), + sep = "\t", quote = FALSE, row.names = FALSE + ) + + col_data <- as.data.frame(colData(x)) + + col_data[] <- lapply( + col_data, function(col) if( is.factor(col) ) as.character(col) else col + ) + + col_types <- apply( + col_data, 2L, function(col) switch( + type(col), character = "categorical", integer = , double = "numeric") + ) + + col_data <- rbind(col_types, col_data) + col_data <- cbind(`sample-id` = c("#q2:types", colnames(x)), col_data) + + write.table( + col_data, paste0(dpath, "metadata.tsv"), + sep = "\t", quote = FALSE, row.names = FALSE + ) + + sel_assay <- data.frame(rownames(x), assay(x, assay.type), row.names = NULL) + colnames(sel_assay)[1L] <- "#OTU ID" + + write.table( + sel_assay, paste0(dpath, assay.type, ".tsv"), + sep = "\t", quote = FALSE, row.names = FALSE + ) + + row_tree <- rowTree(x, tree.name) + + if( !is.null(row_tree) ){ + write.tree(row_tree, paste0(dpath, "tree.nwk")) + } + + if( !is.null(referenceSeq(x)) ){ + write.FASTA(referenceSeq(x), paste0(dpath, "seqs.fna")) + } + invisible(NULL) +}) + + +#' @rdname export-methods +#' @importFrom ape write.tree write.FASTA +setMethod("exportMothur", signature = c(x = "TreeSummarizedExperiment"), + function(x, dpath, assay.type = "counts", tree.name = "phylo", + group.var = NULL){ + + if( !endsWith(dpath, "/") ) dpath <- paste0(dpath, "/") + if( !dir.exists(dpath) ) dir.create(dpath) + + rownames(x) <- gsub("-", "_", rownames(x), fixed = TRUE) + colnames(x) <- gsub("-", "_", colnames(x), fixed = TRUE) + + if( !is.null(group.var) ){ + + group <- rowData(x)[group.var] + group[[group.var]] <- gsub("-", "_", group[[group.var]]) + + write.table( + group, paste0(dpath, group.var, ".group"), + sep = "\t", quote = FALSE, col.names = FALSE + ) + } + + row_data <- apply(rowData(x)[taxonomyRanks(x)], 1L, paste, collapse = ";") + row_data <- gsub("(;|;NA)+$", "", row_data) + row_data <- gsub("-", "_", row_data, fixed = TRUE) + + sel_assay <- assay(x, assay.type) + row_sums <- rowSums(sel_assay) + + row_data <- data.frame( + names(row_data), row_sums, row_data, row.names = NULL + ) + colnames(row_data) <- c("OTU", "Size", "Taxonomy") + + write.table( + row_data, paste0(dpath, "taxonomy.tsv"), + sep = "\t", quote = FALSE, row.names = FALSE + ) + + col_data <- data.frame(group = colnames(x), colData(x), row.names = NULL) + + write.table( + col_data, paste0(dpath, "metadata.tsv"), + sep = "\t", quote = FALSE, row.names = FALSE + ) + + sel_assay <- data.frame( + rownames(sel_assay), total = row_sums, sel_assay, row.names = NULL + ) + colnames(sel_assay)[1L] <- "Representative_Sequence" + + write.table( + sel_assay, paste0(dpath, assay.type, ".tsv"), + sep = "\t", quote = FALSE, row.names = FALSE + ) + + row_tree <- rowTree(x, tree.name) + + if( !is.null(row_tree) ){ + row_tree$tip.label <- gsub("-", "_", row_tree$tip.label) + write.tree(row_tree, paste0(dpath, "tree.nwk")) + } + + if( !is.null(referenceSeq(x)) ){ + write.FASTA(referenceSeq(x), paste0(dpath, "seqs.fna")) + } + invisible(NULL) +}) diff --git a/man/export-methods.Rd b/man/export-methods.Rd new file mode 100644 index 000000000..9936c8e15 --- /dev/null +++ b/man/export-methods.Rd @@ -0,0 +1,110 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/exporters.R +\name{exportRaw} +\alias{exportRaw} +\alias{exportQIIME2} +\alias{exportMothur} +\alias{export-methods} +\alias{exportRaw,TreeSummarizedExperiment-method} +\alias{exportQIIME2,TreeSummarizedExperiment-method} +\alias{exportMothur,TreeSummarizedExperiment-method} +\title{Exporters to common formats for microbiome data outside of R} +\usage{ +exportRaw(x, ...) + +exportQIIME2(x, ...) + +exportMothur(x, ...) + +\S4method{exportRaw}{TreeSummarizedExperiment}( + x, + dpath, + rowdata.file = "rowdata", + coldata.file = "coldata", + assay.dir = "assays", + rowtree.dir = "row_trees", + coltree.dir = "col_trees", + dimred.dir = "dim_reds", + altexp.dir = "alt_exps" +) + +\S4method{exportQIIME2}{TreeSummarizedExperiment}(x, dpath, assay.type = "counts", tree.name = "phylo") + +\S4method{exportMothur}{TreeSummarizedExperiment}( + x, + dpath, + assay.type = "counts", + tree.name = "phylo", + group.var = NULL +) +} +\arguments{ +\item{x}{a \code{\link[TreeSummarizedExperiment]{TreeSummarizedExperiment}} +object.} + +\item{...}{Unused.} + +\item{dpath}{\code{Character scalar}.} + +\item{rowdata.file}{\code{Character scalar}. (Default: \code{"rowdata"})} + +\item{coldata.file}{\code{Character scalar}. (Default: \code{"coldata"})} + +\item{assay.dir}{\code{Character scalar}. (Default: \code{"assays"})} + +\item{rowtree.dir}{\code{Character scalar}. (Default: \code{"row_trees"})} + +\item{coltree.dir}{\code{Character scalar}. (Default: \code{"col_trees"})} + +\item{dimred.dir}{\code{Character scalar}. (Default: \code{"dim_reds"})} + +\item{altexp.dir}{\code{Character scalar}. (Default: \code{"alt_exps"})} + +\item{assay.type}{\code{Character scalar}. (Default: \code{"counts"})} + +\item{tree.name}{\code{Character scalar}. (Default: \code{"phylo"})} + +\item{group.var}{\code{Character scalar}. (Default: \code{NULL})} +} +\value{ +Directory at \code{dpath} with components of \code{x} each stored +as a file in the proper format. +} +\description{ +There are a few very popular external tools for microbiome analysis, +including QIIME2 and mothur. However, R does not currently provide any class +to accommodate those data formats. When exporting data from mia to external +tools, the best approach is therefore to break a data container into its +building blocks (assays, side information, trees, etc.). + +Thanks to \code{exportRaw}, \code{exportQIIME2} and \code{exportMothur}, +it is now possible to export a +\code{\link[TreeSummarizedExperiment]{TreeSummarizedExperiment}} object as +raw elements or near-ready QIIME2 and mothur formats, respectively. This way, +migrating from mia to an external system is still a bad idea, but at least it +is fairly straightforward. +} +\details{ +The output directory contains the elements of \code{x}. For +\code{exportQIIME2} and \code{exportMothur}, data will need some more +processing using the target tool. For some tips, check the rbiom package +vignettes on converting data: +\url{https://cmmr.github.io/rbiom/articles/convert.html} +} +\examples{ +library(TreeSummarizedExperiment) + +tse <- makeTSE() +assayNames(tse) <- "counts" +names(rowData(tse))[1] <- "Genus" + +# Export raw TreeSE components in custom directory +exportRaw(tse, "out") + +# Export TreeSE components in near-ready QIIME2 format +exportQIIME2(tse, "qiime2_dir") + +# Export TreeSE components in near-ready mothur format +exportMothur(tse, "mothur_dir") + +} diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml index d4f4fe208..d9d991444 100644 --- a/pkgdown/_pkgdown.yml +++ b/pkgdown/_pkgdown.yml @@ -51,7 +51,7 @@ reference: - summary - getDominant - getAbundant -- title: Data loading +- title: Import/export/convert data - contents: - importBIOM - importQIIME2 @@ -61,6 +61,9 @@ reference: - importTaxpasta - convertFromDADA2 - convertFromPhyloseq + - exportRaw + - exportMothur + - exportQIIME2 - title: Diversity - subtitle: Alpha Diversity diff --git a/src/.gitignore b/src/.gitignore index 9d22eb46a..db62a4ece 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -1,2 +1,2 @@ *.o -*.so +*.so \ No newline at end of file diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 08ec88198..ec36b8f9b 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -34,10 +34,24 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// unifrac_cpp +Rcpp::NumericVector unifrac_cpp(const Rcpp::NumericMatrix& assay, const Rcpp::List& rowTree, bool weighted); +RcppExport SEXP _mia_unifrac_cpp(SEXP assaySEXP, SEXP rowTreeSEXP, SEXP weightedSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const Rcpp::NumericMatrix& >::type assay(assaySEXP); + Rcpp::traits::input_parameter< const Rcpp::List& >::type rowTree(rowTreeSEXP); + Rcpp::traits::input_parameter< bool >::type weighted(weightedSEXP); + rcpp_result_gen = Rcpp::wrap(unifrac_cpp(assay, rowTree, weighted)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_mia_faith_cpp", (DL_FUNC) &_mia_faith_cpp, 2}, {"_mia_apply_transformation_difference_or_division", (DL_FUNC) &_mia_apply_transformation_difference_or_division, 2}, + {"_mia_unifrac_cpp", (DL_FUNC) &_mia_unifrac_cpp, 3}, {NULL, NULL, 0} }; diff --git a/src/assay.cpp b/src/assay.cpp index c98dd0efe..d396fc8c3 100644 --- a/src/assay.cpp +++ b/src/assay.cpp @@ -18,8 +18,8 @@ using namespace su; -Assay::Assay(const Rcpp::NumericMatrix & assay){ - table = assay; +Assay::Assay(const Rcpp::NumericMatrix & assay): + table(assay) { sample_ids = std::vector(); obs_ids = std::vector(); @@ -27,6 +27,9 @@ Assay::Assay(const Rcpp::NumericMatrix & assay){ Rcpp::StringVector rownames = Rcpp::rownames(table); obs_ids = Rcpp::as>(rownames); + Rcpp::StringVector colnames = Rcpp::colnames(table); + sample_ids = Rcpp::as>(colnames); + n_samples = table.ncol(); n_obs = obs_ids.size(); @@ -38,8 +41,7 @@ Assay::Assay(const Rcpp::NumericMatrix & assay){ sample_counts = get_sample_counts(); } -Assay::~Assay(){ -} +Assay::~Assay(){} void Assay::create_id_index(std::vector &ids, std::unordered_map &map){ @@ -59,6 +61,20 @@ std::vector Assay::get_obs_data(const std::string &id) const { return out; } +std::vector Assay::get_obs_data_range(const std::string &id, unsigned int start, unsigned int end, bool normalize) const { + std::vector out = std::vector(); + uint32_t idx = obs_id_index.at(id); + for(unsigned int i = start; i < end; i++) { + if (normalize) { + out.push_back(table(idx, i)/sample_counts[i]); + } + else { + out.push_back(table(idx, i)); + } + } + return out; +} + std::vector Assay::get_sample_counts(){ std::vector sample_counts = std::vector(); diff --git a/src/assay.h b/src/assay.h index 89b50c262..0aff086c3 100644 --- a/src/assay.h +++ b/src/assay.h @@ -41,15 +41,24 @@ class Assay { /* Get a dense vector of observation data * * @param id The observation ID to fetch - * @param out An allocated array of at least size n_samples. - * Values of an index position [0, n_samples) which do not - * have data will be zero'd. */ std::vector get_obs_data(const std::string &id) const; - private: - Rcpp::NumericMatrix table; // Access to raw sample counts in R's memory + /* get a dense vector of a range of observation data + * + * @param id The observation ID to fetch + * @param start Initial index + * @param end First index past the end + * @param normalize If set, divide by sample_counts + */ + std::vector get_obs_data_range(const std::string &id, + unsigned int start, + unsigned int end, + bool normalize) const; + private: + const Rcpp::NumericMatrix & table; // Access to raw sample counts in R's memory + std::vector get_sample_counts(); /* At construction, lookups mapping IDs -> index position within an @@ -58,11 +67,11 @@ class Assay { std::unordered_map obs_id_index; /* Create an index mapping an ID to its corresponding index - * position. - * - * @param ids A vector of IDs to index - * @param map A hash table to populate - */ + * position. + * + * @param ids A vector of IDs to index + * @param map A hash table to populate + */ void create_id_index(std::vector &ids, std::unordered_map &map); diff --git a/src/propmap.cpp b/src/propmap.cpp index ce02cfa1e..a177b101e 100644 --- a/src/propmap.cpp +++ b/src/propmap.cpp @@ -11,34 +11,25 @@ #include "assay.h" #include "propmap.h" -#include -#include -#include -#include -#include -#include -#include - #include using namespace su; -PropMap::PropMap(uint32_t vecsize) +PropMap::PropMap(uint32_t vecsize) : prop_map() , defaultsize(vecsize) { prop_map.reserve(1000); } -PropMap::~PropMap() { -} +PropMap::~PropMap(){} std::vector PropMap::get(uint32_t i){ if( prop_map.count(i) > 0 ){ return prop_map.at(i); } else { - return(std::vector()); - } + return(std::vector()); + } } void PropMap::clear(uint32_t i){ @@ -49,12 +40,38 @@ void PropMap::update(uint32_t node, std::vector vec){ prop_map[node] = vec; } -std::vector su::set_proportions(const BPTree &tree, +PropMapMulti::PropMapMulti(uint32_t _vecsize) + : vecsize(_vecsize) + , multi(get_num_stacks(), PropMap(DEF_VEC_SIZE)) {} + +PropMapMulti::~PropMapMulti(){} + +// Number of stacks = number of def_sizes that go in vecsize +// Rounding up ensures that there are always enough stacks for full vecsize +uint32_t PropMapMulti::get_num_stacks() const { + return (vecsize + (DEF_VEC_SIZE-1)) / DEF_VEC_SIZE; +} + +// get_start and get_end are used only for passing the values to set_prop_range +// and embed_prop_range +uint32_t PropMapMulti::get_start(uint32_t idx) const { + return idx*DEF_VEC_SIZE; +} + +uint32_t PropMapMulti::get_end(uint32_t idx) const { + return std::min((idx+1)*DEF_VEC_SIZE, vecsize); +} + +PropMap & PropMapMulti::get_prop_map(uint32_t idx){ + return multi[idx]; +} + +std::vector su::set_proportions(const BPTree & tree, uint32_t node, - const Assay &table, - PropMap &ps, + const Assay & table, + PropMap & pm, bool normalize){ - std::vector props = std::vector(); + std::vector props = std::vector(table.n_samples, 0.0); if( tree.isleaf(node) ){ std::string leaf = tree.names[node]; props = table.get_obs_data(leaf); // get row for the specified node @@ -66,23 +83,51 @@ std::vector su::set_proportions(const BPTree &tree, } else { unsigned int current = tree.leftchild(node); unsigned int right = tree.rightchild(node); - - for( unsigned int i = 0; i < table.n_samples; i++ ){ - props.push_back(0); - } - + while( current <= right && current != 0 ){ - std::vector vec = ps.get(current); // Pull from prop map - ps.clear(current); // Remove from prop map - + std::vector vec = pm.get(current); // Pull from prop map + pm.clear(current); // Remove from prop map + for( unsigned int i = 0; i < table.n_samples; i++ ){ props[i] = props[i] + vec[i]; } - + current = tree.rightsibling(current); } } - - ps.update(node, props); + + pm.update(node, props); return(props); } + +std::vector su::set_proportions_range(const su::BPTree & tree, + uint32_t node, + const su::Assay & table, + unsigned int start, + unsigned int end, + PropMap & pm, + bool normalize){ + const unsigned int els = end-start; + std::vector props = std::vector(els, 0.0); + if(tree.isleaf(node)) { + std::string leaf = tree.names[node]; + props = table.get_obs_data_range(leaf, start, end, normalize); + } else { + unsigned int current = tree.leftchild(node); + unsigned int right = tree.rightchild(node); + + while(current <= right && current != 0) { + std::vector vec = pm.get(current); // pull from prop map + pm.clear(current); // remove from prop map, place back on stack + + for(unsigned int i = 0; i < els; i++){ + props[i] = props[i] + vec[i]; + } + + current = tree.rightsibling(current); + } + } + + pm.update(node, props); + return props; +} diff --git a/src/propmap.h b/src/propmap.h index ec936d194..62dbca527 100644 --- a/src/propmap.h +++ b/src/propmap.h @@ -7,8 +7,8 @@ * See LICENSE file for more details */ -#ifndef __FAITH_PROPMAP -#define __FAITH_PROPMAP 1 +#ifndef __FAITH_PROPMAP_H +#define __FAITH_PROPMAP_H 1 #include #include @@ -18,6 +18,7 @@ #include "assay.h" namespace su { + class PropMap { public: PropMap(uint32_t vecsize); @@ -31,10 +32,37 @@ class PropMap { uint32_t defaultsize; }; -std::vector set_proportions(const BPTree &tree, uint32_t node, - const Assay &table, - PropMap &ps, +// Helper class that splits the full proportions vector into smaller chunks of +// pre-defined size +class PropMapMulti { + public: + PropMapMulti(uint32_t _vecsize); + ~PropMapMulti(); + + uint32_t get_num_stacks() const; + uint32_t get_start(uint32_t idx) const; + uint32_t get_end(uint32_t idx) const; + PropMap & get_prop_map(uint32_t idx); + + private: + const uint32_t vecsize; // Size of the full vector, equal to n_samples + static const uint32_t DEF_VEC_SIZE = 1024; // size of the sub-vectors + std::vector multi; +}; + +std::vector set_proportions(const BPTree & tree, uint32_t node, + const Assay & table, + PropMap & pm, bool normalize = true); + +std::vector set_proportions_range(const su::BPTree & tree, + uint32_t node, + const su::Assay & table, + unsigned int start, + unsigned int end, + PropMap & pm, + bool normalize = true); + } -#endif /* __FAITH_PROPMAP */ +#endif /* __FAITH_PROPMAP_H */ diff --git a/src/stripemap.cpp b/src/stripemap.cpp new file mode 100644 index 000000000..a478857d2 --- /dev/null +++ b/src/stripemap.cpp @@ -0,0 +1,44 @@ +/* +* BSD 3-Clause License +* +* Copyright (c) 2016-2021, UniFrac development team. +* All rights reserved. +* +* See LICENSE file for more details +*/ + +#include "stripemap.h" + +using namespace su; + +StripeMap::StripeMap(uint32_t _n_samples) + : stripe_map() + , n_samples(_n_samples) +{ + n_stripes = (n_samples + 1) / 2; + for( unsigned int i = 0; i < n_stripes; i++ ){ + this->update(i, std::vector(n_samples, 0.0)); + } +} + +StripeMap::~StripeMap(){} + +std::vector StripeMap::get(uint32_t i){ + if( stripe_map.count(i) > 0 ){ + return stripe_map.at(i); + } else { + return(std::vector()); + } +} + +void StripeMap::clear(uint32_t i){ + stripe_map[i] = std::vector(); +} + +void StripeMap::update(uint32_t node, std::vector vec){ + stripe_map[node] = vec; +} + +bool StripeMap::is_empty(uint32_t i){ + return get(i).empty(); +} \ No newline at end of file diff --git a/src/stripemap.h b/src/stripemap.h new file mode 100644 index 000000000..879e8e203 --- /dev/null +++ b/src/stripemap.h @@ -0,0 +1,38 @@ +/* +* BSD 3-Clause License +* +* Copyright (c) 2016-2021, UniFrac development team. +* All rights reserved. +* +* See LICENSE file for more details +*/ + +#ifndef __UNIFRAC_STRIPEMAP_H +#define __UNIFRAC_STRIPEMAP_H 1 + +#include +#include +#include +#include + +namespace su { + +class StripeMap { + public: + StripeMap(uint32_t n_samples); + virtual ~StripeMap(); + + void clear(uint32_t i); + void update(uint32_t i, std::vector vec); + std::vector get(uint32_t i); + bool is_empty(uint32_t i); + + private: + std::unordered_map> stripe_map; + uint32_t n_samples; + uint32_t n_stripes; +}; + +} + +#endif /* __UNIFRAC_STRIPEMAP_H */ diff --git a/src/tree.cpp b/src/tree.cpp index 480d92cfd..a7a042390 100644 --- a/src/tree.cpp +++ b/src/tree.cpp @@ -167,8 +167,7 @@ BPTree BPTree::collapse() { return this->mask(collapsemask, new_lengths); } -BPTree::~BPTree(){ -} +BPTree::~BPTree(){} void BPTree::index_and_cache(){ // Should probably do the open/close in here too diff --git a/src/tree.h b/src/tree.h index 05b116e61..0ba293bb2 100644 --- a/src/tree.h +++ b/src/tree.h @@ -19,6 +19,7 @@ #include namespace su { + class BPTree { public: /* Tracked attributes */ @@ -134,7 +135,8 @@ class BPTree { int32_t bwd(uint32_t i, int32_t d) const; int32_t enclose(uint32_t i) const; - }; +}; + } #endif /* __FAITH_TREE_H */ diff --git a/src/unifrac.cpp b/src/unifrac.cpp new file mode 100644 index 000000000..0db135aa1 --- /dev/null +++ b/src/unifrac.cpp @@ -0,0 +1,245 @@ +/* + * BSD 3-Clause License + * + * Copyright (c) 2016-2021, UniFrac development team. + * All rights reserved. + * + * See LICENSE file for more details + */ + +#include "unifrac.h" +#include "propmap.h" +#include "stripemap.h" +#include "tree.h" + +su::mat_t su::one_off(const su::Assay & table, + const su::BPTree & tree, + bool weighted, + bool bypass_tips){ + + //Number of stripes to be used, basically half of samples + const unsigned int stripe_stop = (table.n_samples + 1) / 2; + + su::StripeMap dm_stripes(table.n_samples); + su::StripeMap dm_stripes_total(table.n_samples); + + su::task_parameters task; + + //thread id - currently single-threaded + task.tid = 0; + //Stripes to start and stop on - single task, so the entire thing + task.start = 0; + task.stop = stripe_stop; + task.bypass_tips = bypass_tips; + + task.n_samples = table.n_samples; + + + su::unifrac(std::ref(table), + std::ref(tree), + std::ref(dm_stripes), + std::ref(dm_stripes_total), + weighted, + task); + + + su::mat_t result; + result.n_samples = table.n_samples; + result.cf_size = su::comb_2(table.n_samples); + result.sample_ids = table.sample_ids; + result.is_upper_triangle = true; + result.condensed_form = su::stripes_to_condensed_form(dm_stripes, + table.n_samples, + task.start, + task.stop); + + return result; +} + +void su::unifrac(const su::Assay &table, + const su::BPTree &tree, + su::StripeMap & dm_stripes, + su::StripeMap & dm_stripes_total, + bool weighted, + const su::task_parameters task_p){ + // unweighted + if (weighted == false){ + unifracTT( + table, tree, true, dm_stripes, dm_stripes_total, + task_p ); + } + // weighted unnormalized + else { + unifracTT( + table, tree, false, dm_stripes, dm_stripes_total, + task_p ); + } +} + +template +inline void su::unifracTT(const su::Assay & table, + const su::BPTree & tree, + const bool want_total, + su::StripeMap & dm_stripes, + su::StripeMap & dm_stripes_total, + const su::task_parameters & task_p){ + const unsigned int n_samples = task_p.n_samples; + const uint64_t n_samples_r = ((n_samples + UNIFRAC_BLOCK-1) / + UNIFRAC_BLOCK)*UNIFRAC_BLOCK; // round up + + su::PropMapMulti propmap_multi(table.n_samples); + + const unsigned int max_emb = TaskT::RECOMMENDED_MAX_EMBS; + + TaskT taskObj(dm_stripes, dm_stripes_total, max_emb, task_p); + + std::vector lengths = std::vector(max_emb); + + /* + * The values in the example vectors correspond to index positions of an + * element in the resulting distance matrix. So, in the example below, + * the following can be interpreted: + * + * [0 1 2] + * [1 2 3] + * + * As comparing the sample for row 0 against the sample for col 1, the + * sample for row 1 against the sample for col 2, the sample for row 2 + * against the sample for col 3. + * + * In other words, we're computing stripes of a distance matrix. In the + * following example, we're computing over 6 samples requiring 3 + * stripes. + * + * A; stripe == 0 + * [0 1 2 3 4 5] + * [1 2 3 4 5 0] + * + * B; stripe == 1 + * [0 1 2 3 4 5] + * [2 3 4 5 0 1] + * + * C; stripe == 2 + * [0 1 2 3 4 5] + * [3 4 5 0 1 2] + * + * The stripes end up computing the following positions in the distance + * matrix. + * + * x A B C x x + * x x A B C x + * x x x A B C + * C x x x A B + * B C x x x A + * A B C x x x + * + * However, we store those stripes as vectors, ie + * [ A A A A A A ] + * + * We end up performing N / 2 redundant calculations on the last stripe + * (see C) but that is small over large N. + */ + + unsigned int k = 0; // index in tree + const unsigned int max_k = (tree.nparens / 2) - 1; + + const unsigned int num_prop_chunks = propmap_multi.get_num_stacks(); + // num_prop_chunks = 1 + + while(k node_proportions + = su::set_proportions_range(tree, + node, + table, + tstart, + tend, + propmap); + + if (task_p.bypass_tips && tree.isleaf(node)){ + continue; + } + + // they all do the same thing, so enough for the first to update + // the global state + if (ck==0){ + lengths[filled_emb] = tree.lengths[node]; + filled_emb++; + } + + taskObj.embed_proportions_range(node_proportions, + tstart, + tend, + my_filled_emb); + my_filled_emb++; + } + + // they all do the same thing, so enough for the first to update the + // global state + if (ck==0){ + k=my_k; + } + } + taskObj._run(filled_emb,lengths); + filled_emb=0; + } + + if(want_total){ + const uint64_t start_idx = task_p.start; + const uint64_t stop_idx = task_p.stop; + + for(uint64_t i = start_idx; i < stop_idx; i++){ + for(uint64_t j = 0; j < n_samples; j++) { + uint64_t idx = ((i-start_idx)*n_samples_r)+j; + taskObj.dm_stripes.buf[idx] = taskObj.dm_stripes.buf[idx] / + taskObj.dm_stripes_total.buf[idx]; + } + } + } +} + +std::vector su::stripes_to_condensed_form(su::StripeMap & stripes, + uint32_t n, + unsigned int start, + unsigned int stop) { + // n must be >= 2, but that should be enforced upstream as that would imply + // computing unifrac on a single sample. + + uint64_t comb_N = comb_2(n); + std::vector cf = std::vector(comb_N, 0.0); + + for(unsigned int stripe = start; stripe < stop; stripe++){ + std::vector dm_stripe = stripes.get(stripe); + // compute the (i, j) position of each element in each stripe + uint64_t i = 0; + uint64_t j = stripe + 1; + for(uint64_t k = 0; k < n; k++, i++, j++){ + if(j == n){ + i = 0; + j = n - (stripe + 1); + } + // determine the position in the condensed form vector for a given + // (i, j) + // based on https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.squareform.html + uint64_t comb_N_minus_i = comb_2(n - i); + cf[comb_N - comb_N_minus_i + (j - i - 1)] = dm_stripe[k]; + } + } + return cf; +} \ No newline at end of file diff --git a/src/unifrac.h b/src/unifrac.h new file mode 100644 index 000000000..c9deb04da --- /dev/null +++ b/src/unifrac.h @@ -0,0 +1,87 @@ +/* + * BSD 3-Clause License + * + * Copyright (c) 2016-2021, UniFrac development team. + * All rights reserved. + * + * See LICENSE file for more details + */ + +#ifndef __UNIFRAC_H +#define __UNIFRAC_H 1 + +#include +#include +#include +#include + +#include "assay.h" +#include "tree.h" +#include "propmap.h" +#include "unifrac_task.h" + +namespace su { + +typedef struct mat { + unsigned int n_samples; + unsigned int cf_size; + bool is_upper_triangle; + std::vector condensed_form; + std::vector sample_ids; +} mat_t; + +su::mat_t one_off(const su::Assay & table, + const su::BPTree & tree, + bool weighted, + bool bypass_tips); + +// Chooses the right task for the job and constructs a unifracTT +void unifrac(const su::Assay &table, + const su::BPTree &tree, + su::StripeMap & dm_stripes, + su::StripeMap & dm_stripes_total, + bool weighted, + const su::task_parameters task_p); + +// Works the vectors +template +inline void unifracTT(const su::Assay & table, + const su::BPTree & tree, + const bool want_total, + su::StripeMap & dm_stripes, + su::StripeMap & dm_stripes_total, + const su::task_parameters & task_p); + +inline uint64_t comb_2(uint64_t N) { + // based off of _comb_int_long + // https://github.com/scipy/scipy/blob/v0.19.1/scipy/special/_comb.pyx + + // Compute binom(N, k) for integers. + // + // we're disregarding overflow as that practically should not + // happen unless the number of samples processed is in excess + // of 4 billion + uint64_t val, j, M, nterms; + uint64_t k = 2; + + M = N + 1; + nterms = k < (N - k) ? k : N - k; + + val = 1; + + for(j = 1; j < nterms + 1; j++) { + val *= M - j; + val /= j; + } + return val; +} + +// Stripes to condensed form for the results +std::vector stripes_to_condensed_form(su::StripeMap & stripes, + uint32_t n, + unsigned int start, + unsigned int stop); + +} + +#endif /* __UNIFRAC_H */ diff --git a/src/unifrac_R.cpp b/src/unifrac_R.cpp new file mode 100644 index 000000000..5b882e155 --- /dev/null +++ b/src/unifrac_R.cpp @@ -0,0 +1,87 @@ +/* + * BSD 3-Clause License + * + * Copyright (c) 2016-2021, UniFrac development team. + * All rights reserved. + * + * See LICENSE file for more details + */ + +#include +#include + +#include + +#include "assay.h" +#include "tree.h" +#include "propmap.h" +#include "stripemap.h" + +#include "unifrac.h" + +// Calculate Unifrac +// +// This function calculates Unifrac distances for a given assay and rowTree, +// using a C++ implementation of the Striped Unifrac algorithm. +// +// @details +// This function makes several assumptions about the contents of +// \code{assay} and \code{rowTree}, namely that: +// \itemize{ +// \item \code{assay} and \code{rowTree} are both non-empty. +// \item \code{assay} has row and column names. +// \item \code{rowTree}'s nodes are arranged in cladewise order. +// } +// These checks should all be handled in the surrounding R code. +// +// The C++ code was adapted from an implementation by the Unifrac team +// (Armstrong et al. 2021), which is licensed under the BSD 3-Clause license. +// +// @param assay An R numeric matrix containing the assay of a \code{TreeSE} +// object. +// @param rowTree An \code{ape::phylo} object containing the rowTree of a +// \code{TreeSE} object. +// @param weighted Boolean: Whether to calculate unweighted or weighted Unifrac. +// @param bypass_tips Boolean: Whether to bypass tips during calculations. This +// speeds up calculations considerably, and does not seem to have a noticeable +// effect on the results. +// @return A vector containing Unifrac distances. +// +// @keywords internal +// [[Rcpp::export(.unifrac_cpp)]] +Rcpp::NumericVector unifrac_cpp(const Rcpp::NumericMatrix & assay, + const Rcpp::List & rowTree, + bool weighted){ + + su::BPTree tree = su::BPTree(rowTree); + su::Assay table = su::Assay(assay); + + std::unordered_set to_keep(table.obs_ids.begin(), + table.obs_ids.end()); + + su::BPTree tree_sheared = tree.shear(to_keep).collapse(); + + su::mat_t results = su::one_off(table, tree_sheared, weighted, false); + + unsigned int n = results.condensed_form.size(); + Rcpp::NumericVector unifrac = Rcpp::NumericVector(n); + // + // Rcpp::List::create(Rcpp::Named("n_samples") = result->n_samples, + // Rcpp::Named("is_upper_triangle") = result->is_upper_triangle, + // Rcpp::Named("cf_size") = result->cf_size, + // Rcpp::Named("c_form") = cf) + // + for(unsigned int i = 0; i < n; i++){ + unifrac[i] = results.condensed_form[i]; + } + + unifrac.attr("class") = "dist"; + Rcpp::StringVector labels(table.n_samples); + labels = table.sample_ids; + unifrac.attr("Labels") = labels; + unifrac.attr("Size") = table.n_samples; + unifrac.attr("Diag") = false; + unifrac.attr("Upper") = false; + + return unifrac; +} diff --git a/src/unifrac_task.cpp b/src/unifrac_task.cpp new file mode 100644 index 000000000..c35d32a51 --- /dev/null +++ b/src/unifrac_task.cpp @@ -0,0 +1,354 @@ +/* + * BSD 3-Clause License + * + * Copyright (c) 2016-2021, UniFrac development team. + * All rights reserved. + * + * See LICENSE file for more details + */ + +#include +#include +#include +#include +#include + +#include "tree.h" +#include "unifrac_task.h" + +using namespace su; + +UnifracTaskVector::UnifracTaskVector(su::StripeMap & _dm_stripes, + const su::task_parameters _task_p) + : dm_stripes(_dm_stripes) + , start_idx(_task_p.start) + , n_samples(_task_p.n_samples) + , n_samples_r(((n_samples + UNIFRAC_BLOCK-1)/UNIFRAC_BLOCK)*UNIFRAC_BLOCK) + , task_p(_task_p) +{ + + // The buffer is only needed if the stripes are non-empty + buf = (dm_stripes.is_empty(start_idx)) ? std::vector() : + std::vector(n_samples_r* + (task_p.stop-start_idx), 0.0); + // Copy stripe values to buffer + if (!buf.empty()){ + for(unsigned int stripe=start_idx; stripe < task_p.stop; stripe++) { + std::vector dm_stripe = dm_stripes.get(stripe); + for(unsigned int k=0; k < dm_stripe.size(); k++) { + buf[ (stripe-start_idx)*n_samples_r + k ] = dm_stripe[k]; + } + } + } +} + +UnifracTaskVector::~UnifracTaskVector(){ + // If the buffer isn't empty, copy its values back to the stripes + if (!buf.empty()){ + for(unsigned int stripe=start_idx; stripe < task_p.stop; stripe++) { + std::vector dm_stripe = dm_stripes.get(stripe); + for(unsigned int k=0; k < dm_stripe.size(); k++) { + dm_stripe[k] = buf[ (stripe-start_idx)*n_samples_r + k ]; + } + dm_stripes.update(stripe, dm_stripe); + } + } +} + +template +UnifracTaskBase::UnifracTaskBase(su::StripeMap & _dm_stripes, + su::StripeMap & _dm_stripes_total, + unsigned int _max_embs, + su::task_parameters _task_p) + : dm_stripes(_dm_stripes,_task_p) + , dm_stripes_total(_dm_stripes_total,_task_p) + , task_p(_task_p) + , max_embs(_max_embs) +{ + embedded_proportions = initialize_embedded(dm_stripes.n_samples_r, + _max_embs); +} + +template +UnifracTaskBase::~UnifracTaskBase(){} + +template +std::vector UnifracTaskBase::initialize_embedded( + const uint64_t n_samples_r, + unsigned int max_embs){ + uint64_t bsize = n_samples_r * get_emb_els(max_embs); + return std::vector(bsize); +} + +template +void UnifracTaskBase::embed_proportions(const std::vector & in, + unsigned int emb){ + embed_proportions_range(in,0,dm_stripes.n_samples,emb); +} + + +template +UnifracTask::UnifracTask(su::StripeMap & _dm_stripes, + su::StripeMap & _dm_stripes_total, + unsigned int _max_embs, + su::task_parameters _task_p) + : UnifracTaskBase(_dm_stripes, + _dm_stripes_total, + _max_embs, + _task_p){} + +template +UnifracTask::~UnifracTask(){} + +UnifracUnweightedTask::UnifracUnweightedTask(su::StripeMap & _dm_stripes, + su::StripeMap & _dm_stripes_total, + unsigned int _max_embs, + su::task_parameters _task_p) + : UnifracTask(_dm_stripes,_dm_stripes_total,_max_embs,_task_p) +{ + const unsigned int bsize = _max_embs*32; + sums = std::vector(bsize, 0.0); +} + +UnifracUnweightedTask::~UnifracUnweightedTask(){} + +void UnifracUnweightedTask::run(unsigned int filled_embs, + const std::vector & lengths){ + _run(filled_embs, lengths); +} + +void UnifracUnweightedTask::_run(unsigned int filled_embs, + const std::vector & lengths){ + const uint64_t start_idx = this->task_p.start; + const uint64_t stop_idx = this->task_p.stop; + const uint64_t n_samples = this->task_p.n_samples; + const uint64_t n_samples_r = this->dm_stripes.n_samples_r; + + const uint64_t step_size = UnifracUnweightedTask::step_size; + const uint64_t sample_steps = (n_samples+(step_size-1))/step_size; + + const uint64_t filled_embs_els = filled_embs/64; + const uint64_t filled_embs_rem = filled_embs%64; + + const uint64_t filled_embs_els_round = (filled_embs+63)/64; + + // pre-compute sums of length elements, since they are likely to be accessed + // many times + // We will use a 8-bit map, to keep it small enough to keep in L1 cache + for(uint64_t emb_el=0; emb_el pl = std::vector(8); + + + uint64_t len_off = emb8*8; + + // compute all the combinations for this block (8-bits total) + // psum[0] = 0.0 // +0*pl[0]+0*pl[1]+0*pl[2]+... + // psum[1] = pl[0] // +0*pl[1]+0*pl[2]+... + // psum[2] = pl[1] // +0*pl[0]+0*pl[2]+ + // psum[2] = pl[0] + pl[1] + // ... + // psum[255] = pl[1] +.. + pl[7] // + 0*pl[0] + // psum[255] = pl[0] +pl[1] +.. + pl[7] + for(uint64_t b8_i=0; b8_i<0x100; b8_i++){ + sums[(emb8<<8) + b8_i] = + (((b8_i >> 0) & 1) * lengths[len_off + 0]) + + (((b8_i >> 1) & 1) * lengths[len_off + 1]) + + (((b8_i >> 2) & 1) * lengths[len_off + 2]) + + (((b8_i >> 3) & 1) * lengths[len_off + 3]) + + (((b8_i >> 4) & 1) * lengths[len_off + 4]) + + (((b8_i >> 5) & 1) * lengths[len_off + 5]) + + (((b8_i >> 6) & 1) * lengths[len_off + 6]) + + (((b8_i >> 7) & 1) * lengths[len_off + 7]); + } + } + } + + if (filled_embs_rem>0){ // add also the overflow elements + const uint64_t emb_el=filled_embs_els; + for(uint64_t sub8=0; sub8<8; sub8++){ + // we are summing we have enough buffer in sums + const uint64_t emb8 = emb_el*8+sub8; + + // compute all the combinations for this block, set to 0 any past + // the limit as above + for(uint64_t b8_i=0; b8_i<0x100; b8_i++){ + double val= 0; + for(uint64_t li=(emb8*8); li> (li-(emb8*8))) & 1) * lengths[li]; + } + sums[(emb8<<8) + b8_i] = val; + } + + } + } + + // point of thread + for(uint64_t sk = 0; sk < sample_steps ; sk++){ + for(uint64_t stripe = start_idx; stripe < stop_idx; stripe++){ + for(uint64_t ik = 0; ik < step_size ; ik++){ + // within-stripe index (0:n_samples-1) + const uint64_t k = sk*step_size + ik; + //buffer index + const uint64_t idx = (stripe-start_idx) * n_samples_r; + + if (k>=n_samples) continue; // past the limit + + const uint64_t l1 = (k + stripe + 1)%n_samples; // wraparound + + bool did_update = false; + double my_stripe = 0.0; + double my_stripe_total = 0.0; + + //Main calculation phase + for(uint64_t emb_el=0; emb_el> 8) & 0xff)] + + sums[sums_off + 0x200+((x1 >> 16) & 0xff)] + + sums[sums_off + 0x300+((x1 >> 24) & 0xff)] + + sums[sums_off + 0x400+((x1 >> 32) & 0xff)] + + sums[sums_off + 0x500+((x1 >> 40) & 0xff)] + + sums[sums_off + 0x600+((x1 >> 48) & 0xff)] + + sums[sums_off + 0x700+((x1 >> 56) )]; + my_stripe_total += sums[sums_off + (o1 & 0xff)] + + sums[sums_off + 0x100+((o1 >> 8) & 0xff)] + + sums[sums_off + 0x200+((o1 >> 16) & 0xff)] + + sums[sums_off + 0x300+((o1 >> 24) & 0xff)] + + sums[sums_off + 0x400+((o1 >> 32) & 0xff)] + + sums[sums_off + 0x500+((o1 >> 40) & 0xff)] + + sums[sums_off + 0x600+((o1 >> 48) & 0xff)] + + sums[sums_off + 0x700+((o1 >> 56) )]; + } + } + + if (did_update){ + dm_stripes.buf[idx + k] += my_stripe; + dm_stripes_total.buf[idx + k] += my_stripe_total; + } + } + } + } +} + +UnifracUnnormalizedWeightedTask::UnifracUnnormalizedWeightedTask( + su::StripeMap & _dm_stripes, + su::StripeMap & _dm_stripes_total, + unsigned int _max_embs, + su::task_parameters _task_p) + : UnifracTask(_dm_stripes,_dm_stripes_total,_max_embs,_task_p) +{ + const unsigned int n_samples = this->task_p.n_samples; + zcheck = std::vector(n_samples, 0); + sums = std::vector(n_samples, 0.0); +} + +UnifracUnnormalizedWeightedTask::~UnifracUnnormalizedWeightedTask(){} + +void UnifracUnnormalizedWeightedTask::run(unsigned int filled_embs, + const std::vector & lengths){ + _run(filled_embs, lengths); +} + +void UnifracUnnormalizedWeightedTask::_run(unsigned int filled_embs, + const std::vector & lengths){ + const uint64_t start_idx = this->task_p.start; + const uint64_t stop_idx = this->task_p.stop; + const uint64_t n_samples = this->task_p.n_samples; + const uint64_t n_samples_r = this->dm_stripes.n_samples_r; + + const uint64_t step_size = UnifracUnnormalizedWeightedTask::step_size; + const uint64_t sample_steps = (n_samples+(step_size-1))/step_size; + + // check for zero values and pre-compute single column sums + for(uint64_t k=0; k=n_samples) continue; // past the limit + + const uint64_t l1 = (k + stripe + 1)%n_samples; // wraparound + + const bool allzero_k = zcheck[k]; + const bool allzero_l1 = zcheck[l1]; + + if (allzero_k && allzero_l1) { + // nothing to do, would have to add 0 + } + else { + double my_stripe; + + if (allzero_k || allzero_l1){ + // one side has all zeros + // we can use the distributed property, and use the + // pre-computed values + + const uint64_t ridx = (allzero_k) ? l1 : k; + // if (nonzero_l1) ridx=l1 // fabs(k-l1), with k==0 + // if (nonzero_k) ridx=k // fabs(k-l1), with l1==0 + + my_stripe = sums[ridx]; + + } + else { + // both sides non zero, use the explicit but slow + // approach + my_stripe = 0.0; + for(uint64_t emb=0; emb the number of samples being processed +* start the first stripe to process +* stop the last stripe to process +* tid the thread identifier +* bypass_tips ignore tips on compute, reduces compute by ~50% +* g_unifrac_alpha an alpha value for generalized unifrac +*/ + +struct task_parameters { + uint32_t n_samples; // number of samples + unsigned int start; // starting stripe + unsigned int stop; // stopping stripe + unsigned int tid; // thread ID + bool bypass_tips; // avoid compute at tips +}; + +// Helper class that manages stripes +class UnifracTaskVector { + public: + su::StripeMap & dm_stripes; + const unsigned int start_idx; + const unsigned int n_samples; + const uint64_t n_samples_r; + std::vector buf; + + UnifracTaskVector(su::StripeMap & _dm_stripes, + const su::task_parameters _task_p); + + //Destructor copies the buffer values back into dm_stripes + ~UnifracTaskVector(); + + private: + const su::task_parameters task_p; +}; + +// Base task class to be shared by all tasks +// Templated to allow proportions to be embedded as either doubles (weighted) or +// packed bools (unweighted) +template +class UnifracTaskBase { + public: + UnifracTaskVector dm_stripes; + UnifracTaskVector dm_stripes_total; + + su::task_parameters task_p; + + const unsigned int max_embs; + //Continuous vector - each stripe has n_samples_r elements + //Has at most max_embs stripes - when filled, results stored by the + //task's _run() function and embeds are cleared for the next batch + std::vector embedded_proportions; + + UnifracTaskBase(su::StripeMap & _dm_stripes, + su::StripeMap & _dm_stripes_total, + unsigned int _max_embs, + su::task_parameters _task_p); + + virtual ~UnifracTaskBase(); + + // Templated function used when initializing embeds + static unsigned int get_emb_els(unsigned int max_embs); + + static std::vector initialize_embedded( + const uint64_t n_samples_r, + unsigned int max_embs); + + // Store proportions from in into embedded_proportions + void embed_proportions(const std::vector & in, + unsigned int emb); + + void embed_proportions_range(const std::vector & in, + unsigned int start, + unsigned int end, + unsigned int emb); + + protected: + void embed_proportions_range_straight(std::vector & out, + const std::vector & in, + unsigned int start, + unsigned int end, + unsigned int emb) const { + const unsigned int n_samples = dm_stripes.n_samples; + const uint64_t n_samples_r = dm_stripes.n_samples_r; + const uint64_t offset = emb * n_samples_r; + + //Copy to stripe indicated by emb + //Stripes are all contained in in/out in one mass + //Start/end aren't necessarily the whole stripe? + for(unsigned int i = start; i < end; i++){ + out[offset + i] = in[i-start]; + } + + if (end==n_samples){ + // avoid NaNs + for(unsigned int i = n_samples; i < n_samples_r; i++){ + out[offset + i] = 0.0; + } + } + } + + // packed bool + // Compute (in[:]>0) on each element, and store only the boolean bit. + // The output values are stored in a multi-byte format, one bit per emb + // index, so it will likely take multiple passes to store all the values + // Note: assumes we are processing emb in increasing order, starting + // from 0 + std::vector embed_proportions_range_bool( + std::vector out, + const std::vector & in, + unsigned int start, + unsigned int end, + unsigned int emb) const { + const unsigned int n_packed = sizeof(uint64_t)*8; + const unsigned int n_samples = dm_stripes.n_samples; + const uint64_t n_samples_r = dm_stripes.n_samples_r; + + // The output values are stored in a multi-byte format, one bit per + // emb index + // Compute the element to store the bit into, as well as which bit + // in that element + unsigned int emb_block = emb/n_packed; // beginning of block + unsigned int emb_bit = emb%n_packed; // bit inside the elements + const uint64_t offset = emb_block * n_samples_r; + + if (emb_bit == 0){ + // assign for emb_bit==0, so it clears the other bits + // assumes we processing emb in increasing order starting from 0 + for(unsigned int i = start; i < end; i++){ + out[offset + i] = (in[i - start] > 0); + } + + if (end == n_samples){ + // avoid NaNs + for(unsigned int i = n_samples; i < n_samples_r; i++) { + out[offset + i] = 0; + } + } + } + else { + // just update my bit + for(unsigned int i = start; i < end; i++){ + out[offset + i] |= (uint64_t(in[i-start] > 0) << emb_bit); + } + + // the rest of the els are already OK + } + return out; + } +}; + + +template<> inline void UnifracTaskBase::embed_proportions_range( + const std::vector & in, + unsigned int start, + unsigned int end, + unsigned int emb){ + embed_proportions_range_straight(embedded_proportions,in,start,end,emb); +} + +template<> inline unsigned int UnifracTaskBase::get_emb_els( + unsigned int max_embs){ + return max_embs; +} + +template<> inline void UnifracTaskBase::embed_proportions_range( + const std::vector & in, + unsigned int start, + unsigned int end, + unsigned int emb){ + embedded_proportions = embed_proportions_range_bool(embedded_proportions, + in, + start, + end, + emb); +} + +template<> inline unsigned int UnifracTaskBase::get_emb_els( + unsigned int max_embs){ + return (max_embs+63)/64; +} + + +/* void unifrac tasks +* +* all methods utilize the same function signature. that signature is as follows: +* +* dm_stripes vector the stripes of the distance matrix being accumulated +* into for unique branch length +* dm_stripes vector the stripes of the distance matrix being accumulated +* into for total branch length (e.g., to normalize unweighted unifrac) +* embedded_proportions the proportions vector for a sample, or rather +* the counts vector normalized to 1. this vector is embedded as it is +* duplicated: if A, B and C are proportions for features A, B, and C, the +* vector will look like [A B C A B C]. +* length the branch length of the current node to its parent. +* task_p task specific parameters. +*/ + +template +class UnifracTask : public UnifracTaskBase { + public: + UnifracTask(su::StripeMap & _dm_stripes, + su::StripeMap & _dm_stripes_total, + unsigned int _max_embs, + su::task_parameters _task_p); + + virtual ~UnifracTask(); + + virtual void run(unsigned int filled_embs, + const std::vector & lengths) = 0; + + protected: + // Controls the size of inner loops in the calculation phase + static const unsigned int step_size = 4; + + // Max embs are theoretically optimized for cache performance + static const unsigned int RECOMMENDED_MAX_EMBS_STRAIGHT = 64-16; + static const unsigned int RECOMMENDED_MAX_EMBS_BOOL = 64*32; + +}; + + + +class UnifracUnweightedTask : public UnifracTask { + public: + static const unsigned int RECOMMENDED_MAX_EMBS + = UnifracTask::RECOMMENDED_MAX_EMBS_BOOL; + + // Note: _max_emb MUST be multiple of 64 + UnifracUnweightedTask(su::StripeMap & _dm_stripes, + su::StripeMap & _dm_stripes_total, + unsigned int _max_embs, + su::task_parameters _task_p); + + virtual ~UnifracUnweightedTask(); + + virtual void run(unsigned int filled_embs, + const std::vector & lengths); + + void _run(unsigned int filled_embs, + const std::vector & lengths); + private: + std::vector sums; // temp buffer +}; + + + +class UnifracUnnormalizedWeightedTask : public UnifracTask { + public: + static const unsigned int RECOMMENDED_MAX_EMBS + = UnifracTask::RECOMMENDED_MAX_EMBS_STRAIGHT; + + UnifracUnnormalizedWeightedTask(su::StripeMap & _dm_stripes, + su::StripeMap & _dm_stripes_total, + unsigned int _max_embs, + su::task_parameters _task_p); + + virtual ~UnifracUnnormalizedWeightedTask(); + + virtual void run(unsigned int filled_embs, + const std::vector & lengths); + + void _run(unsigned int filled_embs, + const std::vector & lengths); + + protected: + // temp buffers + std::vector zcheck; + std::vector sums; +}; + +} + +#endif /* __UNIFRAC_TASK_H */ diff --git a/tests/testthat/test-5Unifrac.R b/tests/testthat/test-5Unifrac.R index 74006b6a4..7ea9b7ba0 100644 --- a/tests/testthat/test-5Unifrac.R +++ b/tests/testthat/test-5Unifrac.R @@ -65,6 +65,8 @@ test_that("Unifrac beta diversity", { weighted = FALSE, tree.name = "tree2") ) unifrac_mia <- as.matrix(unifrac_mia) + + skip_if_not(requireNamespace("ecodive", quietly = TRUE)) unifrac_ecodive <- as.matrix(ecodive::unweighted_unifrac(t(assay(tse_ref)), rowTree(tse_ref))) expect_equal(unifrac_mia, unifrac_ecodive)