changeset 32:87bf067cfc53 draft

Uploaded

--- a/edgeR_Differential_Gene_Expression.xml	Tue May 20 05:28:50 2014 -0400
+++ b/edgeR_Differential_Gene_Expression.xml	Wed May 21 04:55:20 2014 -0400
@@ -247,34 +247,36 @@
   write.table(file=output_raw_counts,dge\$counts,sep="\t",row.names=TRUE,col.names=NA)
   
   
-  if(output_MAplot != "/dev/null") {
-    write("Creating MA plot...",stdout())
-    
+  if(output_MAplot != "/dev/null" or output_PValue_distribution_plot != "/dev/null") {
     etable <- topTags(lrt, n=nrow(dge))\$table
     etable <- etable[order(etable\$FDR), ]
-    pdf(output_MAplot)
-    with(etable, plot(logCPM, logFC, pch=20, main="edgeR: Fold change vs abundance"))
-    with(subset(etable, FDR < fdr), points(logCPM, logFC, pch=20, col="red"))
-    abline(h=c(-1,1), col="blue")
-    dev.off()
-  }
+    
+    if(output_MAplot != "/dev/null") {
+      write("Creating MA plot...",stdout())
+      pdf(output_MAplot)
+      with(etable, plot(logCPM, logFC, pch=20, main="edgeR: Fold change vs abundance"))
+      with(subset(etable, FDR < fdr), points(logCPM, logFC, pch=20, col="red"))
+      abline(h=c(-1,1), col="blue")
+      dev.off()
+    }
   
-  if(output_PValue_distribution_plot != "/dev/null") {
-    write("Creating P-value distribution plot...",stdout())
-    pdf(output_PValue_distribution_plot)
-    expressed_genes <- subset(etable, PValue < 0.99)
-    h <- hist(expressed_genes\$PValue,breaks=nrow(expressed_genes)/15,main="Binned P-Values (< 0.99)")
-    center <- sum(h\$counts) / length(h\$counts)
-    lines(c(0,1),c(center,center),lty=2,col="red",lwd=2)
-    k <- ksmooth(h\$mid, h\$counts)
-    lines(k\$x,k\$y,col="red",lwd=2)
-    rmsd <- (h\$counts) - center
-    rmsd <- rmsd^2
-    rmsd <- sum(rmsd)
-    rmsd <- sqrt(rmsd)
-    text(0,max(h\$counts),paste("e=",round(rmsd,2),sep=""),pos=4,col="blue")
-    ## change e into epsilon somehow
-    dev.off()
+    if(output_PValue_distribution_plot != "/dev/null") {
+      write("Creating P-value distribution plot...",stdout())
+      pdf(output_PValue_distribution_plot)
+      expressed_genes <- subset(etable, PValue < 0.99)
+      h <- hist(expressed_genes\$PValue,breaks=nrow(expressed_genes)/15,main="Binned P-Values (< 0.99)")
+      center <- sum(h\$counts) / length(h\$counts)
+      lines(c(0,1),c(center,center),lty=2,col="red",lwd=2)
+      k <- ksmooth(h\$mid, h\$counts)
+      lines(k\$x,k\$y,col="red",lwd=2)
+      rmsd <- (h\$counts) - center
+      rmsd <- rmsd^2
+      rmsd <- sum(rmsd)
+      rmsd <- sqrt(rmsd)
+      text(0,max(h\$counts),paste("e=",round(rmsd,2),sep=""),pos=4,col="blue")
+      ## change e into epsilon somehow
+      dev.off()
+    }
   }
   
   ##output_hierarchical_clustering_plot = args[13]