RNAseq analysis

analysis.1412/DEanalysis.R

@@ -0,0 +1,48 @@
+# Loading libraries
+library(DESeq2)
+library(data.table)
+
+
+# Creation of sample table of counts files for DESeq2 ------------------------
+
+countfiles = character(4)
+samplename = character(4)
+for (i in 1:4) {
+  countfiles[i] = sprintf("ERR9905%ds.counts", i + 56)
+  samplename[i] = sprintf("ERR9905%d", i + 56)
+}
+condition = factor(c("a", "a", "b", "b"))
+
+sampTable = data.table(sampName = samplename, file = countfiles, condition = condition)
+
+
+# Differential gene expression analysis -------------------------------------
+
+DEdataset = DESeqDataSetFromHTSeqCount(sampleTable = sampTable, design = ~condition)
+dds = DESeq(DEdataset)
+res = results(dds)
+
+
+# Plots and table of differential expressions -------------------------------
+
+png(filename="./graphs/MA-plot.png")
+plotMA(dds, ylim=c(-8, 8))
+dev.off()
+
+png(filename="./graphs/DispersionEst-plot.png")
+plotDispEsts(dds)
+dev.off()
+
+png(filename="./graphs/PCA-plot.png")
+plotPCA(DESeqTransform(dds))
+dev.off()
+
+restable = as.data.frame(res)
+setDT(restable, keep.rownames = "geneName")
+
+
+# Extraction of significantly (padj<0.01) differentially expressed genes ----
+
+res_padj = restable[padj < 0.01][order(padj)]
+
+fwrite(res_padj[,"geneName"], "GeneList.csv", col.names = F)