Mercurial > repos > fabio > iwtomics

comparison loadandplot.R @ 18:565118df598a draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded 20170531
author fabio
date Wed, 31 May 2017 11:07:22 -0400
parents 1e677d6b1aaf
children 97b7f7bd7f43
comparison
equal deleted inserted replaced
17:81faa1567178 18:565118df598a
1 if (require("IWTomics",character.only = TRUE,quietly = FALSE)) {
2 require(tools,quietly = FALSE)
3 args=commandArgs(TRUE)
4
5 for (i in seq_along(args)) {
6 message(args[[i]])
7 }
8
9 # get args names and values
10 args_values=strsplit(args,'=')
11 args_names=unlist(lapply(args_values,function(arg) arg[1]))
12 names(args_values)=args_names
13 args_values=lapply(args_values,function(arg) arg[2])
14 # read filenames
15 outrdata=args_values$outrdata
16 outregions=args_values$outregions
17 outfeatures=args_values$outfeatures
18 outpdf=args_values$outpdf
19 regionspaths=unlist(strsplit(args_values$regionspaths,'\\|'))
20 if("regionsheaderfile" %in% args_names){
21 # the file regionsheaderfile must contain as first column the (unique) regionsfilenames,
22 # as second column the corresponding ids and as third column the names
23 tryCatch({
24 regionsheader=read.delim(args_values$regionsheaderfile,header=FALSE,stringsAsFactors=FALSE,row.names=1,sep="\t")
25 regionsfilenames=unlist(strsplit(args_values$regionsfilenames,'\\|'))
26 if(length(setdiff(regionsfilenames,row.names(regionsheader)))) {
27 quit(save="no", status=11)
28 stop('Not all regionsfilenames are present in the first column of regionsheader.')
29 }
30 id_regions=regionsheader[regionsfilenames,1]
31 name_regions=regionsheader[regionsfilenames,2]
32 }, error = function(err) {
33 quit(save="no", status=10) #error on header file
34 stop(err)
35 })
36 }else{
37 eval(parse(text=args[[which(args_names=='regionsgalaxyids')]]))
38 id_regions=paste0('data_',regionsgalaxyids)
39 name_regions=paste0('data_',regionsgalaxyids)
40 }
41 featurespaths=unlist(strsplit(args_values$featurespaths,'\\|'))
42 if("featuresheaderfile" %in% args_names){
43 # the file featuresheaderfile must contain as first column the (unique) featuresfilenames,
44 # as second column the corresponding ids and as third column the names
45 tryCatch({
46 featuresheader=read.delim(args_values$featuresheaderfile,header=FALSE,stringsAsFactors=FALSE,row.names=1,sep="\t")
47 featuresfilenames=unlist(strsplit(args_values$featuresfilenames,'\\|'))
48 if(length(setdiff(featuresfilenames,row.names(featuresheader)))) {
49 quit(save="no", status=21)
50 stop('Not all featuresfilenames are present in the first column of featuresheader.')
51 }
52 id_features=featuresheader[featuresfilenames,1]
53 name_features=featuresheader[featuresfilenames,2]
54 }, error = function(err) {
55 quit(save="no", status=20) #error on header file
56 stop(err)
57 })
58 }else{
59 eval(parse(text=args[[which(args_names=='featuresgalaxyids')]]))
60 id_features=paste0('data_',featuresgalaxyids)
61 name_features=paste0('data_',featuresgalaxyids)
62 }
63 # read parameters (from smoothing on)
64 i_smoothing=which(args_names=='smoothing')
65 for(i in i_smoothing:length(args)){
66 eval(parse(text=args[[i]]))
67 }
68
69 # load data
70 tryCatch({
71 regionsFeatures=IWTomicsData(regionspaths,featurespaths,alignment,
72 id_regions,name_regions,id_features,name_features,start.are.0based=start.are.0based)
73 }, error = function(err) {
74 if(grepl('invalid format',err$message)){
75 quit(save="no", status=31) # error, not enough columns in input file
76 }else if(grepl('duplicated regions',err$message)){
77 quit(save="no", status=32) # error, duplicated regions in region file
78 }else if(grepl('duplicated windows',err$message)){
79 quit(save="no", status=33) # error, duplicated windows in feature file
80 }else if(grepl('overlapping windows',err$message)){
81 quit(save="no", status=34) # error, overlapping windows in feature file
82 }else if(grepl('not all regions in datasets',err$message)){
83 quit(save="no", status=35) # error, windows in feature files do not cover all regions in region files
84 }else if(grepl('ifferent size windows',err$message)){
85 quit(save="no", status=36) # error, all windows in a feature files must have the same size
86 }
87 #error loading data
88
89 stop(err)
90
91 quit(save="no", status=30)
92
93 })
94
95 # smooth data
96 if(smoothing!='no'){
97 tryCatch({
98 if(smoothing=='locpoly'){
99 dist_knots=10
100 }else if(smoothing=='kernel'){
101 degree=3
102 dist_knots=10
103 }else if(smoothing=='splines'){
104 bandwidth=5
105 }
106 if(alignment=='scale'){
107 if(scale==0){
108 regionsFeatures=smooth(regionsFeatures,type=smoothing,fill_gaps=fill_gaps,
109 bandwidth=bandwidth,degree=degree,dist_knots=dist_knots)
110 }else{
111 regionsFeatures=smooth(regionsFeatures,type=smoothing,fill_gaps=fill_gaps,
112 bandwidth=bandwidth,degree=degree,dist_knots=dist_knots,scale_grid=scale)
113 }
114 }else{
115 regionsFeatures=smooth(regionsFeatures,type=smoothing,fill_gaps=fill_gaps,
116 bandwidth=bandwidth,degree=degree,dist_knots=dist_knots)
117 }
118 }, error = function(err) {
119 quit(save="no", status=40) #error on smoothing
120 stop(err)
121 })
122 }
123
124 # plot data
125 pdf(outpdf,width=10,height=8)
126 if(plottype=='boxplot'){
127 # fix repeated probs
128 probs=sort(unique(probs))
129 }else{
130 probs=c(0.25,0.5,0.75)
131 }
132 plot(regionsFeatures,type=plottype,probs=probs,average=average,size=size,ask=FALSE)
133 dev.off()
134
135 # create output
136 #write.table(cbind(unlist(strsplit(args_values$regionsfilenames,'\\|')),idRegions(regionsFeatures),nameRegions(regionsFeatures)),
137 #file=outregions,quote=FALSE,sep='\t',row.names=FALSE,col.names=FALSE)
138 write.table(as.data.frame(t(idRegions(regionsFeatures))),file=outregions,quote=FALSE,sep='\t',row.names=FALSE,col.names=FALSE)
139 #write.table(cbind(unlist(strsplit(args_values$featuresfilenames,'\\|')),idFeatures(regionsFeatures),nameFeatures(regionsFeatures)),
140 #file=outfeatures,quote=FALSE,sep='\t',row.names=FALSE,col.names=FALSE)
141 write.table(as.data.frame(t(idFeatures(regionsFeatures))),file=outfeatures,quote=FALSE,sep='\t',row.names=FALSE,col.names=FALSE)
142 save(regionsFeatures,file=outrdata)
143 }else{
144 quit(save="no", status=255)
145 stop("Missing IWTomics package")
146 }