diff segmentation.R @ 4:3fcbb8030fcc draft

"planemo upload for repository https://github.com/sblanck/MPAgenomics4Galaxy/tree/master/mpagenomics_wrappers commit 40eda5ea3551e8b3bae32d0a8f405fe90ef22646-dirty"
author sblanck
date Mon, 12 Apr 2021 14:47:09 +0000
parents 4d539083cf7f
children
line wrap: on
line diff
--- a/segmentation.R	Tue Jun 16 04:34:09 2020 -0400
+++ b/segmentation.R	Mon Apr 12 14:47:09 2021 +0000
@@ -3,7 +3,7 @@
 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
 
 # we need that to not crash galaxy with an UTF8 error on German LC settings.
-loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
+# loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
 
 library("optparse")
 
@@ -58,9 +58,13 @@
 #signalType=args[8]
 
 library(MPAgenomics)
-workdir=file.path(tmp_dir,"mpagenomics",userId)
+workdir=file.path(tmp_dir)
+if (!dir.exists(workdir))
+  dir.create(workdir, showWarnings = TRUE, recursive = TRUE)
 setwd(workdir)
 
+workdir
+
 if (outputlog){
 	sinklog <- file(log, open = "wt")
 	sink(sinklog ,type = "output")
@@ -92,6 +96,21 @@
     	callobj= callingObject(copynumber=currentSeg$signal, segmented=currentSeg$segmented,chromosome=rep(chr,length(currentSeg$signal)), position=currentSeg$startPos,sampleNames=sample)
 	    currentCall=callingProcess(callobj,nclass=nbcall,cellularity=cellularity,verbose=TRUE)
 	    currentResult=currentCall$segment
+	    if(outputgraph)
+  		{    
+    	
+		currentPos=unlist(currentPositions)
+		figName <- sprintf("%s,%s", sample, chr);
+          	pathname <- file.path(sprintf("%s.png", figName));
+          	png(filename = pathname, width = 1280, height = 480)
+          	plot(NA,xlim=c(min(currentPos),max(currentPos)), ylim=c(0,6),xlab="Position", main=figName,ylab="CN", pch=".")
+		points(currentPos, unlist(currentSignal), pch=".");
+		for(i in 1:nrow(currentResult))
+            		lines(c(currentResult$chromStart[i],currentResult$chromEnd[i]),rep(currentResult$means[i],2),col="red",lwd=3)
+          	dev.off()	
+			
+ 		 }	
+	    
 	    currentResult["sampleNames"]=c(rep(sample,length(currentCall$segment$chrom)))
 	    result=rbind(result,currentResult)
 	}
@@ -119,17 +138,38 @@
 				currentResult["chrom"]=c(rep(chr,length(currentSeg$segment$means)))
 				currentResult["sampleNames"]=c(rep(sample,length(currentSeg$segment$means)))
 				result=rbind(result,currentResult)
-				
+				if(outputgraph)
+	                	{
+
+                		currentPos=unlist(currentPositions) 
+                		figName <- sprintf("%s,%s", sample, chr);
+                		pathname <- file.path(sprintf("%s.png", figName));
+                		png(filename = pathname, width = 1280, height = 480)
+                		plot(NA,xlim=c(min(currentPos),max(currentPos)), ylim=c(0,1),xlab="Position", main=figName,ylab="CN", pch=".")
+                		points(currentPos, unlist(currentSignal), pch=".");
+				print(currentResult)
+				for(i in 1:nrow(currentResult))
+                       			lines(c(currentResult$start[i],currentResult$end[i]),rep(currentResult$means[i],2),col="red",lwd=3)
+                		dev.off()
+                        
+                 	}
+
+	
+	
 			}
 			cat(paste0("OK\n"))
 		}
 	}
 	finalResult=data.frame(sampleNames=result["sampleNames"],chrom=result["chrom"],chromStart=result["start"],chromEnd=result["end"],probes=result["points"],means=result["means"],stringsAsFactors=FALSE)
+	colnames(finalResult)=c("sampleNames","chrom","chromStart","chromEnd","probes","means")
 	write.table(finalResult,output,row.names = FALSE, quote=FALSE, sep = "\t")
 }
 
 if (outputgraph){
-	file.rename(file.path(tmp_dir,"mpagenomics",userId,"Rplots.pdf"), graph)
+	library(zip)
+        files2zip <- dir(pattern=".png")
+        zipr(graph, files = files2zip)
+
 }
 
 if (outputlog){