diff segmentation.R @ 5:b7f3854e08f8

Adding all tools

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/segmentation.R	Wed Apr 29 09:55:41 2015 +0200
@@ -0,0 +1,81 @@
+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, "mpagenomics",userId)
+setwd(workdir)
+
+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
+	    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)
+sink(output)
+print(format(finalResult))
+sink()
+#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)
+				
+			}
+			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)
+	sink(output)
+	print(format(finalResult))
+	sink()
+	#write.table(finalResult,output,row.names = FALSE, quote=FALSE, sep = "\t")
+}
+