Page Comparison

#set the working directory
setwd("/path/to/input/files")

#install the following packages if not yet installed
library("ggrepel")
library("DESeq2")
library("ggplot2")
library("EnhancedVolcano")
library("reshape2")
library("plyr")
require(reshape2)

##STEP1: import gene counts
countData <- read.csv('t1_melanocytes_control__vs__melanocytes_treatment__RawCounts_integers_GeneSymbols.csv', header = TRUE, sep = ",")
head(countData)

##STEP2: import metaData
metaData <- read.csv('metadata_t1.v2.csv', header = TRUE, sep = ",")
metaData

##STEP3: Construct DESEQDataSet Object
dds <- DESeqDataSetFromMatrix(countData=countData,                               
                              colData=metaData,                               
                              design=~celltype, tidy = TRUE)
#define the reference group
dds$celltype <- factor(dds$celltype, levels = c("control","cancer"))

##STEP4: Run DESEQ function
dds <- DESeq(dds)

##STPE5: Take a look at the results table
res <- results(dds)
head(results(dds, tidy=TRUE)) #let's look at the results table

##STEP6: Summary of differential gene expression
summary(res) 

##STEP7: Sort summary list by p-value
res <- res[order(res$padj),]
head(res)

##STEP8: save DEGs to file
write.table(res, file = "t1_melanocytes_control__vs__melanocytes_treatment__DESeq2_toptags_DEGs.txt")

##STEP9: save normalised counts to file
foo <- counts(dds, normalized = TRUE)
write.csv(foo, file="t1_melanocytes_control__vs__melanocytes_treatment__DESeq2_toptags_DEGs__Normalised_Counts.csv")

##STEP10: PCA
#version 1
pdf(file="Figure_t1_melanocytes_control__vs__melanocytes_treatment__PCA.pdf")
vsdata <- vst(dds, blind=TRUE, nsub = 1000)
plotPCA(vsdata, intgroup="dex") 
dev.off()

#version 2:
pdf(file="Figure_t1_melanocytes_control__vs__melanocytes_treatment_label.pdf")
p <- plotPCA(vsdata, intgroup="dex") 
p + geom_text_repel(aes(label=substr(name, start = 1, stop = 5)),size = 2, min.segment.length = 0, seed = 42, box.padding = 0.5, max.overlaps = 20)
dev.off()

#STEP11: VOLCANO PLOTS
#basic volcano plot 
EnhancedVolcano(res,
                lab = rownames(res),
                x = 'log2FoldChange',
                y = 'pvalue')

#save enhanced volcano plot to file
pdf(file="Figure_t1_melanocytes_control__vs__melanocytes_treatment__enhanced_volcano_plot_2.0FC_p1e-5.pdf")
EnhancedVolcano(res,
                lab = rownames(res),
                x = 'log2FoldChange',
                y = 'pvalue',
                title = 'melanocytes control vs melanocytes treatment',
                pCutoff = 1e-5,
                FCcutoff = 2.0,
                pointSize = 1.0,
                labSize = 3.0)
dev.off()

#save as high resolution jpeg
library(Cairo)
Cairo(file="Figure_t1_melanocytes_control__vs__melanocytes_treatment__enhanced_volcano_plot_2.0FC_p1e-5.jpeg",
      type="png",
      units="px",
      width=2500,
      height=2500,
      pointsize=12,
      dpi=300)
EnhancedVolcano(res,
      lab = rownames(res),
      x = 'log2FoldChange',
      y = 'pvalue',
      title = 'melanocytes control vs melanocytes treatment',
      pCutoff = 1e-5,
      FCcutoff = 2.0,
      pointSize = 1.0,
      labSize = 3.0)
dev.off()