# TODO: Add comment # # Author: E.Korsching Jul 16, 2012 ############################################################################### perm.testing.1 <- function(x, refset, testset, method="Pearson", subfolder="path", fname="test01", ylim=c(-1,1)){ ## comprehensive test function: of a given grouping and of the best solution in this grouping # x: data.frame # refset: c(4,5,6,7,8,10) (e.g. CK1 CK10 CK19 CK8.18 CK5 CK14) # testset: c(1:3,9,11,14,17:20) (e.g. Vimentin p53 ERBB2 BCL2 CyclinD1 AR ER PR Mib1 EMA) # method: Pearson, Spearman or L2 # subfolder: path, fname: file name acronym of test approach # ylim: y axis range in regression plot # ini anzSamplings <- 1000 len.ref <- length(refset) len.test <- length(testset) x <- x[,c(testset,refset),drop=F] #order and filter x.names <- names(x) #new names testset <- c(1:len.test) #new test refset <- c((len.test+1):(len.test+len.ref)) #new ref # optimal order # 2011: 3s (6/10 in total approx. 5h) , check var=0 ! x.hy1.f90 <- enumeration.search.dep.boots.f90(x=x, anzSamplings=0, refset=refset, testset=testset, method=method, sampleonly=F, output="boot.raw.col", ssqflag=1) # 2011: 28min x.hy1.b1000f90 <- enumeration.search.dep.boots.f90(x=x, anzSamplings=anzSamplings, refset=refset, testset=testset, method=method, sampleonly=F, output="boot.raw.col", ssqflag=1) #extrema ssq.given.situation.min <- x.hy1.b1000f90[[2]][1,1] sample.min.row <- which.min(x.hy1.b1000f90[[2]][,1]) #642 sample.min.value <- min(x.hy1.b1000f90[[2]][,1], na.rm=T) #151.7538 sample.max.value <- max(x.hy1.b1000f90[[2]][,1], na.rm=T) ssq.sample <- x.hy1.b1000f90[[2]][ ,1] number.better <- sum(ssq.given.situation.min>ssq.sample, na.rm=T) sample.p.value <- number.better/anzSamplings #45/1000 = 0.045 #regression pdf(file=paste(getwd(),"/",subfolder,"/",fname,".regression.",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""), width=11, height=7, onefile=T, title=paste(fname,"regression"), pointsize=12) plot.sample.search.dependencies(x=list( list(x.hy1.f90[[3]][1:len.test,,drop=F], x.hy1.f90[[1]][1,2:(len.test+1)], "black") #[1, (len.test+1):2] ), ylim=ylim, y.at=NULL, ylab="correlation", text=T, description="best", cex.t=0.8, cex.rt=1, mfrow=c(2,3), ssqflag=1 ) plot.sample.search.dependencies(x=list( list(x.hy1.b1000f90[[3]][((sample.min.row*len.test)-(len.test-1)):(sample.min.row*len.test),,drop=F], x.hy1.b1000f90[[1]][sample.min.row,2:(len.test+1)], "black") #11:2 ), ylim=ylim, y.at=NULL, ylab="correlation", text=T, description="bootstrap best", cex.t=0.8, cex.rt=1, mfrow=c(2,3), ssqflag=1 ) dev.off() #base distribution # 2011: 1s x.hy1.f90dist <- enumeration.search.dep.boots.f90.base.dist(x=x, refset=refset, testset=testset, method=method, zz=NULL, prox.only=F, output="boot.raw.col", ssqflag=1) #all optimal orders - 2011: 3h55min x.hy1.extf90 <- enumeration.search.dep.boots.GOnly.extf90(x=x, refset=refset, testset=testset, method=method, ssqflag=1) #extrema ssq.given.situation.max <- max(x.hy1.f90dist[[1]], na.rm=T) comb.min.row <- which.min(x.hy1.extf90[,1]) #406 comb.min.value <- min(x.hy1.extf90[,1], na.rm=T) #210.2306 comb.max.value <- max(x.hy1.extf90[,1], na.rm=T) orders.min <- x.hy1.extf90[comb.min.row, ] #order: best ssq, best test order of 'selected test', selected test, selected ref #hist pdf(file=paste(getwd(),"/",subfolder,"/",fname,".histo.",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""), width=11, height=7, onefile=T, title=paste(fname,"histo"), pointsize=12) hist.plot.color(x=x.hy1.f90dist[[1]], bin.num=30, x.range=NULL, norm=F, n.factor=1, lty=3, lwd=1, xlab="ssqg range", ylab="", xlim=NULL, y.max=NULL, x.at=NULL, y.at=NULL, ylog=F, bar.width=0.7, offset=0, digits=1, col.grad=c("red","blue"), col.f="blue", col.sh=NULL, col.b=NA, cex=1, h.title=paste(fname,"distribution of all ssqg in the selected set"), add=F, rC=F) hist.plot.color(x=x.hy1.b1000f90[[2]][,1], bin.num=30, x.range=NULL, norm=F, n.factor=1, lty=3, lwd=1, xlab="", ylab="", xlim=NULL, y.max=NULL, x.at=NULL, y.at=NULL, ylog=F, bar.width=0.7, offset=0, digits=1, col.grad=c("red","blue"), col.f="blue", col.sh=NULL, col.b=NA, cex=1, h.title=paste(fname,"distribution of bootstrap ssqg in the selected set, p=", round( sum(x.hy1.b1000f90[[2]][,1]ssq2.sample, na.rm=T) sample2.p.value <- number2.better/anzSamplings #regression pdf(file=paste(getwd(),"/",subfolder,"/",fname,".regression.best.",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""), width=11, height=7, onefile=T, title=paste(fname,"regression.best"), pointsize=12) plot.sample.search.dependencies(x=list( list(x.hy1a.f90[[3]][1:len.test,,drop=F], x.hy1a.f90[[1]][1,2:(len.test+1)], "black") #[1, (len.test+1):2] ), ylim=ylim, y.at=NULL, ylab="correlation", text=T, description="best of best combination", cex.t=0.8, cex.rt=1, mfrow=c(2,3), ssqflag=1 ) plot.sample.search.dependencies(x=list( list(x.hy1a.b1000f90[[3]][((sample2.min.row*len.test)-(len.test-1)):(sample2.min.row*len.test),,drop=F], x.hy1a.b1000f90[[1]][sample2.min.row,2:(len.test+1)], "black") #[1, (len.test+1):2] ), ylim=ylim, y.at=NULL, ylab="correlation", text=T, description="bootstrap best of best combination", cex.t=0.8, cex.rt=1, mfrow=c(2,3), ssqflag=1 ) dev.off() #return all information cat("\n\n") cat("\n #### data / approach: ", fname, " , date: ", date(), " ####") cat("\n number of test permutations: ", factorial(len.test), " , number of ref combinations: ", choose(len.ref+len.test, len.ref)) cat("\n number of samplings: ", anzSamplings) cat("\n 1) given grouping") cat("\n test: ", x.names[testset]," , ref: ", x.names[refset]) cat("\n min ssqg given grouping: ", ssq.given.situation.min, " , max (permutation) ssqg: ", ssq.given.situation.max) cat("\n min sampling ssqg: ", sample.min.value, " , max sampling ssqg: ", sample.max.value, " - row with min: ", sample.min.row) cat("\n smaller solutions (", number.better, " / ", anzSamplings, ") , sample p value ", sample.p.value) cat("\n 2) min ssqg of all combinations: ", comb.min.value, " , max ssqg: ", comb.max.value, " - row with min: ", comb.min.row) cat("\n 3) selected min combination") cat("\n test: ", x.names[test.2]," , ref: ", x.names[ref.2]) cat("\n min ssqg selected combination: ", ssq2.given.situation.min, " , max (permutation) ssqg: ", ssq2.given.situation.max) cat("\n min sampling ssqg: ", sample2.min.value, " , max sampling ssqg: ", sample2.max.value, " - row with min: ", sample2.min.row) cat("\n smaller solutions (", number2.better, " / ", anzSamplings, ") , sample p value: ", sample2.p.value) cat("\n\n") #return results - given grouping + best in data grouping return(list(given.f90=x.hy1.f90, bestInData.f90=x.hy1a.f90)) } # 04-2024: the whole run time is: 11s for: #bc4.x <- perm.testing.1(aa1, refset=c(1,2,3,5,9,10), testset=c(4,6,7,8,11,12,13), method="Pearson", subfolder="tech_tests", fname="He589366c2_x_", ylim=c(-1,1)) #horst.test1 <- perm.testing.1(x=horstLCA.2,refset=c(1,2,8,9),testset=c(3,4,5,6,7,10,11),subfolder="result",fname="horstLCA.2")