knitr::opts_chunk$set(echo = TRUE)
library(knitr)

Load Example Results

library(HelpersforDESeq2)
library(DESeq2)

data_dir <- system.file("extdata/", package = "HelpersforDESeq2")


load(file.path(data_dir, "read.counts.rda"))
kable(read.counts[1:10,1:6])


load(file.path(data_dir, "dds.rda"))
dds


my_conditions <- colData(dds)$Sample
my_conditions

my_conditions <- factor(my_conditions, levels = unique(my_conditions)[c(1,3,2,5,6,4)])
colData(dds)$Sample <- my_conditions

my_conditions

Read Stats

par(mfrow=c(2,2), mar = c(4,4,2,2), oma = c(2,1,1,1), mgp = c(2,1,0))

my_colors <- c("#8D9093", "#C199CE", "#DA0EA1", "#5026D9", "#90E7F8", "#1971A9")

plotReadStats(read.counts = read.counts,
              conditions = my_conditions, 
              bar_colors = my_colors)

plotLegend(conditions = my_conditions,
           legend_colors = my_colors, 
           legend_size = 1.0)

# Note multi-mappers are zero here, because they were excluded by the STAR aligner settings

Batch correction

library(sva)

log2_counts_uncor <- log2(counts(dds, normalized = TRUE)+1)

batchVar <- colData(dds)$Batch
modcombat <- model.matrix(~Sample, data = colData(dds))

log2_counts_bcor <- ComBat(dat = log2_counts_uncor, 
                           batch = batchVar, 
                           mod = modcombat, 
                           par.prior = TRUE, prior.plots = FALSE)

PCA

par(mfrow=c(1,2), mar = c(4,4,2,2), oma = c(1,1,1,1), mgp = c(2,1,0))


plottingPCA(log2_counts_uncor,
            xcomp = 1,
            ycomp = 2,
            conditions = my_conditions,
            pca_colors = my_colors,
            main_title = "Before Batch Correction",
            quantiles = c(0,1),
            show_labels = FALSE,
            point_size = 1.1,
            my_xlimits = c(-100,100),
            my_ylimits = c(-100,100))


plottingPCA(log2_counts_bcor,
            xcomp = 1,
            ycomp = 2,
            conditions = my_conditions,
            pca_colors = my_colors,
            main_title = "After Batch Correction",
            quantiles = c(0,1),
            show_labels = FALSE,
            point_size = 1.1,
            my_xlimits = c(-100,100),
            my_ylimits = c(-100,100))

plotLegend(conditions = my_conditions,
           legend_colors = my_colors, 
           legend_size = 0.7)

Soon ...



tschauer/HelpersforDESeq2 documentation built on Nov. 11, 2021, 3:10 a.m.