Mercurial > repos > sblanck > mpagenomics

comparison segmentation.R @ 4:3fcbb8030fcc draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"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
comparison
equal deleted inserted replaced
3:94fc6ed13946 4:3fcbb8030fcc
1 #!/usr/bin/env Rscript 1 #!/usr/bin/env Rscript
2 # setup R error handling to go to stderr 2 # setup R error handling to go to stderr
3 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } ) 3 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
4 4
5 # we need that to not crash galaxy with an UTF8 error on German LC settings. 5 # we need that to not crash galaxy with an UTF8 error on German LC settings.
6 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") 6 # loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
7 7
8 library("optparse") 8 library("optparse")
9 9
10 ##### Read options 10 ##### Read options
11 option_list=list( 11 option_list=list(
56 #method=args[6] 56 #method=args[6]
57 #userId=args[7] 57 #userId=args[7]
58 #signalType=args[8] 58 #signalType=args[8]
59 59
60 library(MPAgenomics) 60 library(MPAgenomics)
61 workdir=file.path(tmp_dir,"mpagenomics",userId) 61 workdir=file.path(tmp_dir)
62 if (!dir.exists(workdir))
63 dir.create(workdir, showWarnings = TRUE, recursive = TRUE)
62 setwd(workdir) 64 setwd(workdir)
65
66 workdir
63 67
64 if (outputlog){ 68 if (outputlog){
65 sinklog <- file(log, open = "wt") 69 sinklog <- file(log, open = "wt")
66 sink(sinklog ,type = "output") 70 sink(sinklog ,type = "output")
67 sink(sinklog, type = "message") 71 sink(sinklog, type = "message")
90 { 94 {
91 currentSeg=segmentation(signal=unlist(currentSignal),position=unlist(currentPositions),method=method) 95 currentSeg=segmentation(signal=unlist(currentSignal),position=unlist(currentPositions),method=method)
92 callobj= callingObject(copynumber=currentSeg$signal, segmented=currentSeg$segmented,chromosome=rep(chr,length(currentSeg$signal)), position=currentSeg$startPos,sampleNames=sample) 96 callobj= callingObject(copynumber=currentSeg$signal, segmented=currentSeg$segmented,chromosome=rep(chr,length(currentSeg$signal)), position=currentSeg$startPos,sampleNames=sample)
93 currentCall=callingProcess(callobj,nclass=nbcall,cellularity=cellularity,verbose=TRUE) 97 currentCall=callingProcess(callobj,nclass=nbcall,cellularity=cellularity,verbose=TRUE)
94 currentResult=currentCall$segment 98 currentResult=currentCall$segment
99 if(outputgraph)
100 {
101
102 currentPos=unlist(currentPositions)
103 figName <- sprintf("%s,%s", sample, chr);
104 pathname <- file.path(sprintf("%s.png", figName));
105 png(filename = pathname, width = 1280, height = 480)
106 plot(NA,xlim=c(min(currentPos),max(currentPos)), ylim=c(0,6),xlab="Position", main=figName,ylab="CN", pch=".")
107 points(currentPos, unlist(currentSignal), pch=".");
108 for(i in 1:nrow(currentResult))
109 lines(c(currentResult$chromStart[i],currentResult$chromEnd[i]),rep(currentResult$means[i],2),col="red",lwd=3)
110 dev.off()
111
112 }
113
95 currentResult["sampleNames"]=c(rep(sample,length(currentCall$segment$chrom))) 114 currentResult["sampleNames"]=c(rep(sample,length(currentCall$segment$chrom)))
96 result=rbind(result,currentResult) 115 result=rbind(result,currentResult)
97 } 116 }
98 } 117 }
99 } 118 }
117 currentSeg=segmentation(signal=unlist(currentSignal),position=unlist(currentPositions),method=method) 136 currentSeg=segmentation(signal=unlist(currentSignal),position=unlist(currentPositions),method=method)
118 currentResult=currentSeg$segment 137 currentResult=currentSeg$segment
119 currentResult["chrom"]=c(rep(chr,length(currentSeg$segment$means))) 138 currentResult["chrom"]=c(rep(chr,length(currentSeg$segment$means)))
120 currentResult["sampleNames"]=c(rep(sample,length(currentSeg$segment$means))) 139 currentResult["sampleNames"]=c(rep(sample,length(currentSeg$segment$means)))
121 result=rbind(result,currentResult) 140 result=rbind(result,currentResult)
122 141 if(outputgraph)
142 {
143
144 currentPos=unlist(currentPositions)
145 figName <- sprintf("%s,%s", sample, chr);
146 pathname <- file.path(sprintf("%s.png", figName));
147 png(filename = pathname, width = 1280, height = 480)
148 plot(NA,xlim=c(min(currentPos),max(currentPos)), ylim=c(0,1),xlab="Position", main=figName,ylab="CN", pch=".")
149 points(currentPos, unlist(currentSignal), pch=".");
150 print(currentResult)
151 for(i in 1:nrow(currentResult))
152 lines(c(currentResult$start[i],currentResult$end[i]),rep(currentResult$means[i],2),col="red",lwd=3)
153 dev.off()
154
155 }
156
157
158
123 } 159 }
124 cat(paste0("OK\n")) 160 cat(paste0("OK\n"))
125 } 161 }
126 } 162 }
127 finalResult=data.frame(sampleNames=result["sampleNames"],chrom=result["chrom"],chromStart=result["start"],chromEnd=result["end"],probes=result["points"],means=result["means"],stringsAsFactors=FALSE) 163 finalResult=data.frame(sampleNames=result["sampleNames"],chrom=result["chrom"],chromStart=result["start"],chromEnd=result["end"],probes=result["points"],means=result["means"],stringsAsFactors=FALSE)
164 colnames(finalResult)=c("sampleNames","chrom","chromStart","chromEnd","probes","means")
128 write.table(finalResult,output,row.names = FALSE, quote=FALSE, sep = "\t") 165 write.table(finalResult,output,row.names = FALSE, quote=FALSE, sep = "\t")
129 } 166 }
130 167
131 if (outputgraph){ 168 if (outputgraph){
132 file.rename(file.path(tmp_dir,"mpagenomics",userId,"Rplots.pdf"), graph) 169 library(zip)
170 files2zip <- dir(pattern=".png")
171 zipr(graph, files = files2zip)
172
133 } 173 }
134 174
135 if (outputlog){ 175 if (outputlog){
136 sink(type="output") 176 sink(type="output")
137 sink(type="message") 177 sink(type="message")