Mercurial > repos > cristian > rbgoa

comparison gomwu.functions.R @ 2:5acf9dfdfa27 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 66a856bcce69986d9a6f1a39820dd9b3f4f6b0db
author cristian
date Wed, 09 Nov 2022 08:57:54 +0000
parents f7287f82602f
children
comparison
equal deleted inserted replaced
1:f7287f82602f 2:5acf9dfdfa27
1 clusteringGOs=function(gen2go,div,cutHeight) { 1 clusteringGOs <- function(gen2go, div, cutHeight) {
2 nn <- strsplit(gen2go, "[/.]") 2 nn <- strsplit(gen2go, "[/.]")
3 if (length(nn[[1]]) == 3) { 3 if (length(nn[[1]]) == 3) {
4 dir <- nn[[1]][1] 4 dir <- nn[[1]][1]
5 name <- nn[[1]][2] 5 name <- nn[[1]][2]
6 ext <- nn[[1]][3] 6 ext <- nn[[1]][3]
7 } else if (length(nn[[1]]) == 2) { 7 } else if (length(nn[[1]]) == 2) {
8 dir <- "." 8 dir <- "."
9 name <- nn[[1]][1] 9 name <- nn[[1]][1]
10 ext <- nn[[1]][2] 10 ext <- nn[[1]][2]
11 } 11 }
12 inname=paste(dir,"/","dissim0_",div,"_",name,".",ext,sep="") 12 inname <- paste(dir, "/", "dissim0_", div, "_", name, ".", ext, sep = "")
13 13
14 outname=paste(dir,"/","cl_dissim0_",div,"_",name,".",ext,sep="") 14 outname <- paste(dir, "/", "cl_dissim0_", div, "_", name, ".", ext, sep = "")
15 if (!file.exists(outname)) { 15 if (!file.exists(outname)) {
16 diss=read.table(inname,sep="\t",header=T,check.names=F) 16 diss <- read.table(inname, sep = "\t", header = T, check.names = F)
17 row.names(diss)=names(diss) 17 row.names(diss) <- names(diss)
18 hc=hclust(as.dist(diss),method="complete") 18 hc <- hclust(as.dist(diss), method = "complete")
19 cc=cutree(hc,h=cutHeight) 19 cc <- cutree(hc, h = cutHeight)
20 write.csv(cc,file=outname,quote=F) 20 write.csv(cc, file = outname, quote = F)
21 } 21 }
22 } 22 }
23 23
24 24
25 25
26 #--------------- 26 #---------------
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){ 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, "[/.]") 28 nn <- strsplit(input, "[/.]")
30 if (length(nn[[1]]) == 3) { 29 if (length(nn[[1]]) == 3) {
31 dir <- nn[[1]][1] 30 dir <- nn[[1]][1]
32 name <- nn[[1]][2] 31 name <- nn[[1]][2]
33 ext <- nn[[1]][3] 32 ext <- nn[[1]][3]
35 dir <- "." 34 dir <- "."
36 name <- nn[[1]][1] 35 name <- nn[[1]][1]
37 ext <- nn[[1]][2] 36 ext <- nn[[1]][2]
38 } 37 }
39 38
40 extraOptions=paste("largest=",largest," smallest=",smallest," cutHeight=",clusterCutHeight,sep="") 39 extraOptions <- paste("largest=", largest, " smallest=", smallest, " cutHeight=", clusterCutHeight, sep = "")
41 if (Module==TRUE) { adjust.multcomp="shuffle" } 40 if (Module == TRUE) {
42 gomwu_a = paste(scriptdir, "gomwu_a.pl", sep="/") 41 adjust.multcomp <- "shuffle"
43 gomwu_b = paste(scriptdir, "gomwu_b.pl", sep="/") 42 }
44 system(paste(perlPath, gomwu_a,goDatabase,goAnnotations,input,goDivision,extraOptions)) 43 gomwu_a <- paste(scriptdir, "gomwu_a.pl", sep = "/")
45 clusteringGOs(goAnnotations,goDivision,clusterCutHeight) 44 gomwu_b <- paste(scriptdir, "gomwu_b.pl", sep = "/")
46 system(paste(perlPath, gomwu_b,goAnnotations,input,goDivision)) 45 system(paste(perlPath, gomwu_a, goDatabase, goAnnotations, input, goDivision, extraOptions))
47 46 clusteringGOs(goAnnotations, goDivision, clusterCutHeight)
48 inname=paste(dir,"/",name,"_",goDivision,".tsv",sep="") 47 system(paste(perlPath, gomwu_b, goAnnotations, input, goDivision))
49 rsq=read.table(inname,sep="\t",header=T) 48
50 rsq$term=as.factor(rsq$term) 49 inname <- paste(dir, "/", name, "_", goDivision, ".tsv", sep = "")
51 50 rsq <- read.table(inname, sep = "\t", header = T)
52 mwut.t=TRUE 51 rsq$term <- as.factor(rsq$term)
53 if (length(levels(as.factor(rsq$value)))==2) { 52
54 cat("Binary classification detected; will perform Fisher's test\n"); 53 mwut.t <- TRUE
55 mwut.t=F 54 if (length(levels(as.factor(rsq$value))) == 2) {
56 rr=fisherTest(rsq) 55 cat("Binary classification detected; will perform Fisher's test\n")
57 } else { 56 mwut.t <- F
58 if (Module==TRUE) { 57 rr <- fisherTest(rsq)
59 rsq.f=rsq 58 } else {
60 rsq.f$value=as.numeric(rsq.f$value>0) 59 if (Module == TRUE) {
61 rf=fisherTest(rsq.f) 60 rsq.f <- rsq
62 rsq.m=rsq[rsq$value>0,] 61 rsq.f$value <- as.numeric(rsq.f$value > 0)
63 rsq.m$term=factor(rsq.m$term,levels=unique(rsq.m$term)) 62 rf <- fisherTest(rsq.f)
64 rrm=mwuTest(rsq.m,"g") 63 rsq.m <- rsq[rsq$value > 0, ]
65 rr0=rf[rf$term %in% rrm$term,] 64 rsq.m$term <- factor(rsq.m$term, levels = unique(rsq.m$term))
66 rr1=rf[!(rf$term %in% rrm$term),] 65 rrm <- mwuTest(rsq.m, "g")
67 rr0=rr0[order(rr0$term),] 66 rr0 <- rf[rf$term %in% rrm$term, ]
68 rrm=rrm[order(rrm$term),] 67 rr1 <- rf[!(rf$term %in% rrm$term), ]
69 rr0$pval=rr0$pval*rrm$pval 68 rr0 <- rr0[order(rr0$term), ]
70 rr=rbind(rr0,rr1) 69 rrm <- rrm[order(rrm$term), ]
71 } else { 70 rr0$pval <- rr0$pval * rrm$pval
72 cat("Continuous measure of interest: will perform MWU test\n"); 71 rr <- rbind(rr0, rr1)
73 rr=mwuTest(rsq,Alternative) 72 } else {
74 } 73 cat("Continuous measure of interest: will perform MWU test\n")
75 } 74 rr <- mwuTest(rsq, Alternative)
76 75 }
77 if (adjust.multcomp=="shuffle"){ 76 }
78 message("shuffling values to calculate FDR, ",shuffle.reps," reps") 77
79 reps=shuffle.reps 78 if (adjust.multcomp == "shuffle") {
80 spv=c() 79 message("shuffling values to calculate FDR, ", shuffle.reps, " reps")
81 for (i in 1:reps) { 80 reps <- shuffle.reps
82 message("replicate ",i) 81 spv <- c()
83 rsqq=rsq 82 for (i in 1:reps) {
84 rsqq$value=sample(rsq$value) 83 message("replicate ", i)
85 if (Module==TRUE) { 84 rsqq <- rsq
86 rsq.f=rsqq 85 rsqq$value <- sample(rsq$value)
87 rsq.f$value=as.numeric(rsq.f$value>0) 86 if (Module == TRUE) {
88 rf=fisherTest(rsq.f) 87 rsq.f <- rsqq
89 rsq.m=rsqq[rsqq$value>0,] 88 rsq.f$value <- as.numeric(rsq.f$value > 0)
90 rsq.m$term=factor(rsq.m$term,levels=unique(rsq.m$term)) 89 rf <- fisherTest(rsq.f)
91 rrm=mwuTest(rsq.m,"g") 90 rsq.m <- rsqq[rsqq$value > 0, ]
92 rr0=rf[rf$term %in% rrm$term,] 91 rsq.m$term <- factor(rsq.m$term, levels = unique(rsq.m$term))
93 rr1=rf[!(rf$term %in% rrm$term),] 92 rrm <- mwuTest(rsq.m, "g")
94 rr0=rr0[order(rr0$term),] 93 rr0 <- rf[rf$term %in% rrm$term, ]
95 rrm=rrm[order(rrm$term),] 94 rr1 <- rf[!(rf$term %in% rrm$term), ]
96 rr0$pval=rr0$pval*rrm$pval 95 rr0 <- rr0[order(rr0$term), ]
97 rs=rbind(rr0,rr1) 96 rrm <- rrm[order(rrm$term), ]
98 } else { 97 rr0$pval <- rr0$pval * rrm$pval
99 if (mwut.t==TRUE) { rs=mwuTest(rsqq,Alternative) } else { rs=fisherTest(rsqq) } 98 rs <- rbind(rr0, rr1)
100 } 99 } else {
101 spv=append(spv,rs$pval) 100 if (mwut.t == TRUE) {
102 } 101 rs <- mwuTest(rsqq, Alternative)
103 fdr=c() 102 } else {
104 for (p in rr$pval){ 103 rs <- fisherTest(rsqq)
105 fdr=append(fdr,(sum(spv<=p)/reps)/sum(rr$pval<=p)) 104 }
106 } 105 }
107 fdr[fdr>1]=1 106 spv <- append(spv, rs$pval)
108 } else { 107 }
109 fdr=p.adjust(rr$pval,method=adjust.multcomp) 108 fdr <- c()
110 } 109 for (p in rr$pval) {
111 110 fdr <- append(fdr, (sum(spv <= p) / reps) / sum(rr$pval <= p))
112 message(sum(fdr<0.1)," GO terms at 10% FDR") 111 }
113 rr$p.adj=fdr 112 fdr[fdr > 1] <- 1
114 fname=paste(dir,"/","MWU_",goDivision,"_",name,".",ext,sep="") 113 } else {
115 write.table(rr,fname,row.names=F) 114 fdr <- p.adjust(rr$pval, method = adjust.multcomp)
115 }
116
117 message(sum(fdr < 0.1), " GO terms at 10% FDR")
118 rr$p.adj <- fdr
119 fname <- paste(dir, "/", "MWU_", goDivision, "_", name, ".", ext, sep = "")
120 write.table(rr, fname, sep = "\t", row.names = F)
116 } 121 }
117 122
118 #--------------------- 123 #---------------------
119 mwuTest=function(gotable,Alternative) { 124 mwuTest <- function(gotable, Alternative) {
120 gos=gotable 125 gos <- gotable
121 terms=levels(gos$term) 126 terms <- levels(gos$term)
122 gos$seq=as.character(gos$seq) 127 gos$seq <- as.character(gos$seq)
123 nrg=gos[!duplicated(gos$seq),5] 128 nrg <- gos[!duplicated(gos$seq), 5]
124 names(nrg)=gos[!duplicated(gos$seq),4] 129 names(nrg) <- gos[!duplicated(gos$seq), 4]
125 # nrg=nrg+rnorm(nrg,sd=0.01) # to be able to do exact wilcox test 130 # nrg=nrg+rnorm(nrg,sd=0.01) # to be able to do exact wilcox test
126 rnk=rank(nrg) 131 rnk <- rank(nrg)
127 names(rnk)=names(nrg) 132 names(rnk) <- names(nrg)
128 pvals=c();drs=c();nams=c();levs=c();nseqs=c() 133 pvals <- c()
129 for (t in terms){ 134 drs <- c()
130 got=gos[gos$term==t,] 135 nams <- c()
131 got=got[!duplicated(got$seq),] 136 levs <- c()
132 ngot=gos[gos$term!=t,] 137 nseqs <- c()
133 ngot=ngot[!duplicated(ngot$seq),] 138 for (t in terms) {
134 ngot=ngot[!(ngot$seq %in% got$seq),] 139 got <- gos[gos$term == t, ]
135 sgo.yes=got$seq 140 got <- got[!duplicated(got$seq), ]
136 n1=length(sgo.yes) 141 ngot <- gos[gos$term != t, ]
137 sgo.no=ngot$seq 142 ngot <- ngot[!duplicated(ngot$seq), ]
138 n2=length(sgo.no) 143 ngot <- ngot[!(ngot$seq %in% got$seq), ]
139 #message(t," wilcox:",n1," ",n2) 144 sgo.yes <- got$seq
140 wi=wilcox.test(nrg[sgo.yes],nrg[sgo.no],alternative=Alternative) # removed correct=FALSE 145 n1 <- length(sgo.yes)
141 r1=sum(rnk[sgo.yes])/n1 146 sgo.no <- ngot$seq
142 r0=sum(rnk[sgo.no])/n2 147 n2 <- length(sgo.no)
143 dr=r1-r0 148 # message(t," wilcox:",n1," ",n2)
144 drs=append(drs,round(dr,0)) 149 wi <- wilcox.test(nrg[sgo.yes], nrg[sgo.no], alternative = Alternative) # removed correct=FALSE
145 levs=append(levs,got$lev[1]) 150 r1 <- sum(rnk[sgo.yes]) / n1
146 nams=append(nams,as.character(got$name[1])) 151 r0 <- sum(rnk[sgo.no]) / n2
147 pvals=append(pvals,wi$p.value) 152 dr <- r1 - r0
148 nseqs=append(nseqs,n1) 153 drs <- append(drs, round(dr, 0))
149 } 154 levs <- append(levs, got$lev[1])
150 res=data.frame(cbind("delta.rank"=drs,"pval"=pvals,"level"=levs,nseqs)) 155 nams <- append(nams, as.character(got$name[1]))
151 res=cbind(res,"term"=as.character(terms),"name"=nams) 156 pvals <- append(pvals, wi$p.value)
152 res$pval=as.numeric(as.character(res$pval)) 157 nseqs <- append(nseqs, n1)
153 res$delta.rank=as.numeric(as.character(res$delta.rank)) 158 }
154 res$level=as.numeric(as.character(res$level)) 159 res <- data.frame(cbind("delta.rank" = drs, "pval" = pvals, "level" = levs, nseqs))
155 res$nseqs=as.numeric(as.character(res$nseqs)) 160 res <- cbind(res, "term" = as.character(terms), "name" = nams)
156 return(res) 161 res$pval <- as.numeric(as.character(res$pval))
162 res$delta.rank <- as.numeric(as.character(res$delta.rank))
163 res$level <- as.numeric(as.character(res$level))
164 res$nseqs <- as.numeric(as.character(res$nseqs))
165 return(res)
157 } 166 }
158 #------------------------ 167 #------------------------
159 fisherTest=function(gotable) { 168 fisherTest <- function(gotable) {
160 gos=gotable 169 gos <- gotable
161 terms=levels(gos$term) 170 terms <- levels(gos$term)
162 gos$seq=as.character(gos$seq) 171 gos$seq <- as.character(gos$seq)
163 pft=c();nam=c();lev=c();nseqs=c() 172 pft <- c()
164 for (t in terms) { 173 nam <- c()
165 got=gos[gos$term==t,] 174 lev <- c()
166 got=got[!duplicated(got$seq),] 175 nseqs <- c()
167 ngot=gos[gos$term!=t,] 176 for (t in terms) {
168 ngot=ngot[!duplicated(ngot$seq),] 177 got <- gos[gos$term == t, ]
169 ngot=ngot[!(ngot$seq %in% got$seq),] 178 got <- got[!duplicated(got$seq), ]
170 go.sig=sum(got$value) 179 ngot <- gos[gos$term != t, ]
171 go.ns=length(got[,1])-go.sig 180 ngot <- ngot[!duplicated(ngot$seq), ]
172 ngo.sig=sum(ngot$value) 181 ngot <- ngot[!(ngot$seq %in% got$seq), ]
173 ngo.ns=length(ngot[,1])-ngo.sig 182 go.sig <- sum(got$value)
174 sig=c(go.sig,ngo.sig) # number of significant genes belonging and not belonging to the tested GO category 183 go.ns <- length(got[, 1]) - go.sig
175 ns=c(go.ns,ngo.ns) # number of not-significant genes belonging and not belonging to the tested GO category 184 ngo.sig <- sum(ngot$value)
176 mm=matrix(c(sig,ns),nrow=2,dimnames=list(ns=c("go","notgo"),sig=c("go","notgo"))) 185 ngo.ns <- length(ngot[, 1]) - ngo.sig
177 ff=fisher.test(mm,alternative="greater") 186 sig <- c(go.sig, ngo.sig) # number of significant genes belonging and not belonging to the tested GO category
178 pft=append(pft,ff$p.value) 187 ns <- c(go.ns, ngo.ns) # number of not-significant genes belonging and not belonging to the tested GO category
179 nam=append(nam,as.character(got$name[1])) 188 mm <- matrix(c(sig, ns), nrow = 2, dimnames = list(ns = c("go", "notgo"), sig = c("go", "notgo")))
180 lev=append(lev,got$lev[1]) 189 ff <- fisher.test(mm, alternative = "greater")
181 nseqs=append(nseqs,length(got[,1])) 190 pft <- append(pft, ff$p.value)
182 } 191 nam <- append(nam, as.character(got$name[1]))
183 res=data.frame(cbind("delta.rank"=rep(0),"pval"=pft,"level"=lev,nseqs,"term"=terms,"name"=nam)) 192 lev <- append(lev, got$lev[1])
184 res[,1]=as.numeric(as.character(res[,1])) 193 nseqs <- append(nseqs, length(got[, 1]))
185 res[,2]=as.numeric(as.character(res[,2])) 194 }
186 res[,3]=as.numeric(as.character(res[,3])) 195 res <- data.frame(cbind("delta.rank" = rep(0), "pval" = pft, "level" = lev, nseqs, "term" = terms, "name" = nam))
187 res$nseqs=as.numeric(as.character(res$nseqs)) 196 res[, 1] <- as.numeric(as.character(res[, 1]))
188 return(res) 197 res[, 2] <- as.numeric(as.character(res[, 2]))
198 res[, 3] <- as.numeric(as.character(res[, 3]))
199 res$nseqs <- as.numeric(as.character(res$nseqs))
200 return(res)
189 } 201 }
190 202
191 #------------------------- 203 #-------------------------
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) { 204 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) {
193 require(ape) 205 require(ape)
194 206
195 nn <- strsplit(inFile, "[/.]") 207 nn <- strsplit(inFile, "[/.]")
196 if (length(nn[[1]]) == 3) { 208 if (length(nn[[1]]) == 3) {
197 dir <- nn[[1]][1] 209 dir <- nn[[1]][1]
198 name <- nn[[1]][2] 210 name <- nn[[1]][2]
199 ext <- nn[[1]][3] 211 ext <- nn[[1]][3]
210 } else if (length(aa[[1]]) == 2) { 222 } else if (length(aa[[1]]) == 2) {
211 adir <- "." 223 adir <- "."
212 aname <- aa[[1]][1] 224 aname <- aa[[1]][1]
213 aext <- aa[[1]][2] 225 aext <- aa[[1]][2]
214 } 226 }
215 # input=inFile 227 # input=inFile
216 in.mwu=paste(dir,"/",paste("MWU",goDivision,name,sep="_"), ".", ext,sep="") 228 in.mwu <- paste(dir, "/", paste("MWU", goDivision, name, sep = "_"), ".", ext, sep = "")
217 in.dissim=paste(dir, "/", paste("dissim",goDivision,name,aname,sep="_"), ".", aext, sep="") 229 in.dissim <- paste(dir, "/", paste("dissim", goDivision, name, aname, sep = "_"), ".", aext, sep = "")
218 230
219 cutoff=-log(level1,10) 231 cutoff <- -log(level1, 10)
220 pv=read.table(in.mwu,header=T) 232 pv <- read.table(in.mwu, header = T)
221 row.names(pv)=pv$term 233 row.names(pv) <- pv$term
222 in.raw=paste(dir,"/",paste(name,goDivision,sep="_"), ".tsv", sep="") 234 in.raw <- paste(dir, "/", paste(name, goDivision, sep = "_"), ".tsv", sep = "")
223 rsq=read.table(in.raw,sep="\t",header=T) 235 rsq <- read.table(in.raw, sep = "\t", header = T)
224 rsq$term=as.factor(rsq$term) 236 rsq$term <- as.factor(rsq$term)
225 237
226 if (adjusted==TRUE) { pvals=pv$p.adj } else { pvals=pv$pval } 238 if (adjusted == TRUE) {
227 heat=data.frame(cbind("pval"=pvals)) 239 pvals <- pv$p.adj
228 row.names(heat)=pv$term 240 } else {
229 heat$pval=-log(heat$pval+1e-15,10) 241 pvals <- pv$pval
230 heat$direction=0 242 }
231 heat$direction[pv$delta.rank>0]=1 243 heat <- data.frame(cbind("pval" = pvals))
232 if (cutoff>0) { 244 row.names(heat) <- pv$term
233 goods=subset(heat,pval>=cutoff) 245 heat$pval <- -log(heat$pval + 1e-15, 10)
234 } else { 246 heat$direction <- 0
235 goods.names=unique(rsq$term[abs(rsq$value)>=absValue]) 247 heat$direction[pv$delta.rank > 0] <- 1
236 goods=heat[row.names(heat) %in% goods.names,] 248 if (cutoff > 0) {
237 } 249 goods <- subset(heat, pval >= cutoff)
238 250 } else {
239 if (is.null(colors) | length(colors)<4 ) { 251 goods.names <- unique(rsq$term[abs(rsq$value) >= absValue])
240 colors=c("dodgerblue2","firebrick1","skyblue2","lightcoral") 252 goods <- heat[row.names(heat) %in% goods.names, ]
241 if (sum(goods$direction)==nrow(goods) | sum(goods$direction)==0) { 253 }
242 colors=c("black","black","grey50","grey50") 254
243 } 255 if (is.null(colors) | length(colors) < 4) {
244 } 256 colors <- c("dodgerblue2", "firebrick1", "skyblue2", "lightcoral")
245 goods.names=row.names(goods) 257 if (sum(goods$direction) == nrow(goods) | sum(goods$direction) == 0) {
246 258 colors <- c("black", "black", "grey50", "grey50")
247 # reading and subsetting dissimilarity matrix 259 }
248 diss=read.table(in.dissim,sep="\t",header=T,check.names=F) 260 }
249 row.names(diss)=names(diss) 261 goods.names <- row.names(goods)
250 diss.goods=diss[goods.names,goods.names] 262
251 263 # reading and subsetting dissimilarity matrix
252 # how many genes out of what we started with we account for with our best categories? 264 diss <- read.table(in.dissim, sep = "\t", header = T, check.names = F)
253 good.len=c();good.genes=c() 265 row.names(diss) <- names(diss)
254 for (g in goods.names) { 266 diss.goods <- diss[goods.names, goods.names]
255 sel=rsq[rsq$term==g,] 267
256 pass=abs(sel$value)>=absValue 268 # how many genes out of what we started with we account for with our best categories?
257 sel=sel[pass,] 269 good.len <- c()
258 good.genes=append(good.genes,as.character(sel$seq)) 270 good.genes <- c()
259 good.len=append(good.len,nrow(sel)) 271 for (g in goods.names) {
260 } 272 sel <- rsq[rsq$term == g, ]
261 ngenes=length(unique(good.genes)) 273 pass <- abs(sel$value) >= absValue
262 274 sel <- sel[pass, ]
263 #hist(rsq$value) 275 good.genes <- append(good.genes, as.character(sel$seq))
264 totSum=length(unique(rsq$seq[abs(rsq$value)>=absValue])) 276 good.len <- append(good.len, nrow(sel))
265 row.names(goods)=paste(good.len,"/",pv[pv$term %in% goods.names,]$nseqs," ",pv[pv$term %in% goods.names,]$name,sep="") 277 }
266 row.names(heat)=paste(good.len,"/",pv$nseqs," ",pv$name,sep="") 278 ngenes <- length(unique(good.genes))
267 row.names(diss.goods)=paste(good.len,"/",pv[pv$term %in% goods.names,]$nseqs," ",pv[pv$term %in% goods.names,]$name,sep="") 279
268 280 # hist(rsq$value)
269 # clustering terms better than cutoff 281 totSum <- length(unique(rsq$seq[abs(rsq$value) >= absValue]))
270 GO.categories=as.dist(diss.goods) 282 row.names(goods) <- paste(good.len, "/", pv[pv$term %in% goods.names, ]$nseqs, " ", pv[pv$term %in% goods.names, ]$name, sep = "")
271 cl.goods=hclust(GO.categories,method="average") 283 row.names(heat) <- paste(good.len, "/", pv$nseqs, " ", pv$name, sep = "")
272 labs=cl.goods$labels[cl.goods$order] # saving the labels to order the plot 284 row.names(diss.goods) <- paste(good.len, "/", pv[pv$term %in% goods.names, ]$nseqs, " ", pv[pv$term %in% goods.names, ]$name, sep = "")
273 goods=goods[labs,] 285
274 labs=sub(" activity","",labs) 286 # clustering terms better than cutoff
275 287 GO.categories <- as.dist(diss.goods)
276 old.par <- par( no.readonly = TRUE ) 288 cl.goods <- hclust(GO.categories, method = "average")
277 289 labs <- cl.goods$labels[cl.goods$order] # saving the labels to order the plot
278 plots=layout(matrix(c(1,2,3),1,3,byrow=T),c(treeHeight,3,1),TRUE) 290 goods <- goods[labs, ]
279 291 labs <- sub(" activity", "", labs)
280 par(mar = c(2,2,0.85,0)) 292
281 plot(as.phylo(cl.goods),show.tip.label=FALSE,cex=0.0000001) 293 old.par <- par(no.readonly = TRUE)
282 step=100 294
283 left=1 295 plots <- layout(matrix(c(1, 2, 3), 1, 3, byrow = T), c(treeHeight, 3, 1), TRUE)
284 top=step*(2+length(labs)) 296
285 297 par(mar = c(2, 2, 0.85, 0))
286 par(mar = c(0,0,0.3,0)) 298 plot(as.phylo(cl.goods), show.tip.label = FALSE, cex = 0.0000001)
287 plot(c(1:top)~c(1:top),type="n",axes=F,xlab="",ylab="") 299
288 ii=1 300 par(mar = c(2, 0, 1, 0))
289 goods$color=1 301 step <- dev.size("px")[2] / (length(goods.names) - 1)
290 goods$color[goods$direction==1 & goods$pval>cutoff]=colors[4] 302 left <- 1
291 goods$color[goods$direction==0 & goods$pval>cutoff]=colors[3] 303 top <- step * (length(labs) - 1)
292 goods$color[goods$direction==1 & goods$pval>(-log(level2,10))]=colors[2] 304 # print(paste("Size:", dev.size("px")))
293 goods$color[goods$direction==0 & goods$pval>(-log(level2,10))]=colors[1] 305 # print(paste("Good:", length(goods.names), "Step:", step, "Top:", top))
294 goods$color[goods$direction==1 & goods$pval>(-log(level3,10))]=colors[2] 306
295 goods$color[goods$direction==0 & goods$pval>(-log(level3,10))]=colors[1] 307 plot(c(1:top), c(1:top), type = "n", axes = F, xlab = "", ylab = "")
296 for (i in length(labs):1) { 308 ii <- 0
297 ypos=top-step*ii 309 goods$color <- 1
298 ii=ii+1 310 goods$color[goods$direction == 1 & goods$pval > cutoff] <- colors[4]
299 if (goods$pval[i]> -log(level3,10)) { 311 goods$color[goods$direction == 0 & goods$pval > cutoff] <- colors[3]
300 text(left,ypos,labs[i],font=2,cex=1*txtsize,col=goods$color[i],adj=c(0,0),family=font.family) 312 goods$color[goods$direction == 1 & goods$pval > (-log(level2, 10))] <- colors[2]
301 } else { 313 goods$color[goods$direction == 0 & goods$pval > (-log(level2, 10))] <- colors[1]
302 if (goods$pval[i]>-log(level2,10)) { 314 goods$color[goods$direction == 1 & goods$pval > (-log(level3, 10))] <- colors[2]
303 text(left,ypos,labs[i],font=1,cex=0.8* txtsize,col=goods$color[i],adj=c(0,0),family=font.family) 315 goods$color[goods$direction == 0 & goods$pval > (-log(level3, 10))] <- colors[1]
304 } else { 316 for (i in length(labs):1) {
305 # if (goods$pval[i]>cutoff) { 317 ypos <- round(top - step * ii)
306 # text(left,ypos,labs[i],font=3,cex=0.8* txtsize,col=goods$color[i],adj=c(0,0),family=font.family) 318 # print(paste("Ypos:", ypos))
307 # } else { 319 ii <- ii + 1
308 text(left,ypos,labs[i],font=3,cex=0.8* txtsize,col=goods$color[i],adj=c(0,0),family=font.family) 320 if (goods$pval[i] > -log(level3, 10)) {
309 #} 321 text(left, ypos, labs[i], font = 2, cex = 1 * txtsize, col = goods$color[i], adj = c(0, 0), family = font.family)
310 } 322 } else {
311 } 323 if (goods$pval[i] > -log(level2, 10)) {
312 } 324 text(left, ypos, labs[i], font = 1, cex = 0.8 * txtsize, col = goods$color[i], adj = c(0, 0), family = font.family)
313 325 } else {
314 par(mar = c(3,1,1,0)) 326 # if (goods$pval[i]>cutoff) {
315 327 # text(left,ypos,labs[i],font=3,cex=0.8* txtsize,col=goods$color[i],adj=c(0,0),family=font.family)
316 plot(c(1:top)~c(1:top),type="n",axes=F,xlab="",ylab="") 328 # } else {
317 text(left,top-step*2,paste("p < ",level3,sep=""),font=2,cex=1* txtsize,adj=c(0,0),family=font.family) 329 text(left, ypos, labs[i], font = 3, cex = 0.8 * txtsize, col = goods$color[i], adj = c(0, 0), family = font.family)
318 text(left,top-step*3,paste("p < ",level2,sep=""),font=1,cex=0.8* txtsize,adj=c(0,0),family=font.family) 330 # }
319 text(left,top-step*4,paste("p < ",10^(-cutoff),sep=""),font=3,col="grey50",cex=0.8* txtsize,adj=c(0,0),family=font.family) 331 }
320 332 }
321 message("GO terms dispayed: ",length(goods.names)) 333 }
322 message("\"Good genes\" accounted for: ", ngenes," out of ",totSum, " ( ",round(100*ngenes/totSum,0), "% )") 334
323 par(old.par) 335 par(mar = c(3, 1, 1, 0))
324 goods$pval=10^(-1*goods$pval) 336
325 return(list(goods,cl.goods)) 337 plot(c(1:top), c(1:top), type = "n", axes = F, xlab = "", ylab = "")
338 text(left, top - step * 2, paste("p < ", level3, sep = ""), font = 2, cex = 1 * txtsize, adj = c(0, 0), family = font.family)
339 text(left, top - step * 3, paste("p < ", level2, sep = ""), font = 1, cex = 0.8 * txtsize, adj = c(0, 0), family = font.family)
340 text(left, top - step * 4, paste("p < ", 10^(-cutoff), sep = ""), font = 3, col = "grey50", cex = 0.8 * txtsize, adj = c(0, 0), family = font.family)
341
342 message("GO terms dispayed: ", length(goods.names))
343 message("\"Good genes\" accounted for: ", ngenes, " out of ", totSum, " ( ", round(100 * ngenes / totSum, 0), "% )")
344 par(old.par)
345 goods$pval <- 10^(-1 * goods$pval)
346 return(list(goods, cl.goods))
326 } 347 }
327 348
328 #------------------ 349 #------------------
329 # returns non-overlapping GO categories based on dissimilarity table 350 # returns non-overlapping GO categories based on dissimilarity table
330 indepGO=function(dissim.table,min.similarity=1) { 351 indepGO <- function(dissim.table, min.similarity = 1) {
331 tt=read.table(dissim.table,sep="\t", header=TRUE) 352 tt <- read.table(dissim.table, sep = "\t", header = TRUE)
332 tt=as.matrix(tt) 353 tt <- as.matrix(tt)
333 diag(tt)=1 354 diag(tt) <- 1
334 for (i in 1:ncol(tt)) { 355 for (i in 1:ncol(tt)) {
335 mins=apply(tt,2,min) 356 mins <- apply(tt, 2, min)
336 if (min(mins)>=min.similarity) { break } 357 if (min(mins) >= min.similarity) {
337 sums=apply(tt,2,sum) 358 break
338 worst=which(sums==min(sums))[1] 359 }
339 # cat("\n",worsts,"\n") 360 sums <- apply(tt, 2, sum)
340 # gw=c() 361 worst <- which(sums == min(sums))[1]
341 # for(w in worsts) { gw=append(gw,sum(tt[,w]==1)) } 362 # cat("\n",worsts,"\n")
342 # cat(gw) 363 # gw=c()
343 # worst=worsts[gw==min(gw)] 364 # for(w in worsts) { gw=append(gw,sum(tt[,w]==1)) }
344 # cat("\n",i," removing",worst,"; sum =",sums[worst]) 365 # cat(gw)
345 tt=tt[-worst,-worst] 366 # worst=worsts[gw==min(gw)]
346 mins=mins[-worst] 367 # cat("\n",i," removing",worst,"; sum =",sums[worst])
347 # cat(" new min =",min(mins)) 368 tt <- tt[-worst, -worst]
348 } 369 mins <- mins[-worst]
349 goods=colnames(tt) 370 # cat(" new min =",min(mins))
350 goods=gsub("GO\\.","GO:",goods) 371 }
351 goods=gsub("\\.GO",";GO",goods) 372 goods <- colnames(tt)
352 } 373 goods <- gsub("GO\\.", "GO:", goods)
353 374 goods <- gsub("\\.GO", ";GO", goods)
375 }