Two-factor analysis with DESeq2

Two_factor_anova_DESeq2.R

library(RColorBrewer)
library(ggplot2)
library(ggrepel)

setwd("<path/to/dir>")

# Load DESeq2 object
load("expression_data/DESeqOBJ.RData")
dds

head(colData(dds))
# Create 2 additional variables: Melatonin and Insulation
newvars <- strsplit(as.character(colData(dds)$Group), "_")
newvars[[1]][1]

melatonin <- unlist(lapply(newvars, function(x) { x[[1]] }))
insulation <- unlist(lapply(newvars, function(x) { x[[2]] }))

colData(dds)$melatonin <- as.factor(melatonin)
colData(dds)$insulation <- as.factor(insulation)

write.csv(as.data.frame(colData(dds)), 
          file =  "sample_layout_interaction.csv")

## Create separate DESeq2 objects for
## testis and epididymis
## Filter out genes with low expressions
dds_tes <- dds[,colData(dds)$Organ == "Testis"]
colData(dds_tes)$Organ
keep <- which(rowSums(counts(dds_tes) >= 10) > 2)
dds_tes <- dds_tes[keep,]
dim(counts(dds_tes))

dds_epi <- dds[,colData(dds)$Organ == "Epididymis"]
colData(dds_epi)$Organ
keep <- which(rowSums(counts(dds_epi) >= 10) > 2)
dds_epi <- dds_epi[keep,]
dim(counts(dds_epi))

## Load annotation 
annot <- read.csv("./docs/gene_annotation.csv", header = T)
head(annot)[1:2,]


######################################################################
## Investigate main effects of insulation and melatonin, as well as 
## their interactions

## ------------------------------------------------------------------ ##
## TESTIS ##
dds_tes$melatonin
dds_tes$insulation <- relevel(dds_tes$insulation, "non-insulated")
dds_tes$insulation

design(dds_tes) <- ~ melatonin + insulation + melatonin:insulation
dds_tes <- DESeq(dds_tes) 
resultsNames(dds_tes)
# "Intercept"                              "melatonin_Melatonin_vs_Control"        
# "insulation_insulated_vs_non.insulated"  "melatoninMelatonin.insulationinsulated"

## The effect of insulation in Controls - Main effect
res <- results(dds_tes, contrast=c("insulation","insulated","non-insulated"))
ix <- which.min(res$padj) # most significant
res <- res[order(res$padj),] # sort
head(res[1:5,-(3:4)])
barplot(assay(dds_tes)[ix,],las=2, main=rownames(dds_tes)[ ix  ]  )

norm_counts <- counts(dds_tes, normalized = T)
colnames(norm_counts) <- dds_tes$SampleID
res_counts <- merge(res, norm_counts, by=0)
head(res_counts)
names(res_counts)
res_annot <- merge(res_counts, annot, by.x= "Row.names", 
                   by.y = "ensembl_gene_id") 
write.csv(res_annot, file="Testis_insulated_vs_non-insulated_in_control.csv",
          row.names = F)

#############################################################################
## The effect of melatonin in non-insulated group
res <- results(dds_tes, contrast=c("melatonin","Melatonin", "Control"))
head(res)

norm_counts <- counts(dds_tes, normalized = T)
colnames(norm_counts) <- dds_tes$SampleID
res_counts <- merge(res, norm_counts, by=0)
head(res_counts)
names(res_counts)
res_annot <- merge(res_counts, annot, by.x= "Row.names", 
                   by.y = "ensembl_gene_id")
write.csv(res_annot, file="Testis_Melatonin_vs_Control_in_non_insulated.csv",
          row.names = F)


#############################################################################
## Interaction between Insulation and Melatonin
res <- results(dds_tes, name="melatoninMelatonin.insulationinsulated")
res <- res[order(res$padj),] # sort
head(res[1:5,-(3:4)])

norm_counts <- counts(dds_tes, normalized = T)
colnames(norm_counts) <- dds_tes$SampleID
res_counts <- merge(res, norm_counts, by=0)
head(res_counts)
names(res_counts)
res_annot <- merge(res_counts, annot, by.x= "Row.names", 
                   by.y = "ensembl_gene_id")
write.csv(res_annot, file="Testis_Melatonin_Insulation_interaction.csv",
          row.names = F)

## Visualize the results of interaction analysis
## Variance stabilizing transformation
vsd_tes <- varianceStabilizingTransformation(dds_tes)
head(assay(vsd_tes))
colData(vsd_tes)
exp_tes <- assay(vsd_tes)
colnames(exp_tes) <- colData(vsd_tes)$SampleID
# Select the genes with significant interaction
res_annot <- res_annot[!is.na(res_annot$padj),]
sig_genes <- res_annot[res_annot$padj < 0.05,]$Row.names
exp_tes <- exp_tes[rownames(exp_tes) %in% sig_genes,] 
exp_tes <- t(exp_tes)
exp_tes <- as.data.frame(exp_tes)
exp_tes <- data.frame(SampleID = rownames(exp_tes),
               melatonin = colData(vsd_tes)$melatonin,
               insulation = colData(vsd_tes)$insulation, exp_tes)
head(exp_tes)
dim(exp_tes)

for (idx in c(4 : 15)) {
  gene_id <- names(exp_tes[,c(1:3, idx)])[4]
  theme_set(theme_bw(16))
  ggplot(exp_tes[,c(1:3, idx)], aes(x=melatonin, y=exp_tes[,c(1:3, idx)][,gene_id])) + 
    geom_boxplot() + stat_summary(fun.y=mean, geom="point", shape=23, size=4) +
    stat_summary(fun=mean, geom="line", aes(group=1)) + ylab(gene_id) +
    facet_grid(~ insulation)
  ggsave( paste(gene_id, "_interaction_boxplot.pdf", sep=""), device = "pdf",
          units = "in", width = 5, height = 5)
}