add corrected volcano plot function

vs_ripchip/scripts/process_mogene.1.Rmd

@@ -124,6 +124,75 @@ PCAplot.sd.eset <- function(eset=NULL,  title=NULL){
 fmt <- function(){ 
   function(x) format(x,nsmall = 1,scientific = FALSE)
 }
+
+vcplot <- function(stats, title="Volcano Plot with Marginal Distributions", pval.cutoff=0.05, lfc.cutoff=1, shade.colour="green", shade.alpha=0.25, point.colour="gray", point.alpha=0.75, point.outline.colour="darkgray", line.colour="gray") {
+  # get range of log fold change and p-value values to setup plot borders
+  range.lfc <- c(floor(min(stats$logFC)), ceiling(max(stats$logFC)))
+  range.pval <- c(floor(min(-log10(stats$adj.P.Val))), ceiling(max(-log10(stats$adj.P.Val))))
+  
+  #make top plot - density plot with fold changes
+  lfcd <- as.data.frame(cbind(density(stats$logFC)$x, density(stats$logFC)$y))
+  hist_top <- ggplot(data=stats, aes(x=logFC))+
+    geom_density(color=line.colour)+
+    geom_ribbon(data=subset(lfcd, V1>lfc.cutoff),aes(x=V1,ymax=V2),ymin=0,fill=shade.colour, alpha=shade.alpha)+
+    theme_bw()+ 
+    theme(axis.title.x=element_blank())+
+    theme(plot.margin=unit(c(3,-5.5,4,3), "mm") )+
+    scale_x_continuous(limits = range.lfc, breaks = range.lfc[1]:range.lfc[2], expand = c(.05,.05))+
+    scale_y_continuous(labels=fmt())
+  
+  # make blank plot
+  empty <- ggplot()+geom_point(aes(1,1), colour="white")+
+    theme(panel.grid=element_blank(),
+          axis.ticks=element_blank(), 
+          panel.background=element_blank(), 
+          axis.text.x=element_blank(),  
+          axis.text.y=element_blank(),          
+          axis.title.x=element_blank(), 
+          axis.title.y=element_blank()
+    )
+  
+  #make scatter volcano plot
+  scat.poly.up <- with(stats, data.frame(x=as.numeric(c(lfc.cutoff,  lfc.cutoff, max(range.lfc),max(range.lfc))), y=as.numeric(c(-log10(pval.cutoff), max(range.pval), max(range.pval),-log10(pval.cutoff)))))
+  
+  scatter <- ggplot(data=stats, aes(x=logFC, y=-log10(adj.P.Val))) +
+    geom_point(alpha=point.alpha, pch=21, fill=point.colour, color=point.outline.colour) +
+    geom_polygon(data=scat.poly.up, aes(x=x,y=y), fill=shade.colour, alpha=shade.alpha) +
+    xlab("log2 fold change") + ylab("-log10(adjusted p-value)") +
+    theme_bw()+
+    theme(legend.position="none") +
+    theme(plot.margin=unit(c(3,-5.5,4,3), "mm") )+
+    scale_x_continuous(limits = range.lfc, breaks = range.lfc[1]:range.lfc[2], expand = c(.05,.05))+
+    scale_y_continuous(labels=fmt(), limits = range.pval)
+  
+  # make right plot - density plot of adjusted pvalues
+  pvald <- as.data.frame(cbind(density(-log10(stats$adj.P.Val))$x, density(-log10(stats$adj.P.Val))$y))
+  hist_right <- ggplot(data=stats, aes(x=-log10(adj.P.Val)))+
+    geom_density(color=line.colour)+
+    geom_ribbon(data=subset(pvald, V1>-log10(pval.cutoff)),aes(x=V1,ymax=V2),ymin=0,fill=shade.colour, alpha=shade.alpha)+
+    theme_bw()+coord_flip()+
+    scale_x_continuous(limits = range.pval)+
+    theme(axis.title.y=element_blank())+ 
+    theme(plot.margin=unit(c(3,-5.5,4,3), "mm"))
+  
+  # plot all plots
+  pp.logfc <- ggplotGrob(hist_top)
+  pp.empty <- ggplotGrob(empty)
+  pp.volc <- ggplotGrob(scatter)
+  pp.pval  <- ggplotGrob(hist_right)
+  grid.arrange(top=textGrob(title),
+               arrangeGrob(pp.logfc,
+                           pp.volc, 
+                           heights=c(1,3),
+                           ncol=1),
+               arrangeGrob(pp.empty,
+                           pp.pval,  
+                           heights=c(1,3),
+                           ncol=1), 
+               ncol=2, 
+               widths=c(3,1))
+} 
+
 ```
 
 ---
@@ -465,8 +534,8 @@ Here we can visulize the relationship between the fold changes in expression obs
 3) Lower right - a density plot (smoothed histogram) of the adjusted pvalued observed for the contrast, the part of the distribution above `r pvalue.cutoff` is highlighted under the curve in `r highlight.color`. Note that for this plot, this highlight also included genes enriched in the input samples.
 
 ```{r ggplotexps, out.width='100%', dev="png"}
-volcano_density_plot(stats=all.results[[1]]$stats.eset, title="Stau2 whole RNA pulldown vs. input", lfc.cutoff = lfc.cutoff, pval.cutoff = pvalue.cutoff, shade.colour=highlight.color )          
-volcano_density_plot(stats=all.results[[3]]$stats.eset, title="WT1 whole RNA pulldown vs. input" , lfc.cutoff = lfc.cutoff, pval.cutoff = pvalue.cutoff, shade.colour=highlight.color)          
+vcplot(stats=all.results[[1]]$stats.eset, title="Stau2 whole RNA pulldown vs. input", lfc.cutoff = lfc.cutoff, pval.cutoff = pvalue.cutoff, shade.colour=highlight.color )          
+vcplot(stats=all.results[[3]]$stats.eset, title="WT1 whole RNA pulldown vs. input" , lfc.cutoff = lfc.cutoff, pval.cutoff = pvalue.cutoff, shade.colour=highlight.color)          
 ```
 
 Using these pvalue and log2 fold change cutoffs we can identify which genes are showing enrichment in the two pulldowns. The cutoffs I have picked here (pvalue<`r pvalue.cutoff` and log2foldchange>`r lfc.cutoff`) are within accepted range, if a bit stringent, but are arbitrary.