# Differential gene expression analysis
# Comparison involved: LL vs ENL

library("DESeq2")
library(ggplot2)

# Loading data
count_data <- read.csv("count_data_GSE16844.csv, header = TRUE, sep = ",")
head(count_data)

# Loading metadata
meta_data <- read.csv("metadata_GSE16844.csv", header = TRUE, sep = ",")
meta_data

# Constructing DESeq dataset object
dds <- DESeqDataSetFromMatrix(countData = count_data, colData = meta_data, design=~group, tidy = TRUE)
dds

# Running DESeq 
dds <- DESeq(dds)

# Storing results
res <- results(dds)
head(results(dds, tidy=TRUE))
summary(res)
res <- res[order(res$padj),]
head(res)

# Making volcano plots
with(res, plot(log2FoldChange, -log10(pvalue), pch=20, main="Volcano plot", xlim=c(-8,8)))

# Adding coloured points based on a threshold
with(subset(res, padj < 0.05 && abs (log2FoldChange) < 1), points(log2FoldChange, -log10(pvalue), pch=20, col="green"))
with(subset(res, padj < 0.05 && abs(log2FoldChange) > 1), points(log2FoldChange, -log10(pvalue), pch=20, col="red"))

# Making PCA plot
X <- vst(dds, blind=FALSE)
plotPCA(X, intgroup="group")