example differential expression with limma voom

voom_de.R

library(limma)
library(tidyverse)

dge <- DGEList(counts(sce_de))
dge <- calcNormFactors(dge)

design <- model.matrix(~ (dbz_cluster_str == "Unknown"), colData(sce_de)) # Your design matrix here

v <- voom(dge, design, plot = TRUE)

fit <- lmFit(v, design)
fit <- eBayes(fit)
res <- decideTests(fit)

# Correct the p-values and put into a data frame
qvals <- p.adjust(fit$p.value[,2], method = 'BH')

df_limma <- data_frame(log2foldchange = fit$coefficients[,2], 
                       pval = fit$p.value[,2],
                       qval = qvals,
                       ensembl_id = rownames(sce_de),
                       gene_symbol = rowData(sce_de)$symbol,
                       is_significant = qval < 0.05) 
                       
# Nice volcano plot (will need customized)

sig_cols <- c("FALSE" = "grey80", "TRUE" = "darkred")
theme_set(cowplot::theme_cowplot(font_size = 11))

df_text <- dplyr::filter(df_limma, -log10(qval) > 160)

ggplot(df_limma, aes(x = log2foldchange, y = -log10(qval), color = is_significant)) +
  geom_point(alpha = 0.7) +
  scale_color_manual(values = sig_cols, name = "Significantly differentially expressed") +
  geom_text_repel(data = df_text, aes(label = gene_symbol), color = 'black', size = 2) +
  labs(x = expression(log[2]~"(fold change)"),
       y = expression(-log[10]~"(q-value)"),
       subtitle = "Differential expression on non-CNV genes") +
  theme(legend.position = "bottom",
        legend.box.background = element_rect(linetype = 1, size = 1)) +
  annotate("label", x = 2, y = 0, label = "Upregulated in Unknown") +
  annotate("label", x = -1.2, y = 0, label = "Upregulated in Ciliated cells")