Mercurial > repos > sblanck > mpagenomics

comparison segmentation.R @ 0:4d539083cf7f 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 689d0d8dc899a683ee18700ef385753559850233-dirty
author sblanck
date Tue, 12 May 2020 10:40:36 -0400
parents
children 3fcbb8030fcc
comparison
equal deleted inserted replaced
-1:000000000000 0:4d539083cf7f
1 #!/usr/bin/env Rscript
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 ) } )
4
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")
7
8 library("optparse")
9
10 ##### Read options
11 option_list=list(
12 make_option("--input",type="character",default=NULL, dest="input"),
13 make_option("--output",type="character",default=NULL, dest="output"),
14 make_option("--new_file_path",type="character",default=NULL, dest="new_file_path"),
15 make_option("--nbcall",type="character",default=NULL, dest="nbcall"),
16 make_option("--outputgraph",type="character",default=NULL, dest="outputgraph"),
17 make_option("--graph",type="character",default=NULL, dest="graph"),
18 make_option("--signalType",type="character",default=NULL, dest="signalType"),
19 make_option("--cellularity",type="character",default=NULL, dest="cellularity"),
20 make_option("--outputlog",type="character",default=NULL, dest="outputlog"),
21 make_option("--log",type="character",default=NULL, dest="log"),
22 make_option("--user_id",type="character",default=NULL, dest="userid"),
23 make_option("--method",type="character",default=NULL, dest="method")
24 );
25
26 opt_parser = OptionParser(option_list=option_list);
27 opt = parse_args(opt_parser);
28
29 if(is.null(opt$input)){
30 print_help(opt_parser)
31 stop("input required.", call.=FALSE)
32 }
33
34 #loading libraries
35
36 input=opt$input
37 output=opt$output
38 tmp_dir=opt$new_file_path
39 nbcall=as.numeric(opt$nbcall)
40 outputgraph=type.convert(opt$outputgraph)
41 cellularity=as.numeric(opt$cellularity)
42 userId=opt$userid
43 method=opt$method
44 log=opt$log
45 outputlog=opt$outputlog
46 graph=opt$graph
47 signalType=opt$signalType
48
49 #args<-commandArgs(TRUE)
50 #
51 #input=args[1]
52 #tmp_dir=args[2]
53 #nbcall=as.numeric(args[3])
54 #cellularity=as.numeric(args[4])
55 #output=args[5]
56 #method=args[6]
57 #userId=args[7]
58 #signalType=args[8]
59
60 library(MPAgenomics)
61 workdir=file.path(tmp_dir,"mpagenomics",userId)
62 setwd(workdir)
63
64 if (outputlog){
65 sinklog <- file(log, open = "wt")
66 sink(sinklog ,type = "output")
67 sink(sinklog, type = "message")
68 }
69
70 CN=read.table(input,header=TRUE)
71 uniqueChr=unique(CN$chromosome)
72 drops=c("chromosome","position","probeName")
73 CNsignal=CN[,!(names(CN)%in% drops),drop=FALSE]
74
75 samples=names(CNsignal)
76
77 if (signalType=="CN")
78 {
79
80 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)
81
82 for (chr in uniqueChr)
83 {
84 currentSubset=subset(CN, chromosome==chr)
85 currentPositions=currentSubset["position"]
86 for (sample in samples)
87 {
88 currentSignal=currentSubset[sample]
89 if (length(which(!is.na(unlist(currentSignal))))>1)
90 {
91 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)
93 currentCall=callingProcess(callobj,nclass=nbcall,cellularity=cellularity,verbose=TRUE)
94 currentResult=currentCall$segment
95 currentResult["sampleNames"]=c(rep(sample,length(currentCall$segment$chrom)))
96 result=rbind(result,currentResult)
97 }
98 }
99 }
100 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)
101
102 write.table(finalResult,output,row.names = FALSE, quote=FALSE, sep = "\t")
103 } else {
104 result=data.frame(sampleNames=character(0),chrom=character(0),start=numeric(0),end=numeric(0),points=numeric(0),means=numeric(0),stringsAsFactors=FALSE)
105
106 for (chr in uniqueChr)
107 {
108 cat(paste0("chromosome ",chr,"\n"))
109 currentSubset=subset(CN, chromosome==chr)
110 currentPositions=currentSubset["position"]
111 for (sample in samples)
112 {
113 cat(paste0(" sample ",sample,"..."))
114 currentSignal=currentSubset[sample]
115 if (length(which(!is.na(unlist(currentSignal))))>1)
116 {
117 currentSeg=segmentation(signal=unlist(currentSignal),position=unlist(currentPositions),method=method)
118 currentResult=currentSeg$segment
119 currentResult["chrom"]=c(rep(chr,length(currentSeg$segment$means)))
120 currentResult["sampleNames"]=c(rep(sample,length(currentSeg$segment$means)))
121 result=rbind(result,currentResult)
122
123 }
124 cat(paste0("OK\n"))
125 }
126 }
127 finalResult=data.frame(sampleNames=result["sampleNames"],chrom=result["chrom"],chromStart=result["start"],chromEnd=result["end"],probes=result["points"],means=result["means"],stringsAsFactors=FALSE)
128 write.table(finalResult,output,row.names = FALSE, quote=FALSE, sep = "\t")
129 }
130
131 if (outputgraph){
132 file.rename(file.path(tmp_dir,"mpagenomics",userId,"Rplots.pdf"), graph)
133 }
134
135 if (outputlog){
136 sink(type="output")
137 sink(type="message")
138 close(sinklog)
139 }