## Read in input file (output from DESeq2)

igg <- read.delim("201704_HNadimpalli_DESeq2.txt", sep="\t", header=T, as.is=T)
rownames(igg) <- igg$external_gene_name
# select columns relevant for Volcano plots
igg2 <- igg[,c("log2FoldChange_wt_dcr2_vs_wt_igg", "padj_wt_dcr2_vs_wt_igg")]
colnames(igg2) <- c("log2FC", "padj")
# remove NA padj
igg2 <- igg2[!is.na(igg2$padj),]
# convert padj to -log10(padj)
igg2$padj <- -log10(igg2$padj)
# add column with significance
igg2$significant <- "no"
# add 2 lines: 2 levels of significance:
igg2$significant[igg2$padj > -log10(0.05) & igg2$log2FC > log2(1.5)] <- "medium"
igg2$significant[igg2$padj > -log10(0.05) & igg2$log2FC > 1] <- "yes"


mygenes1 <- c("bcd", "Tl")
mygenes2 <- c("bcd", "Tl", "r2d2", "pgc")

# Open pdf file
pdf("20170920_Volcano_igg.pdf")

for(myg in list(mygenes1, mygenes2)){
igg2$names <- NA
igg2$names[rownames(igg2)%in%myg] <- rownames(igg2)[rownames(igg2)%in%myg]

igg2$significant_name <- igg2$significant
igg2$significant_name[!is.na(igg2$names)] <- "name"
igg2$significant_name <- factor(igg2$significant_name, levels=c("no", "medium", "yes", "name"), ordered=T)
p <- ggplot(data=igg2%>%arrange(significant_name), aes(x=log2FC, y=padj, color=significant_name, label=names, size=significant_name, shape=significant_name)) + geom_point() + 
        geom_hline(yintercept=-log10(0.05), col="black") + 
        scale_color_manual(values=c("gray68", "darkred", "red", "blue")) + 
#       geom_vline(xintercept=c(-1, -log2(1.5), log2(1.5), 1), col="black") + 
        geom_vline(xintercept=c(log2(1.5), 1), col="black") + 
        ggtitle("(wt_ip-wt_input) vs (mutant_ip-mutant_input)", subtitle=paste("adjusted P-value < 0.05, positive log2 Fold Change > 1 or >", round(log2(1.5), 2))) + 
        geom_text(aes(label=names), col="blue", size=4, hjust=1.3) + 
        scale_size_manual(values=c(rep(0.6,3), 2)) + 
        scale_shape_manual(values=c(rep(19,3), 17)) + 
        theme_bw() + 
        theme(legend.position="none") +
        scale_y_continuous(name="-log10(adjusted P-value)")
p2 <- p + scale_x_continuous(name="log2(Fold Change)", limits=c(-6, 6))
print(p2)

# zoomed
p3 <- p2 + scale_y_continuous(name="-log10(adjusted P-value)", limits=c(-5, 120))
print(p3)
}
dev.off()