|
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 adjustedpvaluematrix=args_values$adjustedpvaluematrix
|
|
|
11 iwtomicsrespdf=args_values$iwtomicsrespdf
|
|
|
12 iwtomicssumpdf=args_values$iwtomicssumpdf
|
|
|
13 regionids=args_values$regionids
|
|
|
14 featureids=args_values$featureids
|
|
|
15 rdatafile=args_values$rdatafile
|
|
|
16 iwtomicsrdata=args_values$iwtomicsrdata
|
|
|
17 iwtomicstests=args_values$iwtomicstests
|
|
|
18 iwtomicsselectedfeatures=args_values$iwtomicsselectedfeatures
|
|
|
19 # read parameters (from region1 on)
|
|
|
20 i_region1=which(args_names=='region1')
|
|
|
21 for(i in i_region1:length(args)){
|
|
|
22 eval(parse(text=args[[i]]))
|
|
|
23 }
|
|
|
24
|
|
|
25 # load RData
|
|
|
26 load(rdatafile)
|
|
|
27 # read regionids and featureids
|
|
|
28 regionids=as.character(read.delim(regionids,header=FALSE,sep='\t',stringsAsFactors=FALSE))
|
|
|
29 featureids=as.character(read.delim(featureids,header=FALSE,sep='\t',stringsAsFactors=FALSE))
|
|
|
30 # retrieve region1, region2 and features_subset ids and check they are in the RData
|
|
|
31 id_region1=regionids[region1]
|
|
|
32 id_region2=regionids[region2]
|
|
|
33 id_features_subset=featureids[features_subset]
|
|
|
34 if(length(setdiff(c(id_region1,id_region2),idRegions(regionsFeatures)))!=0){
|
|
|
35 quit(save="no", status=10)
|
|
|
36 stop('Wrong region ids')
|
|
|
37 }
|
|
|
38 if(length(setdiff(id_features_subset,idFeatures(regionsFeatures)))!=0){
|
|
|
39 quit(save="no", status=20)
|
|
|
40 stop('Wrong feature ids')
|
|
|
41 }
|
|
|
42 if(sum(duplicated(paste0(id_region1,id_region2)))){
|
|
|
43 quit(save="no", status=30)
|
|
|
44 stop('Same test repeated multiple times.')
|
|
|
45 }
|
|
|
46
|
|
|
47 # perform test
|
|
|
48 tryCatch({
|
|
|
49 # fix repeated probs
|
|
|
50 if(statistics=='quantile'){
|
|
|
51 # fix repeated probs
|
|
|
52 testprobs=sort(unique(testprobs))
|
|
|
53 }else{
|
|
|
54 testprobs=0.5
|
|
|
55 }
|
|
|
56 regionsFeatures_test=IWTomicsTest(regionsFeatures,id_region1,id_region2,id_features_subset,
|
|
|
57 statistics=statistics,probs=testprobs,B=B)
|
|
|
58 # create adjustedvaluematrix output
|
|
|
59 for(test in seq_along(id_region1)){
|
|
|
60 for(id_feature in id_features_subset){
|
|
|
61 write(paste0('Test: ',id_region1[test],' vs ',id_region2[test],', on feature ',id_feature),
|
|
|
62 file=adjustedpvaluematrix,append=TRUE,sep='\t')
|
|
|
63 pval=regionsFeatures_test@test$result[[test]][[id_feature]]$adjusted_pval_matrix
|
|
|
64 row.names(pval)=paste('Scale',rev(seq_len(nrow(pval))))
|
|
|
65 write.table(pval,file=adjustedpvaluematrix,append=TRUE,sep='\t',quote=FALSE,row.names=TRUE,col.names=FALSE)
|
|
|
66 write('',file=adjustedpvaluematrix,append=TRUE,sep='\t')
|
|
|
67 }
|
|
|
68 }
|
|
|
69 }, error = function(err) {
|
|
|
70 quit(save="no", status=40) #error testing
|
|
|
71 stop(err)
|
|
|
72 })
|
|
|
73
|
|
|
74 # plot test results
|
|
|
75 pdf(iwtomicsrespdf,width=5,height=7)
|
|
|
76 if(plottype=='boxplot'){
|
|
|
77 # fix repeated probs
|
|
|
78 probs=sort(unique(probs))
|
|
|
79 }else{
|
|
|
80 probs=c(0.25,0.5,0.75)
|
|
|
81 }
|
|
|
82 plotTest(regionsFeatures_test,alpha=testalpha,type=plottype,probs=probs,average=average,size=size,ask=FALSE)
|
|
|
83 dev.off()
|
|
|
84
|
|
|
85 # plot summary results
|
|
|
86 if(groupby!='none'){
|
|
|
87 pdf(iwtomicssumpdf,width=15,height=10)
|
|
|
88 plotSummary(regionsFeatures_test,alpha=summaryalpha,only_significant=only_significant,groupby=groupby,ask=FALSE,append=TRUE)
|
|
|
89 dev.off()
|
|
|
90 }
|
|
|
91
|
|
|
92 # create output
|
|
|
93 write.table(as.data.frame(t(paste(id_region1,'vs',id_region2))),file=iwtomicstests,quote=FALSE,sep='\t',row.names=FALSE,col.names=FALSE)
|
|
|
94 write.table(as.data.frame(t(idFeatures(regionsFeatures_test))),file=iwtomicsselectedfeatures,quote=FALSE,sep='\t',row.names=FALSE,col.names=FALSE)
|
|
|
95 save(regionsFeatures_test,file=iwtomicsrdata)
|
|
|
96 }else{
|
|
|
97 quit(save="no", status=255)
|
|
|
98 stop("Missing IWTomics package")
|
|
|
99 } |
