In steinbaugh/DESeqAnalysis: Acid Genomics DESeq2 Analysis Utilities

## nolint start
suppressPackageStartupMessages({
    library(basejump)
    library(DESeqAnalysis)
})
prepareTemplate()
source("_setup.R")
## nolint end
## This is required for tabset R Markdown headers to work in `lapply()` call.
knitr::opts_chunk[["set"]](message = FALSE)
initDir(params[["output_dir"]])

Load DESeqAnalysis object

if (is.character(params[["object"]])) {
    file <- params[["object"]]
    object <- import(file)
    name <- params[["name"]]
    if (!isTRUE(nzchar(name))) {
        name <- basenameSansExt(file)
    }
    rm(file)
} else {
    object <- params[["object"]]
    name <- params[["name"]]
}
stopifnot(
    is(object, "DESeqAnalysis"),
    isTRUE(nzchar(name))
)
if (!is.null(params[["alpha_threshold"]])) {
    alphaThreshold(object) <- params[["alpha_threshold"]]
}
if (!is.null(params[["lfc_shrink"]])) {
    lfcShrink(object) <- params[["lfc_shrink"]]
}
if (!is.null(params[["lfc_threshold"]])) {
    lfcThreshold(object) <- params[["lfc_threshold"]]
}
if (!is.null(params[["base_mean_threshold"]])) {
    baseMeanThreshold(object) <- params[["base_mean_threshold"]]
}
resultsNames <- resultsNames(object)
invisible(validObject(object))
print(object)

Sample data

sampleData(object)

DEG summary

plotDegStackedBar(object)

pdf(
    file = file.path(
        params[["output_dir"]],
        "plot-deg-upset.pdf"
    ),
    width = 15L,
    height = 10L
)
plotDegUpset(object)
dev.off()

MA plot {.tabset}

An MA plot compares transformed counts on M (log ratio) and A (mean average) scales [@Yang2002-sx].

invisible(Map(
    f = function(object, i) {
        markdownHeader(i, level = 2L, asis = TRUE)
        print(plotMa(object = object, i = i))
    },
    i = resultsNames,
    MoreArgs = list("object" = object)
))

Volcano plot {.tabset}

A volcano plot compares significance (BH-adjusted P value) against fold change (log2) [@Cui2003-rn; @Li2014-ll]. Genes in the green box with text labels have an adjusted P value are likely to be the top candidate genes of interest.

invisible(Map(
    f = function(object, i) {
        markdownHeader(i, level = 2L, asis = TRUE)
        print(plotVolcano(object = object, i = i))
    },
    i = resultsNames,
    MoreArgs = list("object" = object)
))

DEG PCA {.tabset}

invisible(Map(
    f = function(object, i, ...) {
        markdownHeader(i, level = 2L, asis = TRUE)
        print(plotDegPca(object = object, i = i, ...))
    },
    i = resultsNames,
    MoreArgs = list(
        "object" = object,
        "contrastSamples" = params[["contrast_samples"]]
    )
))

DEG Heatmap {.tabset}

This plot shows only differentially expressed genes on a per-sample basis. We have scaled the data by row (z-score) and used the ward.D2 method for clustering [@Ward1963-xf].

invisible(Map(
    f = function(object, i, ...) {
        markdownHeader(i, level = 2L, asis = TRUE)
        plotDegHeatmap(object = object, i = i, ...)
    },
    i = resultsNames,
    MoreArgs = list(
        "object" = object,
        "contrastSamples" = params[["contrast_samples"]]
    )
))

Top tables {.tabset}

Only the top up- and down-regulated genes (arranged by adjusted P value) are shown.

invisible(Map(
    f = function(object, i, n) {
        markdownHeader(text = i, level = 2L, asis = TRUE)
        markdownTables(object = object, i = i, n = n)
    },
    i = resultsNames,
    MoreArgs = list(
        "object" = object,
        "n" = 25L
    )
))

Export

Here we are exporting the data stored inside our DESeqDataSet, DESeqTransform, and DESeqResults objects (results = unshrunken LFC; lfcShrink = shrunken LFC), along with the differentially expressed gene (DEG) results tables (resultsTables).

The DEG results tables are sorted by BH-adjusted P value (except for the "all" file, which is unmodified), and contain the following columns:

• baseMean: Mean of the normalized counts per gene for all samples.
• log2FoldChange: log2 fold change.
• lfcSE: log2 standard error.
• stat: Wald statistic.
• pvalue: Walt test P value.
• padj: BH adjusted Wald test P value (corrected for multiple comparisons; aka FDR).

These files also contain genome annotation information and the size-factor adjusted normalized counts stored inside the DESeqDataSet used to perform the differential expression analysis.

Note that colData files contain annotations corresponding to the columns (i.e. samples) and the rowData files contain annotations corresponding to the rows (i.e. genes). You can obtain gene-to-symbol mappings using this rowData file, for example.

export(
    object = object,
    con = params[["output_dir"]],
    compress = FALSE
)

steinbaugh/DESeqAnalysis documentation built on April 1, 2024, 8:30 a.m.