#################################################################################################
##
#################################################################################################
getSkipFromGPR <- function(file){
        con <- file(file,"r")
        twoLines <- readLines(con,n=2)
    close(con)
        if(substring(twoLines,1,3)[1]=="ATF"){
                skip <- substring(twoLines,1,2)
        }
        else{
                skip=c(0,-2)
        }
        as.numeric(skip[2])+2
}

###########################################################################################
## creates a new GAL file from the information from the gpr files
###########################################################################################
createGALFromGPR <- function(file,galfile="dummy.gal",colnames=c("Block","Column","Row","Name","ID")){
        table <- read.table(file,as.is=TRUE,comment.char="",check.names=FALSE,header=TRUE,skip=getSkipFromGPR(file),sep="\t")
        names <- colnames(table)
        result <- data.frame(matrix(ncol=5,nrow=length(table[,1])))
        colnames(result) <- colnames
        for(i in 1:length(result[1,])){
                result[,i] <- table[,colnames[i]]
        }
        #result[,1:5] <- table[,1:5]
        #colnames(result) <- names[1:5]
        write.table(file=galfile,result,sep="\t",row.name=FALSE,quote=FALSE)
}

plotHistogram <- function(R,G,r.col=2,g.col=3,breaks=NULL,log.transform=FALSE,add=FALSE,weights=NULL,return.breaks=FALSE,ylim=NULL,xlim=NULL,title=NULL,relative=FALSE,...){
        if(is.null(weights)){
                weights <- matrix(ncol=1,nrow=length(R))
                weights[,] <- 1
        }

        if(log.transform==TRUE){
                r <- log2(R)
                g <- log2(G)
                by <- 0.1
                toadd <- 0.1
        }
        else{
                r <- R
                g <- G
                by <- 50
                toadd <- 50
        }
        r <- r[weights==1]
        g <- g[weights==1]
        r <- r[!is.na(r)]
        g <- g[!is.na(g)]
        r <- r[!is.infinite(r)]
        g <- g[!is.infinite(g)]

        if(is.null(breaks)){
                breaks <- seq(min(0,(min(r,g))),(max(max(r),max(g))+toadd),by=by)
        }
        
        hr <- hist(r,plot=FALSE,breaks=breaks)
        hg <- hist(g,plot=FALSE,breaks=breaks)
        if(is.null(xlim)){
                xlim <- range(breaks)
        }

        r.counts <- hr$counts
        g.counts <- hg$counts
        ylabb <- "Counts"
        if(relative){
                r.counts <- r.counts / length(r)
                g.counts <- g.counts / length(g)
                ylabb <- "relative Counts"
                yy <- -0.001
        }
        else{
                yy <- -10
        }

        if(is.null(ylim)){
                ylim <- c(0,max(r.counts,g.counts))
        }
                
        if(!add){
                plot(hr$mids,r.counts,col=r.col,xlab="Intensity",ylab=ylabb,type="l",ylim=ylim,xlim=xlim,main=title,...)
        }
        else{
                points(hr$mids,r.counts,col=r.col,type="l",...)
        }
        points(hg$mids,g.counts,col=g.col,type="l",...)
        points(medianWithoutNA(r),yy,col=r.col,pch=3)
        points(medianWithoutNA(g),yy,col=g.col,pch=3)
        
        if(return.breaks){
                breaks
        }
#       hr <- density(r,na.rm=TRUE)
#       hg <- density(g,na.rm=TRUE)
#       xlim <- c(min(hr$x,hg$x),max(hr$x,hg$x))
#       ylim <- c(min(hr$y,hg$y),max(hr$y,hg$y))
#       plot(hr$x,hr$y,col=2,type="l",ylim=ylim,...)
#       points(hg$x,hg$y,col=3,type="l")
#       hr <- hist(r,breaks=breaks,plot=FALSE)
#       hg <- hist(g,breaks=breaks,plot=FALSE)
}


#############################################################################
##
#############################################################################
plotHistogramFromSlide <- function(lObject,slide=1,use.weights=TRUE,log.transform=FALSE,breaks=NULL,bg.lty=3,...){
        if(slide=="all"){
        
        
        
        }
        if(is.null(lObject$R)){
                R <- MAtoR(lObject$M[,slide],lObject$A[,slide])
                G <- MAtoG(lObject$M[,slide],lObject$A[,slide])
                Rb <- NULL
                Gb <- NULL
                draw.bg <- FALSE
                names <- colnames(lObject$M)
        }
        else{
                R <- lObject$R[,slide]
                G <- lObject$G[,slide]
                Rb <- lObject$Rb[,slide]
                Gb <- lObject$Gb[,slide]
                draw.bg <- !is.null(Rb)
                names <- colnames(lObject$R)
        }
        weights <- matrix(ncol=1,nrow=length(R))
        weights[,] <- 1
        if(is.null(lObject$weights)){
        }
        else{
                if(use.weights){
                        weights <- lObject$weights[,slide]
                }
        }
#       r <- R[weights==1]
#       g <- G[weights==1]
#       rb <- Rb[weights==1]
#       gb <- Gb[weights==1]
        
        breaks <- plotHistogram(R=R,G=G,log.transform=log.transform,weights=weights,breaks=breaks,return.breaks=TRUE,title=names[slide],...)
        if(draw.bg){
                plotHistogram(R=Rb,G=Gb,log.transform=log.transform,breaks=NULL,weights=weights,lty=bg.lty,add=TRUE,...)
        }
}




checkExistantFile <- function(file,path=".",no.subdir=TRUE){
        file <- gsub(pattern="_",replacement=".",file)
        file <- gsub(pattern="/",replacement="-",file)
        counter <- 0
        suff <- NULL
        ending <- NULL
        ffile <- file
        splitted <- unlist(strsplit(file,split="\\."))
        if(length(splitted)>1){
                ending <- paste(".",splitted[length(splitted)],sep="")
                ffile <- paste(splitted[1:(length(splitted)-1)],collapse=".")
                #cat(ffile,"_",ending,"\n")
        }
        while(TRUE){
                newfile <- paste(ffile,suff,ending,sep="")
                if(file.exists(newfile)){
                        #file does exist...
                        counter <- counter +1
                        suff <- paste("(",counter,")",sep="")
                }
                else{
                        break
                }
        }
        newfile
}




saveMAList <- function(object,filename,remove.flagged=TRUE,na.value="NA",subset=NULL, allvalues=TRUE){
        genes <- NULL
        weights <- NULL
        g.length <- 0
        col.names <- NULL
        if(class(object)=="MAList"){
                if(is.null(dim(object$M))){
                        object$M <- as.matrix(object$M)
                        object$A <- as.matrix(object$A)
                        if(!is.null(object$weights)){
                                object$weights <- as.matrix(object$weights)
                        }
                }
                if(is.null(subset)){
                        subset <- rep(TRUE,nrow(object$M))
                }
                M <- as.matrix(object$M[subset,])
                CN <- colnames(object$M)
                if(is.null(CN)){
                        CN <- object$targets[1,1]
                }
                colnames(M) <- CN
#               colnames(M) <- colnames(object$M)
                A <- as.matrix(object$A[subset,])
                colnames(A) <- CN
#               colnames(A) <- colnames(object$A)
                if(!is.null(object$weights)){
                        weights <- as.matrix(object$weights[subset,])
                }
                if(!is.null(object$genes)){
                        genes <- object$genes[subset,]
                        if(!is.null(genes)){
                                g.length <- ncol(genes)
                        }
                }
        }
        else if(class(object)=="EexprSet" | class(object)=="MadbSet"){
                if(object@type=="TwoColor"){
                        if(is.null(subset)){
                                subset <- rep(TRUE,nrow(exprs(object)))
                        }
                        M <- getM(object)
                        Colnames <- colnames(M)
                        M <- as.matrix(M[subset,])
                        colnames(M) <- Colnames
#                       M <- as.matrix(getM(object)[subset,])
#                       cat(class(M),colnames(M),"\n")
#                       A <- as.matrix(getA(object)[subset,])
                        A <- getA(object)
                        A <- as.matrix(A[subset,])
                        colnames(A) <- Colnames
                        #cat("1")
                        weights <- as.matrix(getWeights(object)[subset,])
                        #cat("2")
                        genes <- as.matrix(object@genes[subset,])
                        if(!is.null(genes)){
                                g.length <- ncol(genes)
                                if(g.length==1){
                                        colnames(genes) <- "ID"
                                }
                        }
                        #cat("g.length",g.length,"dim M",dim(M),"\n")
                        #cat("3")
                }
        }
        else if(class(object)=="RGList"){
                if(is.null(subset)){
                        subset <- rep(TRUE,nrow(object$R))
                }
                if(!.is.log(object$R)){
                        object$R <- log2(object$R)
                }
                if(!.is.log(object$G)){
                        object$G <- log2(object$G)
                }
                M <- object$R[subset,] - object$G[subset,]
                A <- 0.5*(object$R[subset,] + object$G[subset,])
                if(!is.null(object$weights)){
                        weights <- object$weights[subset,]
                }
                if(!is.null(object$genes)){
                        genes <- object$genes[subset,]
                        if(!is.null(genes)){
                                g.length <- ncol(genes)
                        }
                }
        }
        else{
                stop("Objects from the class ",class(object)," can not be handled by this function\n")
        }
        #cat("ncol=",g.length+(ncol(M)*2),"nrow=",nrow=nrow(M),"\n")
        Table <- matrix(ncol=(g.length+(ncol(M)*2)),nrow=nrow(M))
        if(!is.null(genes)){
                Table[,1:ncol(genes)] <- as.matrix(genes)
                col.names <- colnames(genes)
        }

        #cat(colnames(M),"\n")

        for(i in 1:ncol(M)){
                col.names <- c(col.names,paste(colnames(M)[i],"M"),paste(colnames(A)[i],"A"))
                Table[,((i*2)-1+g.length)] <- M[,i]
                if(!is.null(weights) & remove.flagged){
                        Table[weights[,i]==0,((i*2)-1+g.length)] <- NA
                }
                Table[is.na(Table[,((i*2)-1+g.length)]),((i*2)-1+g.length)] <- NA
                
                Table[,((i*2)+g.length)] <- A[,i]
                if(!is.null(weights) & remove.flagged){
                        Table[weights[,i]==0,((i*2)+g.length)] <- NA
                }
                Table[is.na(Table[,((i*2)+g.length)]),((i*2)+g.length)] <- NA
        }
        colnames(Table) <- col.names
        
        if(allvalues){
                ## calculate the red and green intensities from the M and A values...
                R <- MAtoR(M, A)
                G <- MAtoG(M, A)
                if(is.null(ncol(R))){
                        R <- matrix(R, ncol=1)
                        G <- matrix(G, ncol=1)
                }
                colnames(R) <- paste(colnames(M), "R")
                colnames(G) <- paste(colnames(M), "G")
                Table <- cbind(Table, R, G)
        }
        
        write.table(Table,file=filename,sep="\t",na=na.value,row.names=FALSE,quote=FALSE)
}

## closes all connections of gpr files...
closeMyConnections <- function(pattern="gpr"){
    All <- showConnections(TRUE)
    what <- grep(All[,"description"],pattern=pattern)
    what <- what-1
    for(theCon in what){
        close(getConnection(theCon))
    }
}

checkGPRs <- function(files=NULL,overwrite=TRUE){
        ##
        # check for correct axon format...
        requiredcolumns <- c("Block","Column","Row","Name","ID","F635 Median","F635 Mean","B635 Median","B635 Mean","F532 Median","F532 Mean","B532 Median","B532 Mean","Flags")
        for(i in 1:length(files)){
                con <- file(files[i],"r")
                twoLines <- readLines(con,n=2)
                close(con)
#        twoLines <- scan(file=files[i], what="character",nlines=2, quiet=TRUE)
#        closeMyConnections()
                if(substring(twoLines,1,3)[1]!="ATF"){
                        # is not in correct Axon file format
                        cat("The file",files[i],"is not in correct Axon file format!\n")
                        header <- read.table(files[i],header=TRUE,nrows=1,sep="\t",check.names=FALSE)
                        if(sum(colnames(header) %in% requiredcolumns)==length(requiredcolumns)){
                                BigFile <- readLines(con=files[i])
                                writeLines(c("ATF\t1.0",paste("1\t",ncol(header),sep=""),"Provider=CARMAweb"),con=files[i])
                                con <- file(files[i],"at")
                                writeLines(con=con,BigFile)
                flush(con)
                                close(con)
                                cat("Successfully adjusted file",files[i],".\n")
#                closeMyConnections()
                        }
                        else{
                                stop("The file",files[i],"is not in correct Axon file format! Was this file scanned with an Axon GenePix Scanner?\n")
                        }
                }
        }

        skips <- NULL
        for(i in 1:length(files)){
                skips <- c(skips,getSkipFromGPR(files[i]))
        }
        if(length(unique(skips))!=1){
                cat("GPR files have different number of header rows!\nAdjusting GPR files\n")
                for(i in 1:length(files)){
                        dummy <- read.table(files[i],skip=skips[i],header=TRUE,comment.char="",sep="\t",check.names=FALSE)
                        #cat(colnames(dummy),"\n")
                        write.table(dummy,file=files[i],sep="\t",quote=FALSE,row.names=FALSE)
                }
        checkGPRs(files)    # just re-run the check once again...
        }
}

createM <- function(object,red.channels,green.channels){
        M <- matrix(ncol=length(red.channels),nrow=nrow(exprs(object)))
        weights <- getWeights(object)
        col.names <- NULL
        for(i in 1:ncol(M)){
                r.idx <- myGrep(colnames(exprs(object)),pattern=red.channels[i],exact.match=TRUE)
                g.idx <- myGrep(colnames(exprs(object)),pattern=green.channels[i],exact.match=TRUE)
                M[,i] <- getM(object,r=r.idx,g=g.idx)
                M[weights[,r.idx]==0 | weights[,g.idx]==0,i] <- NA
                col.names <- c(col.names,paste(red.channels[i]," vs ",green.channels[i],sep=""))
        }
        colnames(M) <- col.names
        M
}


createA <- function(object,red.channels,green.channels){
        A <- matrix(ncol=length(red.channels),nrow=nrow(exprs(object)))
        weights <- getWeights(object)
        col.names <- NULL
        for(i in 1:ncol(A)){
                r.idx <- myGrep(colnames(exprs(object)),pattern=red.channels[i],exact.match=TRUE)
                g.idx <- myGrep(colnames(exprs(object)),pattern=green.channels[i],exact.match=TRUE)
                A[,i] <- getA(object,r=r.idx,g=g.idx)
                A[weights[,r.idx]==0 | weights[,g.idx]==0,i] <- NA
                col.names <- c(col.names,paste(red.channels[i]," vs ",green.channels[i],sep=""))
        }
        colnames(A) <- col.names
        A
}


##
# to remove all those genes that have flags in 
excludeFromTest <- function(data,classlabels,weights=NULL,nr.present=2){
        if(!is.null(weights)){
                if(sum(dim(data)==dim(weights))==2){
                        data <- .removeFlagged(data,weights)
                }
        }
        result <- apply(data,MARGIN=1,.dummyFun,cl=classlabels,np=nr.present)
        result
}

.removeFlagged <- function(data,weights,flagged.value=NA){
        for(i in 1:ncol(data)){
                data[weights[,i]==0,i] <- flagged.value
        }
        data
}

.dummyFun <- function(x,cl,np){
        res <- sum(!is.na(x[cl==0])) < np | sum(!is.na(x[cl==1])) < np
        res
}

tukey.biweight.mod <- function(x,c=5,epsilon=1e-4,na.rm=TRUE){
        x <- x[!is.na(x)]
        ret <- tukey.biweight(x,c=c,epsilon=epsilon)
}

if( !isGeneric("filterOnVariance") )
        setGeneric("filterOnVariance", function(object,variance=0,array.names=NULL,remove.flagged=TRUE,v=TRUE,...)
        standardGeneric("filterOnVariance"))

setMethod("filterOnVariance","MadbSet",
        function(object,variance=0,array.names=NULL,remove.flagged=TRUE,v=TRUE,...){
                if(is.null(array.names)){
                        array.names <- colnames(exprs(object))
                }
                E <- exprs(object)[,array.names]
                if(remove.flagged){
                        weights <- getWeights(object)[,array.names]
                        for(i in 1:ncol(weights)){
                                E[weights[,i]==0,i] <- NA
                        }
                }
                sds <- apply(E,MARGIN=1,sd,na.rm=TRUE)
                quant <- quantile(sds,variance,na.rm=TRUE)
                sel <- sds >= quant
                sel[is.na(sel)] <- FALSE
                if(v){
                        cat(sum(sel)," features out of ",length(sel)," have a sd bigger than ",quant,"\n")
                }
                ret <- object
                exprs(ret) <- exprs(object)[sel,]
                ret@weights <- getWeights(object)[sel,]
                if(nrow(exprs(object))==nrow(object@genes)){
                        if(ncol(object@genes)==1){
                                G <- matrix(object@genes[sel,],ncol=1)
                                rownames(G) <- rownames(exprs(ret))
                                colnames(G) <- colnames(object@genes)
                                ret@genes <- G
                        }
                        else{
                                ret@genes <- object@genes[sel,]
                        }
                }
                ret
})


##
# this function is a generalization of Robert Gentlemens GOHyperG funciton from the GOstats package.
# it allows to submit all EntrezGenes present on the used microarray.
# in principal this function performs the same steps as the GOHyperG function, nevertheless the results
# are sligtly different (maybe because of the different annotation packages used:
#  mapping to GO terms using the affymetrix probe sets -> lokuslink -> GO terms (GOstats) or mapping the locuslink ids directly
#  to GO terms using the GO package like it is done here.)
calculateGOHyperG <- function(x,allLLs=NULL,chip=NULL,what="MF",use.all=FALSE,v=TRUE){
        require(GO)
        require(GOstats)
        if (chip == "fly.db0.db"){chip = "fly.db0"}
        if(is.null(allLLs) & is.null(chip) & !use.all){
                stop("You have to submit either all EntrezGene IDs present on the array used, or the Affymetrix GeneChip used, or select use.all=TRUE to use all EntrezGenes available in GO!\n\n")
        }
        if(use.all){
                ##
                # just that the code below throws no exception...
                allLLs <- x
        }
        ##
        # if chip was submitted retrieve all LLs from the specific chip
        if(!is.null(chip)){
                require(chip,character.only=TRUE) || stop("need data package",chip)
                allLLs <- as.character(unlist(as.list(get(paste(chip,"LOCUSID",sep=""),mode="environment"))))
        }
        ##
        # prepare the x (myLLs) and allLLs. remove duplicates, NAs, use only those x that are present in allLLs
        # remove possible LL. or LL strings (old way of writing LocusLink IDs)
        x <- sub("LL.","",x)
        x <- sub("LL","",x)
        allLLs <- sub("LL.","",allLLs)
        allLLs <- sub("LL","",allLLs)
        allLLs <- unique(allLLs)
        allLLs <- allLLs[!is.na(allLLs) & allLLs!=""]
        x <- unique(x)
        x <- x[!is.na(x) & x!=""]
        ##
        # just interested in those x that are present in allLLs
        x <- allLLs[match(x,allLLs)]
        x <- x[!is.na(x) & x!=""]
        
        ##
        # do some checking...
        if(length(x)==0 || is.na(x)){
                stop("None of the submitted interesting genes are present in the list of all EntrezGene (LocusLink) IDs submitted!\nPossible causes:\n-wrong list of genes (wrong IDs) submitted\n-EntrezGene IDs of the interesting genes and list of EntrezGene IDs of all genes on the array are not present on the same microarray\n")
        }
        
        
        
        ##
        # ok, now map x to GO, retrieve all GO terms where one of the x EntrezGene IDs is present (using GOALLLOCUSID from the GO package)
        GOLL <- as.list(get("GOALLLOCUSID",mode="environment"))     ## warning, GOALLLOCUSID is deprecated in newer versions, use GOALLENTREZID
        GOLL <- GOLL[!is.na(GOLL)] # just removing all the GO ids that are not mapped to any LL
        PresentGO <- sapply(GOLL,function(z){
                if(is.na(z) || length(z)==0)
                        return(FALSE)
                any(x %in% z)
                }
        )
        ##
        # another checking...
        if(max(PresentGO)==0){
                stop("The submitted list of EntrezGene IDs of the interesting genes can not be mapped to any of the GO terms.\nMaybe the list contained no EntrezGene (LocusLink) IDs?\n")
        }
        if(v){
                cat("The",length(x),"submitted EntrezGene IDs map to",sum(PresentGO),"GO terms in all tree categories.\n")
        }
        GOLL <- GOLL[PresentGO]
        ##
        # next restrict the GO terms to the selected category (specified by "what")
        goCat <- unlist(getGOOntology(names(GOLL)))
        goodGO <- goCat == what
        goodGO[is.na(goodGO)] <- FALSE
        if(v){
                cat(sum(goodGO),"of them are in the GO category we are interested in.\n")
        }
        mm <- match(names(goCat),names(GOLL))
        mm <- mm[goodGO]
        GOLL <- GOLL[mm]
        
        ##
        # next we have to remove all the EntrezGene IDs from each GO term, that was not present on the microarray used.
        if(!use.all){
                myGOLL <- sapply(GOLL,function(z){
                        if(length(z)==0 || is.na(z))
                                return(NA)
                        mylls <- unique(z[match(allLLs,z)])
                        mylls <- mylls[!is.na(mylls)]
                        mylls
                }
                )
                ##
                # removing all those GO terms that have no entries left (just in the case there are some...)
                empties <- sapply(myGOLL,function(z){length(z)==1 && is.na(z)})
                myGOLL <- myGOLL[!empties]
                GOLL <- GOLL[!empties]
        }
        else{
                ##
                # if use.all was selected, the proportion of the submitted EntrezGene IDs per GO term is compared to all
                # EntrezGenes mapped to that GO term.
                myGOLL <- GOLL
        }
        

        ##
        # the rest is the same as in the GOHyperG function...
        cLLs <- unique(unlist(myGOLL))
        nLL <- length(cLLs)                    # nr of unique LL IDs present in all GO terms
        goCounts <- sapply(myGOLL,length)      # nr of LL IDs present on the array per GO term
        #allCounts <- sapply(GOLL,length)       # nr of all LL IDs per GO term
        ##
        # some LL ids are represented more than once in a GO term...
        uniqueAllCounts <- sapply(GOLL,
                function(z){
                        if(length(z)==0 || is.na(z))
                                return(0)
                        lls <- unique(z)
                        lls <- lls[!is.na(lls)]
                        return(length(lls))
                })
        whGood <- x[x %in% cLLs]               # those LLs from the interesting ones that are present in the GO terms selected
        nInt <- length(whGood)
        if(nInt==0){
                warning("No interesting genes left")
        }
        
        MappingGOLL <- lapply(myGOLL,function(z){whGood[whGood %in% z]})
        useCts <- sapply(myGOLL,function(z){sum(whGood %in% z)})   # nr of interesting LLs per GO term
        ##
        # calculate the hypergeometric p values...
        pvs <- phyper(useCts-1,nInt,nLL-nInt,goCounts,lower.tail=FALSE)
        ord <- order(pvs)

        return(
                list(pvalues=pvs[ord],goCounts=goCounts[ord],intCounts=useCts[ord],allCounts=uniqueAllCounts[ord],GOLL=myGOLL[ord],intMapping=MappingGOLL[ord],numLL=nLL,numInt=nInt,intLLs=x)
        )
        
}


################
## get EntrezGenes for the specified GO term
getEntrezIDsFromGO <- function(GO, package="GO", entrezgenes){
        match.arg(package, c("GO", "humanLLMappings", "mouseLLMappings", "ratLLMappings"))
        require(package, character.only=TRUE)
        if(package=="GO"){
                ALLIDs <- mget(GO, env=GOALLENTREZID, ifnotfound=NA)
                if(!missing(entrezgenes)){
                        ## subset to the specified entrezgene ids.
                        ALLIDs <- lapply(ALLIDs, function(x){ x[ x %in% entrezgenes] })
                }
                return(ALLIDs)
        }
        if(package=="humanLLMappings"){
                ALLIDs <- mget(GO, env=humanLLMappingsGO2LL, ifnotfound=NA)
                if(!missing(entrezgenes)){
                        ## subset to the specified entrezgene ids.
                        ALLIDs <- lapply(ALLIDs, function(x){ x[ x %in% entrezgenes] })
                }
                return(ALLIDs)          
        }
        if(package=="mouseLLMappings"){
                ALLIDs <- mget(GO, env=mouseLLMappingsGO2LL, ifnotfound=NA)
                if(!missing(entrezgenes)){
                        ## subset to the specified entrezgene ids.
                        ALLIDs <- lapply(ALLIDs, function(x){ x[ x %in% entrezgenes] })
                }
                return(ALLIDs)          
        }
        if(package=="ratLLMappings"){
                ALLIDs <- mget(GO, env=ratLLMappingsGO2LL, ifnotfound=NA)
                if(!missing(entrezgenes)){
                        ## subset to the specified entrezgene ids.
                        ALLIDs <- lapply(ALLIDs, function(x){ x[ x %in% entrezgenes] })
                }
                return(ALLIDs)          
        }
}


#getLLGOmapping <- function(x,what="MF"){
#       gos <- mget(x, env=GOLOCUSID2GO,ifnotfound=NA)
#       if(length(gos)==1)
#               bd = is.na(gos[[1]])
#       else
#               bd = is.na(gos)
#       gos <- gos[!bd]
#       mymapping <- lapply(gos,getOntology,ontology=what)
#       emptygos <- sapply(mymapping,function(x) length(x)==0 )
#       mymapping <- mymapping[!emptygos]
#       mymapping
#}

#getGOFromLL <- function(x,Ontology="MF"){
#       require(GO) || stop("no GO library")
#       match.arg(Ontology, c("MF", "BP", "CC"))
#       goterms <- mget(x, env = GOLOCUSID2GO, ifnotfound=NA)
#       if(length(goterms)==1){
#               goterms <- goterms[!is.na(goterms[[1]])]
#       }
#       else{
#               goterms <- goterms[!is.na(goterms)]
#       }
#       goterms <- lapply(goterms, function(x) x[sapply(x,
#                      function(x) {if (is.na(x$Ontology) )
#                                          return(FALSE)
#                                      else
#                                          x$Ontology == Ontology})])
#       goids <- NULL
#       if(length(goterms)==0){
#               stop("The submitted list of EntrezGene IDs of the interesting genes can not be mapped to any of the GO terms.\nMaybe the list contained no EntrezGene (LocusLink) IDs?\n")
#       }
#       for(i in 1:length(goterms)){
#               goids <- c(goids,unique(unlist(sapply(goterms[[i]], function(x) x$GOID))))
#       }
#       goids
#}
#

##
# check if package is available
isAnnotationAvailable <- function(chip,develOK=FALSE){
        ok = FALSE
        if (chip == "fly.db0.db"){chip  = "fly.db0"}
       if(!require("reposTools",character.only=TRUE)){
               warning("reposTools package is not installed.\n")
               if(!require(chip,character.only=TRUE)){
                       ok = FALSE
               }
               else{
                       ok = TRUE
               }
       }
       else{
                if(!require(chip,character.only=TRUE)){
             source("http://www.bioconductor.org/getBioC.R")
                      library(reposTools)
                      # have to install the package
                      tryCatch(install.packages2(pkg=chip,develOK=develOK,getAllDeps=TRUE),error=return(FALSE))
                      try(getBioC(chip))
                      if(!require(chip,character.only=TRUE)){
                          ok=FALSE
                      }
                       else{
                           ok=TRUE
                       }
                }
                else{
                        ok = TRUE
                }
       }
        ok
}

saveSAMResults <- function(SAM,table,delta,filename){
        dummy <- print(SAM)
        for(i in 1:nrow(dummy)){
                saveSAMResult(SAM,table,delta=dummy[i,"Delta"],filename=filename)
        }
}

saveSAMResult <- function(SAM, table,delta,filename){
        sum.SAM <- (summary(SAM,delta))@mat.sig
        if(!is.null(sum.SAM) & nrow(sum.SAM)>0){
                table.sub <- table[sum.SAM[,"Row"],]
                if(class(table)=="aafTable"){
                        for(i in 2:(ncol(sum.SAM))){
                                table.sub <- merge(table.sub,aafTable(sum.SAM[,i],colnames=colnames(sum.SAM)[i]))
                        }
                        Filename <- checkExistantFile(paste(filename,"-delta-",delta,".txt",sep=""))
                        saveText(table.sub,filename=Filename)
                }
                else{
                        if(nrow(sum.SAM)==1){
                                table.sub <- matrix(table.sub,nrow=1)
                                colnames(table.sub) <- colnames(table)
                        }
                        table.sub <- cbind(as.matrix(table.sub),as.matrix(sum.SAM[,2:ncol(sum.SAM)]))
                        Filename <- checkExistantFile(paste(filename,"-delta-",delta,".txt",sep=""))
                        write.table(table.sub,file=Filename,sep="\t",quote=FALSE,row.names=FALSE)
                }
                cat(nrow(table.sub),"significant genes (delta=",delta,") were saved to the file:",Filename,"\n")
        }
        else{
                cat("No significant genes with a delta of",delta,"\n")
        }
}

###
# fills missing spots
fillMissingSpots <- function(object){
        ret <- object
        if(class(ret)=="RGList"){
                if(!is.null(ret$printer)){
                        layout <- object$printer
                        nr.row <- layout$ngrid.r * layout$ngrid.c * layout$nspot.r * layout$nspot.c
                        missing <- (nr.row) - nrow(object$R)
                        if(missing!=0){
                                cat("Have to fill in",missing,"missing spots\n")
                                index <- as.numeric(apply(object$genes,MARGIN=1,FUN=.getIndexFromBCR,nspot.r=layout$nspot.r,nspot.c=layout$nspot.c))
                                R <- matrix(ncol=ncol(object$R),nrow=nr.row)
                                colnames(R) <- colnames(object$R)
                                rownames(R) <- 1:nrow(R)

                                ## check if length of index and nrow of R are the same (<- problem with custom GAL files...)
                                if(length(index)!=nrow(object$R)){
                                        stop("ERROR: number of data values does not match number of annotations!\nPlease try again without specifying a GAL file!")
                                }
                                R[index,] <- object$R
                                ret$R <- R
                                R[index,] <- object$G
                                ret$G <- R
                                if(!is.null(object$Rb)){
                                        R[index,] <- object$Rb
                                        ret$Rb <- R
                                        R[index,] <- object$Gb
                                        ret$Gb <- R
                                }
                                if(!is.null(object$weights)){
                                        R[,] <- 0
                                        R[index,] <- object$weights
                                        ret$weights <- R
                                }
                                else{
                                        R[,] <- 0
                                        R[index,] <- 1
                                        ret$weights <- R
                                }
                                rm(R)
                                g <- gc()
                                Genes = data.frame(matrix(ncol=ncol(object$genes),nrow=nr.row))
                                colnames(Genes) <- colnames(object$genes)
                                rownames(Genes) <- 1:nrow(Genes)
#                               cat("nrow Genes",nrow(Genes),"ncol Genes",ncol(Genes),"length index",length(index),"\n")
                                Genes[index,1:ncol(Genes)] = object$genes[,1:ncol(Genes)]
#                               cat("genes filled\n","ncol= ",ncol(Genes)," nrow: ",nrow(Genes),"\n")
                                ret$genes = NULL
                                ret$genes = Genes
#                               ret$genes = Genes
#                               cat("genes set\n")
                                rm(Genes)
                                g <- gc()
                        }
                }
        }
        ret
}

.getIndexFromBCR <- function(x,nspot.c,nspot.r){
        ((as.numeric(x["Block"])-1)*nspot.c*nspot.r) + ((as.numeric(x["Row"])-1)*nspot.c) + as.numeric(x["Column"])
}

if( !isGeneric("drawBoxplot") )
        setGeneric("drawBoxplot", function(x,new.plot=TRUE,xlim=NULL,at=0,exclude.flagged=FALSE,...)
        standardGeneric("drawBoxplot"))

setMethod("drawBoxplot","MadbSet",
        function(x,new.plot=TRUE,xlim=NULL,at=0,exclude.flagged=FALSE,...){
                if(exclude.flagged){
                        W <- getWeights(x)
                        for(i in 1:ncol(exprs(x))){
                                exprs(x)[W[,i]==0,i] <- NA
                        }
                }
                boxplots(x=exprs(x),new.plot=new.plot,xlim=xlim,at=at,...)
        }
)

boxplots <- function(x,new.plot=TRUE,xlim=NULL,at=0,col=NULL,log2.transform=FALSE,ylim=NULL,...){
        data <- x
#       if(class(x)=="exprSet" | class(x)=="EexprSet"){
#               data <- exprs(x)
#       }
        if(log2.transform){
                data <- log2(data)
        }
        if(is.null(ncol(data))){
                data <- as.matrix(data)
        }
        ## remove those nasty infinite values
        for(i in 1:ncol(data)){
                data[is.infinite(data[,i]),i] <- NA
        }

        CN <- colnames(data)
        if(is.null(CN)){
                CN <- 1:ncol(data)
        }
        if(is.null(col)){
                col=rep(0,ncol(data))
        }
        if(length(col)!=ncol(data)){
                col=rep(col[1],ncol(data))
        }
        if(new.plot){
        if(is.null(xlim)){
                xlim=c(0.5,(ncol(data)+0.5))
        }
        par(omi=0.7*c(min(2,max(strwidth(CN,units="i"))),0,0,0),cex.axis=0.7)
        if(is.null(ylim)){
                ylim=c(min(data,na.rm=TRUE),max(data,na.rm=TRUE))
        }
        plot(1,1,pch=NA,ylab="expression",xaxt="n",ylim=ylim,xlim=xlim,xlab=NA)
        }
        axis(side=1,labels=CN,at=(1+at):(ncol(data)+at),las=3)
        for(i in 1:ncol(data)){
                boxplot(data[,i],at=(i+at),add=TRUE,col=col[i],range=0,...)
        }
}


plotRankedVariance <- function(data,subset=NULL,cut){
        if(!is.null(subset)){
                data <- data[,subset]
        }
        if(!.is.log(data)){
                data <- log2(data)
        }
        SD <- apply(data,MARGIN=1,sd,na.rm=TRUE)
        CUT <- quantile(SD,cut,na.rm=TRUE)
        plot(sort(SD),pch=16,cex=0.33,ylab="sorted variance")
        abline(h=CUT,col=2)
        text((length(SD)/20),(CUT+0.2),label=format(CUT,digits=4),col=2)
}

plotSortedPValues <- function(x,top=NULL,supported.cols=c("rawp","adjp","minP","maxT","BH","BY","Bonferroni","Hochberg","Holm","SidakSS","SidakSD"),...){
        if(is.null(nrow(x))){
                stop("Please submit a matrix, data.frame or instance of aafTable!\n")
        }
        if(!is.null(top)){
                top <- min(c(nrow(x),top))
        }
        else{
                top <- nrow(x)
        }
        Columns <- supported.cols[supported.cols %in% colnames(x)]
        if(length(Columns)==0){
                stop("The submitted matrix must contain at least one column with p-values that is named according to one of the supported column names!")
        }
        if(class(x)=="aafTable"){
                require(annaffy)
                #have to do some transformations...
                x.new <- matrix(ncol=length(Columns),nrow=nrow(x))
                colnames(x.new) <- Columns
                for(i in Columns){
                        x.new[,i] <- as.numeric(getText(x[[i]]))
                }
                x <- x.new
        }
        #ok, now i can start drawing...
        x <- apply(x[,Columns],MARGIN=2,sort,na.last=TRUE)
        Max <- max(x[1:top,],na.rm=TRUE)

        plot(x[1:top,Columns[1]],xlab="number of significant genes",ylab="sorted p-value",col=1,ylim=c(0,Max),main=paste("Top",top,"sorted p-values"),...)
        if(length(Columns)>1){
                for(i in 2:length(Columns)){
                        points(x[1:top,Columns[i]],col=i,...)
                }
        }
        legend(top-top/4,Max,Columns,col=1:length(Columns),pch=16)
}

.is.log <- function(data,maxi=100){
        data <- data[!is.na(data)]
        data <- data[!is.infinite(data)]
#       cat("Try to check wheter or not the values are in log scale...")
        if(max(data>maxi)){
                log.values = FALSE
#               cat("signals are not in log2 scale, transforming...\n")
        }
        else{
                log.values = TRUE
#               cat("signals are already in log2 scale \n")
        }
        log.values
}

toNumeric <- function(x,Sub=","){
        if(!is.null(dim(x))){
                RN <- rownames(x)
                if(class(x[1,])=="numeric"){
                        return(x)
                }
                if(!is.null(Sub)){
                        for(s in Sub){
                                x <- apply(x,MARGIN=2,sub,pattern=s,replacement="")
                        }
                }
                x <- tryCatch(apply(x,MARGIN=2,as.numeric),error=function(error){stop("ERROR: Can not perform calculations on non-numeric data! Please submit only numeric values.")},warning=function(warning){stop("ERROR: Can not perform calculations on non-numeric data! Please submit only numeric values.")})
                rownames(x) <- RN
                x
        }
}

####
## function similar to PWamat function from the annotate package,
## generates a matrix with rows unique LL ids, columns KEGG pathways
## containing 1 if gene can be associated with the pathway, 0 otherwise.
makeLL2KEGGMappingMatrix <- function(x,nrGenes=0,verbose=TRUE){
        if(is.null(x)){
                stop("No EntrezGene (LocusLink) IDs have been submitted!\n")
        }
        require(KEGG)
        x <- unique(as.character(x))
        paths <- mget(x,KEGGEXTID2PATHID,ifnotfound=NA)
        not_found <- sapply(paths,function(x){if(length(x)==1 && is.na(x)) TRUE else FALSE})
        paths <- paths[!not_found]
        if(verbose){
                cat(length(paths),"EntrezGene IDs from the inputlist can be associated with at least one KEGG pathway.\n",sum(not_found),"IDs can not be associated with KEGG pathways.\n")
        }
        unique_paths <- unique(unlist(paths))
        Test <- sapply(paths,function(x){
                mtch <- match(x,unique_paths)
                res <- rep(0,length(unique_paths))
                res[mtch] <- 1
                res
        })
        Test <- t(Test)
        rownames(Test) <- x
        colnames(Test) <- unique_paths
        if(nrGenes>0){
                to_exclude <- colSums(Test) < nrGenes
                if(verbose){
                        cat(sum(to_exclude),"Pathways are excluded, because less than",nrGenes,"genes can be mapped to them.\n")
                }
                Test <- Test[,!to_exclude]
        }
        return(Test)
}