R/write_to_biom.R

In compbiomed/animalcules: Interactive microbiome analysis toolkit

#' Output biom
#'
#' @param MAE A multi-assay experiment object
#' @param path_to_output The folder to output biom file
#' @return A message
#'
#' @import MultiAssayExperiment
#' @importFrom Matrix Matrix
#'
#' @export

write_to_biom <- function(MAE, path_to_output) {
    # Extract data
    if (!requireNamespace("biomformat", quietly = TRUE)) {
        stop("Package 'biomformat' required to continue.",
            "Please install before proceeding.")
    }
    microbe <- MAE[["MicrobeGenetics"]] # double bracket subsetting is easier
    
    tax_table <- as.data.frame(SummarizedExperiment::rowData(microbe))
    sam_table <- as.data.frame(SummarizedExperiment::colData(microbe))
    counts_table <- as.data.frame(SummarizedExperiment::assays(microbe))[
        ,
        rownames(sam_table)
    ] # organism x sample
    
    counts_table_dge <- as(Matrix::Matrix(as.matrix(counts_table)), "unpackedMatrix")
    
    y <- biomformat::make_biom(counts_table, sam_table, tax_table)
    
    
    metadata <- list()
    tax_order_7 <- c("k__", "p__", "c__", "o__", "f__", "g__", "s__")
    
    for (i in seq_len(nrow(tax_table))) {
        metadata[[i]] <- as.character(tax_table[i, ])
        if (length(metadata[[i]]) == 7) {
            for (j in seq_len(7)) {
                metadata[[i]][j] <- paste0(tax_order_7[j], metadata[[i]][j])
            }
        } else if (length(metadata[[i]]) == 6) {
            for (j in seq_len(7)) {
                if (j < 7) {
                    metadata[[i]][j] <- paste0(tax_order_7[j], metadata[[i]][j])
                } else {
                    metadata[[i]] <- c(metadata[[i]], tax_order_7[7])
                }
            }
        } else {
            warning(
                "Please check why the number of taxonomy levels is less than 6."
                )
        }
    }
    
    data_dir <- system.file("extdata/biom_example.rds", package = "animalcules")
    xx <- readRDS(data_dir)
    
    y_rows_new <- NULL
    for (i in seq_len(nrow(tax_table))) {
        mm <- xx$rows[1, ]
        mm[1, 1] <- rownames(tax_table)[i]
        mm$metadata$taxonomy[[1]] <- metadata[[i]]
        rownames(mm) <- as.character(i)
        rownames(mm$metadata) <- as.character(i)
        y_rows_new <- rbind(y_rows_new, mm)
    }
    
    y$rows <- y_rows_new
    
    
    y_cols_new <- NULL
    for (i in seq_len(nrow(sam_table))) {
        mm <- xx$columns[1, ]
        mm[1, 1] <- rownames(sam_table)[i]
        mm$metadata <- sam_table[i, ]
        rownames(mm) <- as.character(i)
        rownames(mm$metadata) <- as.character(i)
        y_cols_new <- rbind(y_cols_new, mm)
    }
    
    y$columns <- y_cols_new
    
    biomformat::write_biom(y, path_to_output)
    return("Successfully output the biom file!")
}