R/quality_control.R

Defines functions controlDensityPerPlate controlDensityPerScreen controlDensity compareReplicateSDPerScreen compareReplicateSD compareReplicaPlates compareReplicates spatialDistrib makeBoxplot4PlateType makeBoxplotPerPlate makeBoxplotPerScreen makeBoxplotControlsPerScreen makeBoxplotControlsPerPlate makeBoxplotControls replicatesCV replicatesSpearmancor plotQQperscreen plotQQperplate plotQQ plotHistoPerscreen plotHistoPerplate plotHisto plotControlHistoPerscreen plotControlHistoPerplate plotControlHisto DRQualControl SNRQualControl ZPRIMEQualControl percCellQualControl numCellQualControl dynamicRange ZPrime ylabCompute ZScorePlotTwo plotBar plotBarGeruest channelPlot discardWells discardLabtek

Documented in channelPlot compareReplicaPlates compareReplicates compareReplicateSD compareReplicateSDPerScreen controlDensity controlDensityPerPlate controlDensityPerScreen discardLabtek discardWells DRQualControl makeBoxplot4PlateType makeBoxplotControls makeBoxplotControlsPerPlate makeBoxplotControlsPerScreen makeBoxplotPerPlate makeBoxplotPerScreen numCellQualControl percCellQualControl plotBar plotControlHisto plotControlHistoPerplate plotControlHistoPerscreen plotHisto plotHistoPerplate plotHistoPerscreen plotQQ plotQQperplate plotQQperscreen replicatesCV replicatesSpearmancor SNRQualControl spatialDistrib ZPRIMEQualControl ZScorePlotTwo

discardLabtek<-function(data, screenNr, labtekNr){
    subset<-createSubset(data, data$ScreenNb, screenNr)
    IX<-indexSubset(data$ScreenNb, screenNr)
    subset$SpotType[which(subset$LabtekNb == labtekNr)]<-(-1)
    data[IX, ]<-subset
    invisible(data)
}


discardWells<-function(data, screenNr, labtekNr, vecPositions){
    subset<-createSubset(data, data$ScreenNb, screenNr)
    IX<-indexSubset(data$ScreenNb, screenNr)

    subsubset<-createSubset(subset, subset$LabtekNb, labtekNr)
    IXX<-indexSubset(subset$LabtekNb, labtekNr)

    subsubset$SpotType[vecPositions]<-(-1)

    subset[IXX, ]<-subsubset
    data[IX, ]<-subset
    
    invisible(data)
}


channelPlot <- function(header, dataset, vecOfChannels, flag, plotTitle,
showPlot, smSpan=2/3){

  dataset<-dataset[which(dataset$SpotType!=-1), ]

  if (flag == 0){

    headerTemp<-strsplit(header[1], ",")
    plotName<-paste(headerTemp[[1]][2], plotTitle, sep="_")

    combiCounter<-0
    for (j in 1:(length(vecOfChannels)-1)){
        for (k in (j+1):length(vecOfChannels)){
            combiCounter<-combiCounter+1
            ch1<-vecOfChannels[j]
            ch2<-vecOfChannels[k]

            dataset<-dataset[which(!is.na(dataset[[get("ch1")]])), ]
            dataset<-dataset[which(!is.na(dataset[[get("ch2")]])), ]

            tt<-paste(plotTitle, " - Channel ", j, " vs. Channel ", k, sep="")
            x<-dataset[[get("ch1")]]
            y<-dataset[[get("ch2")]]

            if (showPlot == 1){
                if(interactive()){
                    x11()
                    plot(x, y, main=tt, xlab="Ch1", ylab="Ch2", cex.main=0.7)
                    lines(lowess(x, y, f=smSpan), col="red")
                }
            }
            pdf(paste(plotName, "(", combiCounter, ").pdf", sep=""))
                plot(x, y, main=tt, xlab="Ch1", ylab="Ch2", cex.main=0.7)
                lines(lowess(x, y, f=smSpan), col="red")
            dev.off()
            png(paste(plotName, "(", combiCounter, ").png", sep=""), width=300, height=300)
                plot(x, y, main=tt, xlab="Ch1", ylab="Ch2", cex.main=0.7)
                lines(lowess(x, y, f=smSpan), col="red")
            dev.off()
        }
    }
    return(plotName)
  }

  if (flag == 1){

    minOfScreens<-min(dataset$ScreenNb)
    numOfScreens<-max(dataset$ScreenNb)

    headerTemp<-strsplit(header[1], ",")
    plotName<-paste(headerTemp[[1]][2], plotTitle, sep="_")

    for (i in minOfScreens:numOfScreens){
        if (length(which(dataset$ScreenNb == i))>0){
        
            subset<-dataset[which(dataset$ScreenNb == i), ]

            combiCounter<-0
            for (j in 1:(length(vecOfChannels)-1)){
                for (k in (j+1):length(vecOfChannels)){
                    combiCounter<-combiCounter+1
                    ch1<-vecOfChannels[j]
                    ch2<-vecOfChannels[k]

                    subset<-subset[which(!is.na(subset[[get("ch1")]])), ]
                    subset<-subset[which(!is.na(subset[[get("ch2")]])), ]
                    
                    tt<-paste(plotTitle, " - Channel ", j, " vs. Channel ", k, sep="")
		    x<-subset[[get("ch1")]]
                    y<-subset[[get("ch2")]]

                    if (showPlot == 1){
                        if(interactive()){
                            x11()
                            plot(x, y, main=tt, xlab="Ch1", ylab="Ch2", cex.main=0.7)
                            lines(lowess(x, y, f=smSpan), col="red")
                        }
                    }
                    pdf(paste(plotName, "(", combiCounter, ") (Exp. ", i, ").pdf", sep=""))
                        plot(x, y, main=tt, xlab="Ch1", ylab="Ch2", cex.main=0.7)
                        lines(lowess(x, y, f=smSpan), col="red")
                    dev.off()
                    png(paste(plotName, "(", combiCounter, ") (Exp. ", i, ").png", sep=""), 
                    width=300, height=300)
                        plot(x, y, main=tt, xlab="Ch1", ylab="Ch2", cex.main=0.7)
                        lines(lowess(x, y, f=smSpan), col="red")
                    dev.off()
                }
            }
        }
    }
    return(list(plotName, minOfScreens, numOfScreens))
  }

  if (flag == 2){

    numOfScreens<-max(dataset$ScreenNb)
    minOfScreens<-min(dataset$ScreenNb)

    for (j in minOfScreens:numOfScreens){
        if (length(which(dataset$ScreenNb == j))>0){
        
            subset<-dataset[which(dataset$ScreenNb == j), ]
            
            minOfPlates<-min(subset$LabtekNb)
            numOfPlates<-max(subset$LabtekNb)

            headerTemp<-strsplit(header[1], ",")
            plotName<-paste(headerTemp[[1]][2], plotTitle, sep="_")

            for (i in minOfPlates:numOfPlates){
                if (length(which(subset$LabtekNb == i))>0){
                    
                    subsubset<-subset[which(subset$LabtekNb == i), ]
                    
                    combiCounter<-0
                    for (o in 1:(length(vecOfChannels)-1)){
                        for (p in (o+1):length(vecOfChannels)){
                            combiCounter<-combiCounter+1
                            ch1<-vecOfChannels[o]
                            ch2<-vecOfChannels[p]

                            subsubset<-subsubset[which(!is.na(subsubset[[get("ch1")]])), ]
                            subsubset<-subsubset[which(!is.na(subsubset[[get("ch2")]])), ]
                            
                            tt<-paste(plotTitle, " - Channel ", o, " vs. Channel ", p, sep="")
   		            x<-subsubset[[get("ch1")]]
                            y<-subsubset[[get("ch2")]]

                            if (showPlot == 1){
                                if(interactive()){
                                    x11()
                                    plot(x, y, main=tt, xlab="Ch1", ylab="Ch2", cex.main=0.7)
                                    lines(lowess(x, y, f=smSpan), col="red")
                                }
                            }
                            pdf(paste(plotName, "(", combiCounter, ")_Exp_", j, "_Plate_", i, "_.pdf", 
                            sep=""))
                                plot(x, y, main=tt, xlab="Ch1", ylab="Ch2", cex.main=0.7)
                                lines(lowess(x, y, f=smSpan), col="red")
                            dev.off()
                            png(paste(plotName, "(", combiCounter, ")_Exp_", j, "_Plate_", i, "_.png", 
                            sep=""), width=300, height=300)
                                plot(x, y, main=tt, xlab="Ch1", ylab="Ch2", cex.main=0.7)
                                lines(lowess(x, y, f=smSpan), col="red")
                            dev.off()
                        }
                    }
                }
            }
        }
    }
    return(list(plotName, c(minOfScreens, numOfScreens), 
    c(minOfPlates, numOfPlates)))
  }
}


plotBarGeruest <- function(dataset, col4val){

  wells <- nrow(dataset)
  pos<-which(dataset$SpotType == 1)
  neg<-which(dataset$SpotType == 0)
  colors <- rep("black", wells)
  colors[pos] <- "green"
  colors[neg] <- "red"
  data<-dataset[[get("col4val")]]
  med <- median(data, na.rm=T)
  dev <- mad(data, na.rm=T)
  if (is.na(med))
    med <- 0
  if (is.na(dev))
    dev <- 1
  minim <- max(min(as.numeric(data), na.rm=TRUE), med-8*dev)
  maxim <- min(max(as.numeric(data), na.rm=TRUE), med+8*dev)
  if (minim  ==  Inf)
    minim <- -5
  if (maxim  ==  -Inf)
    maxim <- -5

  invisible(list(minim, maxim, data, colors, med, dev))
}

plotBar <- function(header, dataset, col4val, flag, plotTitle, showPlot){

  dataset<-dataset[which(dataset$SpotType!=-1), ]


  if (flag == 0){
      plotData<-plotBarGeruest(dataset, col4val)

    minim<-plotData[[1]]
    maxim<-plotData[[2]]
    data<-plotData[[3]]
    colors<-plotData[[4]]
    med<-plotData[[5]]
    dev<-plotData[[6]]
    
    headerTemp<-strsplit(header[1], ",")
    plotName<-paste(headerTemp[[1]][2], plotTitle, sep="_")
    
    if (showPlot == 1){
        if(interactive()){
            boxplot(data.frame(t(data)), border=colors, col=colors, ylim=c(minim, maxim))
            abline(med, 0, col="blue")
            abline(med+dev, 0, col="green")
            abline(med-dev, 0, col="green")
            abline(med+2*dev, 0, col="red")
            abline(med-2*dev, 0, col="red")
        }
    }
    pdf(paste(plotName, ".pdf", sep=""))
        boxplot(data.frame(t(data)), border=colors, col=colors, ylim=c(minim, maxim))
        abline(med, 0, col="blue")
        abline(med+dev, 0, col="green")
        abline(med-dev, 0, col="green")
        abline(med+2*dev, 0, col="red")
        abline(med-2*dev, 0, col="red")
    dev.off()
    png(paste(plotName, ".png", sep=""), width=300, height=300)
        boxplot(data.frame(t(data)), border=colors, col=colors, ylim=c(minim, maxim))
        abline(med, 0, col="blue")
        abline(med+dev, 0, col="green")
        abline(med-dev, 0, col="green")
        abline(med+2*dev, 0, col="red")
        abline(med-2*dev, 0, col="red")
    dev.off()

    return(plotName)
  }

  if (flag == 1){

    minOfScreens<-min(dataset$ScreenNb)
    numOfScreens<-max(dataset$ScreenNb)

    headerTemp<-strsplit(header[1], ",")
    plotName<-paste(headerTemp[[1]][2], plotTitle, sep="_")

    for (i in minOfScreens:numOfScreens){
        if (length(which(dataset$ScreenNb == i))>0){
        
            subset<-dataset[which(dataset$ScreenNb == i), ]

            plotData<-plotBarGeruest(subset, col4val)
            minim<-plotData[[1]]
            maxim<-plotData[[2]]
            data<-plotData[[3]]
            colors<-plotData[[4]]
            med<-plotData[[5]]
            dev<-plotData[[6]]

            if (showPlot == 1){
                if(interactive()){
                    x11()
                    boxplot(data.frame(t(data)), border=colors, col=colors, ylim=c(minim, maxim))
                    abline(med, 0, col="blue")
                    abline(med+dev, 0, col="green")
                    abline(med-dev, 0, col="green")
                    abline(med+2*dev, 0, col="red")
                    abline(med-2*dev, 0, col="red")
                }
            }

            pdf(paste(plotName, "(Exp. ", i, ").pdf", sep=""))
                boxplot(data.frame(t(data)), border=colors, col=colors, ylim=c(minim, maxim))
                abline(med, 0, col="blue")
                abline(med+dev, 0, col="green")
                abline(med-dev, 0, col="green")
                abline(med+2*dev, 0, col="red")
                abline(med-2*dev, 0, col="red")
            dev.off()
            png(paste(plotName, "(Exp. ", i, ").png", sep=""), width=300, height=300)
                boxplot(data.frame(t(data)), border=colors, col=colors, ylim=c(minim, maxim))
                abline(med, 0, col="blue")
                abline(med+dev, 0, col="green")
                abline(med-dev, 0, col="green")
                abline(med+2*dev, 0, col="red")
                abline(med-2*dev, 0, col="red")
            dev.off()
        }
    }
    return(list(plotName, minOfScreens, numOfScreens))
  }

  if (flag == 2){

    numOfScreens<-max(dataset$ScreenNb)
    minOfScreens<-min(dataset$ScreenNb)

    for (j in minOfScreens:numOfScreens){
        if (length(which(dataset$ScreenNb == j))>0){
        
            subset<-dataset[which(dataset$ScreenNb == j), ]
            
            minOfPlates<-min(subset$LabtekNb)
            numOfPlates<-max(subset$LabtekNb)

            headerTemp<-strsplit(header[1], ",")
            plotName<-paste(headerTemp[[1]][2], plotTitle, sep="_")

            for (i in minOfPlates:numOfPlates){
                if (length(which(subset$LabtekNb == i))>0){
                    
                    subsubset<-subset[which(subset$LabtekNb == i), ]
                    plotData<-plotBarGeruest(subsubset, col4val)
                    minim<-plotData[[1]]
                    maxim<-plotData[[2]]
                    data<-plotData[[3]]
                    colors<-plotData[[4]]
                    med<-plotData[[5]]
                    dev<-plotData[[6]]

                    if (showPlot == 1){
                        if(interactive()){
                            x11()
                            boxplot(data.frame(t(data)), border=colors, col=colors, ylim=c(minim, maxim))
                            abline(med, 0, col="blue")
                            abline(med+dev, 0, col="green")
                            abline(med-dev, 0, col="green")
                            abline(med+2*dev, 0, col="red")
                            abline(med-2*dev, 0, col="red")
                        }
                    }

                    pdf(paste(plotName, "Exp", j, "Plate", i, ".pdf", sep="_"))
                        boxplot(data.frame(t(data)), border=colors, col=colors, ylim=c(minim, maxim))
                        abline(med, 0, col="blue")
                        abline(med+dev, 0, col="green")
                        abline(med-dev, 0, col="green")
                        abline(med+2*dev, 0, col="red")
                        abline(med-2*dev, 0, col="red")
                    dev.off()
                    png(paste(plotName, "Exp", j, "Plate", i, ".png", sep="_"), width=300, height=300)
                        boxplot(data.frame(t(data)), border=colors, col=colors, ylim=c(minim, maxim))
                        abline(med, 0, col="blue")
                        abline(med+dev, 0, col="green")
                        abline(med-dev, 0, col="green")
                        abline(med+2*dev, 0, col="red")
                        abline(med-2*dev, 0, col="red")
                    dev.off()
                }
            }
        }
    }
    return(list(plotName, c(minOfScreens, numOfScreens), c(minOfPlates, 
    numOfPlates)))
  }
}




ZScorePlotTwo<-function(header, dataset, flag, flag2, col4plot, col4anno, 
plotTitle, showPlot){

    dataset<-dataset[which(dataset$SpotType!=-1), ]
    
    if (flag == 0){

        yy<-ylabCompute(dataset, col4plot, flag2)
        yymin<-yy[[1]]
        yymax<-yy[[2]]
        currData<-dataset[[get("col4plot")]]
        
        if (flag2 == 0){
            
            sd2a<-mean(currData, na.rm=T)+2*sd(currData, na.rm=T)
            sd2b<-mean(currData, na.rm=T)-2*sd(currData, na.rm=T)
            sd1a<-mean(currData, na.rm=T)+sd(currData, na.rm=T)
            sd1b<-mean(currData, na.rm=T)-sd(currData, na.rm=T)
        }
        if (flag2 == 1){
            
            sd2a<-median(currData, na.rm=T)+2*mad(currData, na.rm=T)
            sd2b<-median(currData, na.rm=T)-2*mad(currData, na.rm=T)
            sd1a<-median(currData, na.rm=T)+mad(currData, na.rm=T)
            sd1b<-median(currData, na.rm=T)-mad(currData, na.rm=T)
        }
        
    
        if (showPlot == 1){
            if(interactive()){
                plot(currData, xaxt='n', ylim=c(yymin, yymax), ylab=col4plot, 
                main=plotTitle, cex.main=0.8)
            
                if (flag2 == 0){
                    abline(mean(currData, na.rm=T),0,col="blue")
                }

                if (flag2 == 1){
                    abline(median(currData, na.rm=T),0,col="blue")
                }
            
                abline(sd2a,0,col="red")
                abline(sd2b,0,col="red")
                abline(sd1a,0,col="green")
                abline(sd1b,0,col="green")
                identify(dataset[[get("col4plot")]], labels=dataset[[get("col4anno")]])
            }
        }
        
        headerTemp<-strsplit(header[1], ",")
        plotName<-paste(headerTemp[[1]][2], plotTitle, sep="_")
        pdf(paste(plotName, ".pdf", sep=""))
            plot(dataset[[get("col4plot")]], xaxt='n', ylim=c(yymin, yymax), ylab=col4plot, 
            main=plotTitle, cex.main=0.8)
            
            if (flag2 == 0){
                abline(mean(currData, na.rm=T),0,col="blue")
            }
            if (flag2 == 1){
                abline(median(currData, na.rm=T),0,col="blue")
            }
            
            abline(sd2a,0,col="red")
            abline(sd2b,0,col="red")
            abline(sd1a,0,col="green")
            abline(sd1b,0,col="green")
        dev.off()
        
        png(paste(plotName, ".png", sep=""), width=300, height=300)
            plot(dataset[[get("col4plot")]], xaxt='n', ylim=c(yymin, yymax), ylab=col4plot, 
            main=plotTitle, cex.main=0.8)
            
            if (flag2 == 0){
                abline(mean(currData, na.rm=T),0,col="blue")
            }
            if (flag2 == 1){
                abline(median(currData, na.rm=T),0,col="blue")
            }
            
            abline(sd2a,0,col="red")
            abline(sd2b,0,col="red")
            abline(sd1a,0,col="green")
            abline(sd1b,0,col="green")
        dev.off()
        return(plotName)
    }


    if (flag == 1){
    
        minOfScreens<-min(dataset$ScreenNb)
        numOfScreens<-max(dataset$ScreenNb)
    
        headerTemp<-strsplit(header[1], ",")
        plotName<-paste(headerTemp[[1]][2], plotTitle, sep="_")
    
        for (i in minOfScreens:numOfScreens){
            if (length(which(dataset$ScreenNb == i))>0){
            
                subset<-dataset[which(dataset$ScreenNb == i), ]
                currData<-subset[[get("col4plot")]]
                
                if(sum(!is.na(currData))>0){
                
                    yy<-ylabCompute(subset, col4plot, flag2)
                    yymin<-yy[[1]]
                    yymax<-yy[[2]]
                
                    if (flag2 == 0){
                    
                        sd2a<-mean(currData, na.rm=T)+2*sd(currData, na.rm=T)
       	                sd2b<-mean(currData, na.rm=T)-2*sd(currData, na.rm=T)
                        sd1a<-mean(currData, na.rm=T)+sd(currData, na.rm=T)
                        sd1b<-mean(currData, na.rm=T)-sd(currData, na.rm=T)
                    }
                    if (flag2 == 1){
                    
                        sd2a<-median(currData, na.rm=T)+2*mad(currData, na.rm=T)
                        sd2b<-median(currData, na.rm=T)-2*mad(currData, na.rm=T)
                        sd1a<-median(currData, na.rm=T)+mad(currData, na.rm=T)
                        sd1b<-median(currData, na.rm=T)-mad(currData, na.rm=T)
                    }
                
            
                    if (showPlot == 1){
                        if(interactive()){
                            plot(currData, xaxt='n', ylim=c(yymin, yymax), ylab=col4plot, 
                            main=plotTitle, cex.main=0.8)
                    
                            if (flag2 == 0){                        
                                abline(mean(currData, na.rm=T),0,col="blue")
                            }
        
                            if (flag2 == 1){
                                abline(median(currData, na.rm=T),0,col="blue")
                            }
                    
                            abline(sd2a,0,col="red")
                            abline(sd2b,0,col="red")
                            abline(sd1a,0,col="green")
                            abline(sd1b,0,col="green")
                            identify(subset[[get("col4plot")]], labels=subset[[get("col4anno")]])
                        }
                    }
                }else{
                    plot.new()
                    text(0.5, 0.75, paste("Cannot generate plot for exp ", i, sep=""))
		    text(0.5, 0.25, "- only NAs available")
                }
                
                headerTemp<-strsplit(header[1], ", ")
                plotName<-paste(headerTemp[[1]][2], plotTitle, sep="_")
                
                if(sum(!is.na(currData))>0){
                    pdf(paste(plotName, "(Exp. ", i, ").pdf", sep=""))
                        plot(currData, xaxt='n', ylim=c(yymin, yymax), ylab=col4plot, 
                        main=plotTitle, cex.main=0.8)
                    
                        if (flag2 == 0){
                            abline(mean(currData, na.rm=T),0,col="blue")
                        }
                        if (flag2 == 1){
                            abline(median(currData, na.rm=T),0,col="blue")
                        }
                        abline(sd2a,0,col="red")
		        abline(sd2b,0,col="red")
		        abline(sd1a,0,col="green")
                        abline(sd1b,0,col="green")
                    dev.off()
                
                    png(paste(plotName, "(Exp. ", i, ").png", sep=""), width=300, height=300)
                        plot(currData, xaxt='n', ylim=c(yymin, yymax), ylab=col4plot, 
                        main=plotTitle, cex.main=0.8)
                    
                        if (flag2 == 0){
                            abline(mean(currData, na.rm=T),0,col="blue")
                        }
                        if (flag2 == 1){
                            abline(median(currData, na.rm=T),0,col="blue")
                        }
                        abline(sd2a,0,col="red")
		        abline(sd2b,0,col="red")
		        abline(sd1a,0,col="green")
                        abline(sd1b,0,col="green")
                    dev.off()
                }else{
                    pdf(paste(plotName, "(Exp. ", i, ").pdf", sep=""))
                        plot.new()
                        text(0.5, 0.75, paste("Cannot generate plot for exp ", i, sep=""))
		        text(0.5, 0.25, "- only NAs available")
                    dev.off()
                    png(paste(plotName, "(Exp. ", i, ").png", sep=""), width=300, height=300)
                        plot.new()
                        text(0.5, 0.75, paste("Cannot generate plot for exp ", i, sep=""))
		        text(0.5, 0.25, "- only NAs available")
                    dev.off()
                }
            }
        }
        return(list(plotName, minOfScreens, numOfScreens))
    }
    
    if (flag == 2){
    
        numOfScreens<-max(dataset$ScreenNb)
        minOfScreens<-min(dataset$ScreenNb)
    
        for (j in minOfScreens:numOfScreens){
            if (length(which(dataset$ScreenNb == j))>0){
            
                subset<-dataset[which(dataset$ScreenNb == j), ]
                
                minOfPlates<-min(subset$LabtekNb)
                numOfPlates<-max(subset$LabtekNb)
    
                headerTemp<-strsplit(header[1], ",")
                plotName<-paste(headerTemp[[1]][2], plotTitle, sep="_")
    
                for (i in minOfPlates:numOfPlates){
                    if (length(which(subset$LabtekNb == i))>0){
                        
                        subsubset<-subset[which(subset$LabtekNb == i), ]
                        currData<-subsubset[[get("col4plot")]]

                        if(sum(!is.na(currData))>0){

                            yy<-ylabCompute(subsubset, col4plot, flag2)
                            yymin<-yy[[1]]
                            yymax<-yy[[2]]
                        
                            if (flag2 == 0){
                                sd2a<-mean(currData, na.rm=T)+2*sd(currData, na.rm=T)
		                sd2b<-mean(currData, na.rm=T)-2*sd(currData, na.rm=T)
		                sd1a<-mean(currData, na.rm=T)+sd(currData, na.rm=T)
                                sd1b<-mean(currData, na.rm=T)-sd(currData, na.rm=T)
                            }
                            if (flag2 == 1){
                                sd2a<-median(currData, na.rm=T)+2*mad(currData, na.rm=T)
 		                sd2b<-median(currData, na.rm=T)-2*mad(currData, na.rm=T)
 		                sd1a<-median(currData, na.rm=T)+mad(currData, na.rm=T)
 		                sd1b<-median(currData, na.rm=T)-mad(currData, na.rm=T)
 		            }
                    
                            if (showPlot == 1){
                                if(interactive()){
                                    plot(currData, xaxt='n', ylim=c(yymin, yymax), ylab=col4plot, 
                                    main=plotTitle, cex.main=0.8)
                            
                                    if (flag2 == 0){
                                        abline(mean(currData, na.rm=T),0,col="blue")
                                    }
                
                                    if (flag2 == 1){
                                        abline(median(currData, na.rm=T),0,col="blue")
                                    }
                                    abline(sd2a,0,col="red")
	        	            abline(sd2b,0,col="red")
        		            abline(sd1a,0,col="green")
                                    abline(sd1b,0,col="green")
                                    identify(subsubset[[get("col4plot")]], labels=subsubset[[get("col4anno")]])
                                }
                            }
                        }else{
                            plot.new()
                            text(0.5, 0.75, "Cannot generate plot for")
                            text(0.5, 0.25, paste("plate ", i, " in exp ", j, sep=""))
                            text(0.5, 0.05, "- only NAs available")
                        }
                        
                        headerTemp<-strsplit(header[1], ", ")
                        plotName<-paste(headerTemp[[1]][2], plotTitle, sep="_")
                        
                        if(sum(!is.na(currData))>0){
                            pdf(paste(plotName, "Exp", j, "Plate", i, ".pdf", sep="_"))

                                plot(currData, xaxt='n', ylim=c(yymin, yymax), ylab=col4plot, 
                                main=plotTitle, cex.main=0.8)
                            
                                if (flag2 == 0){
                                    abline(mean(currData, na.rm=T),0,col="blue")
                                }
                                if (flag2 == 1){
                                    abline(median(currData, na.rm=T),0,col="blue")
                                }
                                abline(sd2a,0,col="red")
		                abline(sd2b,0,col="red")
    		                abline(sd1a,0,col="green")
                                abline(sd1b,0,col="green")
                            dev.off()
                        
                            png(paste(plotName, "Exp", j, "Plate", i, ".png", sep="_"), 
                            width=300, height=300)
                        
                                plot(currData, xaxt='n', ylim=c(yymin, yymax), ylab=col4plot, 
                                main=plotTitle, cex.main=0.8)
                            
                                if (flag2 == 0){
                                    abline(mean(currData, na.rm=T),0,col="blue")
                                }
                                if (flag2 == 1){
                                    abline(median(currData, na.rm=T),0,col="blue")
                                }
                                abline(sd2a,0,col="red")
		                abline(sd2b,0,col="red")
    		                abline(sd1a,0,col="green")
                                abline(sd1b,0,col="green")
                            dev.off()
                        }else{
                            pdf(paste(plotName, "Exp", j, "Plate", i, ".pdf", sep="_"))
                                plot.new()
                                text(0.5, 0.75, "Cannot generate plot for")
                                text(0.5, 0.25, paste("plate ", i, " in exp ", j, sep=""))
                                text(0.5, 0.05, "- only NAs available")
                            dev.off()
                            png(paste(plotName, "Exp", j, "Plate", i, ".png", sep="_"), 
                            width=300, height=300)
                                plot.new()
                                text(0.5, 0.75, "Cannot generate plot for")
                                text(0.5, 0.25, paste("plate ", i, " in exp ", j, sep=""))
                                text(0.5, 0.05, "- only NAs available")
                            dev.off()
                        }
                    }
                }
            }
        }
        return(list(plotName, c(minOfScreens, numOfScreens), 
        c(minOfPlates, numOfPlates)))
    }
}


ylabCompute<-function(dataset, col4plot, flag2){

    currData<-dataset[[get("col4plot")]]
    
    if(sum(!is.na(currData))>0){

        if (flag2 == 0){
            
	    dev1<-mean(currData, na.rm=T)-2*sd(currData, na.rm=T)
	    dev2<-mean(currData, na.rm=T)+2*sd(currData, na.rm=T)
	
            if (dev1<min(currData, na.rm=T)){
                yymin<-mean(currData, na.rm=T)-2*sd(currData, na.rm=T)
            }else{
                yymin<-min(currData, na.rm=T)
            }
        
            if (dev2>max(currData, na.rm=T)){
                yymax<-mean(currData, na.rm=T)+2*sd(currData, na.rm=T)
            }else{
                yymax<-max(currData, na.rm=T)
            }
        }

        if (flag2 == 1){
    
            dev1<-median(currData, na.rm=T)-2*mad(currData, na.rm=T)
            dev2<-median(currData, na.rm=T)+2*mad(currData, na.rm=T)

            if (dev1<min(currData, na.rm=T)){
                yymin<-median(currData, na.rm=T)-2*mad(currData, na.rm=T)
            }else{
                yymin<-min(currData, na.rm=T)
            }
            if (dev2>max(currData, na.rm=T)){
                yymax<-median(currData, na.rm=T)+2*mad(currData, na.rm=T)
            }else{
                yymax<-max(currData, na.rm=T)
            }
        }
    }else{
        yymin<-NA
        yymax<-NA
    }

    invisible(list(yymin, yymax))
}




ZPrime<-function(posControls, negControls){

    ZPrimeFactor<-1-3*((mad(posControls, na.rm=TRUE)+mad(negControls, na.rm=TRUE))/
    (abs(median(posControls, na.rm=TRUE)-median(negControls, na.rm=TRUE))))

}

dynamicRange<-function(dataset, channel){
    
    dataset<-dataset[which(dataset$SpotType!=-1), ]

    numOfScreens<-max(dataset$ScreenNb)
    minOfScreens<-min(dataset$ScreenNb)

    totNbPlates4allScreens<-0
    for (j in minOfScreens:numOfScreens){
        if (length(which(dataset$ScreenNb == j))>0){
            subset<-createSubset(dataset, dataset$ScreenNb, j)
            totNbPlates4allScreens<-totNbPlates4allScreens+length(unique(subset$LabtekNb))
        }
    }

    DR<-rep(0, totNbPlates4allScreens)

    count<-0
    for (j in minOfScreens:numOfScreens){
        if (length(which(dataset$ScreenNb == j))>0){
            partSet<-dataset[which(dataset$ScreenNb == j), ]
            minOfPlates<-min(partSet$LabtekNb)
            numOfPlates<-max(partSet$LabtekNb)
            
            for (i in minOfPlates:numOfPlates){
                if (length(which(partSet$LabtekNb == i))>0){
                    count<-count+1    
                    subset<-partSet[which(partSet$LabtekNb == i), ]
                    DR[count]<-mean(subset[[get("channel")]][which(subset$SpotType == 0)])/
                    mean(subset[[get("channel")]][which(subset$SpotType == 1)])
                }
            }
        }
    }
    invisible(DR)
}


numCellQualControl<-function(DataSetFile, nbLinesHeader, plotTitle){

    header<-readLines(DataSetFile, nbLinesHeader)
    data<-read.table(DataSetFile, skip=nbLinesHeader, colClasses=c(NA, NA, NA, NA, 
    "factor", NA, NA, NA, NA, NA, NA, NA, NA, NA), stringsAsFactors=FALSE)
    
    threshNumCells<-mean(data$NbCells, na.rm=T)-3*sd(data$NbCells, na.rm=T)
    upperThreshNumCells<-mean(data$NbCells, na.rm=T)+3*sd(data$NbCells, na.rm=T)

    linesUnderThreshNumCells<-data[which(data$NbCells<threshNumCells), ]
    lineNbsUnderThreshNumCells<-which(data$NbCells<threshNumCells)

    linesOverUpperThreshNumCells<-data[which(data$NbCells>upperThreshNumCells), ]
    lineNbsOverUpperThreshNumCells<-which(data$NbCells>upperThreshNumCells)

    
    s1<-"Number of cells per well; Data for wells under lower threshold:"
    print(paste(s1, threshNumCells, sep=" "), quote=F)
    print("", quote=F)
    
    xlabHist<-"Number of cells per well"
    plotTitleHist<-"Distribution of the number of cells"
    
    if(interactive()){
        showThreshHisto<-hist(data$NbCells, breaks=20, xlab=xlabHist, 
        main=plotTitleHist, cex.main=0.8)
        lines(c(threshNumCells, threshNumCells), c(0, max(showThreshHisto$counts)), 
        col="red")
        lines(c(upperThreshNumCells, upperThreshNumCells), 
        c(0, max(showThreshHisto$counts)), col="red")
    }

    if (nrow(linesUnderThreshNumCells) == 0){
        print("No wells under lower threshold.")
    }else{
        print(cbind(lineNbsUnderThreshNumCells, linesUnderThreshNumCells))
    }
    
    if(interactive()){
        s3<-"\r\nTo set new lower threshold enter value, otherwise \"n\".\r\n\r\n"
    	userInput1<-readline(s3)
    }else{
    	userInput1<-"n"
    }
    
    if (!is.na(as.integer(userInput1))){
    ##set new user-defined threshold:
        linesUnderThreshNumCells<-data[which(data$NbCells<as.double(userInput1)), ]
        lineNbsUnderThreshNumCells<-which(data$NbCells<as.double(userInput1))
        threshNumCells<-as.double(userInput1)

        s4<-"Number of cells per well; Data for wells under new lower threshold:"
        print(paste(s4, threshNumCells, sep=" "), quote=F)
        if (nrow(linesUnderThreshNumCells) == 0){
            print("No wells under lower threshold.")
        }else{
            print(cbind(lineNbsUnderThreshNumCells, linesUnderThreshNumCells))
        }

    }else{
        ##generates warning "NAs introduced by coercion"
        if(is.na(match(userInput1, "n")) | match(userInput1, "n")!=1){
            print("Error. Wrong input.")
        }
    }

    if(interactive()){
        s5<-"\r\nDiscard all wells under lower threshold? [y-n]\r\n\r\n"
        userInput2<-readline(s5)
    }else{
        userInput2<-"y"
    }
    if (!is.na(match(userInput2, "y")) & match(userInput2, "y") == 1){
        
        data$SpotType[lineNbsUnderThreshNumCells]<-(-1)

        ##save list of discarded siRNA values:
        headerTemp<-strsplit(header[1], ",")
        s6<-"numCellQualControl_discarded_lower.txt"
        filnam<-paste(headerTemp[[1]][2], s6, sep="_")
        makeHeader<-""
        for (h in 1:length(colnames(data))){
            if (h == 1){
                makeHeader<-colnames(data)[h]
            }else{
                makeHeader<-paste(makeHeader, colnames(data)[h], sep="\t")
            }
        }
        write.table(makeHeader, file=filnam, quote=F, col.names=F, row.names=F)
        write.table(cbind(lineNbsUnderThreshNumCells, linesUnderThreshNumCells), 
        file=filnam, quote=F, sep="\t", col.names=F, row.names=F, append=T)

    }else{
        if(is.na(match(userInput2, "n")) | match(userInput2, "n")!=1){
            print("Error. Wrong input.")
        }
    }

    s7<-"Number of cells per well; Data for wells over upper threshold:"
    print(paste(s7, upperThreshNumCells, sep=" "), quote=F)
    print("", quote=F)
    
    if (nrow(linesOverUpperThreshNumCells) == 0){
        print("No wells over upper threshold.")
    }else{
        print(cbind(lineNbsOverUpperThreshNumCells, linesOverUpperThreshNumCells))
    }

    if(interactive()){
        s8<-"\r\nTo set new upper threshold enter value, otherwise \"n\".\r\n\r\n"
        userInput3<-readline(s8)
    }else{
        userInput3<-"n"
    }
    
    if (!is.na(as.integer(userInput3))){
        linesOverUpperThreshNumCells<-data[which(data$NbCells>as.double(userInput3)), ]
        lineNbsOverUpperThreshNumCells<-which(data$NbCells>as.double(userInput3))
        upperThreshNumCells<-as.double(userInput3)


        s9<-"Number of cells per well; Data for wells over new upper threshold:"
        print(paste(s9, upperThreshNumCells, sep=" "), quote=F)
        if (nrow(linesOverUpperThreshNumCells) == 0){
            print("No wells over upper threshold.")
        }else{
            print(cbind(lineNbsOverUpperThreshNumCells, linesOverUpperThreshNumCells))
        }

    }else{
        ##here warning "NAs introduced by coercion"?
        if(is.na(match(userInput3, "n")) | match(userInput3, "n")!=1){
            print("Error. Wrong input.")
        }
    }

    if(interactive()){
        s1<-"\r\nDiscard all wells over upper threshold? [y-n]\r\n\r\n"
        userInput4<-readline(s1)
    }else{
        userInput4<-"n"
    }
    if (!is.na(match(userInput4, "y")) & match(userInput4, "y") == 1){

        data$SpotType[lineNbsOverUpperThreshNumCells]<-(-1)

        headerTemp<-strsplit(header[1], ",")
        histoName<-paste(headerTemp[[1]][2], plotTitle, ".pdf", sep="_")
        
        pdf(histoName)
            if (!is.na(match(userInput2, "y")) & match(userInput2, "y") == 1){
                showThreshHisto<-hist(data$NbCells, breaks=20, xlab="Number of cells per well", 
                main=plotTitle, cex.main=0.8)
                
                lines(c(upperThreshNumCells, upperThreshNumCells), 
                c(0, max(showThreshHisto$counts)), col="red")
                lines(c(threshNumCells, threshNumCells), c(0, max(showThreshHisto$counts)), 
                col="red")
            }
            if (!is.na(match(userInput2, "n")) & match(userInput2, "n") == 1){
                showThreshHisto<-hist(data$NbCells, breaks=20, xlab="Number of cells per well", 
                main=plotTitle, cex.main=0.8)
                lines(c(upperThreshNumCells, upperThreshNumCells), 
                c(0, max(showThreshHisto$counts)), col="red")
            }
        dev.off()
        ##save list of discarded siRNA values:
        filnam<-paste(headerTemp[[1]][2], "numCellQualControl_discarded_higher.txt", sep="_")
        makeHeader<-""
        for (h in 1:length(colnames(data))){
            if (h == 1){
                makeHeader<-colnames(data)[h]
            }else{
                makeHeader<-paste(makeHeader, colnames(data)[h], sep="\t")
            }
        }
        write.table(makeHeader, file=filnam, quote=F, col.names=F, row.names=F)
        write.table(cbind(lineNbsOverUpperThreshNumCells, linesOverUpperThreshNumCells), 
        file=filnam, quote=F, sep="\t", col.names=F, row.names=F, append=T)

    }else{

        headerTemp<-strsplit(header[1], ",")
        histoName<-paste(headerTemp[[1]][2], plotTitle, ".pdf", sep="_")
        pdf(histoName)
            if (!is.na(match(userInput2, "y")) & match(userInput2, "y") == 1){
                showThreshHisto<-hist(data$NbCells, breaks=20, xlab="Number of cells per well", 
                main=plotTitle, cex.main=0.8)
                lines(c(threshNumCells, threshNumCells), c(0, max(showThreshHisto$counts)), 
                col="red")
            }
            if (!is.na(match(userInput2, "n")) & match(userInput2, "n") == 1){
                showThreshHisto<-hist(data$NbCells, breaks=20, xlab="Number of cells per well", 
                main=plotTitle, cex.main=0.8)
            }
        dev.off()

        if(is.na(match(userInput4, "n")) | match(userInput4, "n")!=1){
            print("Error. Wrong input.")
        }
    }

    write.table(header, file=DataSetFile, quote=F, col.names=F, row.names=F)
    write.table(data, file=DataSetFile, sep="\t", quote=F, col.names=T, append=T)
    invisible(histoName)
}


percCellQualControl<-function(DataSetFile, nbLinesHeader, plotTitle){

    header<-readLines(DataSetFile,nbLinesHeader)
    data<-read.table(DataSetFile, skip=nbLinesHeader, colClasses=c(NA, NA, NA, NA, 
    "factor", NA, NA, NA, NA, NA, NA, NA, NA, NA), stringsAsFactors=FALSE)
	
    threshPcCells<-mean(data$PercCells, na.rm=T)-3*sd(data$PercCells, na.rm=T)
    upperThreshPcCells<-mean(data$PercCells, na.rm=T)+3*sd(data$PercCells, na.rm=T)
    
    linesUnderthreshPcCells<-data[which(data$PercCells<threshPcCells), ]
    lineNbsUnderthreshPcCells<-which(data$PercCells<threshPcCells)
    
    linesOverUpperthreshPcCells<-data[which(data$PercCells>upperThreshPcCells), ]
    lineNbsOverUpperthreshPcCells<-which(data$PercCells>upperThreshPcCells)

    print("", quote=F)
    s1<-"Percentage of cells per well; Data for wells under lower threshold:"
    print(paste(s1, threshPcCells, sep=" "), quote=F)
    print("", quote=F)
    
    xlabs<-"Percentage of cells per well"
    tit<-"Distribution of the percentage of cells"
    if(interactive()){
        showThreshHisto<-hist(data$PercCells, breaks=20, xlab=xlabs, main=tit, 
        cex.main=0.8)
        lines(c(threshPcCells, threshPcCells), c(0, max(showThreshHisto$counts)), 
        col="red")
    }

    if (nrow(linesUnderthreshPcCells) == 0){
        print("No wells under lower threshold.")
    }else{
        print(cbind(lineNbsUnderthreshPcCells, linesUnderthreshPcCells))
    }

    if(interactive()){
        s2<-"\r\n\r\nTo set new lower threshold enter value, otherwise \"n\"\r\n\r\n"
        userInput1<-readline(s2)
    }else{
        userInput1<-"n"
    }
    
    if (!is.na(as.integer(userInput1))){
        linesUnderthreshPcCells<-data[which(data$PercCells<as.double(userInput1)), ]
        lineNbsUnderthreshPcCells<-which(data$PercCells<as.double(userInput1))
        threshPcCells<-as.double(userInput1)

        s3<-"Percentage of cells per well; Data for wells under new lower threshold:"
        print(paste(s3, threshPcCells, sep=" "), quote=F)
        if (nrow(linesUnderthreshPcCells) == 0){
            print("No wells under lower threshold.")
        }else{
            print(cbind(lineNbsUnderthreshPcCells, linesUnderthreshPcCells))
        }


    }else{
        if(is.na(match(userInput1, "n")) | match(userInput1, "n")!=1){
            print("Error. Wrong input.")
        }
    }

    if(interactive()){
        s4<-"\r\nDiscard all wells under lower threshold? [y-n]\r\n\r\n"
        userInput2<-readline(s4)
    }else{
        userInput2<-"y"
    }
    if (!is.na(match(userInput2, "y")) & match(userInput2, "y") == 1){

        data$SpotType[lineNbsUnderthreshPcCells]<-(-1)

        ##save list of discarded siRNA values:
        headerTemp<-strsplit(header[1], ",")
        s5<-"percCellQualControl_discarded_lower.txt"
        filnam<-paste(headerTemp[[1]][2], s5, sep="_")
        makeHeader<-""
        for (h in 1:length(colnames(data))){
            if (h == 1){
                makeHeader<-colnames(data)[h]
            }else{
                makeHeader<-paste(makeHeader, colnames(data)[h], sep="\t")
            }
        }
        write.table(makeHeader, file=filnam, quote=F, col.names=F, row.names=F)
        write.table(cbind(lineNbsUnderthreshPcCells, linesUnderthreshPcCells), 
        file=filnam, quote=F, sep="\t", col.names=F, row.names=F, append=T)

    }else{
        if(is.na(match(userInput2, "n")) | match(userInput2, "n")!=1){
            print("Error. Wrong input.")
        }
    }

    s6<-"Percentage of cells per well; Data for wells over upper threshold:"
    print(paste(s6, upperThreshPcCells, sep=" "), quote=F)
    print("", quote=F)
    
    if (nrow(linesOverUpperthreshPcCells) == 0){
        print("No wells over upper threshold.")
    }else{
        print(cbind(lineNbsOverUpperthreshPcCells, linesOverUpperthreshPcCells))
    }

    if(interactive()){
        s7<-"\r\nTo set new upper threshold enter value, otherwise \"n\".\r\n\r\n"
        userInput3<-readline(s7)
    }else{
        userInput3<-"n"
    }
    
    if (!is.na(as.integer(userInput3))){
        linesOverUpperthreshPcCells<-data[which(data$PercCells>as.double(userInput3)), ]
        lineNbsOverUpperthreshPcCells<-which(data$PercCells>as.double(userInput3))
        upperThreshPcCells<-as.double(userInput3)

        s8<-"Percentage of cells per well; Data for wells over new upper threshold:"
        print(paste(s8, upperThreshPcCells, sep=" "), quote=F)
        if (nrow(linesOverUpperthreshPcCells) == 0){
            print("No wells over upper threshold.")
        }else{
            print(cbind(lineNbsOverUpperthreshPcCells, linesOverUpperthreshPcCells))
        }

    }else{
        ##generates warning "NAs introduced by coercion"
        if(is.na(match(userInput3, "n")) | match(userInput3, "n")!=1){
            print("Error. Wrong input.")
        }
    }

    if(interactive()){
        s9<-"\r\nDiscard all wells over upper threshold? [y-n]\r\n\r\n"
        userInput4<-readline(s9)
    }else{
        userInput4<-"n"
    }
    if (!is.na(match(userInput4, "y")) & match(userInput4, "y") == 1){
        data$SpotType[lineNbsOverUpperthreshPcCells]<-(-1)

        headerTemp<-strsplit(header[1], ",")
        histoName<-paste(headerTemp[[1]][2], plotTitle, ".pdf", sep="_")
        pdf(histoName)
            if (!is.na(match(userInput2, "y")) & match(userInput2, "y") == 1){
            
                showThreshHisto<-hist(data$PercCells, breaks=20, 
                xlab="Percentage of cells per well", main=plotTitle, cex.main=0.8)
                
                lines(c(upperThreshPcCells, upperThreshPcCells), 
                c(0, max(showThreshHisto$counts)), col="red")
                lines(c(threshPcCells, threshPcCells), c(0, max(showThreshHisto$counts)), 
                col="red")
            }
            if (!is.na(match(userInput2, "n")) & match(userInput2, "n") == 1){
                showThreshHisto<-hist(data$PercCells, breaks=20, 
                xlab="Percentage of cells per well", main=plotTitle, cex.main=0.8)
                
                lines(c(upperThreshPcCells, upperThreshPcCells), 
                c(0, max(showThreshHisto$counts)), col="red")
            }
        dev.off()
        ##save list of discarded siRNA values:
        s1<-"percCellQualControl_discarded_higher.txt"
        filnam<-paste(headerTemp[[1]][2], s1, sep="_")
        makeHeader<-""
        for (h in 1:length(colnames(data))){
            if (h == 1){
                makeHeader<-colnames(data)[h]
            }else{
                makeHeader<-paste(makeHeader, colnames(data)[h], sep="\t")
            }
        }
        write.table(makeHeader, file=filnam, quote=F, col.names=F, row.names=F)
        write.table(cbind(lineNbsOverUpperthreshPcCells, linesOverUpperthreshPcCells), 
        file=filnam, quote=F, sep="\t", col.names=F, row.names=F, append=T)

    }else{

        headerTemp<-strsplit(header[1], ",")
        histoName<-paste(headerTemp[[1]][2], plotTitle, ".pdf", sep="_")
        pdf(histoName)
            if (!is.na(match(userInput2, "y")) & match(userInput2, "y") == 1){
                showThreshHisto<-hist(data$PercCells, breaks=20, 
                xlab="Percentage of cells per well", main=plotTitle, cex.main=0.8)
                
                lines(c(threshPcCells, threshPcCells), c(0, max(showThreshHisto$counts)), 
                col="red")
            }
            if (!is.na(match(userInput2, "n")) & match(userInput2, "n") == 1){
                showThreshHisto<-hist(data$PercCells, breaks=20, 
                xlab="Percentage of cells per well", main=plotTitle, cex.main=0.8)
            }
        dev.off()

        if(is.na(match(userInput4, "n")) | match(userInput4, "n")!=1){
            print("Error. Wrong input.")
        }
    }

    write.table(header, file=DataSetFile, quote=F, col.names=F, row.names=F)
    write.table(data, file=DataSetFile, sep="\t", quote=F, col.names=T, append=T)
    invisible(histoName)
}


ZPRIMEQualControl<-function(header, data, channel, plotTitle, showPlot){

    data<-data[which(data$SpotType!=-1), ]

    numOfScreens<-max(data$ScreenNb)
    minOfScreens<-min(data$ScreenNb)

    totalNbOfPlates4allScreens<-0
    for (j in minOfScreens:numOfScreens){
        if (length(which(data$ScreenNb == j))>0){
            subset<-createSubset(data, data$ScreenNb, j)
            lenData<-length(unique(subset$LabtekNb))
            totalNbOfPlates4allScreens<-totalNbOfPlates4allScreens+lenData
        }
    }

    ZPRIME<-rep(0, totalNbOfPlates4allScreens)

    count<-0
    for (j in minOfScreens:numOfScreens){
        if (length(which(data$ScreenNb == j))>0){
        
            smallSet<-data[which(data$ScreenNb == j), ]
            
            minOfPlates<-min(smallSet$LabtekNb)
            numOfPlates<-max(smallSet$LabtekNb)
        
            for (i in minOfPlates:numOfPlates){
                if (length(which(smallSet$LabtekNb == i))>0){
                    count<-count+1
                
                    subset<-smallSet[which(smallSet$LabtekNb == i), ]
                    posContr<-subset[[get("channel")]][which(subset$SpotType == 1)]
                    negContr<-subset[[get("channel")]][which(subset$SpotType == 0)]
                    ZPRIME[count]<-ZPrime(posContr, negContr)
                }
            }
        }
    }

    headerTemp<-strsplit(header[1], ",")
    filnam<-paste(headerTemp[[1]][2], "Z'Scores.txt", sep="_")

    if (showPlot == 1){
        print("", quote=F)
    }
    vekForPlot<-rep(0, length(ZPRIME))
    count<-0
    Exp<-rep(0, length(ZPRIME))
    Pla<-rep(0, length(ZPRIME))
    for (j in minOfScreens:numOfScreens){
        if (length(which(data$ScreenNb == j))>0){
    
            smallSet<-data[which(data$ScreenNb == j), ]
                
            minOfPlates<-min(smallSet$LabtekNb)
            numOfPlates<-max(smallSet$LabtekNb)
    
            for (i in minOfPlates:numOfPlates){
                if (length(which(smallSet$LabtekNb == i))>0){
            
                    count<-count+1
                    s2<-"The Z\'Score for Experiment"
                    ausgabe1<-paste(s2, j, "Plate", i, "is", ZPRIME[count], sep=" ")
                    Exp[count]<-j
                    Pla[count]<-i
                    if (showPlot == 1){
                        print(ausgabe1, quote=F)
                    }
                    vekForPlot[count]<-paste(j, i, sep="_")
                }
            }
        }
    }
    if (showPlot == 1){
        print("", quote=F)
    }
    
    ausgabe2<-data.frame(Exp, Pla, ZPRIME)
    colnames(ausgabe2)<-c("Experiment", "Plate", "Z_Prime_Score")
    write.table(ausgabe2, file=filnam, quote=F, row.names=F)
    
    if (showPlot == 1){
        if(interactive()){
            plot(ZPRIME, xaxt='n', main=plotTitle, cex.main=0.8)
            axis(1, 1:length(vekForPlot), labels=vekForPlot, xlab="Exp_Plate", cex.axis=0.6, 
            las=2)
        }
    }

    headerTemp<-strsplit(header[1], ",")
    plotName<-paste(headerTemp[[1]][2], plotTitle, sep="_")
    pdf(paste(plotName, ".pdf", sep=""))
        plot(ZPRIME, xaxt='n', main=plotTitle, cex.main=0.8)
        axis(1, 1:length(vekForPlot), labels=vekForPlot, xlab="Exp_Plate", cex.axis=0.6, 
        las=2)
    dev.off()
    png(paste(plotName, ".png", sep=""), width=300, height=300)
        plot(ZPRIME, xaxt='n', main=plotTitle, cex.main=0.8, cex.lab=0.8, cex.axis=0.8)
        axis(1, 1:length(vekForPlot), labels=vekForPlot, xlab="Exp_Plate", cex.axis=0.6, 
        las=2)
    dev.off()

    ZPrimeTabelle<-read.table(filnam)
    invisible(list(plotName, ZPrimeTabelle))
}


SNRQualControl<-function(dataSetFile, nbLinesHeader, channel, noise, plotTitle, showPlot){

    header<-readLines(dataSetFile, nbLinesHeader)
    data<-read.table(dataSetFile, skip=nbLinesHeader, colClasses=c(NA, NA, NA, NA, 
    "factor", NA, NA, NA, NA, NA, NA, NA, NA, NA), stringsAsFactors=FALSE)

    data<-data[which(data$SpotType!=-1), ]

    SNR<-data[[get("channel")]]/data[[get("noise")]]

    if (showPlot==1){
        if(interactive()){
            hist(SNR, breaks=20, xlab="SNR per well")
        }
    }

    numOfScreens<-max(data$ScreenNb)
    minOfScreens<-min(data$ScreenNb)

    c1<-ceiling(sqrt(length(unique(data$ScreenNb))))
    c2<-ceiling(length(unique(data$ScreenNb))/c1)
    
    if (showPlot==1){
        if(interactive()){
            x11()
            par(mfrow=c(c1, c2))
        }
    }
    for (i in minOfScreens:numOfScreens){
        if (length(which(data$ScreenNb == i))>0){
            subsetSnr<-SNR[which(data$ScreenNb == i)]
            
            if (showPlot==1){
                if(interactive()){
                    hist(subsetSnr, breaks=20, xlab=paste("SNR per well", sep=" "), 
                    main=paste(plotTitle, "for Exp.", i, sep=" "))
                }
            }
        }
    }

    for (j in minOfScreens:numOfScreens){
        if (length(which(data$ScreenNb == j))>0){
        
            subset<-data[which(data$ScreenNb == j), ]
            
            minOfPlates<-min(subset$LabtekNb)
            numOfPlates<-max(subset$LabtekNb)
            
            if (showPlot==1){
                if(interactive()){
                    x11()
                    c3<-ceiling(sqrt(length(unique(subset$LabtekNb))))
                    c4<-ceiling(length(unique(subset$LabtekNb))/c3)
                    par(mfrow=c(c3, c4))
                }
            }
            
            for (i in minOfPlates:numOfPlates){
                if (length(which(subset$LabtekNb == i))>0){
                
                    subsetSnr<-SNR[which(subset$LabtekNb == i)]
                    
                    if (showPlot==1){
                        if(interactive()){
                            hist(subsetSnr, breaks=20, xlab=paste("SNR per well", sep=" "), 
                            main=paste(plotTitle, "for Exp.", j, "Plate", i, sep=" "))
                        }
                    }
                }
            }
        }
    }

    headerTemp<-strsplit(header[1], ",")
    histoName<-paste(headerTemp[[1]][2], plotTitle, ".pdf", sep="_")
    pdf(histoName)
        hist(SNR, breaks=20, xlab="SNR per well")
    dev.off()

    histoName<-paste(headerTemp[[1]][2], plotTitle, "PerExp.pdf", sep="_")
    pdf(histoName)
        par(mfrow=c(c1, c2))
        for (i in minOfScreens:numOfScreens){
            if (length(which(data$ScreenNb == i))>0){
            
                subsetSnr<-SNR[which(data$ScreenNb == i)]
                
                hist(subsetSnr, breaks=20, xlab=paste("SNR per well", sep=" "), 
                main=paste(plotTitle, "for Exp.", i, sep=" "), cex.main=0.8)
            }
        }
    dev.off()

    for (j in minOfScreens:numOfScreens){
        if (length(which(data$ScreenNb == j))>0){
            subset<-data[which(data$ScreenNb == j), ]
            
            minOfPlates<-min(subset$LabtekNb)
            numOfPlates<-max(subset$LabtekNb)
            
            histoName<-paste(headerTemp[[1]][2], plotTitle, "PerLabtek_Exp", j, ".pdf", sep="_")
            
            pdf(histoName)
                c3<-ceiling(sqrt(length(unique(subset$LabtekNb))))
                c4<-ceiling(length(unique(subset$LabtekNb))/c3)
                par(mfrow=c(c3, c4))
                for (i in minOfPlates:numOfPlates){
                    if (length(which(subset$LabtekNb == i))>0){
                
                        subsetSnr<-SNR[which(subset$LabtekNb == i)]
                        
                        hist(subsetSnr, breaks=20, xlab=paste("SNR per well", sep=" "), 
                        main=paste(plotTitle, "for Exp.", j, "Plate", i, sep=" "), cex.main=0.8)
                    }
                }
            dev.off()
        }
    }
}

DRQualControl<-function(header, data, nbLinesHeader, channel, plotTitle, showPlot){

    data<-data[which(data$SpotType!=-1), ]

    numOfScreens<-max(data$ScreenNb)
    minOfScreens<-min(data$ScreenNb)

    DR<-dynamicRange(data, channel)

    headerTemp<-strsplit(header[1], ",")
    filnam<-paste(headerTemp[[1]][2], "DR.txt", sep="_")
    write.table(paste("Experiment", "Plate", "DR", sep="\t"), file=filnam, quote=F, 
    col.names=F, row.names=F)

    print("", quote=F)
    vekForPlot<-rep(0, length(DR))
    count<-0
    for (j in minOfScreens:numOfScreens){
        if (length(which(data$ScreenNb == j))>0){

            subset<-createSubset(data, data$ScreenNb, j)
            minOfPlates<-min(subset$LabtekNb)
            numOfPlates<-max(subset$LabtekNb)
        
            for (i in minOfPlates:numOfPlates){
                if (length(which(subset$LabtekNb == i))>0){
            
                    count<-count+1
                    s1<-"The dynamic range for Experiment"
                    ausgabe1<-paste(s1, j, "Plate", i, "is", DR[count], sep=" ")
                    ausgabe2<-paste(j, i, DR[count], sep="\t")
                    print(ausgabe1, quote=F)
                    vekForPlot[count]<-paste(j, i, sep="_")
                    write.table(ausgabe2, file=filnam, quote=F, col.names=F, row.names=F, append=T)
                }
            }
        }
    }
    print("", quote=F)
    
    if (showPlot==1){
        if(interactive()){
            plot(DR, xaxt='n', main=plotTitle, cex.main=0.8)
            axis(1, 1:length(vekForPlot), labels=vekForPlot, xlab="Exp_Plate", 
            cex.axis=0.6, las=2)
        }
    }
    
    headerTemp<-strsplit(header[1], ",")
    plotName<-paste(headerTemp[[1]][2], plotTitle, ".pdf", sep="_")
    pdf(plotName)
        plot(DR, xaxt='n', main=plotTitle, cex.main=0.8)
        axis(1, 1:length(vekForPlot), labels=vekForPlot, xlab="Exp_Plate", 
        cex.axis=0.6, las=2)
    dev.off()
}




plotControlHisto<-function(header, dataset, channel, plotTitle, showPlot){

    All<-dataset[[get("channel")]][which(dataset$SpotType!=-1)]

    Neg<-dataset[[get("channel")]][which(dataset$SpotType == 0)]
    Pos<-dataset[[get("channel")]][which(dataset$SpotType == 1)]
    Other<-dataset[[get("channel")]][which(dataset$SpotType == 2)]

    if (sum(!is.na(All)) == 0){
        stop("Cannot plot histogram (only NAs in dataset)")
    }
#######
#    if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#        a<-min(All, na.rm=T)
#        b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#        d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#        computeHistoBreaks=seq(a, b, d)
#    }else{
#        a<-min(All, na.rm=T)
#        b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
#        d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#        computeHistoBreaks=seq(a, b, d)
#    }
######### Changes Will Greenwood (22.07.2009):
    if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
        a<-min(All, na.rm=T)
        b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
        d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
        computeHistoBreaks=seq(a, b, length=d)
    }else{
        a<-min(All, na.rm=T)
        b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
        d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
        computeHistoBreaks=seq(a, b, length=d)
    }
######### End of changes Will Greenwood (22.07.2009)

    all<-hist(Other, plot=F, breaks=computeHistoBreaks)
    pos<-hist(Pos, plot=F, breaks=computeHistoBreaks)
    neg<-hist(Neg, plot=F, breaks=computeHistoBreaks)

    if (showPlot == 1){
        if(interactive()){
            plot(all, main=plotTitle, xlab=channel, cex.main=0.7)
            par(fg="green")
            lines(pos, angle=45, density=20)
            par(fg="red")
            lines(neg, angle=-45, density=20)
            par(fg="black")
            legend("topleft", c("Data", "Positive Controls", "Negative Controls"), 
            angle=c(0, 45, -45), density=c(0, 20, 20), bty="n", cex=0.8)
        }
    }

    headerTemp<-strsplit(header[1], ",")
    histoName<-paste(headerTemp[[1]][2], plotTitle, sep="_")
    pdf(paste(histoName, ".pdf", sep=""))
        plot(all, main=plotTitle, xlab=channel, cex.main=0.7)
        par(fg="green")
        lines(pos, angle=45, density=20)
        par(fg="red")
        lines(neg, angle=-45, density=20)
        par(fg="black")
        legend("topleft", c("Data", "Positive Controls", "Negative Controls"), 
        angle=c(0, 45, -45), density=c(0, 20, 20), bty="n", cex=0.8)
    dev.off()
    png(paste(histoName, ".png", sep=""), width=300, height=300)
        plot(all, main=plotTitle, xlab=channel, cex.main=0.8)
        par(fg="green")
        lines(pos, angle=45, density=20)
        par(fg="red")
        lines(neg, angle=-45, density=20)
        par(fg="black")
        legend("topleft", c("Data", "Positive Controls", "Negative Controls"), 
        angle=c(0, 45, -45), density=c(0, 20, 20), bty="n", cex=0.7)
    dev.off()
    invisible(histoName)
}


plotControlHistoPerplate<-function(header, dataset, channel, plotTitle, 
plotDesign, showPlot){

    dataset<-dataset[which(dataset$SpotType!=-1), ]

    numOfScreens<-max(dataset$ScreenNb)
    minOfScreens<-min(dataset$ScreenNb)

    if (showPlot == 1){
        if(interactive()){
            for (j in minOfScreens:numOfScreens){
                if (length(which(dataset$ScreenNb == j))>0){
        
                    subset<-dataset[which(dataset$ScreenNb == j), ]
            
                    minOfPlates<-min(subset$LabtekNb)
                    numOfPlates<-max(subset$LabtekNb)
            
                    x11()
                    if (plotDesign == 1){
                        c1<-ceiling(sqrt(length(unique(subset$LabtekNb))))
                        c2<-ceiling(length(unique(subset$LabtekNb))/c1)
                        par(mfrow=c(c1, c2), oma=c(0, 0, 2, 0))
                    }

                    for (i in minOfPlates:numOfPlates){
                        if (length(which(subset$LabtekNb == i))>0){
                
                            subsubset<-subset[which(subset$LabtekNb == i), ]
                            All<-subsubset[[get("channel")]][which(subsubset$SpotType!=-1)]
                            Neg<-subsubset[[get("channel")]][which(subsubset$SpotType == 0)]
                            Pos<-subsubset[[get("channel")]][which(subsubset$SpotType == 1)]
                            Other<-subsubset[[get("channel")]][which(subsubset$SpotType == 2)]

                            if (sum(!is.na(All)) == 0){
                                plot.new()
                                text(0.5, 0.75, "Cannot plot histogram for")
				text(0.5, 0.25, paste("exp ", j, " plate ", i, sep=""))
                                text(0.5, 0.05, "- only NAs available")
                            }else{

#######
#                                if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                                    a<-min(All, na.rm=T)
#                                    b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                                    d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                                    computeHistoBreaks=seq(a, b, d)
#                                }else{
#                                    a<-min(All, na.rm=T)
#                                    b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
#                                    d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                                    computeHistoBreaks=seq(a, b, d)
#                                }
######### Changes Will Greenwood (22.07.2009):
                                if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                                    a<-min(All, na.rm=T)
                                    b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    computeHistoBreaks=seq(a, b, length=d)
                                }else{
                                    a<-min(All, na.rm=T)
                                    b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                                    d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    computeHistoBreaks=seq(a, b, length=d)
                                }
######### End of changes Will Greenwood (22.07.2009)

                                all<-hist(Other, plot=F, breaks=computeHistoBreaks)
                                pos<-hist(Pos, plot=F, breaks=computeHistoBreaks)
                                neg<-hist(Neg, plot=F, breaks=computeHistoBreaks)

                                plot(all, main=paste(plotTitle, "for plate", i, sep=" "), xlab=channel, 
                                cex.main=0.6)
                                par(fg="green")
                                lines(pos, angle=45, density=20)
                                par(fg="red")
                                lines(neg, angle=-45, density=20)
                                par(fg="black")
                                legend("topleft", c("Data", "Positive Controls", "Negative Controls"), 
                                angle=c(0, 45, -45), density=c(0, 20, 20), bty="n", cex=0.8)
                                mtext(paste(plotTitle, "for Experiment", j, sep=" "), side=3, outer=T, cex=0.8)
                            }
                        }
                    }
                }
            }
        }
    }

    for (j in minOfScreens:numOfScreens){
        if (length(which(dataset$ScreenNb == j))>0){

            subset<-dataset[which(dataset$ScreenNb == j), ]
            
            minOfPlates<-min(subset$LabtekNb)
            numOfPlates<-max(subset$LabtekNb)

            headerTemp<-strsplit(header[1], ",")
            histoName<-paste(headerTemp[[1]][2], plotTitle, sep="_")

            if (plotDesign == 1){
            
                pdf(paste(histoName, "Exp", j, "PerPlate", ".pdf", sep="_"))
                    c1<-ceiling(sqrt(length(unique(subset$LabtekNb))))
                    c2<-ceiling(length(unique(subset$LabtekNb))/c1)
                    par(mfrow=c(c1, c2), oma=c(0, 0, 2, 0))

    
                    for (i in minOfPlates:numOfPlates){
                        if (length(which(subset$LabtekNb == i))>0){
                    
                            subsubset<-subset[which(subset$LabtekNb == i), ]
         
                            All<-subsubset[[get("channel")]][which(subsubset$SpotType!=-1)]    
                            Neg<-subsubset[[get("channel")]][which(subsubset$SpotType == 0)]
                            Pos<-subsubset[[get("channel")]][which(subsubset$SpotType == 1)]
                            Other<-subsubset[[get("channel")]][which(subsubset$SpotType == 2)]

                            if (sum(!is.na(All))!=0){   
                            
#######
#                                if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                                    a<-min(All, na.rm=T)
#                                    b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                                    d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                                    computeHistoBreaks=seq(a, b, d)
#                                }else{
#                                    a<-min(All, na.rm=T)
#                                    b<-max(All, na.rm=T)+abs(a+0.5)+((max(All, na.rm=T)-a+1)/20)
#                                    d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                                    computeHistoBreaks=seq(a, b, d)
#                                }
######### Changes Will Greenwood (22.07.2009):
                                if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                                    a<-min(All, na.rm=T)
                                    b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    computeHistoBreaks=seq(a, b, length=d)
                                }else{
                                    a<-min(All, na.rm=T)
                                    b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                                    d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    computeHistoBreaks=seq(a, b, length=d)
                                }
######### End of changes Will Greenwood (22.07.2009)

                                all<-hist(Other, plot=F, breaks=computeHistoBreaks)
                                pos<-hist(Pos, plot=F, breaks=computeHistoBreaks)
                                neg<-hist(Neg, plot=F, breaks=computeHistoBreaks)

                                plot(all, main=paste(plotTitle, "for plate", i, sep=" "), xlab=channel, 
                                cex.main=0.7)
                                par(fg="green")
                                lines(pos, angle=45, density=20)
                                par(fg="red")
                                lines(neg, angle=-45, density=20)
                                par(fg="black")
                                legend("topleft", c("Data", "Positive Controls", "Negative Controls"), 
                                angle=c(0, 45, -45), density=c(0, 20, 20), bty="n", cex=0.5)
                                mtext(paste(plotTitle, "for Experiment", j, sep=" "), side=3, outer=T)
                            }else{
                                plot.new()
                                text(0.5, 0.75, "Cannot plot histogram for")
				text(0.5, 0.25, paste("exp ", j, " plate ", i, sep=""))
                                text(0.5, 0.05, "- only NAs available") 
                            }
                        }
                    }
                dev.off()
                
                png(paste(histoName, "Exp", j, "PerPlate", ".png", sep="_"))
                    c1<-ceiling(sqrt(length(unique(subset$LabtekNb))))
                    c2<-ceiling(length(unique(subset$LabtekNb))/c1)
                    par(mfrow=c(c1, c2), oma=c(0, 0, 2, 0))

    
                    for (i in minOfPlates:numOfPlates){
                        if (length(which(subset$LabtekNb == i))>0){
                    
                            subsubset<-subset[which(subset$LabtekNb == i), ]
                            All<-subsubset[[get("channel")]][which(subsubset$SpotType!=-1)]    
                            Neg<-subsubset[[get("channel")]][which(subsubset$SpotType == 0)]
                            Pos<-subsubset[[get("channel")]][which(subsubset$SpotType == 1)]
                            Other<-subsubset[[get("channel")]][which(subsubset$SpotType == 2)]

                            if (sum(!is.na(All))!=0){    
#######
#                                if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                                    a<-min(All, na.rm=T)
#                                    b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                                    d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                                    computeHistoBreaks=seq(a, b, d)
#                                }else{
#                                    a<-min(All, na.rm=T)
#                                    b<-max(All, na.rm=T)+abs(a+0.5)+((max(All, na.rm=T)-a+1)/20)
#                                    d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                                    computeHistoBreaks=seq(a, b, d)
#                                }
######### Changes Will Greenwood (22.07.2009):
                                if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                                    a<-min(All, na.rm=T)
                                    b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    computeHistoBreaks=seq(a, b, length=d)
                                }else{
                                    a<-min(All, na.rm=T)
                                    b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                                    d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    computeHistoBreaks=seq(a, b, length=d)
                                }
######### End of changes Will Greenwood (22.07.2009)

                                all<-hist(Other, plot=F, breaks=computeHistoBreaks)
                                pos<-hist(Pos, plot=F, breaks=computeHistoBreaks)
                                neg<-hist(Neg, plot=F, breaks=computeHistoBreaks)

                                plot(all, main=paste(plotTitle, "for plate", i, sep=" "), xlab=channel, 
                                cex.main=0.7)
                                par(fg="green")
                                lines(pos, angle=45, density=20)
                                par(fg="red")
                                lines(neg, angle=-45, density=20)
                                par(fg="black")
                                legend("topleft", c("Data", "Positive Controls", "Negative Controls"), 
                                angle=c(0, 45, -45), density=c(0, 20, 20), bty="n", cex=0.5)
                                mtext(paste(plotTitle, "for Experiment", j, sep=" "), side=3, outer=T)
                            }else{
                                plot.new()
                                text(0.5, 0.75, "Cannot plot histogram for")
				text(0.5, 0.25, paste("exp ", j, " plate ", i, sep=""))
                                text(0.5, 0.05, "- only NAs available")
                            }
                        }
                    }
                dev.off()
                
            }else{
                    
                for (i in minOfPlates:numOfPlates){
                    if (length(which(subset$LabtekNb == i))>0){
                    
                        subsubset<-subset[which(subset$LabtekNb == i), ]         
                        All<-subsubset[[get("channel")]][which(subsubset$SpotType!=-1)]
                        Neg<-subsubset[[get("channel")]][which(subsubset$SpotType == 0)]
                        Pos<-subsubset[[get("channel")]][which(subsubset$SpotType == 1)]
                        Other<-subsubset[[get("channel")]][which(subsubset$SpotType == 2)]
                        
                        if (sum(!is.na(All))!=0){
#######    
#                            if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                                a<-min(All, na.rm=T)
#                                b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                                d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                                computeHistoBreaks=seq(a, b, d)
#                            }else{
#                                a<-min(All, na.rm=T)
#                                b<-max(All, na.rm=T)+abs(a+0.5)+((max(All, na.rm=T)-a+1)/20)
#                                d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                                computeHistoBreaks=seq(a, b, d)
#                            }
######### Changes Will Greenwood (22.07.2009):
                            if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                                a<-min(All, na.rm=T)
                                b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                computeHistoBreaks=seq(a, b, length=d)
                            }else{
                                a<-min(All, na.rm=T)
                                b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                                d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                computeHistoBreaks=seq(a, b, length=d)
                            }
######### End of changes Will Greenwood (22.07.2009)

                            all<-hist(Other, plot=F, breaks=computeHistoBreaks)
                            pos<-hist(Pos, plot=F, breaks=computeHistoBreaks)
                            neg<-hist(Neg, plot=F, breaks=computeHistoBreaks)

                            pdf(paste(histoName, "Exp", j, "PerPlate", i, ".pdf", sep="_"))

                                plot(all, main=paste(plotTitle, "for plate", i, sep=" "), xlab=channel, 
                                cex.main=0.7)
                                par(fg="green")
                                lines(pos, angle=45, density=20)
                                par(fg="red")
                                lines(neg, angle=-45, density=20)
                                par(fg="black")
                                legend("topleft", c("Data", "Positive Controls", "Negative Controls"), 
                                angle=c(0, 45, -45), density=c(0, 20, 20), bty="n", cex=0.5)
                                mtext(paste(plotTitle, "for Experiment", j, sep=" "), side=3, outer=T)
                            dev.off()
                            
                            png(paste(histoName, "Exp", j, "PerPlate", i, ".png", sep="_"), 
                            width=300, height=300)

                                plot(all, main=paste(plotTitle, "for plate", i, sep=" "), xlab=channel, 
                                cex.main=0.5)
                                par(fg="green")
                                lines(pos, angle=45, density=20)
                                par(fg="red")
                                lines(neg, angle=-45, density=20)
                                par(fg="black")
                                legend("topleft", c("Data", "Positive Controls", "Negative Controls"), 
                                angle=c(0, 45, -45), density=c(0, 20, 20), bty="n", cex=0.5)
                                mtext(paste(plotTitle, "for Experiment", j, sep=" "), side=3, outer=T)
                            dev.off()
                        }else{
                            pdf(paste(histoName, "Exp", j, "PerPlate", i, ".pdf", sep="_"))
                                plot.new()
                                text(0.5, 0.75, "Cannot plot histogram for")
				text(0.5, 0.25, paste("exp ", j, " plate ", i, sep=""))
                                text(0.5, 0.05, "- only NAs available")
                            dev.off()
                            
                            png(paste(histoName, "Exp", j, "PerPlate", i, ".png", sep="_"), 
                            width=300, height=300)
                                plot.new()
                                text(0.5, 0.75, "Cannot plot histogram for")
				text(0.5, 0.25, paste("exp ", j, " plate ", i, sep=""))
                                text(0.5, 0.05, "- only NAs available")
                            dev.off()
                        }
                    }
                }
            }
        }
    }
    invisible(list(histoName, c(minOfScreens, numOfScreens), 
    c(minOfPlates, numOfPlates)))
}

plotControlHistoPerscreen<-function(header, dataset, channel, plotTitle, 
plotDesign, showPlot){

    dataset<-dataset[which(dataset$SpotType!=-1), ]
    
    minOfScreens<-min(dataset$ScreenNb)
    numOfScreens<-max(dataset$ScreenNb)
    
    if (showPlot == 1){
        if(interactive()){
    
            if (plotDesign == 1){
                c1<-ceiling(sqrt(length(unique(dataset$ScreenNb))))
                c2<-ceiling((length(unique(dataset$ScreenNb)))/c1)
                par(mfrow=c(c1, c2))
            }

            for (i in minOfScreens:numOfScreens){
                if (length(which(dataset$ScreenNb == i))>0){
    
                    subset<-dataset[which(dataset$ScreenNb == i), ]     
                    All<-subset[[get("channel")]][which(subset$SpotType!=-1)]
                    Neg<-subset[[get("channel")]][which(subset$SpotType == 0)]
                    Pos<-subset[[get("channel")]][which(subset$SpotType == 1)]
                    Other<-subset[[get("channel")]][which(subset$SpotType == 2)]

                    if (sum(!is.na(All)) == 0){
                        s1<-"Cannot plot histogram for Exp."
                        print(paste(s1, i, "(only NAs in dataset)", sep=" "))
                    }else{
#######
#                        if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                            a<-min(All, na.rm=T)
#                            b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                            d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                            computeHistoBreaks=seq(a, b, d)
#                        }else{
#                            a<-min(All, na.rm=T)
#                            b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
#                            d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                            computeHistoBreaks=seq(a, b, d)
#                        }
######### Changes Will Greenwood (22.07.2009):
                        if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                            a<-min(All, na.rm=T)
                            b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            computeHistoBreaks=seq(a, b, length=d)
                        }else{
                            a<-min(All, na.rm=T)
                            b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                            d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            computeHistoBreaks=seq(a, b, length=d)
                        }
######### End of changes Will Greenwood (22.07.2009)    
    
                        all<-hist(Other, plot=F, breaks=computeHistoBreaks)
                        pos<-hist(Pos, plot=F, breaks=computeHistoBreaks)
                        neg<-hist(Neg, plot=F, breaks=computeHistoBreaks)

                        plot(all, main=paste(plotTitle, "for Exp.", i, sep=" "), xlab=channel, 
                        cex.main=0.7)
                        par(fg="green")
                        lines(pos, angle=45, density=20)
                        par(fg="red")
                        lines(neg, angle=-45, density=20)
                        par(fg="black")
                        legend("topleft", c("Data", "Positive Controls", "Negative Controls"), 
                        angle=c(0, 45, -45), density=c(0, 20, 20), bty="n", cex=0.8)
                    }
                }
            }
        }
    }

    headerTemp<-strsplit(header[1], ",")
    histoName<-paste(headerTemp[[1]][2], plotTitle, sep="_")
    
    if (plotDesign == 1){
        pdf(paste(histoName, ".pdf", sep=""))
            c1<-ceiling(sqrt(length(unique(dataset$ScreenNb))))
            c2<-ceiling((length(unique(dataset$ScreenNb)))/c1)
            par(mfrow=c(c1, c2))
            for (i in minOfScreens:numOfScreens){
                if (length(which(dataset$ScreenNb == i))>0){
        
                    subset<-dataset[which(dataset$ScreenNb == i), ]
 
                    All<-subset[[get("channel")]][which(subset$SpotType!=-1)]
                    Neg<-subset[[get("channel")]][which(subset$SpotType == 0)]
                    Pos<-subset[[get("channel")]][which(subset$SpotType == 1)]
                    Other<-subset[[get("channel")]][which(subset$SpotType == 2)]

                    if (sum(!is.na(All))!=0){

#######
#                        if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                            a<-min(All, na.rm=T)
#                            b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                            d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                            computeHistoBreaks=seq(a, b, d)
#                        }else{
#                            a<-min(All, na.rm=T)
#                            b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
#                            d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                            computeHistoBreaks=seq(a, b, d)
#                        }
######### Changes Will Greenwood (22.07.2009):
                        if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                            a<-min(All, na.rm=T)
                            b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            computeHistoBreaks=seq(a, b, length=d)
                        }else{
                            a<-min(All, na.rm=T)
                            b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                            d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            computeHistoBreaks=seq(a, b, length=d)
                        }
######### End of changes Will Greenwood (22.07.2009) 
    
                        all<-hist(Other, plot=F, breaks=computeHistoBreaks)
                        pos<-hist(Pos, plot=F, breaks=computeHistoBreaks)
                        neg<-hist(Neg, plot=F, breaks=computeHistoBreaks)

                        plot(all, main=paste(plotTitle, "for Exp.", i, sep=" "), xlab=channel, 
                        cex.main=0.5)
                        par(fg="green")
                        lines(pos, angle=45, density=20)
                        par(fg="red")
                        lines(neg, angle=-45, density=20)
                        par(fg="black")
                        legend("topleft", c("Data", "Positive Controls", "Negative Controls"), 
                        angle=c(0, 45, -45), density=c(0, 20, 20), bty="n", cex=0.5)
                    }
                }
            }
        dev.off()
        
        png(paste(histoName, ".png", sep=""))

            c1<-ceiling(sqrt(length(unique(dataset$ScreenNb))))
	    c2<-ceiling((length(unique(dataset$ScreenNb)))/c1)
            par(mfrow=c(c1, c2))
            
            for (i in minOfScreens:numOfScreens){
                if (length(which(dataset$ScreenNb == i))>0){
        
                    subset<-dataset[which(dataset$ScreenNb == i), ] 
                    All<-subset[[get("channel")]][which(subset$SpotType!=-1)]
                    Neg<-subset[[get("channel")]][which(subset$SpotType == 0)]
                    Pos<-subset[[get("channel")]][which(subset$SpotType == 1)]
                    Other<-subset[[get("channel")]][which(subset$SpotType == 2)]

                    if (sum(!is.na(All))!=0){

#######
#                        if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                            a<-min(All, na.rm=T)
#                            b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                            d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                            computeHistoBreaks=seq(a, b, d)
#                        }else{
#                            a<-min(All, na.rm=T)
#                            b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
#                            d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                            computeHistoBreaks=seq(a, b, d)
#                        }
######### Changes Will Greenwood (22.07.2009):
                        if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                            a<-min(All, na.rm=T)
                            b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            computeHistoBreaks=seq(a, b, length=d)
                        }else{
                            a<-min(All, na.rm=T)
                            b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                            d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            computeHistoBreaks=seq(a, b, length=d)
                        }
######### End of changes Will Greenwood (22.07.2009) 
    
                        all<-hist(Other, plot=F, breaks=computeHistoBreaks)
                        pos<-hist(Pos, plot=F, breaks=computeHistoBreaks)
                        neg<-hist(Neg, plot=F, breaks=computeHistoBreaks)

                        plot(all, main=paste(plotTitle, "for Exp.", i, sep=" "), xlab=channel, 
                        cex.main=0.5)
                        par(fg="green")
                        lines(pos, angle=45, density=20)
                        par(fg="red")
                        lines(neg, angle=-45, density=20)
                        par(fg="black")
                        legend("topleft", c("Data", "Positive Controls", "Negative Controls"), 
                        angle=c(0, 45, -45), density=c(0, 20, 20), bty="n", cex=0.5)
                    }
                }
            }
        dev.off()

    }else{

        for (i in minOfScreens:numOfScreens){
            if (length(which(dataset$ScreenNb == i))>0){
        
                subset<-dataset[which(dataset$ScreenNb == i), ]

                All<-subset[[get("channel")]][which(subset$SpotType!=-1)]
                Neg<-subset[[get("channel")]][which(subset$SpotType == 0)]
                Pos<-subset[[get("channel")]][which(subset$SpotType == 1)]
                Other<-subset[[get("channel")]][which(subset$SpotType == 2)]

                if (sum(!is.na(All))!=0){
#######
#                    if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                        a<-min(All, na.rm=T)
#                        b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                        d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                        computeHistoBreaks=seq(a, b, d)
#                    }else{
#                        a<-min(All, na.rm=T)
#                        b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
#                        d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                        computeHistoBreaks=seq(a, b, d)
#                    }
######### Changes Will Greenwood (22.07.2009):
                    if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                        a<-min(All, na.rm=T)
                        b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                        d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                        computeHistoBreaks=seq(a, b, length=d)
                    }else{
                        a<-min(All, na.rm=T)
                        b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                        d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                        computeHistoBreaks=seq(a, b, length=d)
                    }
######### End of changes Will Greenwood (22.07.2009) 
    
                    all<-hist(Other, plot=F, breaks=computeHistoBreaks)
                    pos<-hist(Pos, plot=F, breaks=computeHistoBreaks)
                    neg<-hist(Neg, plot=F, breaks=computeHistoBreaks)

                    pdf(paste(histoName, "(Exp. ", i, ").pdf", sep=""))
                        plot(all, main=paste(plotTitle, "for Exp.", i, sep=" "), xlab=channel, 
                        cex.main=0.5)
                        par(fg="green")
                        lines(pos, angle=45, density=20)
                        par(fg="red")
                        lines(neg, angle=-45, density=20)
                        par(fg="black")
                        legend("topleft", c("Data", "Positive Controls", "Negative Controls"), 
                        angle=c(0, 45, -45), density=c(0, 20, 20), bty="n", cex=0.5)
                    dev.off()
                    
                    png(paste(histoName, "(Exp. ", i, ").png", sep=""), width=300, height=300)
                        plot(all, main=paste(plotTitle, "for Exp.", i, sep=" "), xlab=channel, 
                        cex.main=0.5)
                        par(fg="green")
                        lines(pos, angle=45, density=20)
                        par(fg="red")
                        lines(neg, angle=-45, density=20)
                        par(fg="black")
                        legend("topleft", c("Data", "Positive Controls", "Negative Controls"), 
                        angle=c(0, 45, -45), density=c(0, 20, 20), bty="n", cex=0.5)
                    dev.off()
                    
                }
            }
        }
    }
    invisible(list(histoName, minOfScreens, numOfScreens))
}






plotHisto<-function(header, dataset, channel, plotTitle, showPlot){

    All<-dataset[[get("channel")]][which(dataset$SpotType!=-1)]

    if (sum(!is.na(All)) == 0){
        stop("Cannot plot histogram (only NAs in dataset)")
    }

#######
#    if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#        a<-min(All, na.rm=T)
#        b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#        d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#        computeHistoBreaks=seq(a, b, d)
#    }else{
#        a<-min(All, na.rm=T)
#        b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
#        d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#        computeHistoBreaks=seq(a, b, d)
#    }
######### Changes Will Greenwood (22.07.2009):
    if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
        a<-min(All, na.rm=T)
        b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
        d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
        computeHistoBreaks=seq(a, b, length=d)
    }else{
        a<-min(All, na.rm=T)
        b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
        d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
        computeHistoBreaks=seq(a, b, length=d)
    }
######### End of changes Will Greenwood (22.07.2009) 

    if (showPlot == 1){
        if(interactive()){
            hist(All, breaks=computeHistoBreaks, main=plotTitle, xlab=channel, cex.main=0.8)
        }
    }
    
    headerTemp<-strsplit(header[1], ",")
    histoName<-paste(headerTemp[[1]][2], plotTitle, sep="_")
    pdf(paste(histoName, ".pdf", sep=""))
        hist(All, breaks=computeHistoBreaks, main=plotTitle, xlab=channel, cex.main=0.8)
    dev.off()
    
    png(paste(histoName, ".png", sep=""), width=300, height=300)
        hist(All, breaks=computeHistoBreaks, main=plotTitle, xlab=channel, cex.main=0.8)
    dev.off()
    invisible(histoName)
}

plotHistoPerplate<-function(header, dataset, channel, plotTitle, plotDesign, 
showPlot){

    dataset<-dataset[which(dataset$SpotType!=-1), ]
    
    numOfScreens<-max(dataset$ScreenNb)
    minOfScreens<-min(dataset$ScreenNb)

    if (showPlot == 1){
        if(interactive()){
            for (j in minOfScreens:numOfScreens){
                if (length(which(dataset$ScreenNb == j))>0){
        
                    subset<-dataset[which(dataset$ScreenNb == j), ]
                
                    minOfPlates<-min(subset$LabtekNb)
                    numOfPlates<-max(subset$LabtekNb)
            
                    if (plotDesign == 1){
                        x11()
                        c1<-ceiling(sqrt(length(unique(subset$LabtekNb))))
                        c2<-ceiling(length(unique(subset$LabtekNb))/c1)
                        par(mfrow=c(c1, c2), oma=c(0, 0, 2, 0))
                    }

                    for (i in minOfPlates:numOfPlates){
                        if (length(which(subset$LabtekNb == i))>0){
                
                            subsubset<-subset[which(subset$LabtekNb == i), ]     
                            All<-subsubset[[get("channel")]][which(subsubset$SpotType!=-1)]

                            if (sum(!is.na(All)) == 0){
                                s1<-"Cannot plot histogram for Exp."
                                print(paste(s1, j, "plate", i, "(only NAs in dataset)", sep=" "))
                            }else{
#######
#                                if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                                    a<-min(All, na.rm=T)
#                                    b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                                    d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                                    computeHistoBreaks=seq(a, b, d)
#                                }else{
#                                    a<-min(All, na.rm=T)
#                                    b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
#                                    d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                                    computeHistoBreaks=seq(a, b, d)
#                                }
######### Changes Will Greenwood (22.07.2009):
                                if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                                    a<-min(All, na.rm=T)
                                    b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    computeHistoBreaks=seq(a, b, length=d)
                                }else{
                                    a<-min(All, na.rm=T)
                                    b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                                    d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    computeHistoBreaks=seq(a, b, length=d)
                                }
######### End of changes Will Greenwood (22.07.2009) 
                        
                                hist(All, breaks=computeHistoBreaks, 
                                main=paste(plotTitle, "for plate", i, sep=" "), xlab=channel, cex.main=0.6)
                                mtext(paste(plotTitle, "for Experiment", j, sep=" "), side=3, outer=T, cex=0.8)
                            }
                        }
                    }
                }
            }
        }
    }

    for (j in minOfScreens:numOfScreens){
        if (length(which(dataset$ScreenNb == j))>0){

            subset<-dataset[which(dataset$ScreenNb == j), ]
            
            minOfPlates<-min(subset$LabtekNb)
            numOfPlates<-max(subset$LabtekNb)

            headerTemp<-strsplit(header[1], ",")
            histoName<-paste(headerTemp[[1]][2], plotTitle, sep="_")

            if (plotDesign == 1){

                pdf(paste(histoName, "Exp", j, "PerPlate", ".pdf", sep="_"))
                    c1<-ceiling(sqrt(length(unique(subset$LabtekNb))))
                    c2<-ceiling(length(unique(subset$LabtekNb))/c1)
                    par(mfrow=c(c1, c2), oma=c(0, 0, 2, 0))
    
                    for (i in minOfPlates:numOfPlates){
                        if (length(which(subset$LabtekNb == i))>0){
                    
                            subsubset<-subset[which(subset$LabtekNb == i), ]
                            All<-subsubset[[get("channel")]][which(subsubset$SpotType!=-1)]
    
                            if (sum(!is.na(All))!=0){
#######
#                                if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                                    a<-min(All, na.rm=T)
#                                    b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                                    d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                                    computeHistoBreaks=seq(a, b, d)
#                                }else{
#                                    a<-min(All, na.rm=T)
#                                    b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
#                                    d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                                    computeHistoBreaks=seq(a, b, d)
#                                }
######### Changes Will Greenwood (22.07.2009):
                                if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                                    a<-min(All, na.rm=T)
                                    b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    computeHistoBreaks=seq(a, b, length=d)
                                }else{
                                    a<-min(All, na.rm=T)
                                    b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                                    d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    computeHistoBreaks=seq(a, b, length=d)
                                }
######### End of changes Will Greenwood (22.07.2009) 

                                hist(All, breaks=computeHistoBreaks, 
                                main=paste(plotTitle, "for plate", i, sep=" "), xlab=channel, cex.main=0.7)
                                mtext(paste(plotTitle, "for Experiment", j, sep=" "), side=3, outer=T)
                            }
                        }
                    }
                dev.off()
                
                png(paste(histoName, "Exp", j, "PerPlate", ".png", sep="_"))
                    
                    c1<-ceiling(sqrt(length(unique(subset$LabtekNb))))
                    c2<-ceiling(length(unique(subset$LabtekNb))/c1)
                    par(mfrow=c(c1, c2), oma=c(0, 0, 2, 0))
    
                    for (i in minOfPlates:numOfPlates){
                        if (length(which(subset$LabtekNb == i))>0){
                    
                            subsubset<-subset[which(subset$LabtekNb == i), ]
                            All<-subsubset[[get("channel")]][which(subsubset$SpotType!=-1)]
    
                            if (sum(!is.na(All))!=0){
#######
#                                if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                                    a<-min(All, na.rm=T)
#                                    b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                                    d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                                    computeHistoBreaks=seq(a, b, d)
#                                }else{
#                                    a<-min(All, na.rm=T)
#                                    b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
#                                    d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                                    computeHistoBreaks=seq(a, b, d)
#                                }
######### Changes Will Greenwood (22.07.2009):
                                if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                                    a<-min(All, na.rm=T)
                                    b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    computeHistoBreaks=seq(a, b, length=d)
                                }else{
                                    a<-min(All, na.rm=T)
                                    b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                                    d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                    computeHistoBreaks=seq(a, b, length=d)
                                }
######### End of changes Will Greenwood (22.07.2009) 

                                hist(All, breaks=computeHistoBreaks, 
                                main=paste(plotTitle, "for plate", i, sep=" "), xlab=channel, cex.main=0.7)
                                mtext(paste(plotTitle, "for Experiment", j, sep=" "), side=3, outer=T)
                            }
                        }
                    }
                dev.off()
                
            }else{
                for (i in minOfPlates:numOfPlates){
                    if (length(which(subset$LabtekNb == i))>0){
                    
                        subsubset<-subset[which(subset$LabtekNb == i), ]         
                        All<-subsubset[[get("channel")]][which(subsubset$SpotType!=-1)]
                        
                        if (sum(!is.na(All))!=0){
#######   
#                            if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                                a<-min(All, na.rm=T)
#                                b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                                d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                                computeHistoBreaks=seq(a, b, d)
#                            }else{
#                                a<-min(All, na.rm=T)
#                                b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
#                                d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                                computeHistoBreaks=seq(a, b, d)
#                            }
######### Changes Will Greenwood (22.07.2009):
                            if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                                a<-min(All, na.rm=T)
                                b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                computeHistoBreaks=seq(a, b, length=d)
                            }else{
                                a<-min(All, na.rm=T)
                                b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                                d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                                computeHistoBreaks=seq(a, b, length=d)
                            }
######### End of changes Will Greenwood (22.07.2009) 

                            pdf(paste(histoName, "Exp", j, "PerPlate", i, ".pdf", sep=""))
                                hist(All, breaks=computeHistoBreaks, 
                                main=paste(plotTitle, "for plate", i, sep=" "), xlab=channel, cex.main=0.7)
                                mtext(paste(plotTitle, "for Experiment", j, sep=" "), side=3, outer=T)
                            dev.off()
                            
                            png(paste(histoName, "Exp", j, "PerPlate", i, ".png", sep=""), 
                            width=300, height=300)
                                hist(All, breaks=computeHistoBreaks, 
                                main=paste(plotTitle, "for plate", i, sep=" "), xlab=channel, cex.main=0.7)
                                mtext(paste(plotTitle, "for Experiment", j, sep=" "), side=3, outer=T)
                            dev.off()
                        }
                    }
                }            
            }
        }
    }
    invisible(list(histoName, c(minOfScreens, numOfScreens), 
    c(minOfPlates, numOfPlates)))
}


plotHistoPerscreen<-function(header, dataset, channel, plotTitle, plotDesign, 
showPlot){
    
    dataset<-dataset[which(dataset$SpotType!=-1), ]
    
    minOfScreens<-min(dataset$ScreenNb)
    numOfScreens<-max(dataset$ScreenNb)
    
    if (showPlot == 1){
        if(interactive()){
    
            if (plotDesign == 1){
                c1<-ceiling(sqrt(length(unique(dataset$ScreenNb))))
                c2<-ceiling((length(unique(dataset$ScreenNb)))/c1)
                par(mfrow=c(c1, c2))
            }

            for (i in minOfScreens:numOfScreens){
                if (length(which(dataset$ScreenNb == i))>0){
    
                    subset<-dataset[which(dataset$ScreenNb == i), ]
                    All<-subset[[get("channel")]][which(subset$SpotType!=-1)]

                    if (sum(!is.na(All)) == 0){
                        s1<-"Cannot plot histogram for Exp."
                        print(paste(s1, i, "(only NAs in dataset)", sep=" "))
                    }else{
#######                     
#                        if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                            a<-min(All, na.rm=T)
#                            b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                            d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                            computeHistoBreaks=seq(a, b, d)
#                        }else{
#                            a<-min(All, na.rm=T)
#                            b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
#                            d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                            computeHistoBreaks=seq(a, b, d)
#                        }
######### Changes Will Greenwood (22.07.2009):
                        if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                            a<-min(All, na.rm=T)
                            b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            computeHistoBreaks=seq(a, b, length=d)
                        }else{
                            a<-min(All, na.rm=T)
                            b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                            d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            computeHistoBreaks=seq(a, b, length=d)
                        }
######### End of changes Will Greenwood (22.07.2009) 
    
                        hist(All, breaks=computeHistoBreaks, 
                        main=paste(plotTitle, "for Exp.", i, sep=" "), xlab=channel, cex.main=0.8)
                    }
                }
            }
        }
    }

    headerTemp<-strsplit(header[1], ",")
    histoName<-paste(headerTemp[[1]][2], plotTitle, sep="_")
    
    if (plotDesign == 1){
        pdf(paste(histoName, ".pdf", sep=""))

            c1<-ceiling(sqrt(length(unique(dataset$ScreenNb))))
            c2<-ceiling((length(unique(dataset$ScreenNb)))/c1)
            par(mfrow=c(c1, c2))
            for (i in minOfScreens:numOfScreens){
                if (length(which(dataset$ScreenNb == i))>0){
        
                    subset<-dataset[which(dataset$ScreenNb == i), ]
 
                    All<-subset[[get("channel")]][which(subset$SpotType!=-1)]

                    if (sum(!is.na(All))!=0){
#######                     
#                        if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
#                            a<-min(All, na.rm=T)
#                            b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
#                            d<-round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)
#                            computeHistoBreaks=seq(a, b, d)
#                        }else{
#                            a<-min(All, na.rm=T)
#                            b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
#                            d<-(max(All, na.rm=T)-min(All, na.rm=T)+1)/20
#                            computeHistoBreaks=seq(a, b, d)
#                        }
######### Changes Will Greenwood (22.07.2009):
                        if (round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20)!=0){
                            a<-min(All, na.rm=T)
                            b<-max(All, na.rm=T)+(round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            d<-ceiling((b-a)/round((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            computeHistoBreaks=seq(a, b, length=d)
                        }else{
                            a<-min(All, na.rm=T)
                            b<-max(All, na.rm=T)+abs(min(All, na.rm=T)+0.5)
                            d<-ceiling((b-a)/((max(All, na.rm=T)-min(All, na.rm=T)+1)/20))
                            computeHistoBreaks=seq(a, b, length=d)
                        }
######### End of changes Will Greenwood (22.07.2009) 
    
                        hist(All, breaks=computeHistoBreaks, 
                        main=paste(plotTitle, "for Exp.", i, sep=" "), xlab=channel, cex.main=0.6)
                    }
                }
            }
        dev.off()
        
        png(paste(histoName, ".png", sep=""))

            c1<-ceiling(sqrt(length(unique(dataset$ScreenNb))))
            c2<-