R/bcbio.R

In vbarrera/bcbioLite: bcbioLite

#' bcbreader
#'
#' This function reads bcbio output.
#' @keywords dataLoading
#' @param projectDir path to final bcbio project folder. Normally
#'   named with date format.
#' @export
#' @examples
#' bcbreader()
bcbreader <- function(projectDir, sampleMetadata = NULL){
    projectDir <- normalizePath(projectDir)
    if (is.null(sampleMetadata)){
        sampleMetadata <- read.csv(file.path(projectDir,"metadata.csv"), row.names = 1)
    }else{
        sampleMetadata <- sampleMetadata
    }
    stopifnot(!is.null(sampleMetadata))
    metrics <- import(file.path(projectDir,"multiqc","multiqc_data","multiqc_bcbio_metrics.txt"))
    metrics <- clean_names(metrics,case = "small_camel") %>% remove_empty("cols")
    sampleDirs = file.path(projectDir,"..", rownames(sampleMetadata))
    tx2genes_file <- file.path(projectDir,"tx2gene.csv")
    salmon_files = file.path(sampleDirs, "salmon", "quant.sf")
    sf_files = salmon_files
    names(sf_files) = rownames(sampleMetadata)
    # FIX check files exists and match sampleMetadata
    tx2gene = read.csv(tx2genes_file, sep=",", row.names=NULL, header=FALSE)
    rownames(tx2gene) <- tx2gene$V1
    txi.salmon = tximport(sf_files,
                          type="salmon",
                          tx2gene = tx2gene,
                          ignoreTxVersion = TRUE,
                          countsFromAbundance="lengthScaledTPM")
    rawCounts = round(data.frame(txi.salmon$counts, check.names=FALSE),0) %>%
        as.matrix()
    colData <- sampleMetadata %>% as.data.frame() %>%
        .[colnames(rawCounts),, drop = FALSE] %>%
        rownames_to_column() %>%
        mutate_all(as.factor) %>%
        mutate_all(droplevels) %>%
        column_to_rownames() %>%
        as("data.frame") %>%
        remove_empty("cols")

    metadata <- list(metrics = metrics,
                     countsFromAbundance = txi.salmon$countsFromAbundance)
    se <- SummarizedExperiment(assays = list(raw = rawCounts,
                                             tpm = txi.salmon$abundance,
                                             length = txi.salmon$length),
                               colData = colData,
                               metadata = metadata)
    invisible(se)
}

#' Change colData information
#'
#' This function will change the colData information from a
#' SummarizedExperiment object.
#'
#' @param se SummarizedExperiment object.
#' @examples
#' data(se)
#' bcbcoldata(se, colData(se))
#' @export
bcbcoldata <- function(se, sampleMetadata){
    sampleMetadata <- DataFrame(sampleMetadata)
    common <- intersect(rownames(sampleMetadata), colnames(se))
    if (!(all(common  %in% colnames(se)))){
        message("samples in the new colData is different than")
        message("the samples already in the SE object.")
        message("Only using common samples.")
        message(paste(common, collapse = " "))
    }
    se <- se[, common]
    colData(se) <- sampleMetadata[common,]
    se
}

#' Normalize SummarizedExperiment from bcbreader function.
#'
#' This function will detect if the object contains the slots from
#' salmon quantification or featureCounts and normalize the data
#' using DESeq2. It will update raw, normalized and vst slots in a
#' SummarizedExperiment object.
#'
#' @param se SummarizedExperiment object generated by readBCBIO.
#' @examples
#' data(se)
#' bcbrun(se, design = ~condition)
#' @export
bcbrun <- function(se, design = ~1){

    stopifnot(class(se) == "SummarizedExperiment")

    if (sum(c("length", "tpm")  %in% names(assays(se))) == 2){
        message("Using tximport data.")
        txi <- list(abundance = assays(se)[["tpm"]],
                    counts = assays(se)[["raw"]],
                    length = assays(se)[["length"]],
                    countsFromAbundance = metadata(se)[["countsFromAbundance"]])
        dds <- DESeqDataSetFromTximport(txi, colData(se), design = design) %>%
            DESeq()
        vst <- varianceStabilizingTransformation(dds)

    }else{
        message("Using raw count assay in se object.")
        dds <- DESeqDataSetFromTximport(assay(se),
                                        colData(se),
                                        design = design) %>%
            DESeq()
        vst <- varianceStabilizingTransformation(dds)
    }
    assays(se)[["raw"]] <-counts(dds)
    assays(se)[["normalized"]] <- counts(dds, normalized = TRUE)
    assays(se)[["vst"]] <- assay(vst)
    se
}

#' Add metrics to colData information
#'
#' This function adds the metrics data to the colData information from a
#' SummarizedExperiment object. This helps the subsetting of
#' SummarizedExperiment object while keeping the original values.
#'
#' @param se SummarizedExperiment object.
#' @examples
#' data(se)
#' addMetrics2colData(se)
#' @export
addMetrics2colData <- function(se){
    coldata_original <- colData(se)
    metrics_original <- se@metadata$metrics
    common <- intersect(rownames(coldata_original), metrics_original$sample)
    if (!(all(common  %in% colnames(se)))){
        message("samples in the new colData is different than")
        message("the samples already in the SE object.")
        message("Only using common samples.")
        message(paste(common, collapse = " "))
    }
    se <- se[, common]

    coldata_new <- coldata_original %>%
        as.data.frame() %>%
        rownames_to_column(var = "sample") %>%
        inner_join(metrics_original, by = "sample") %>%
        remove_rownames() %>%
        column_to_rownames(var = "sample") %>%
        as("DFrame")
    colData(se) <- coldata_new
    se
}