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