In tschauer/HelpersforDESeq2: This is a collection of helper functions for DESeq2

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

Load Gene Annotation

library(HelpersforDESeq2)

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

library(rtracklayer)
# Annotation rda file was generated from the gtf annotation
# gtf <- import("dmel-all-r6.17.gtf")
# save(gtf, file = "dmel-all-r6.17.rda")

load(file.path(data_dir, "dmel-all-r6.17.rda"))

gtf[1:5,1:8]

Make Count Table

SampleTable <- read.table(file.path(data_dir, "SampleTable.txt"), 
                           header = T, sep = "\t", stringsAsFactors = FALSE)

SampleTable <- SampleTable[order(SampleTable$ID),]

kable(SampleTable, row.names = F)

count_files <- file.path(data_dir, list.files(path = data_dir, pattern = "ReadsPerGene"))
count_names <- gsub(".*\\/|_[G,A,T,C].*","", count_files)

if(identical(count_names, SampleTable$ID)){

        read.counts <-  makeCountTable(count_names = count_names, 
                               count_files = count_files,
                               stranded = TRUE)
}

save(read.counts, file = paste0(data_dir, "read.counts.rda"))

kable(read.counts[1:10,1:6])

Setup DESeq2

library(DESeq2)

read.counts.genes <- read.counts[-1:-4,]

if(identical(colnames(read.counts), SampleTable$ID)){

        dds <- setupDDS(CountTableName = "read.counts.genes", 
                        SampleTableName = "SampleTable",
                        SampleIdName = "ID", 
                        ConditionName = "Condition", 
                        BatchName = "Batch",
                        n_samples_for_filtering = ncol(read.counts.genes)*0.75,
                        min_number_of_reads = 1)
        }

save(dds, file = paste0(data_dir, "dds.rda"))

dds

Get Results

contrast_list <- list(c("HighMLEpATP", "Input"),
                      c("LowMLEpATP",  "Input"),
                      c("HighMLEpATP", "HighMLEmATP"),
                      c("LowMLEpATP", "LowMLEmATP"))

for(contrast_name in contrast_list){

        res <- getResults(dds = dds,
                          contrast = contrast_name,
                          lfc_cutoff = 0,
                          shrink = FALSE,
                          annotation = "gtf",
                          anno_symbol = "gene_symbol",
                          anno_id = "gene_id")


        res_name <- paste0("res.", contrast_name[1],"-",contrast_name[2])
        assign(res_name, res)

        write.table(res, file = paste0(data_dir, res_name, ".txt"), 
                    quote = F, sep = "\t", row.names = T, col.names = NA) 

}


head(res)

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

plottingMA(res = res,
           main_title = gsub("res.","",res_name),
           selection_ids = c("roX1","roX2"),
           selection_id_type = "gene_symbol",
           selection_point_size = 1,
           selection_text_label = TRUE,
           selection_shadow = FALSE,
           xlims = c(0, 6),
           ylims = c(-10,10),
           x_axis_by = 2,
           padj_cutoff = 0.01,
           show_legend = TRUE)

Interaction Term

SampleTable2 <- SampleTable[grep("ATP",SampleTable$Condition),]

SampleTable2$ATP <- factor(as.integer(grepl("mATP", SampleTable2$Condition)))
SampleTable2$MLE <- factor(as.integer(grepl("High", SampleTable2$Condition)))

kable(SampleTable2, row.names = F)

read.counts.genes.2 <- read.counts.genes[,colnames(read.counts.genes) %in% SampleTable2$ID]

if(identical(colnames(read.counts.genes.2), SampleTable2$ID)){

        ddsI <- setupDDSwithInteraction(CountTableName = "read.counts.genes.2",
                                        SampleTableName = "SampleTable2",
                                        SampleIdName = "ID",
                                        FactorAName = "ATP",
                                        FactorBName = "MLE",
                                        BatchName = "Batch",
                                        n_samples_for_filtering = ncol(read.counts.genes.2)*0.75,
                                        min_number_of_reads = 1)
}

ddsI

resultsNames(ddsI)

resI <- getResults(dds = ddsI,
                  result_name = "FactorA1.FactorB1",
                  lfc_cutoff = 0,
                  annotation = "gtf",
                  anno_symbol = "gene_symbol",
                  anno_id = "gene_id")
head(resI)

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

res1_name = "res.HighMLEpATP-HighMLEmATP"
res2_name = "res.LowMLEpATP-LowMLEmATP"

plotLog2FC(res1 = get(res1_name),
           res2 = get(res2_name),
           main_title = "",
           x_label = paste("log2FC \n", gsub("res.","",res1_name)),
           y_label = paste("log2FC \n", gsub("res.","",res2_name)),
           lims = c(-10,10),
           point_size = 0.33,
           selection_ids = resI$gene_symbol[resI$padj < 0.1 & abs(resI$log2FoldChange) > 2],
           selection_id_type = "gene_symbol",
           selection_point_size = 1,
           selection_legend = NULL,
           selection_text_label = TRUE,
           selection_text_size = 0.6)