annotate loadandplot.R @ 28:97b7f7bd7f43 draft

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