DeSeq2 sample script

################################################################################
# Gaëlle LELANDAIS <[email protected]>
# Stéphane LE CROM <[email protected]>
################################################################################

# library loading
library(DESeq2)

# data reading
countData = read.table("count_data_diffAnalysis.txt", row.names = 1, header = TRUE)

# general information
row.names(countData)
head(countData)

################################################################################

# Loading metadata for the experiment

colData = read.table("design.txt", row.names = 1, header = TRUE)
# DESeqDataSet object creation
dds = DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~condition)
# normalization of counts
# calculation of sizeFactors
dds = estimateSizeFactors(dds)
sizeFactors(dds)
# bar plots
barplot(colSums(counts(dds)))
# boxplots
boxplot(log2(counts(dds)+1))
boxplot(log2(counts(dds,normalized=TRUE)+1))

# variance estimations
dds = estimateDispersions(dds)
plotDispEsts(dds)

# differential analysis
dds = nbinomWaldTest(dds)
res = results(dds)
mcols(res,use.names=TRUE)

# diagnostic plots
# MA plot
plotMA(res)
# Histogram of the p values
hist(res$pvalue,breaks=20,col="grey")
# PCA plots
vsd <- varianceStabilizingTransformation(dds, blind=TRUE)
p = plotPCA(vsd)
p = update(p, panel = function(x, y, ...) {lattice::panel.xyplot(x, y, ...);lattice::ltext(x=x, y=y, labels=rownames(colData(vsd)), pos=1, offset=1, cex=0.8)})
print(p)


# genes are sorted according to adjusted p-values
res = res[order(res$padj),]
write.table(res, "res.txt", row.name=T, quote=F, sep='\t')