Mercurial > repos > cristian > rbgoa
comparison gomwu.functions.R @ 1:f7287f82602f draft
"planemo upload commit 486235d6560c9e95bd42152ad19bf7c3941cdc1b"
author | cristian |
---|---|
date | Tue, 19 Apr 2022 08:28:43 +0000 |
parents | 91261b42c07e |
children | 5acf9dfdfa27 |
comparison
equal
deleted
inserted
replaced
0:91261b42c07e | 1:f7287f82602f |
---|---|
1 clusteringGOs=function(gen2go,div,cutHeight) { | 1 clusteringGOs=function(gen2go,div,cutHeight) { |
2 inname=paste("dissim0_",div,"_",gen2go,sep="") | 2 nn <- strsplit(gen2go, "[/.]") |
3 outname=paste("cl_",inname,sep="") | 3 if (length(nn[[1]]) == 3) { |
4 dir <- nn[[1]][1] | |
5 name <- nn[[1]][2] | |
6 ext <- nn[[1]][3] | |
7 } else if (length(nn[[1]]) == 2) { | |
8 dir <- "." | |
9 name <- nn[[1]][1] | |
10 ext <- nn[[1]][2] | |
11 } | |
12 inname=paste(dir,"/","dissim0_",div,"_",name,".",ext,sep="") | |
13 | |
14 outname=paste(dir,"/","cl_dissim0_",div,"_",name,".",ext,sep="") | |
4 if (!file.exists(outname)) { | 15 if (!file.exists(outname)) { |
5 diss=read.table(inname,sep="\t",header=T,check.names=F) | 16 diss=read.table(inname,sep="\t",header=T,check.names=F) |
6 row.names(diss)=names(diss) | 17 row.names(diss)=names(diss) |
7 hc=hclust(as.dist(diss),method="complete") | 18 hc=hclust(as.dist(diss),method="complete") |
8 cc=cutree(hc,h=cutHeight) | 19 cc=cutree(hc,h=cutHeight) |
12 | 23 |
13 | 24 |
14 | 25 |
15 #--------------- | 26 #--------------- |
16 gomwuStats=function(input,goDatabase,goAnnotations, goDivision, scriptdir, Module=FALSE, Alternative="t", adjust.multcomp="BH", clusterCutHeight=0.25,largest=0.1,smallest=5,perlPath="perl", shuffle.reps=20){ | 27 gomwuStats=function(input,goDatabase,goAnnotations, goDivision, scriptdir, Module=FALSE, Alternative="t", adjust.multcomp="BH", clusterCutHeight=0.25,largest=0.1,smallest=5,perlPath="perl", shuffle.reps=20){ |
28 | |
29 nn <- strsplit(input, "[/.]") | |
30 if (length(nn[[1]]) == 3) { | |
31 dir <- nn[[1]][1] | |
32 name <- nn[[1]][2] | |
33 ext <- nn[[1]][3] | |
34 } else if (length(nn[[1]]) == 2) { | |
35 dir <- "." | |
36 name <- nn[[1]][1] | |
37 ext <- nn[[1]][2] | |
38 } | |
17 | 39 |
18 extraOptions=paste("largest=",largest," smallest=",smallest," cutHeight=",clusterCutHeight,sep="") | 40 extraOptions=paste("largest=",largest," smallest=",smallest," cutHeight=",clusterCutHeight,sep="") |
19 if (Module==TRUE) { adjust.multcomp="shuffle" } | 41 if (Module==TRUE) { adjust.multcomp="shuffle" } |
20 gomwu_a = paste(scriptdir, "gomwu_a.pl", sep="/") | 42 gomwu_a = paste(scriptdir, "gomwu_a.pl", sep="/") |
21 gomwu_b = paste(scriptdir, "gomwu_b.pl", sep="/") | 43 gomwu_b = paste(scriptdir, "gomwu_b.pl", sep="/") |
22 system(paste(perlPath, gomwu_a,goDatabase,goAnnotations,input,goDivision,extraOptions)) | 44 system(paste(perlPath, gomwu_a,goDatabase,goAnnotations,input,goDivision,extraOptions)) |
23 clusteringGOs(goAnnotations,goDivision,clusterCutHeight) | 45 clusteringGOs(goAnnotations,goDivision,clusterCutHeight) |
24 system(paste(perlPath, gomwu_b,goAnnotations,input,goDivision)) | 46 system(paste(perlPath, gomwu_b,goAnnotations,input,goDivision)) |
25 | 47 |
26 inname=paste(goDivision,"_",input,sep="") | 48 inname=paste(dir,"/",name,"_",goDivision,".tsv",sep="") |
27 rsq=read.table(inname,sep="\t",header=T) | 49 rsq=read.table(inname,sep="\t",header=T) |
28 rsq$term=as.factor(rsq$term) | 50 rsq$term=as.factor(rsq$term) |
29 | 51 |
30 mwut.t=TRUE | 52 mwut.t=TRUE |
31 if (length(levels(as.factor(rsq$value)))==2) { | 53 if (length(levels(as.factor(rsq$value)))==2) { |
87 fdr=p.adjust(rr$pval,method=adjust.multcomp) | 109 fdr=p.adjust(rr$pval,method=adjust.multcomp) |
88 } | 110 } |
89 | 111 |
90 message(sum(fdr<0.1)," GO terms at 10% FDR") | 112 message(sum(fdr<0.1)," GO terms at 10% FDR") |
91 rr$p.adj=fdr | 113 rr$p.adj=fdr |
92 fname=paste("MWU_",inname,sep="") | 114 fname=paste(dir,"/","MWU_",goDivision,"_",name,".",ext,sep="") |
93 write.table(rr,fname,row.names=F) | 115 write.table(rr,fname,row.names=F) |
94 } | 116 } |
95 | 117 |
96 #--------------------- | 118 #--------------------- |
97 mwuTest=function(gotable,Alternative) { | 119 mwuTest=function(gotable,Alternative) { |
168 | 190 |
169 #------------------------- | 191 #------------------------- |
170 gomwuPlot=function(inFile,goAnnotations,goDivision,level1=0.1,level2=0.05,level3=0.01,absValue=-log(0.05,10),adjusted=TRUE,txtsize=1,font.family="sans",treeHeight=0.5,colors=NULL) { | 192 gomwuPlot=function(inFile,goAnnotations,goDivision,level1=0.1,level2=0.05,level3=0.01,absValue=-log(0.05,10),adjusted=TRUE,txtsize=1,font.family="sans",treeHeight=0.5,colors=NULL) { |
171 require(ape) | 193 require(ape) |
172 | 194 |
173 input=inFile | 195 nn <- strsplit(inFile, "[/.]") |
174 in.mwu=paste("MWU",goDivision,input,sep="_") | 196 if (length(nn[[1]]) == 3) { |
175 in.dissim=paste("dissim",goDivision,input,goAnnotations,sep="_") | 197 dir <- nn[[1]][1] |
198 name <- nn[[1]][2] | |
199 ext <- nn[[1]][3] | |
200 } else if (length(nn[[1]]) == 2) { | |
201 dir <- "." | |
202 name <- nn[[1]][1] | |
203 ext <- nn[[1]][2] | |
204 } | |
205 aa <- strsplit(goAnnotations, "[/.]") | |
206 if (length(aa[[1]]) == 3) { | |
207 adir <- aa[[1]][1] | |
208 aname <- aa[[1]][2] | |
209 aext <- aa[[1]][3] | |
210 } else if (length(aa[[1]]) == 2) { | |
211 adir <- "." | |
212 aname <- aa[[1]][1] | |
213 aext <- aa[[1]][2] | |
214 } | |
215 # input=inFile | |
216 in.mwu=paste(dir,"/",paste("MWU",goDivision,name,sep="_"), ".", ext,sep="") | |
217 in.dissim=paste(dir, "/", paste("dissim",goDivision,name,aname,sep="_"), ".", aext, sep="") | |
176 | 218 |
177 cutoff=-log(level1,10) | 219 cutoff=-log(level1,10) |
178 pv=read.table(in.mwu,header=T) | 220 pv=read.table(in.mwu,header=T) |
179 row.names(pv)=pv$term | 221 row.names(pv)=pv$term |
180 in.raw=paste(goDivision,input,sep="_") | 222 in.raw=paste(dir,"/",paste(name,goDivision,sep="_"), ".tsv", sep="") |
181 rsq=read.table(in.raw,sep="\t",header=T) | 223 rsq=read.table(in.raw,sep="\t",header=T) |
182 rsq$term=as.factor(rsq$term) | 224 rsq$term=as.factor(rsq$term) |
183 | 225 |
184 if (adjusted==TRUE) { pvals=pv$p.adj } else { pvals=pv$pval } | 226 if (adjusted==TRUE) { pvals=pv$p.adj } else { pvals=pv$pval } |
185 heat=data.frame(cbind("pval"=pvals)) | 227 heat=data.frame(cbind("pval"=pvals)) |