|
0
|
1 if (require("IWTomics",character.only = TRUE,quietly = FALSE)) {
|
|
|
2 args=commandArgs(TRUE)
|
|
|
3
|
|
|
4 # get args names and values
|
|
|
5 args_values=strsplit(args,'=')
|
|
|
6 args_names=unlist(lapply(args_values,function(arg) arg[1]))
|
|
|
7 names(args_values)=args_names
|
|
|
8 args_values=lapply(args_values,function(arg) arg[2])
|
|
|
9 # read filenames
|
|
|
10 adjustedpvalue=args_values$adjustedpvalue
|
|
|
11 iwtomicsrespdf=args_values$iwtomicsrespdf
|
|
|
12 iwtomicssumpdf=args_values$iwtomicssumpdf
|
|
|
13 iwtomicsrdata=args_values$iwtomicsrdata
|
|
|
14 iwtomicstests=args_values$iwtomicstests
|
|
|
15 iwtomicsselectedfeatures=args_values$iwtomicsselectedfeatures
|
|
|
16 test_subset=paste0('c(',strsplit(args_values$test_subset,'\\|')[[1]],')')
|
|
|
17 feature_subset=paste0('c(',strsplit(args_values$feature_subset,'\\|')[[1]],')')
|
|
|
18 # read parameters (from test_subset on)
|
|
|
19 i_scale_subset=which(args_names=='scale_subset')
|
|
|
20 for(i in i_scale_subset:length(args)){
|
|
|
21 eval(parse(text=args[[i]]))
|
|
|
22 }
|
|
|
23
|
|
|
24 # load RData
|
|
|
25 load(iwtomicsrdata)
|
|
|
26 # read testids and featureids and check them
|
|
|
27 unlisted=lapply(seq_along(test_subset),
|
|
|
28 function(i){
|
|
|
29 test_subset_i=eval(parse(text=test_subset[i]))
|
|
|
30 feature_subset_i=eval(parse(text=feature_subset[i]))
|
|
|
31 test_subset_i=rep(test_subset_i,each=length(feature_subset_i))
|
|
|
32 feature_subset_i=rep(feature_subset_i,length.out=length(test_subset_i))
|
|
|
33 scale_subset_i=rep(scale_subset[i],length(test_subset_i))
|
|
|
34 return(list(test_subset=test_subset_i,feature_subset=feature_subset_i,scale_subset=scale_subset_i))
|
|
|
35 })
|
|
|
36 test_subset=unlist(lapply(unlisted,function(l) l$test_subset))
|
|
|
37 feature_subset=unlist(lapply(unlisted,function(l) l$feature_subset))
|
|
|
38 scale_subset=unlist(lapply(unlisted,function(l) l$scale_subset))
|
|
|
39 testids=as.character(read.delim(iwtomicstests,header=FALSE,sep='\t',stringsAsFactors=FALSE))
|
|
|
40 featureids=as.character(read.delim(iwtomicsselectedfeatures,header=FALSE,sep='\t',stringsAsFactors=FALSE))
|
|
|
41 id_features_subset=featureids[feature_subset]
|
|
|
42 if(sum(testids!=paste(testInput(regionsFeatures_test)$id_region1,'vs',testInput(regionsFeatures_test)$id_region2))){
|
|
|
43 quit(save="no", status=10)
|
|
|
44 stop('Wrong test ids')
|
|
|
45 }
|
|
|
46 if(sum(featureids!=idFeatures(regionsFeatures_test))){
|
|
|
47 quit(save="no", status=20)
|
|
|
48 stop('Wrong feature ids')
|
|
|
49 }
|
|
|
50 # retrieve test and features_subset ids
|
|
|
51 id_features_subset=featureids[feature_subset]
|
|
|
52 if(sum(duplicated(paste0(test_subset,id_features_subset)))){
|
|
|
53 quit(save="no", status=30)
|
|
|
54 stop('Two scale thresholds selected for the same test and feature.')
|
|
|
55 }
|
|
|
56 # If scale_subset=0, do not change the threshold
|
|
|
57 default=(scale_subset==0)
|
|
|
58 scale_subset=scale_subset[!default]
|
|
|
59 test_subset=test_subset[!default]
|
|
|
60 id_features_subset=id_features_subset[!default]
|
|
|
61
|
|
|
62 # get scale threshold
|
|
|
63 scale_threshold=lapply(regionsFeatures_test@test$result,
|
|
|
64 function(result) unlist(lapply(result,function(feature) feature$max_scale)))
|
|
|
65 for(i in seq_along(test_subset)){
|
|
|
66 if(scale_threshold[[test_subset[i]]][id_features_subset[i]]<scale_subset[i]){
|
|
|
67 quit(save="no", status=40)
|
|
|
68 stop('Scale threshold too high.')
|
|
|
69 }
|
|
|
70 scale_threshold[[test_subset[i]]][id_features_subset[i]]=scale_subset[i]
|
|
|
71 }
|
|
|
72
|
|
|
73 # create adjustedvalue output
|
|
|
74 pval=adjusted_pval(regionsFeatures_test,scale_threshold=scale_threshold)
|
|
|
75 for(test in seq_along(pval)){
|
|
|
76 for(id_feature in idFeatures(regionsFeatures_test)){
|
|
|
77 write(paste0('Test: ',testids[test],', on feature ',id_feature),
|
|
|
78 file=adjustedpvalue,append=TRUE,sep='\t')
|
|
|
79 pval_i=as.data.frame(t(pval[[test]][[id_feature]]))
|
|
|
80 row.names(pval_i)=paste('Scale',scale_threshold[[test]][[id_feature]])
|
|
|
81 write.table(pval_i,file=adjustedpvalue,append=TRUE,sep='\t',quote=FALSE,row.names=TRUE,col.names=FALSE)
|
|
|
82 write('',file=adjustedpvalue,append=TRUE,sep='\t')
|
|
|
83 }
|
|
|
84 }
|
|
|
85
|
|
|
86
|
|
|
87 # plot test results
|
|
|
88 pdf(iwtomicsrespdf,width=5,height=7)
|
|
|
89 if(plottype=='boxplot'){
|
|
|
90 # fix repeated probs
|
|
|
91 probs=sort(unique(probs))
|
|
|
92 }else{
|
|
|
93 probs=c(0.25,0.5,0.75)
|
|
|
94 }
|
|
|
95 plotTest(regionsFeatures_test,alpha=testalpha,type=plottype,probs=probs,average=average,size=size,scale_threshold=scale_threshold,ask=FALSE)
|
|
|
96 dev.off()
|
|
|
97
|
|
|
98 # plot summary results
|
|
|
99 if(groupby!='none'){
|
|
|
100 tryCatch({
|
|
|
101 pdf(iwtomicssumpdf,width=15,height=10)
|
|
|
102 plotSummary(regionsFeatures_test,alpha=summaryalpha,only_significant=only_significant,groupby=groupby,scale_threshold=scale_threshold,ask=FALSE,append=TRUE)
|
|
|
103 dev.off()
|
|
|
104 }, error = function(err) {
|
|
|
105 if (grepl('selected features with different resolution',err$message)) {
|
|
|
106 quit(save="no", status=50) #error: groupby 'test' but selected features with different resolution.
|
|
|
107 stop(err)
|
|
|
108 }
|
|
|
109 quit(save="no", status=60) #error
|
|
|
110 stop(err)
|
|
|
111 })
|
|
|
112 }
|
|
|
113
|
|
|
114 }else{
|
|
|
115 quit(save="no", status=255)
|
|
|
116 stop("Missing IWTomics package")
|
|
|
117 } |
