R/noise.bench.R

In netbenchmark: Benchmarking of several gene network inference methods

noise.bench <- function(methods="all.fast",datasources.names="all",
    experiments=150,eval="AUPR",no.topedges=20,datasets.num=3,
    local.noise=seq(0,100,len=3),global.noise=0,noiseType="normal",
    sym=TRUE,seed=NULL)
{
    options(warn=1)
    # set random number generator seed if seed is given
    if (!is.null(seed)){
        set.seed(seed)
    }else{
        seed <- as.double(Sys.time())
        set.seed(seed)
    }
    if(all("all.fast" %in% tolower(methods))) {
        methods <- c("aracne.wrap","c3net.wrap","clr.wrap","GeneNet.wrap",
            "mrnetb.wrap","mutrank.wrap","pcit.wrap","zscore.wrap",
            methods[tolower(methods)!="all.fast"])
    }else if(all("all" %in% tolower(methods))) {
        methods <- c("aracne.wrap","c3net.wrap","clr.wrap","GeneNet.wrap",
            "Genie3.wrap","mrnet.wrap","mutrank.wrap","mrnetb.wrap",
            "pcit.wrap","zscore.wrap",methods[tolower(methods)!="all"])
    }
    if(length(datasources.names)==1){
        if (tolower(datasources.names)=="all"){
            datasources.names <- c("rogers1000","syntren1000","syntren300",
                "gnw1565","gnw2000")
        }else{
            if(tolower(datasources.names)=="toy"){
                datasources.names <- "toy"
            }
        }
    }
    if(length(local.noise)!=length(global.noise)){
        if(length(local.noise)==1){
            points <- length(global.noise)
            local.noise <- rep(local.noise,points)
        }else if(length(global.noise)==1){
            points <- length(local.noise)
            global.noise <- rep(global.noise,points)
        }else{
            stop("Error: mismatch in lengths of local.noise and global.noise")
        }
    }
    nmeths <- length(methods)
    ndata <- length(datasources.names)
    results <- as.data.frame(matrix(0,points*ndata,nmeths+4))
    pval <- as.data.frame(matrix(0,points*ndata,nmeths+4))
    Availabledata <- eval(parse(text="Availabledata"))
    seeds <- as.list(round(runif(length(Availabledata),max=1e9)))
    names(seeds) <- Availabledata
    rown <- character()
    if (!all(datasources.names %in% Availabledata)) stop("unknown datasource")
    for(n in seq_len(ndata)){
        message(paste("Datasource:",datasources.names[n],"\n"))
        datasource <- eval(parse(text=paste(datasources.names[n],
            ".data",sep="")))
        true.net <- eval(parse(text=paste(datasources.names[n],".net",
            sep="")))
        s <- dim(datasource)
        ngenes <- s[2] #number of genes in the network
        npos <- sum(true.net) #number of true links in the network
        nlinks <- ngenes^2-ngenes #number of posible links in the network
        if(sym){
            nlinks <- nlinks/2
        }
        no.edges <- round(nlinks*no.topedges/100)
        m <- matrix(0,points,nmeths+1)
        pval.table <- matrix(0,points,nmeths+1)
        tp.local.mat <- matrix(0,no.edges,nmeths+1)
        colnames(tp.local.mat) <- c(methods,"rand")
        l.seed <- eval(parse(text=paste("seeds$",datasources.names[n])))
        set.seed(l.seed)
        spd <- datasource.subsample(datasource,experiments=experiments,
            datasets.num = datasets.num,local.noise = 0,global.noise = 0)
        for(i in seq_len(points)){
            m.local <- matrix(0,datasets.num,nmeths+1)
            rdata <- vector('list',datasets.num)
            for(k in seq_len(datasets.num)){
                if(local.noise[i]!=0){
                    rdata[[k]] <- apply(spd[[k]],2,.cont,
                        noise=local.noise[i],noiseType=noiseType)
                }else{
                    rdata[[k]] <- spd[[k]]
                }
                if(global.noise[i]!=0){
                    sds <- apply(spd[[k]], 2, sd)
                    if(noiseType=="normal"){
                        Gnoise <- matrix(rnorm(s[1]*s[2],mean=0,
                            sd=mean(sds)*global.noise[i]/100), s[1], s[2])
                    }
                    if(noiseType=="lognormal"){
                        Gnoise <- matrix(rlnorm(s[1]*s[2],meanlog=0,
                            sdlog=mean(sds)*global.noise[i]/100), s[1], s[2])
                    }
                    rdata[[k]] <- rdata[[k]]+Gnoise
                }
            }
            for(j in seq_len(nmeths)){
                message(paste(methods[j],"\n"))
                for(k in seq_len(datasets.num)){
                    net <- do.call(methods[j],list(rdata[[k]]))
                    r <- evaluate(net,true.net,extend=no.edges,sym=sym)
                    tp.local.mat[,j] <- tp.local.mat[,j]+r[1:no.edges,"TP"]
                    if(tolower(eval)=="no.truepos"){
                        m.local[k,j] <- mean(r[1:no.edges,"TP"])
                    }else if (tolower(eval)== "aupr"){
                        m.local[k,j] <- aupr(r,no.edges)
                    }else if (tolower(eval)== "auroc"){
                        m.local[k,j] <- auroc(r,no.edges)
                    }else stop("unknown evaluation metric")
                }
                tp.local.mat[,j] <- tp.local.mat[,j]/datasets.num
            }
            m[i,] <- apply(m.local,2,mean)
            M <- which.max(m[i,])
            precision <- tp.local.mat/matrix(rep(1:no.edges,nmeths+1),
                no.edges)
            for(j in seq_len(nmeths)){
                if(j!=M){
                    aux <- wilcox.test(m.local[,j],m.local[,M])
                    pval.table[i,j]=aux[[3]]
                }else{
                    pval.table[i,j]=1
                }
            }
            rand.net <- matrix(runif(ngenes^2),ngenes,ngenes)
            diag(rand.net) <- 0
            colnames(rand.net) <- colnames(net)
            rownames(rand.net) <- colnames(net)
            r <- evaluate(rand.net,true.net,extend=no.edges,sym=sym)
            tp.local.mat[,nmeths+1] <- r[1:no.edges,"TP"]
            precision <- tp.local.mat/matrix(rep(1:no.edges,nmeths+1),
                no.edges)
            if(tolower(eval)=="no.truepos"){
                m[i,nmeths+1]=mean(r[1:no.edges,"TP"])
            }else if (tolower(eval)== "aupr"){
                m[i,nmeths+1]=aupr(r,no.edges)
            }else if (tolower(eval)== "auroc"){
                m[i,nmeths+1]=auroc(r,no.edges)
            }
            aux <- wilcox.test(precision[,nmeths+1],precision[,M])
            pval.table[i,nmeths+1]=aux[[3]]
        }
        rown <- c(rown,rep(datasources.names[n],points))
        results[(1:points)+(n-1)*points,1] <- rep(datasources.names[n],
            points)
        results[(1:points)+(n-1)*points,2] <- local.noise
        results[(1:points)+(n-1)*points,3] <- global.noise
        results[(1:points)+(n-1)*points,4:(nmeths+4)] <- m  
        pval[(1:points)+(n-1)*points,1] <- rep(datasources.names[n],points)
        pval[(1:points)+(n-1)*points,2] <- local.noise
        pval[(1:points)+(n-1)*points,3] <- global.noise
        pval[(1:points)+(n-1)*points,4:(nmeths+4)] <- pval.table
    }
    colnames(results) <- c("Datasource","local.noise","global.noise",
        methods,"rand")
    colnames(pval) <- c("Datasource","local.noise","global.noise",
        methods,"rand")
    list("results"=results,"pval"=pval,"seed"=seed)
}

.cont <- function(x,noise=0,noiseType="normal"){
    s.d <- noise*sd(x)/100
    if(noiseType=="normal"){
        n <- rnorm(length(x),mean=0,sd=s.d)
    }
    if(noiseType=="lognormal"){
        n <- rlnorm(length(x), meanlog = 0, sdlog = s.d)
    }
    return(x+n)
}