Mercurial > repos > cristian > rbgoa

comparison gomwu.functions.R @ 1:f7287f82602f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"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))