outputing empty pca plots when batch set to False

Author: Ubuntu

src/preprocess/batch_removal.R

@@ -1,7 +1,10 @@
 suppressMessages(library(sva))
 suppressMessages(library(limma))
 suppressMessages(library(optparse))
-
+suppressMessages(library(ggplot2))
+suppressMessages(library(ggfortify))
+suppressMessages(library(RColorBrewer))
+suppressMessages(library(ggpubr))
 
 
 # make option list and parse command line
@@ -20,6 +23,22 @@ option_list <- list(
               help="Output files [Required]"))
 )
 
+pca_plot <- function(exprTable, annot,title,Batch) {
+  batch_n <- length(unique(as.character(annot[colnames(exprTable),"Batch"])))
+  print(paste("there are ", batch_n, " batches in your data"))
+  df <- cbind.data.frame(t(exprTable),batch = as.character(annot[colnames(exprTable),"Batch"]))
+  pca_plot <- autoplot(prcomp(t(exprTable)), data = df, col = Batch, size = 1, frame = TRUE, frame.type = 'norm')+
+    labs(title=title)+
+    scale_color_manual(values = brewer.pal(name = "Set1", n = 9)[1:batch_n])+
+    theme_bw()+
+    theme(plot.title = element_text(hjust=0.5,vjust = 0.5, margin = margin(l=100,r=50,t=10,b=10),face = "bold", colour = "black"),
+          axis.text.x=element_blank(),
+          axis.text.y=element_text(size=12,face = "bold",hjust=1),
+          axis.title.x = element_text(size = 12, face = "bold"),
+          axis.title.y = element_text(size = 12, face = "bold"))
+  return(pca_plot)
+}
+
 
 opt_parser <- OptionParser(option_list=option_list);
 opts <- parse_args(opt_parser);
@@ -69,19 +88,32 @@ exprZero <- subset(exprZero, !(rownames(exprZero) %in% names(filt_genes)))
 
 
 if(opts$covariates == "False"){
-        print('No batches used !')
-	expr.limma <- expr.dat
-    } else {
-	print('Running limma for batch removal')
+      print('No batches used !')
+      expr.limma <- expr.dat
+      expr.limma = rbind(expr.limma,exprZero)
+      writeDF(ssgsvaFormat(expr.limma),opts$output_after)
+    
+  } else {
+      samples$Batch <- samples[,opts$covariates] 
+      print(samples)
+      
+      print('Running limma for batch removal')
       expr.limma = tryCatch(
       removeBatchEffect(as.matrix(expr.dat),
-                          samples$opts$covariates),
+                          samples$Batch),
       error = function(e){
       print(e)
       })
-    }
-print(head(expr.limma))
-expr.limma = rbind(expr.limma,exprZero)
-writeDF(ssgsvaFormat(expr.limma),opts$output_after)
+      expr.limma = rbind(expr.limma,exprZero)
+      writeDF(ssgsvaFormat(expr.limma),opts$output_after)
+      
+      #png(file = paste(opts$out_before,'.png',sep=''), res = 300, height = 1200, width = 1500)
+      #pca_plot(expr.dat, annot = samples, title = "PCA plot Before Batch Removal",Batch=Batch)
+      #dev.off()
+     
+      #png(file = paste(opt$out_after,'.png',sep=''), res = 300, height = 1200, width = 1500)
+      #pca_plot(expr.limma, annot = samples, title = "PCA plot Before Batch Removal",Batch=Batch)
+      #dev.off()
 
+    }
 

src/preprocess/pca.R

@@ -58,12 +58,23 @@ pca_plot <- function(exprTable, annot,title, Batch) {
 
 
 ###print pca plot
-png(file = opt$out_before, res = 300, height = 1200, width = 1500)
-pca_plot(exprTable = expr_before, annot = annot, title = "PCA plot Before Batch Removal",Batch=Batch)
-dev.off()
 
-png(file = opt$out_after, res = 300, height = 1200, width = 1500)
-pca_plot(exprTable = expr_after, annot = annot, title = "PCA plot After Batch Removal",Batch = Batch)
-dev.off()
+if (opt$column != "False") {
+	png(file = opt$out_before, res = 300, height = 1200, width = 1500)
+	pca_plot(exprTable = expr_before, annot = annot, title = "PCA plot Before Batch Removal",Batch=Batch)
+	dev.off()
 
+	png(file = opt$out_after, res = 300, height = 1200, width = 1500)
+	pca_plot(exprTable = expr_after, annot = annot, title = "PCA plot After Batch Removal",Batch = Batch)
+	dev.off()
+} else
+{
+	png(file = opt$out_before, res = 300, height = 1200, width = 1500)
+	plot.new()
+	dev.off()
+
+	png(file = opt$out_after, res = 300, height = 1200, width = 1500)
+	plot.new()
+	dev.off()
+}
 ###