deseq2.Rnw

%%% Applies DESeq on a matrix of count data
%%% Inputs: counts matrix and pheno design matrix

\documentclass{article}
\usepackage[margin=1in]{geometry}
\usepackage{here}

\title{DESeq Analysis}
\author{Raymond Lim}

\begin{document}

\maketitle

\SweaveOpts{cache=T, prefix.string=graphics/deSeq2}


<<init, echo = F, cache = F, results = hide>>=
dir.create('graphics', showWarnings = F)
options(width = 100)

includeGraphic <- function(filename, caption = NULL, width = 1) {
    if (is.null(caption)) {
        cat("\\includegraphics[width=", width, "\\linewidth]{", filename, "}\n", sep = "")
    } else {
        cat("\\begin{figure}[h!]\n")
        cat("\\includegraphics[width=", width, "\\linewidth]{", filename, "}\n", sep = "")
        cat("\\caption{", caption, "}\n", sep = "")
        cat("\\end{figure}\n")
    }

}

includeGraphics <- function(filenames, width = 1, caption) {
    cat("\\begin{figure}[h!]\n")
    for (filename in filenames) {
        includeGraphic(filename, width)
    }
    cat("\\caption{", caption, "}\n", sep = "")
    cat("\\end{figure}\n")
}
@


<<libs, echo = F, results = hide>>=
suppressPackageStartupMessages(library("optparse"))
suppressPackageStartupMessages(library("DESeq2"))
suppressPackageStartupMessages(library("GenomicFeatures"))
suppressPackageStartupMessages(library("Rsamtools"))
suppressPackageStartupMessages(library("xtable"))
#library(multicore)
@

<<optParse, echo = F>>=
optList <- list(
                make_option("--condition", default = 'Condition', help = "Factor of interest in pheno file [default %default]"),
                make_option("--refCondition", default = 'Normal', help = "Reference condition [default %default]"),
                make_option("--outFile", default = NULL, help = "Output results to this file [optional]"));

parser <- OptionParser(usage = "%prog [options] [counts file] [pheno file]", option_list = optList);

arguments <- parse_args(parser, positional_arguments = T, arg = arguments);
opt <- arguments$options;

if (length(arguments$args) != 2) {
    cat("Need pheno design file and counts data");
    print_help(parser);
    stop();
} 
@

<<loadData, echo = F>>=
phenoFile <- arguments$args[2];
countsFile <- arguments$args[1];

pheno <- read.table(phenoFile, header = T, sep = '\t', row.names = 1);
pheno[, opt$condition] <- relevel(pheno[, opt$condition], opt$refCondition);

counts <- read.table(countsFile, header = T, sep = '\t', na.strings = "", comment.char = "", stringsAsFactors = F);
counts <- counts[!duplicated(counts[,1]), ]
rownames(counts) <- counts[,1]
counts <- counts[-1]

if (!all(colnames(counts) %in% rownames(pheno))) {
    cat("Design does not match data");
}
counts <- counts[, rownames(pheno)]

cds <- newCountDataSet(counts, pheno[, opt$condition])

dds <- DESeqDataSetFromMatrix(counts, pheno, design = ~ opt$condition)
@

Estimate the effective library size:

<<effectiveLibSize>>=
cds <- estimateSizeFactors(cds)
sizeFactors(cds)
@

Estimate dispersion/variance: 

<<estimateDispersions>>=
cds <- estimateDispersions(cds)

str(fitInfo(cds))
@

<<diffEx>>=
res <- nbinomTest(cds, levels(pData(cds)$condition)[1], levels(pData(cds)$condition)[2])
@

\begin{figure}
<<MAplot, fig = T>>=
plot(res$baseMean, res$log2FoldChange, log = "x", pch = 20, cex = .3, col = ifelse(res$padj < .1, "red", "black"), ylab = 'M', xlab = 'A')
@
    \caption{MA plot, normalised mean vs. log2 fold change}
\end{figure}

\begin{figure}
<<pvalHist, fig = T>>=
hist(res$pval, breaks = 100, col = 'skyblue', border = 'slateblue', main = "", xlab = 'p-value')
@
    \caption{Histogram of p-values}
\end{figure}

<<topGenes, results = tex>>=
capt <- 'Top differentially expressed genes'
print(xtable(head(res[order(res$padj), ], 20), caption = capt))
@

<<writeResults, echo = F>>=
if (!is.null(opt$outFile)) {
    write.table(res, file = opt$outFile, sep = '\t', quote = F)
}
@

\end{document}