DESeq2

Legacy
##STEP1: import gene counts
countData <- read.table('salmon.merge.gene_counts.tsv', header = TRUE, sep = "\t")
head(countData)

##STEP2: import metaData
metaData <- read.csv('metadata.tsv*', header = TRUE, sep = ",")
metaData

##STEP3: Construct DESEQDataSet Object
dds <- DESeqDataSetFromMatrix(countData=countData, 
                              colData=metaData, 
                              design=~celltype, tidy = TRUE)

#define the reference group set as the first mentioned in levels
dds$celltype <- factor(dds$celltype, levels = c("CON","TB"))

##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) #summary of results


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

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


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


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

#version 2:
pdf(file="Figure_NoHeadInjury_CON__vs__HeadInjury_TB_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()

## jpeq
library(Cairo)

Cairo(file="Figure_NoHeadInjury_CON__vs__HeadInjury_TB__PCA_label.jpeg",
      type="png",
      units="px", 
      width=2500, 
      height=2000, 
      pointsize=12, 
      dpi=300)
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_NoHeadInjury_CON__vs__HeadInjury_TB__enhanced_volcano_plot_2.0FC_p1e-2.pdf")
EnhancedVolcano(res,
                lab = rownames(res),
                x = 'log2FoldChange',
                y = 'pvalue',
                title = 'NoHeadInjury_CON__vs__HeadInjury_TB',
                pCutoff = 1e-2,
                FCcutoff = 1.5,
                pointSize = 1.0,
                labSize = 3.5)
dev.off()

pdf(file="Figure_NoHeadInjury_CON__vs__HeadInjury_TB__enhanced_volcano_plot_2.0FC_p10e-3.pdf")
EnhancedVolcano(res,
                      lab = rownames(res),
                      x = "log2FoldChange",
                      y = "pvalue",
                      pCutoff = 10e-3,
                      FCcutoff = 2,
                      ylim = c(0, -log10(10e-12)),
                      pointSize = c(ifelse(res$log2FoldChange>2, 2, 2)),
                      labSize = 3.0,
                      #shape = c(6, 6, 19, 16),
                      title = "DESeq2: NoHeadInjury_CON vs. HeadInjury_TB",
                      subtitle = "Differential expression",
                      caption = bquote(~Log[2]~ "fold change cutoff, 2; p-value cutoff, 10e-3"),
                      #legendPosition = "right",
                      legendLabSize = 8,
                      col = c("grey30", "forestgreen", "royalblue", "red2"),
                      colAlpha = 0.9,
                      drawConnectors = TRUE,
                      hline = c(10e-3),
                      widthConnectors = 0.5)
dev.off()


#save as high resolution jpeg
library(Cairo)

## Initialize the device for saving first plot.
## Plot will have length=500 pixel, width=500 pixel

Cairo(file="Figure_NoHeadInjury_CON__vs__HeadInjury_TB__enhanced_volcano_plot_2.0FC_p1e-3.jpeg",
      type="png",
      units="px", 
      width=2800, 
      height=2500, 
      pointsize=12, 
      dpi=300)
EnhancedVolcano(res,
                lab = rownames(res),
                x = 'log2FoldChange',
                y = 'pvalue',
                title = 'NoHeadInjury_CON__vs__HeadInjury_TB',
                pCutoff = 1e-3,
                FCcutoff = 2.0,
                pointSize = 1.0,
                labSize = 3.0)
dev.off()


#option2:
Cairo(file="Figure_NoHeadInjury_CON__vs__HeadInjury_TB__enhanced_volcano_plot_2.0FC_p10e-3_v2.jpeg",
      type="png",
      units="px", 
      width=2500, 
      height=2200, 
      pointsize=14, 
      dpi=300)
EnhancedVolcano(res,
                lab = rownames(res),
                x = "log2FoldChange",
                y = "pvalue",
                pCutoff = 10e-3,
                FCcutoff = 2,
                ylim = c(0, -log10(10e-12)),
                pointSize = c(ifelse(res$log2FoldChange>2, 2, 2)),
                labSize = 3.5,
                #shape = c(6, 6, 19, 16),
                title = "DESeq2: NoHeadInjury_CON vs. HeadInjury_TB",
                subtitle = "Differential expression",
                caption = bquote(~Log[2]~ "fold change cutoff, 2; p-value cutoff, 10e-3"),
                #legendPosition = "right",
                legendLabSize = 8,
                col = c("grey30", "forestgreen", "royalblue", "red2"),
                colAlpha = 0.9,
                drawConnectors = TRUE,
                hline = c(10e-3),
                widthConnectors = 0.5)
dev.off()




#save as high res tiff

tiff(file="Figure_NoHeadInjury_CON__vs__HeadInjury_TB__enhanced_volcano_plot_2.0FC_p1e-5.tiff",
     width=8, height=8, units="in", res=150)
EnhancedVolcano(res,
                lab = rownames(res),
                x = 'log2FoldChange',
                y = 'pvalue',
                title = 'NoHeadInjury_CON__vs__HeadInjury_TB',
                pCutoff = 1e-5,
                FCcutoff = 2.0,
                pointSize = 1.5,
                labSize = 3.5)

dev.off()