annotate testandplot.R @ 91:2a83c5aafade draft

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