view segmentation.R @ 6:7076911e5c64 draft default tip

"planemo upload for repository https://github.com/sblanck/MPAgenomics4Galaxy/tree/master/mpagenomics_wrappers commit a644ed69951bcc1ac46426c5e6c9a0af1003a9a8-dirty"

#!/usr/bin/env Rscript
# setup R error handling to go to stderr
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")

library("optparse")

##### Read options
option_list=list(
		make_option("--input",type="character",default=NULL, dest="input"),
		make_option("--output",type="character",default=NULL, dest="output"),
		make_option("--new_file_path",type="character",default=NULL, dest="new_file_path"),
		make_option("--nbcall",type="character",default=NULL, dest="nbcall"),
		make_option("--outputgraph",type="character",default=NULL, dest="outputgraph"),
		make_option("--graph",type="character",default=NULL, dest="graph"),
		make_option("--signalType",type="character",default=NULL, dest="signalType"),
		make_option("--cellularity",type="character",default=NULL, dest="cellularity"),
		make_option("--outputlog",type="character",default=NULL, dest="outputlog"),
		make_option("--log",type="character",default=NULL, dest="log"),
		make_option("--user_id",type="character",default=NULL, dest="userid"),
		make_option("--method",type="character",default=NULL, dest="method")
);

opt_parser = OptionParser(option_list=option_list);
opt = parse_args(opt_parser);

if(is.null(opt$input)){
	print_help(opt_parser)
	stop("input required.", call.=FALSE)
}

#loading libraries

input=opt$input
output=opt$output
tmp_dir=opt$new_file_path
nbcall=as.numeric(opt$nbcall)
outputgraph=type.convert(opt$outputgraph)
cellularity=as.numeric(opt$cellularity)
userId=opt$userid
method=opt$method
log=opt$log
outputlog=opt$outputlog
graph=opt$graph
signalType=opt$signalType

#args<-commandArgs(TRUE)
#
#input=args[1]
#tmp_dir=args[2]
#nbcall=as.numeric(args[3])
#cellularity=as.numeric(args[4]) 
#output=args[5]
#method=args[6]
#userId=args[7]
#signalType=args[8]

library(MPAgenomics)
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")
	sink(sinklog, type = "message")
} 

CN=read.table(input,header=TRUE)
uniqueChr=unique(CN$chromosome)
drops=c("chromosome","position","probeName")
CNsignal=CN[,!(names(CN)%in% drops),drop=FALSE]

samples=names(CNsignal)

if (signalType=="CN")
{

result=data.frame(sampleNames=character(0),chrom=character(0),chromStart=numeric(0),chromEnd=numeric(0),probes=numeric(0),means=numeric(0),calls=character(0),stringsAsFactors=FALSE)

for (chr in uniqueChr)
{
currentSubset=subset(CN, chromosome==chr)
currentPositions=currentSubset["position"]
for (sample in samples) 
  {
    currentSignal=currentSubset[sample]
	if (length(which(!is.na(unlist(currentSignal))))>1)
	{
		currentSeg=segmentation(signal=unlist(currentSignal),position=unlist(currentPositions),method=method)
    	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)
	}
  }
}
finalResult=data.frame(sampleNames=result["sampleNames"],chrom=result["chrom"],chromStart=result["chromStart"],chromEnd=result["chromEnd"],probes=result["probes"],means=result["means"],calls=result["calls"],stringsAsFactors=FALSE)

write.table(finalResult,output,row.names = FALSE, quote=FALSE, sep = "\t")
} else {
	result=data.frame(sampleNames=character(0),chrom=character(0),start=numeric(0),end=numeric(0),points=numeric(0),means=numeric(0),stringsAsFactors=FALSE)
	
	for (chr in uniqueChr)
	{
		cat(paste0("chromosome ",chr,"\n"))
		currentSubset=subset(CN, chromosome==chr)
		currentPositions=currentSubset["position"]
		for (sample in samples) 
		{
			cat(paste0("  sample ",sample,"..."))
			currentSignal=currentSubset[sample]
			if (length(which(!is.na(unlist(currentSignal))))>1)
			{
				currentSeg=segmentation(signal=unlist(currentSignal),position=unlist(currentPositions),method=method)
				currentResult=currentSeg$segment
				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){
	library(zip)
        files2zip <- dir(pattern=".png")
        zipr(graph, files = files2zip)

}

if (outputlog){
	sink(type="output")
	sink(type="message")
	close(sinklog)
}