SafeQuant: A Toolbox for the Analysis of Proteomics Data

Documented in plotExpDesign plotVolcano

# TODO: Add comment
#
# Author: erikahrne
###############################################################################

### SET COLOR VECTOR
#' color vector
#' @export
COLORS <- as.character(c(
				"red"
				,"darkgreen"
				,"blue"
				,"darkmagenta"
				,"darkorange"
				,"darkblue"
				,"black"
				,"brown"
				,"darkgrey"
				,"red"
				,"brown"
				,rev(colors())) ### if that's not enough
)


#' @export
#' @import Biobase
#' @importFrom graphics plot par grid lines points legend
.idOverviewPlots <- function(userOptions=userOptions
		,esetNorm=esetNorm
		,fileType=fileType
		,sqaPeptide=sqaPeptide
		,sqaProtein=sqaProtein){

	# HACK TO avoid check error. See function call in safeQuant.R @TODO re-structure this messsy function
	if(is.na(sqaPeptide)){rm(sqaPeptide)}
	if(is.na(sqaProtein)){rm(sqaProtein)}

	######################## OVERVIEW PLOT
	fdr <-userOptions$fdrCutoff
	nbPSM <- sum(!fData(esetNorm)$isFiltered,na.rm=T)
	cex=1
	par(cex.names=1.5, cex.axis= 1.25, cex.lab= 1.25, mar=c(5.1,3.1,4.1,1.1))
	if(userOptions$verbose) cat("IDENTIFICATIONS OVERVIEW PLOT \n")
	plot(0,0,type="n",xlim=c(0,8),ylim=c(-1,1.5), main=paste("\nIDENTIFICATIONS OVERVIEW\n (FDR ",fdr,")",sep=""), axes=F, xlab="", ylab="", cex=cex)
#text(0.5,2,paste("FDR",fdr),cex=cex, pos=4)
	xPos <- -0.5
	yPos <- 1

	if(fileType != "ProgenesisProtein"){
		text(xPos,yPos,paste(nbPSM," PSM"),cex=cex, pos=4)
		yPos <- yPos - 0.5

	}

	if(exists("sqaPeptide")){

		isMod <- nchar(as.character(fData(sqaPeptide$eset)$ptm)) > 0
		nbPeptides <- sum(!fData(sqaPeptide$eset)$isFiltered,na.rm=T)
		nbUnModPeptides <- sum(!fData(sqaPeptide$eset[!isMod,])$isFiltered,na.rm=T)
		nbModPeptides <- sum(!fData(sqaPeptide$eset[isMod,])$isFiltered,na.rm=T)
		nbProteins <-  length(unique(fData(sqaPeptide$eset)[!fData(sqaPeptide$eset)$isFiltered,]$proteinName))

		text(xPos,yPos,paste(nbPeptides, " PEPTIDES"),cex=cex, pos=4)
		yPos <- yPos - 0.25
		text(xPos+0.5,yPos,paste(nbUnModPeptides," UN-MOD."),cex=cex, pos=4)
		yPos <- yPos - 0.25
		text(xPos+0.5,yPos,paste(nbModPeptides," MOD."),cex=cex, pos=4)
		yPos <- yPos - 0.5



	}else if("peptide" %in% names(fData(esetNorm))) { ### TMT export
		text(xPos,yPos,paste(length(unique(fData(esetNorm)$peptide))," PEPTIDES" ),cex=cex, pos=4)
		yPos <- yPos - 0.5
	}

	if(exists("sqaProtein")){
		nbProteins <- sum(!fData(sqaProtein$eset)$isFiltered,na.rm=T)
	}
	text(xPos,yPos,paste(nbProteins," PROTEINS"),cex=cex, pos=4)
	par(mar=c(5.1,4.1,4.1,2.1))
	######################## OVERVIEW PLOT END

	### charge state
	if("charge" %in% names(fData(esetNorm)) & (sum(!is.na(fData(esetNorm)$charge),na.rm=T)>0) ){
		if(userOptions$verbose) cat("CHARGE STATE PLOT \n")
		chargeTable <- table(fData(esetNorm)$charge[!fData(esetNorm)$isFiltered] )
		barplot2(chargeTable, xlab="Charge State", ylab="PSM Counts", col="blue", plot.grid = TRUE, grid.col="lightgrey")

	}
	if(exists("sqaPeptide")){
		if(userOptions$verbose) cat("INFO: NB. MIS-CLEAVAGES PLOT \n")
		### mis-cleavage
		sel <- !fData(esetNorm)$isFiltered
		nbRows <- sum(sel,na.rm=T)
		nbSel <-min(c(nbRows,500))
		nbMCTable <- (table(getNbMisCleavages(fData(esetNorm)$peptide[sel][sample(nbRows,nbSel)] ))/nbSel)*100
		barplot2(nbMCTable, xlab="Nb. Mis-cleavages", ylab="Peptide Counts (%)", col="blue", plot.grid = TRUE, grid.col="lightgrey")

		### don't plot if many (more than 10) NA motifs (NA motid due to wrongly specified fasta)
		if("motifX" %in% names(fData(sqaPeptide$eset)) & (sum(is.na(fData(sqaPeptide$eset)$motifX)) < 10 ) ){

			if(userOptions$verbose) cat("INFO: MOTIF-X PLOT \n")
			motifTable <- table(.getUniquePtmMotifs(sqaPeptide$eset,format=(fileType == "ScaffoldTMT")+1)$ptm)

			if(nrow(motifTable) > 0){ # make sure some non NA motifs were found
				bp <- barplot2(motifTable, ,ylab="Modif. Site Counts", col="blue", plot.grid = TRUE, xaxt="n", grid.col="lightgrey")
				mtext(names(motifTable),side=1,at=bp[,1], line=0.2, las=2,cex=0.6)
			}

		}else{

			if(fileType == "ScaffoldTMT"){

				ptmTag <- as.character(fData(sqaPeptide$eset)$ptm)[!fData(sqaPeptide$eset)$isFiltered]
				ptmTag[(nchar(ptmTag) == 0)] <- "Unmod"
				ptmTag <- gsub("[0-9]","",unlist(strsplit(ptmTag,"\\, ")))
				ptmTable <- table(ptmTag[!fData(sqaPeptide$eset)$isFiltered] )

			}else{
				### ptm PROGENSIS
				ptmTag <- as.character(fData(sqaPeptide$eset)$ptm)[!fData(sqaPeptide$eset)$isFiltered]
				ptmTag[nchar(ptmTag) == 0] <- "Unmod"
				ptmTag <- gsub("\\[[0-9]*\\] {1,}","",unlist(strsplit(ptmTag,"\\|")))
				ptmTable <- table(ptmTag[!fData(sqaPeptide$eset)$isFiltered] )

			}
			if(userOptions$verbose) cat("PTM PLOT \n")
			#barplot2(ptmTable, ylab="Peptide Counts", col="blue", plot.grid = TRUE, las=2, grid.col="lightgrey", cex.names=0.7)
			bp <- barplot2(ptmTable, ylab="Peptide Counts", col="blue", plot.grid = TRUE, xaxt="n", grid.col="lightgrey")
			mtext(names(ptmTable),side=1,at=bp[,1], line=0, cex=0.6, las=2)

		}

		if("nbPtmsPerPeptide"  %in% names(fData(sqaPeptide$eset))){
			if(userOptions$verbose) cat("INFO: NB.PTMS PER PEPTIDE PLOT \n")
			### nb. ptms per peptide
			nbPtmPerPeptideTable <- table(fData(sqaPeptide$eset)$nbPtmsPerPeptide[!fData(sqaPeptide$eset)$isFiltered])
			barplot2(nbPtmPerPeptideTable, xlab="Nb. PTM per Peptide", ylab="Peptide Counts", col="blue", plot.grid = TRUE, grid.col="lightgrey")
		}

		cysteinFreqBarplot(fData(sqaPeptide$eset)$peptide)

	}else if("peptide" %in% names(fData(esetNorm))){
		if(userOptions$verbose) cat("INFO: NB. MIS-CLEAVAGES PLOT 2 \n")
		### mis-cleavage
		#nbMCTable <- table(getNbMisCleavages(fData(esetNorm)$peptide, protease="trypsin")[!fData(esetNorm)$isFiltered] )
		#barplot2(nbMCTable, xlab="Nb. Mis-cleavages", ylab="PSM Counts", col="blue", plot.grid = TRUE, grid.col="lightgrey")
		### mis-cleavage
		sel <- !fData(esetNorm)$isFiltered
		nbRows <- sum(sel,na.rm=T)
		nbSel <-min(c(nbRows,500))
		nbMCTable <- (table(getNbMisCleavages(fData(esetNorm)$peptide[sel][sample(nbRows,nbSel)] ))/nbSel)*100
		barplot2(nbMCTable, xlab="Nb. Mis-cleavages", ylab="Peptide Counts (%)", col="blue", plot.grid = TRUE, grid.col="lightgrey")

	}

	#### peptides per protein
	if(userOptions$verbose) cat("PEPTIDES PER PROTEIN PLOT\n")
	if(exists("sqaProtein")){
		peptidesPerProtein <- fData(sqaProtein$eset)$nbPeptides[!fData(sqaProtein$eset)$isFiltered]
	}else{
		peptidesPerProtein <- fData(sqaPeptide$eset)$nbPeptides[!fData(sqaPeptide$eset)$isFiltered]
		peptidesPerProtein <- peptidesPerProtein[unique(names(peptidesPerProtein))]
	}



	### discard filtered out proteins
	peptidesPerProtein <- peptidesPerProtein[peptidesPerProtein > 0]
	#counts <- max(c(min(peptidesPerProtein,na.rm=T),1),na.rm=T):max(c(max(peptidesPerProtein,na.rm=T),2))
	xPeptides <- min(peptidesPerProtein,na.rm=T):max(c(max(peptidesPerProtein,na.rm=T),10))
	yCount <- unlist(lapply(xPeptides,function(t){
						sum(unlist(peptidesPerProtein) == t,na.rm=T) }))

	yCount[yCount == 0] <- NA
	plot(xPeptides,yCount,type="n", log="x",xlab="Peptides Per Protein", ylab="Protein Counts")
	grid()
	lines(xPeptides,yCount,type="h",col="blue",lwd=2.5)

	### Id's vs Retention Time
	if(exists("sqaPeptide") && ("retentionTime"  %in% names(fData(sqaPeptide$eset)))) plotNbIdentificationsVsRT(sqaPeptide$eset)

}

### some quality control plots
#' @export
.qcPlots <- function(eset,selection=1:7,nbFeatures=500, userOptions=userOptions, ...){

	if(nrow(eset) < nbFeatures) nbFeatures <- nrow(eset)
	sel <- sample(nrow(eset),nbFeatures,replace=F)

	par(mfrow=c(2,2), mar=c(5.1,4.1,4.1,2.1))
	if(1 %in% selection) {
		barplotMSSignal(eset, ...)
		plotMSSignalDistributions(log2(exprs(eset)), col=as.character(.getConditionColors(eset)[pData(eset)$condition,]), lwd=1.5, ...)
		boxplot(getAllCV(eset)*100,  col=as.character(.getConditionColors(eset)[pData(eset)$condition,]))

	}

	par(mfrow=c(1,1), mar=c(5.1,4.1,4.1,2.1))

	if(3 %in% selection)pairsAnnot(log2(exprs(eset)[sel,order(pData(eset)$condition)]), ...)

	if((4 %in% selection) && ("pMassError" %in% names(fData(eset)))){
		plotPrecMassErrorDistrib(eset,pMassTolWindow=userOptions$precursorMassFilter)
	}

	### retention time normalization plot
	if((5 %in% selection) && ("retentionTime" %in% names(fData(eset))) ){
		plotRTNormSummary(getRTNormFactors(eset, minFeaturesPerBin=100))
	}
	#if(6 %in% selection) invisible(hClustHeatMap(eset[sel & !fData(eset)$isFiltered,],main=paste("SELECTED NB. FEATURES:", nbFeatures), ...))
}

### some quality control plots
#' @export
.idPlots <- function(eset,selection=1:7, qvalueThrs=0.01,...){

	if(!is.null(fData(eset)$idScore)){
		isDec <- isDecoy(fData(eset)$proteinName)
		scores <- fData(eset)$idScore
		qvalues <- fData(eset)$idQValue

		if(1 %in% selection) {
			plotScoreDistrib(scores[!isDec],scores[isDec], ...)

#			### score cutoffs abline
#			if("protein level"){
#				cutOff <- min(fData(eset)$idScore[fData(eset)$idQValue < qvalueThrs],na.rm=T )
#				if(is.finite(cutOff)) abline(v=cutOff, col="blue")
#			}else{
#				isMod <- nchar(as.character(fData(eset)$ptm)) > 0
#				modCutOff <- min(fData(eset[isMod,])$idScore[fData(eset[isMod,])$idQValue < qvalueThrs],na.rm=T )
#				noModCutOff <- min(fData(eset[!isMod,])$idScore[fData(eset[!isMod,])$idQValue < qvalueThrs],na.rm=T )
#				if(is.finite(modCutOff)) abline(v=modCutOff, col="blue")
#				if(is.finite(noModCutOff))abline(v=noModCutOff, "darkgreen")
#				legend("right",c("unmod. cutoff","mod. cutoff"),col=c("blue","darkgreen"), lty=1)
#			}

		}

		if(2 %in% selection) plotIdScoreVsFDR(scores,qvalues,qvalueThrs, ...)
		if(3 %in% selection) plotROC(qvalues,qvalueThrs,...)
	}
}

#' @export
.getConditionColors <- function(eset){
	return(data.frame(colors=as.character(COLORS[1:length(levels(pData(eset)$condition))]), row.names=levels(pData(eset)$condition)))
}

### color strip for volcano plot
#' @export
#' @importFrom graphics image mtext axis
.dotColorstrip <- function(colors,minSignal=0, maxSignal = maxSignal,lab= "C.V. (%)"  )
{
	bottom <- 1
	count <- length(colors)
	m <- matrix(1:count, count, 1)
	image(m, col=colors, ylab="", axes=FALSE)

	labels <- round(seq(minSignal,maxSignal,length.out=3))

	### round off to nearest 5
	labels <- as.character( round(labels / 5) *5 )

	at <- seq(0,1,length.out=3)

	axis(bottom,tick=TRUE, labels=labels , at=at )
	mtext(lab, bottom, adj=0.5, line=2)
}

### called from plotVolcano. Creates volcano plot form data.frame input
#' @export
#' @importFrom graphics plot par grid lines points
.plotVolcano <- function(d
		, ratioCutOffAbsLog2=0
		, absLog10pValueCutOff=2
		, xlim = range(d[,1],na.rm=T)
		, ylim = range(abs(        log10(d[,2])[is.finite(  log10(d[,2])  )]    )      ,na.rm=T) # pValue or qValue
		, maxSignal = max(d[,3][is.finite(d[,3])],na.rm=T)
		#, higlightSel = rep(F,nrow(d))
		, controlCondition = "control"
		, caseCondition = "case"
		,...

){

	colorPalette <- rev(rich.colors(32))

	# no p/q-values
	if(!all(is.finite(ylim))){
		ylim <-c(0,1)
	}

	# if d[,3] is all NA
	if(all(!is.finite(d[,3]))){
		dotColors <- rep("darkgrey",nrow(d))
	}else{
		dotColors <- colorPalette[round(1+ d[,3]/maxSignal *31)]
	}

	par(fig=c(0,1,0.18,1))
	plot(0,0
			, xlab= paste("log2(",caseCondition,"/",controlCondition,")", sep="" )
			, ylab= paste("-log10(",names(d)[2],")",sep="")
			, type="n", ylim=ylim,xlim=xlim
			,...
	)

	grid()

	### d[,2] qValues or pValues
	points(d[,1], abs(log10(d[,2])), col=dotColors, pch=20)

	### draw valid squares contianing valid proteins/peptides
	lines(c(-ratioCutOffAbsLog2,-ratioCutOffAbsLog2),c(absLog10pValueCutOff,1000), col="grey")
	lines(c(ratioCutOffAbsLog2,ratioCutOffAbsLog2),c(absLog10pValueCutOff,1000), col="grey")
	lines(c(-1000,-ratioCutOffAbsLog2),c(absLog10pValueCutOff,absLog10pValueCutOff), col="grey")
	lines(c(ratioCutOffAbsLog2,1000),c(absLog10pValueCutOff,absLog10pValueCutOff), col="grey")

	# @TODO remove this feature
	#points(d$ratio[higlightSel],abs(log10(d[,2]))[higlightSel],col="darkgrey")

	### add heat color bar
	if(!all(!is.finite(d[,3]))){
		par(fig=c(0,1,0,0.3), new=TRUE)
		.dotColorstrip(colorPalette, maxSignal=maxSignal*100)
		par(mfrow=c(1,1))
	}
}

#' @export
#' @importFrom graphics arrows
.errorBar <- function(x, y, upper, lower=upper, length=0.1,...){
	if(length(x) != length(y) | length(y) !=length(lower) | length(lower) != length(upper))
		stop("vectors must be same length")
	arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
}

#' @export
#' @importFrom graphics legend plot par
.outerLegend <- function(...) {
	opar <- par(fig=c(0, 1, 0, 1), oma=c(0, 0, 0, 0),
			mar=c(0.5, 0.5, 0.5, 1), new=TRUE)
	on.exit(par(opar))
	plot(0, 0, type='n', bty='n', xaxt='n', yaxt='n')
	legend(...)
}

#' @export
#' @import corrplot
#' @importFrom graphics legend
.correlationPlot <- function(d, textCol="black", labels=colnames(d),... ){

	corrplot(cor(d,use="complete")^2
			, type = "upper"
			, col = colorRampPalette(c("white","white","white","lightblue","blue"))(100)
			, cl.lim = c(0, 1)
			, tl.col = textCol
			#, method="circle"
			, method="pie"
			, addgrid.col="white"
			,...
	)
#	legend("bottomleft"
#			, labels
#			, fill=unique(textCol)
#			#, fill=textCol
#			, box.col="white"
#	)
#	legend("topright"
#			, c("","",as.expression(bquote(R^2*"        ")))
#			#, fill=unique(textCol)
#			#, fill=textCol
#			, box.col="white"
#	)
}

#' @export
.allpValueHist <- function(sqa, col=as.character(.getConditionColors(sqa$eset)[names(sqa$pValue ),])){

	i <- 1
	for(cond in names(sqa$pValue ) ){
		hist(sqa$pValue[,cond],breaks=seq(0,1,length=50), main=cond, xlab="pValue",col=col[i])
		i <- i +1
	}
}

#' Plots volcano, data points colored by max cv of the 2 compared conditions
#' @param obj safeQuantAnalysis object or data.frame
#' @param adjusted TRUE/FALSE plot qValues or pValues on y-axis
#' @param ratioThrs default 1
#' @param pValueThreshold default 0.01
#' @param ... see plot
#' @import Biobase gplots
#' @export
#' @note  No note
#' @details data.frame input object should contain 3 columns (ratio,qValue,cv)
#' @references NA
#' @examples print("No examples")
plotVolcano <- function(obj
		,ratioThrs=1
		,pValueThreshold=0.01
		,adjusted = T

		,...
){

	ratioCutOffAbsLog2 <- abs(log2(ratioThrs))
	absLog10pValueCutOff <- abs(log10(pValueThreshold))

	### plot volcanos for all case control comparisons
	if(class(obj) == "safeQuantAnalysis"){

		### need at least two condiotions
		if(length(unique(pData(obj$eset)$condition)) == 1){
			return(cat("INFO: Only one condition no plotVolcano \n"))
		}

		# ensure the same range on all volcano plots
		xlim <- range(obj$ratio, na.rm=T)

		if(adjusted){
			ylim <- range(abs(log10(obj$qValue)),na.rm=T)
		}else{
			ylim <- range(abs(log10(obj$pValue)),na.rm=T)
		}

		# ensure same scale on color legend
		cvMax <- max(as.vector(unlist(obj$cv)),na.rm=T)

		controlCondition <- .getControlCondition(obj$eset)
		for(caseCondition in colnames(obj$pValue)){

			### create data.frame listing all data points (ratio,pvalue,cv)
			if(adjusted){

				d <- data.frame( ratio=obj$ratio[,caseCondition]
						,qValue=obj$qValue[,caseCondition]
						,cv=apply(obj$cv[,c(caseCondition,controlCondition)],1,max,na.rm=T)
				)
			}else{
				d <- data.frame( ratio=obj$ratio[,caseCondition]
						,pValue=obj$pValue[,caseCondition]
						,cv=apply(obj$cv[,c(caseCondition,controlCondition)],1,max,na.rm=T)
				)
			}

			.plotVolcano(d
					, ratioCutOffAbsLog2=ratioCutOffAbsLog2
					, absLog10pValueCutOff=absLog10pValueCutOff
					, xlim = xlim
					, ylim = ylim # pValue or qValue
					, maxSignal = cvMax
					#, higlightSel = higlightSel
					, controlCondition = controlCondition
					, caseCondition = caseCondition

					,...)
		}
	}else{
		.plotVolcano(obj
				, ratioCutOffAbsLog2=ratioCutOffAbsLog2
				, absLog10pValueCutOff=absLog10pValueCutOff
				#, higlightSel = higlightSel

				,...)
	}

}

#' Display experimental design, high-lighting the control condition
#' @param eset ExpressionSet
#' @param condColors condition colors
#' @param version version number
#' @import Biobase
#' @importFrom graphics text plot
#' @export
#' @note  No note
#' @details No details
#' @references NA
#' @examples print("No examples")
plotExpDesign <- function(eset, condColors=.getConditionColors(eset),  version="X"){

	### plot ctrl at the bottom
	pData(eset) <- rbind(pData(eset)[pData(eset)$isControl,],pData(eset)[!pData(eset)$isControl,])

	conditionNames <- as.character(unique(pData(eset)$condition))
	nbConditions <- length(conditionNames)
	controlCondition = .getControlCondition(eset)

	sampleNames <- rownames(pData(eset))
	nbSamples <- nrow(pData(eset))

	xlim <- c(-1,6)
	ylim <- c(-2,nbSamples+2)

	plot(0,0,type="n", xlim=xlim, ylim=ylim, main="Experimental Design", axes=FALSE, xlab="", ylab="")

	condYPosStep <- (nbSamples+2)/(nbConditions+1)
	sampleNb <- 1

	for(condNb in 1:length(conditionNames)){

		condName <- conditionNames[condNb]
		condCol = as.character(condColors[condName,])

		### control condition in  and underline
		if(condName == controlCondition){
			text(1,(condNb)*condYPosStep,bquote(underline(bold(.(condName)))), col=condCol)
		}else{
			text(1,(condNb)*condYPosStep,condName, col=condCol)

		}

		for(i in sampleNames[as.character(pData(eset)$condition) == condName]){
			text(4,sampleNb,paste(paste(sampleNb,":",sep=""), sampleNames[sampleNb]), col=condCol)
			sampleNb <- sampleNb + 1
		}
	}

	text(-1,-2,paste("SafeQuant v.", version), pos=4)

}

#' Plot lower triangle Pearsons R^2. Diagonal text, upper triangle all against all scatter plots with lm abline
#' @param data data.frame
#' @param textCol text color
#' @param diagText diagnoal text
#' @param col dot col
#' @param isHeatCol heat colors
#' @param cexTxt cex txt
#' @param ... see plot
#' @note  No note
#' @export
#' @importFrom graphics pairs abline par points
#' @importFrom stats cor
#' @details No details
#' @references NA
#' @examples print("No examples")
pairsAnnot<-
		function(data,textCol=rep(1,ncol(data)), diagText=colnames(data),col= rgb(0,100,0,50,maxColorValue=255),isHeatCol=F,cexTxt=2,...) {

	### we need at least two samples
	if(ncol(data.frame(data)) < 2 ){
		return(warning("INFO: Only one sample no pairsAnnot \n"))
	}

	### cex as a function of numbers of columns
	cex <- 0.8
	if(ncol(data) < 6){
		cex <- 2
	}else if(ncol(data) < 9){
		cex <- 1.4
	}
	else if(ncol(data) < 12){
		cex <- 1
	}

	count <- 1

	panel.lm <-
			function (x, y, col = par("col"), bg = NA, pch = par("pch"),
					#	cex = 1, col.lm = "red", lwd=par("lwd"), ...)
					col.lm = "red", lwd=par("lwd"))
	{

		if(isHeatCol){
			df <- data.frame(x,y)
			## Use densCols() output to get density at each point
			densityCol <- densCols(x,y, colramp=colorRampPalette(c("black", "white")))
			df$densityCol <- col2rgb(densityCol)[1,] + 1L
			## Map densities to colors
			cols <-  colorRampPalette(c("#000099", "#00FEFF", "#45FE4F",
							"#FCFF00", "#FF9400", "#FF3100"))(256)

			# HACK to please CRAN CHECK "rollUp: no visible binding for global variable hCol"
			hCol <- NULL
			df$hCol <- cols[df$densityCol]

			points(y~x, data=df[order(df$densityCol),], pch = 20, col =hCol, bg = bg)
			col <<- "black"

		}else{
			points(x, y, pch = 20, col =col, bg = bg)
		}

		ok = is.finite(x) & is.finite(y) & (length(x) > 1)
		if (any(ok)){
			abline(lm(y~x,subset=ok), col = "blue", lwd=1.5)
			abline(coef=c(0,1),lty=2)
		}
	}

	panel.sse <-
			function(y, x, digits=2,...)
	{

		usr = par("usr"); on.exit(par(usr))
		par(usr = c(0, 1, 0, 1))

		### discard non-fininte values
		ok = is.finite(x) & is.finite(y)
		x <- x[ok]
		y <- y[ok]

		model = summary(lm(y~x))
		r2= model$r.squared

		txt = round(r2, digits)
		txt = bquote(R^2*"" == .(txt))

		text(0.5, 0.5, txt,cex=cexTxt)

	}

	#panel.txt <- function(x, y, labels, cex, font,digits=2, ...){
	panel.txt <- function(x, y, labels, font,digits=2,...){

		txt = diagText[count]
		text(0.5, 0.6, txt, cex=cexTxt, col=textCol[count])

		count <<-count+1
	}

	pairs(data,lower.panel=panel.sse,upper.panel=panel.lm, text.panel=panel.txt,col=col,...)


}

#' Plot Percentage of Features with with missing values
#' @param eset ExpressionSet
#' @param col col
#' @param cex.axis cex.axis
#' @param cex.lab cex.lab
#' @param ... see plot
#' @import Biobase
#' @importFrom graphics mtext
#' @export
#' @note  No note
#' @details No details
#' @references NA
#' @examples print("No examples")
missinValueBarplot <- function(eset, col=as.character(.getConditionColors(eset)[pData(eset)$condition,]), cex.axis=1.25, cex.lab=1.25, ...){

	#par(mar=c(7.1,4.1,4.1,2.1))

	eset <- eset[!fData(eset)$isFiltered,]

	d = (colSums(is.na(exprs(eset))) / nrow(eset)) * 100
	ylim <- c(0,range(d)[2])
	if(max(d,na.rm=T) == 0) ylim <- c(0,100)

	bp <- barplot2( d
			, las=2
			, ylab="% Missing Values"
			, col=col
			, cex.axis = cex.axis
			, cex.lab=cex.lab
			, ylim=ylim
			, xaxt="n"
			, plot.grid=T
			, grid.col="lightgrey"
			,...
	)

	mtext(names(d),side=1,at=bp[,1], las=2, line=0.3,cex=0.9)
	#par(mar=c(5.1,4.1,4.1,2.1))

}

###
#' Barplot of ms-signal per column
#' @param eset expressionSet
#' @param col default condition colors
#' @param method c("median","sum","sharedSignal")
#' @param cex.lab default 1.25
#' @param cex.axis default 1.25
#' @param cex.names default 0.9
#' @param labels labels
#' @param ... see plot
#' @import Biobase
#' @export
#' @importFrom graphics mtext
#' @note  No note
#' @details No details
#' @references NA
#' @examples print("No examples")
barplotMSSignal <- function(eset, col = as.character(.getConditionColors(eset)[pData(eset)$condition,])
		,method=c("sum","sharedSignal")
		,cex.lab=1.25
		,cex.axis=1.25
		,cex.names=0.9
		,labels=rownames(pData(eset))
		,...){

	sel <- !fData(eset)$isFiltered

	### only use feature qunatified in all runs for normalization
	if("sharedSignal" %in% method) sel <- sel & (as.vector(apply(is.finite(exprs(eset)),1,sum) == ncol(eset)))

	eset <- eset[sel,]

	if("median" %in% method){
		profile <- apply(exprs(eset),2,median,na.rm=T)
		ylab <- "Median MS-Signal (Scaled)"
	}else{
		profile <- apply(exprs(eset),2,sum,na.rm=T)
		ylab <- "Summed MS-Signal (Scaled)"
	}
	profile <- profile/max(profile)

	bp <- barplot2(profile,las=2, col=col,ylab=ylab, xaxt="n"
			, plot.grid=T
			, grid.col="lightgrey"
			, cex.lab=cex.lab
			, cex.axis=cex.axis
			,...)
	mtext(labels,side=1,at=bp[,1], las=2, line=0.3,cex=cex.names)

}

#' C.V. boxplot
#' @param eset ExpressionSet
#' @param col col
#' @param cex.names default 0.9
#' @param cex.axis default 1.25
#' @param cex.lab default 1.25
#' @param ylab C.V.
#' @param ... see plot
#' @note  No note
#' @export
#' @importFrom graphics boxplot grid par
#' @details No details
#' @references NA
#' @examples print("No examples")
cvBoxplot <- function(eset,col=as.character(.getConditionColors(eset)[unique(pData(eset)$condition),]), cex.names=0.9,cex.axis=1.25,cex.lab=1.25,ylab="C.V. (%)",...){

	eset <- eset[!fData(eset)$isFiltered,]

	cv <- getAllCV(eset)
	### avoid crach when not enough repliocates
	if(sum(!is.na(cv)) > 0){
		boxplot(cv*100,yaxt="n",xaxt="n",...)

		grid(nx=NA, ny=NULL) #grid over boxplot
		par(new=TRUE)

		boxplot(cv*100,  col=col,las=2, ylab=ylab
				, cex.names=cex.names
				, cex.axis=cex.axis
				, cex.lab=cex.lab
				, textcolor=0
				, ...)

	}
}

#' Plot ms.signal distributions
#' @param d matrix of ms-signals
#' @param col color
#' @param ylab default "Frequnecy"
#' @param xlab default "MS-Signal"
#' @param ... see plot
#' @note  No note
#' @export
#' @importFrom graphics plot lines
#' @importFrom graphics plot hist abline legend grid mtext
#' @details No details
#' @references NA
#' @examples print("No examples")
plotMSSignalDistributions <- function(d, col=1:100,ylab="Frequnecy", xlab="MS-Signal",... ){


	### create a sms-signal histogram per column,
	breaks <- seq(min(d,na.rm=T),max(d,na.rm=T),length=40)
	mids <-  hist(d[,1],breaks=breaks,plot=F)$mids
	countPerCol <- data.frame(row.names=mids)
	ncol(countPerCol)

	for(c in colnames(d)){
		count <- hist(d[,c],breaks=breaks,plot=F)$count
		countPerCol <- cbind(countPerCol,count)
	}
	names(countPerCol) <-  colnames(d)

	#par(mar = c(5, 4, 4.1, 7.5))
	### plot histogram trend lines
	plot(mids,mids,ylim=range(countPerCol)
			, type="n"
			,xlab=xlab
			,ylab=ylab
			,...)
	for(i in 1:ncol(countPerCol)){
		lines(mids,countPerCol[,i],col=col[i],...)
	}
	#.outerLegend("right", names(countPerCol),lty=1 , col=col, bg="white")
	#par(mar = c(5, 4.1, 4.1, 2.1))

	#legend("topleft", names(countPerCol),lty=1 , col=col)

}





###
#' Plot Total Number of diffrentially Abundant Features (applying ratio cutoff) vs. qValue/pValue for all conditions
#' @param sqa SafeQuantAnalysis Object
#' @param upRegulated TRUE/FALSE select for upregulated features
#' @param log2RatioCufOff log2 ratio cut-off
#' @param pvalCutOff pValue/qValue cut-off
#' @param isLegend TRUE/FALSE display legend
#' @param isAdjusted TRUE/FALSE qValues/pValue on x-axis
#' @param ylab default  Nb. Features
#' @param xlim see plot
#' @param ylim see plot
#' @param ... see plot
#' @note  No note
#' @export
#' @importFrom graphics plot lines legend grid
#' @details No details
#' @references NA
#' @examples print("No examples")
plotNbValidDeFeaturesPerFDR <- function(sqa,upRegulated=T,log2RatioCufOff=log2(1)
		,pvalCutOff=1, isLegend=T,isAdjusted=T,ylab="Nb. Features", xlim=NA, ylim=NA, ... ){

	### we need at least two conditions
	if(length(unique(pData(sqa$eset)$condition)) == 1){
		return(cat("INFO: Only one condition no plotNbValidDeFeaturesPerFDR \n"))
	}

	if(isAdjusted){
		pvaluesPerCond <- sqa$qValue
		xlab= "False Discovery Rate (qValue)"
	}else{
		pvaluesPerCond <- sqa$pValue
		xlab= "pValue"
	}

	ratiosPerCond <- sqa$ratio
	conditionColors <- .getConditionColors(sqa$eset)

	if(is.na(xlim[1])){
		xlim=c(0,0.3)
	}

	pvalCutOffs <- seq(xlim[1],xlim[2], length.out=10)
	conditions <- names(pvaluesPerCond)

	### create data farme of roc curve per condition
	### store curves
	nbPassingCutOffsPerCond <- data.frame(row.names=pvalCutOffs)

	### iterate over all conditions
	for(cond in conditions){

		pvals <- pvaluesPerCond[cond]
		ratios <- ratiosPerCond[cond]

		#iterate over all cutoffs
		nbPassingCutOffs <- c()
		for(qCutOff in pvalCutOffs){

			if(upRegulated){
				nbPassingCutOffs <- c(nbPassingCutOffs, sum( (pvals < qCutOff) & (ratios >  log2RatioCufOff) ,na.rm=T) )
			}else{
				nbPassingCutOffs <- c(nbPassingCutOffs, sum( (pvals < qCutOff) & (ratios <  -log2RatioCufOff) ,na.rm=T ) )
			}
		}

		nbPassingCutOffsPerCond <- data.frame(nbPassingCutOffsPerCond, nbPassingCutOffs)

	}

	names(nbPassingCutOffsPerCond) <- conditions


	if(is.na(ylim[1])){
		ylim=c(0,max(nbPassingCutOffsPerCond))
	}

	# plot roc curves
	plot(0,0, ylim=ylim, xlim= xlim, type="n",ylab=ylab,xlab=xlab,  ...)
	grid()
	for(cond in conditions){
		lines(pvalCutOffs,nbPassingCutOffsPerCond[,cond], col= as.character(conditionColors[cond ,]), lwd=2)
	}
	abline(v=pvalCutOff,col="grey",lwd=1.5)
	if(isLegend){
		legend("topleft", conditions, fill=as.character(conditionColors[conditions,]), cex=0.6)
	}


}

#' Plot identifications target decoy distribution
#' @param targetScores target Scores
#' @param decoyScores decoy Scores
#' @param xlab default "Identification Score"
#' @param ylab default "Counts"
#' @param ... see plot
#' @note  No note
#' @export
#' @importFrom graphics plot hist points legend grid
#' @details No details
#' @references NA
#' @examples print("No examples")
plotScoreDistrib <-function(targetScores,decoyScores,xlab="Identification Score",ylab="Counts", ...){

	if(length(targetScores) >0 & length(decoyScores) >0){

		breaks <- hist(c(targetScores,decoyScores),plot=FALSE,breaks=100)$breaks

		targetHist = hist(targetScores, breaks=breaks, plot=FALSE)
		decoyHist = hist(decoyScores, breaks=breaks, plot=FALSE)

		ylim = c(0,max(c(decoyHist$counts,targetHist$counts),na.rm=T))

		plot(0,0,type='n', ylim=ylim, xlab=xlab, ylab=ylab, xlim=range(targetScores,decoyScores,na.rm=T),...)
		grid()
		points(targetHist$mids[targetHist$counts > 0], targetHist$counts[targetHist$counts > 0], col=1, type="h", lwd=4)
		points(decoyHist$mids[decoyHist$counts > 0], decoyHist$counts[decoyHist$counts > 0], col=2, type="h", lwd=5)

	}else{

		if(length(targetScores) >0){
			targetHist = hist(targetScores, breaks=100, plot=FALSE)
			plot(targetHist$mids, targetHist$counts, xlab=xlab, ylab=ylab,type="n",... )
			points(targetHist$mids[targetHist$counts > 0], targetHist$counts[targetHist$counts > 0], col=1, type="h", lwd=4)
		}

		cat("plotScoreDistrib: Not enough target or decoy scores\n")
	}

	legend("topright",c("target","decoy"),fill=c(1,2))
}



#' Plot FDR vs. identification score
#' @param idScore vector of identification scores
#' @param qvals vector of q-valres
#' @param qvalueThrs threshold indicated by horizontal line
#' @param xlab default Identification Score
#' @param ylab default False Discovery Rate
#' @param ... see plot
#' @export
#' @importFrom graphics plot grid abline
#' @note  No note
#' @details No details
#' @references NA
#' @examples print("No examples")
plotIdScoreVsFDR <-function(idScore,qvals,qvalueThrs=0.01, ylab="False Discovery Rate", xlab="Identification Score",...){
	plot(sort(idScore),rev(sort(qvals)),type="l",ylab=ylab,xlab=xlab, ... )
	grid()
	abline(h=qvalueThrs,col="grey")
}



### plot identification scores ROC-curve, fdr vs. # valid identifications
#' Plot Number of Identifications vs. FDR
#' @param qvals vector of q-values
#' @param qvalueThrs threshold indicated by vertical line
#' @param breaks see breaks for hist function
#' @param xlab default "False Discovery Rate"
#' @param ylab default "Nb. Valid Identifications"
#' @param xlim default c(0,0.1)
#' @param col default blue
#' @param lwd default 1.5
#' @param ... see plot
#' @export
#' @importFrom graphics plot abline grid
#' @note  No note
#' @details No details
#' @references NA
#' @examples print("No examples")
plotROC <- function(qvals
		,qvalueThrs=0.01
		,xlab="False Discovery Rate"
		,ylab="Nb. Valid Identifications"
		,xlim=c(0,0.1)
		,breaks=100
		,col="blue"
		,lwd=1.5
		,... ){

	if(length(breaks) == 1 ){
		breaks = seq(xlim[1],xlim[2],length=breaks)
	}else{  ### breaks override xlim
		xlim <- c(min(breaks),max(breaks))
	}

	validIds = c()
	for(fdr in breaks){
		validIds = c(validIds,sum(qvals < fdr, na.rm=T))
	}

	plot(breaks,sort(validIds), ylab=ylab,xlab=xlab, xlim=xlim,type="l",col=col,lwd=lwd,...)
	grid()
	abline(v=qvalueThrs,lwd=lwd,col="grey")
}



### Plot Precursor Mass Error Distribution
#' Plot Precursor Mass Error Distribution
#' @param eset ExpressionSet
#' @param pMassTolWindow Precursor Mass Error Tolerance Window
#' @param ... see plot
#' @export
#' @importFrom graphics plot hist abline legend grid mtext axis
#' @note  No note
#' @details No details
#' @references NA
#' @examples print("No examples")
plotPrecMassErrorDistrib <- function(eset,pMassTolWindow=c(-10,10), ...){

	# do nothing if all NA
	if(sum(is.finite(fData(eset)$pMassError)) == 0 ){return(NA)}

#	par(mar = c(5, 4, 4.1, 3.5))
#	plot(rnorm(50), rnorm(50), col=c("steelblue", "indianred"), pch=20)
#	outerLegend("right", legend=c("Foo", "Bar"), pch=20,
#			col=c("steelblue", "indianred"),
#			horiz=F, bty='n', cex=0.8)

	#par(mar=c(5.1,4.1,4.1,4.1))
	isDec <- isDecoy(fData(eset)$proteinName)
	pMassError <- fData(eset)$pMassError

	xRange 	<- range(pMassError,pMassTolWindow,na.rm=T)
	breaks <- seq(xRange[1],xRange[2],length=50)

	### decoy hist
	decoyMdHist 	<- hist(pMassError[isDec],breaks=breaks,plot=F)

	### target hist
	targetMdHist 	<- hist(pMassError[!isDec],breaks=breaks,plot=F)

	yRange <- range(decoyMdHist$counts,targetMdHist$counts)

	### plot histogram contours
	plot(decoyMdHist$mids,decoyMdHist$counts,type="l", col="red",ylim=yRange,xlab="mass diff. (ppm)",ylab="PSM Frequnecy",lwd=1.5, ...)
	grid()
	lines(decoyMdHist$mids,targetMdHist$counts,lwd=1.5)
	abline(v=pMassTolWindow[1],col="grey",lwd=2)
	abline(v=pMassTolWindow[2],col="grey",lwd=2)

	ratio <- (decoyMdHist$counts+1)/(targetMdHist$counts+1)
	ratio[ratio > 1] <- 1
	ratio <- ratio*yRange[2]

	lines(decoyMdHist$mids,ratio,lwd=1.5,col="blue",lty=2)

	### add ratio ratio line
	#axis(4,at=seq(0,yRange[2],length=5), labels=seq(0,1,length=5),col="blue", col.ticks="blue")
	axis(4,col="blue", col.ticks="blue", labels=NA, at=seq(0,yRange[2],length=5))
	mtext(seq(0,1,length=5),side=4,at=seq(0,yRange[2],length=5),col="blue", line=0.5)
	mtext("decoy-target PSM count ratio ",side=4, col="blue",line=1.5)

	legend("left",c("target","decoy"),  fill= c("black","red"))

	#par(mar=c(5.1,4.1,4.1,2.1))

}

###
#' Plot precursorMass error v.s score highlighting decoy and displaying user specified user specified precursor mass filter
#' @param eset ExpressionSet
#' @param pMassTolWindow Precursor Mass Error Tolerance Window
#' @param ... see plot
#' @export
#' @importFrom graphics plot points abline legend grid
#' @note  No note
#' @details No details
#' @references NA
#' @examples print("No examples")
plotPrecMassErrorVsScore <- function(eset, pMassTolWindow=c(-10,10) ,...){

	# do nothing if all NA
	if(sum(is.finite(fData(eset)$pMassError)) == 0 ){return(NA)}

	pMassError <- fData(eset)$pMassError
	idScore <- fData(eset)$idScore
	isDec <- isDecoy(fData(eset)$proteinName)

	withinTol <- (pMassError >= pMassTolWindow[1]) & (pMassError <= pMassTolWindow[2])

	plot(pMassError, idScore
			, type="n"
			#, col = isDecoy+1
			,ylab="score", xlab="mass diff. (ppm)", xlim=range(c(pMassTolWindow,pMassError),na.rm=T)
			,...
	)

	grid()

	points(pMassError[withinTol], idScore[withinTol], col="black", pch=20 )
	points(pMassError[!withinTol], idScore[!withinTol], col="grey", pch=20 )
	points(pMassError[isDec], idScore[isDec], col="red", pch=20 )

	abline(v=pMassTolWindow[1],col="lightgrey",lwd=1)
	abline(v=pMassTolWindow[2],col="lightgrey",lwd=1)

	legend("topright",c("target-kept","decoy","target-discarded"), pch=20, col= c("black","red","grey"))

}

#' Scatter plot with density coloring
#' @param x number vector
#' @param y number vector
#' @param isFitLm fit linear model
#' @param disp  c("abline","R","Rc") display options
#' @param legendPos see legend
#' @param cexLegend = legend txt cex
#' @param pch see plot
#' @param ... see plot
#' @import epiR
#' @importFrom graphics plot lines abline legend
#' @importFrom grDevices col2rgb colorRampPalette colors densCols rgb
#' @importFrom stats lm
#' @note  No note
#' @export
#' @references NA
#' @examples print("No examples")
plotXYDensity <- function(x,y
                          ,isFitLm=T
                          ,legendPos="bottomright"
                          ,cexLegend=2
                          ,disp=c("abline","R","Rc"),pch=20,  ...){

	df <- data.frame(x,y)

	## Use densCols() output to get density at each point
	densityCol <- densCols(x,y, colramp=colorRampPalette(c("black", "white")))
	df$densityCol <- col2rgb(densityCol)[1,] + 1L

	## Map densities to colors
	cols <-  colorRampPalette(c("#000099", "#00FEFF", "#45FE4F",
					"#FCFF00", "#FF9400", "#FF3100"))(256)
	df$col <- cols[df$densityCol]

	## Plot it, reordering rows so that densest points are plotted on top
	plot(y~x, data=df[order(df$densityCol),], pch=pch, col=col, ...)

	### disp linerar model
	if(isFitLm){
		ok <- is.finite(x) & is.finite(y)
		fit <- lm(y[ok] ~ x[ok])

		if("abline" %in% disp){
			abline(coef=c(0,1),lty=2)
			abline(fit)
			#library(MASS)
			#abline(rlm(y[ok] ~ x[ok]))
		}

		if("lowess" %in% disp){
			lines(lowess(x[ok],y[ok],...))
		}
		legd <- c()

		if("R" %in% disp) legd <- c(legd,as.expression(bquote(R^2*"" == .(round(summary(fit)$r.squared,2)))))

		if("Rc" %in% disp) legd <- c(legd,as.expression(bquote(R[c]*"" == .(round(epi.ccc(x,y)$rho.c$est,2)))))

		if(length(legd) > 0 ){
			legend(legendPos
					,legend= legd
					,text.col=1, box.col="transparent", cex=cexLegend, bty="n"
			)
		}
		return(fit)
	}
	return(NA)
}

#' Plot absolut Estimation calibration Curve
#' @param fit simple log-linear model
#' @param dispElements c("formula","lowess","stats")
#' @param cex.lab expansion factor for axis labels
#' @param cex.axis expansion factor for axis
#' @param cex.text expansion factor for legend
#' @param cex.dot expansion factor for plotted dots
#' @param predictorName predictorName
#' @param text add names beside each dot
#' @param xlab xlab
#' @param ylab ylab
#' @param main main
#' @param ... see plot
#' @note  No note
#' @export
#' @importFrom graphics par plot lines mtext abline legend
#' @importFrom stats coef predict median
#' @references NA
#' @examples print("No examples")
plotAbsEstCalibrationCurve <- function(fit
		,dispElements = c("formula","lowess","stats")
		,xlab="Conc. (CPC) "
		,ylab="Pred. Conc. (CPC) "
		,predictorName = paste("log10(",names(coef(fit))[2],")",sep="")
		,text=F
		,cex.lab=1
		,cex.axis=1
		,cex.text=1
		,cex.dot=1
		,main = ""
		,...){
	x <- predict(fit) +  fit$residuals
	y <- predict(fit)


	### some extra margin for axis labels
	par(mar=c(5.5,5.5,4.1,2.1))
	plot(10^x, 10^y,log="xy",xlab="",ylab="",main=main,cex.axis=cex.axis,cex=cex.dot,... )

	if(text){
		text(10^x, 10^y,rownames(fit$qr$qr))
	}


	abline(coef=c(0,1),lty=2)
	### add axis labels
	mtext(side=1,xlab,las=1, line=4,cex=cex.lab, ...)
	mtext(side=2,ylab,las=3, line=4,cex=cex.lab, ...)

	if( "formula" %in% dispElements){
		legend("bottomright"
				,paste("log10(",ylab,")"," = ", signif(coef(fit)[1],2)," + ",signif(coef(fit)[2],2)," * ",predictorName ,sep="")
				,box.lwd=0
				,box.col="white"
				,cex=cex.text
				,...
		)
	}

	if( "lowess" %in% dispElements){
		lws <- lowess(y ~ x)
		lines(10^lws$x, 10^lws$y,col="red",...)
	}

	if( "stats" %in% dispElements){

		df <- data.frame(cpc =  x,signal = y)
		medianFoldError <- median(abs(getLoocvFoldError(df)[,1]),na.rm=T)

		legend("topleft"
				,legend=c(as.expression(bquote(R^2*"" == .(round(summary(fit)$r.squared,2))))
						,paste("Median Fold Error = ",round(medianFoldError,2))
				)
				#,text.col=c(1,2)
				,box.lwd=0
				,box.col="white"
				,cex=cex.text
				,...
		)
	}
	### reset margins
	par(mar=c(5.1,4.1,4.1,2.1))
}


#' Plot all retention time normalization profiles
#' @param eset ExpressionSet
#' @param ... see plot
#' @param col condition colors
#' @note  No note
#' @export
#' @importFrom graphics plot lines abline
#' @details No details
#' @seealso  \code{\link{getRTNormFactors}}
#' @references In Silico Instrumental Response Correction Improves Precision of Label-free Proteomics and Accuracy of Proteomics-based Predictive Models, Lyutvinskiy et al. (2013), \url{http://www.ncbi.nlm.nih.gov/pubmed/23589346}
#' @examples print("No examples")
plotRTNormSummary <- function(eset, col = as.character(.getConditionColors(eset)[pData(eset)$condition,1]),...){

	rtNormFactors <- getRTNormFactors(eset, minFeaturesPerBin=100)
	ylim <- range(rtNormFactors,na.rm=T)
	runNames <- names(rtNormFactors)

	parDef <- par()
	par(mar = c(5, 4, 4.1, 7.5))
	plot(rownames(rtNormFactors),rtNormFactors[,1], type="n", ylab="log2(Ratio)",xlab="Retention Time (min)",ylim=ylim, ...)
	for(i in 1:ncol(rtNormFactors)){
		lines(rownames(rtNormFactors),rtNormFactors[,i], col=col[i], ...   )
	}
	abline(h=0, lty=2)
	.outerLegend("right", runNames, lty=1, col=col, lwd=2, bg="white")
	#.outerLegend("right", runNames, lty=1, col=col, lwd=1.5)
	par(parDef)
}

#' Plot all retention time profile overalying ratios
#' @param rtNormFactors data.frame of normalization factor per r.t bin and sample, obtained by getRTNormFactors
#' @param eset  ExprsssionSet
#' @param samples specify samples (sample numbers) to be plotted
#' @param main main
#' @param ... see plot see plot
#' @note  No note
#' @export
#' @importFrom graphics plot lines abline
#' @details No details
#' @seealso  \code{\link{getRTNormFactors}}
#' @references In Silico Instrumental Response Correction Improves Precision of Label-free Proteomics and Accuracy of Proteomics-based Predictive Models, Lyutvinskiy et al. (2013), \url{http://www.ncbi.nlm.nih.gov/pubmed/23589346}
#' @examples print("No examples")
plotRTNorm <- function(rtNormFactors,eset,samples=1:ncol(rtNormFactors),main="", ...){

	### select for anchor proteins
	eset <- eset[fData(eset)$isNormAnchor,]

	# get all ratios to sample 1
	# @TODO How to select reference run?
	#ratios <- log2(exprs(eset)) - log2(exprs(eset)[,1])
	ratios <- log2(exprs(eset)) - apply(log2(exprs(eset)),1,mean,na.rm=T)

	for(samplesNb in samples){
		plot(fData(eset)$retentionTime,	ratios[,samplesNb]
				#, col="lightgrey"
				, col=rgb(0,100,0,50,maxColorValue=255)
				, xlab="Retention Time (min)"
				, ylab="log2(Ratio)"
				, cex.lab=1.5
				, cex.axis=1.5
				, main=paste(main,names(rtNormFactors)[samplesNb]), ...)
		lines(as.numeric(rownames(rtNormFactors)),as.vector(unlist(rtNormFactors[,samplesNb])),lwd=2, ... )
		abline(h=0,lty=2,...)
	}
}


#' Plot the number of identified Features per Reteintion Time minute.
#' @param eset ExpressionSet
#' @param ... see plot see plot
#' @param cex.axis default 1.25
#' @param cex.lab default 1.25
#' @param col default "blue"
#' @param lwd default 2
#' @note  No note
#' @export
#' @importFrom graphics plot
#' @references NA
#' @examples print("No examples")
plotNbIdentificationsVsRT <- function(eset, cex.axis=1.25,cex.lab=1.25, col="blue", lwd=2, ...){
	rtTable <- table(round(fData(eset)$retentionTime))
	plot(as.numeric(names(rtTable)),rtTable[names(rtTable)], type="h", xlab="Retention Time (min)"
			, ylab="# Identified Features", cex.axis=cex.axis,cex.lab=cex.lab, col=col, lwd=lwd,...)
}


#' Plot qValue vs pValue
#' @param sqa SafeQuantAnalysis Object
#' @param lim x-axis and y-axis range
#' @param ... see plot
#' @note  No note
#' @export
#' @importFrom graphics plot lines abline legend
#' @details No details
#' @references NA
#' @examples print("No examples")
plotQValueVsPValue <- function(sqa, lim=c(0,1), ...){

	### we need at least two conditions
	if(length(unique(pData(sqa$eset)$condition)) == 1){
		return(cat("INFO: Only one condition no plotQValueVsPValue \n"))
	}

	conditionColors <- .getConditionColors(sqa$eset)
	conditions <- names(sqa$pValue)


	plot(0,0,xlab="pValue",ylab="False Discovery Rate (qValue)", type="n", xlim=lim,ylim=lim)
	grid()
	for(cond in conditions){

		pVal <- sqa$pValue[,cond]
		qVal <- sqa$qValue[,cond]

		# discard NA's
		pVal <- pVal[!is.na(qVal)]
		qVal <- qVal[!is.na(qVal)]

		o <- order(pVal)
		lines(pVal[o],qVal[o],col= as.character(conditionColors[cond ,]), lwd=2)
	}

	abline(coef=c(0,1),lty=2)
	legend("bottomright", conditions, fill=as.character(conditionColors[conditions,]), cex=0.6)
}


#' Plot adjusted tmt ratios vs original ratios
#' @param ratio data.frame
#' @param unAdjustedRatio data.frame
#' @note  No note
#' @details plot adjusted tmt ratios vs original ratios
#' @references NA
#' @examples print("No examples")
#' @export
plotAdjustedVsNonAdjustedRatio <- function(ratio,unAdjustedRatio){

	for(i in 1:ncol(ratio)){
		lim <- range(ratio)
		lim <- c(-max(lim),max(lim))
		plot(unAdjustedRatio[,i],ratio[,i]
				, xlim=lim
				, ylim=lim
				, cex.axis=1.25
				, cex.lab=1.5
				, pch=19
				, col="grey"
				, xlab="TMT Ratio"
				, ylab="Adjusted TMT Ratio"
		)
		abline(coef=c(0,1), lty=2,lwd=1.5)
		fit <- lm(ratio[,i] ~ unAdjustedRatio[,i] )
		abline(fit, lwd=1.5)

		# R2
		legd <- c(as.expression(bquote(R^2*"" == .(round(summary(fit)$r.squared,2)))))
		# slope
		legd <- c(legd,as.expression(bquote(K*"" == .(round(coef(fit)[2],2)))))

		legend("bottomright"
				,legend= legd
				,text.col=1, box.col="transparent", cex=1.5
		)
	}
}




### @TODO add unit test
#' @export
#' @importFrom graphics par plot lines mtext abline legend
.plotCalibrationCurve <- function(fit
		,dispElements = c("formula","lowess","stats")
		,xlab="Protein Copies/Cell Measured using SID"
		,ylab="Protein Copies/Cell Estimated using iBAQ"
		,cex=1.5
		,...){
	x <- predict(fit) + fit$residuals
	y <- predict(fit)


	### some extra margin for axis labels
	par(mar=c(6.3,6.3,4.1,2.1))
	plot(10^x, 10^y
			,log="xy"
			,xlab=""
			,ylab=""
#		,yaxt="n"
#		,xaxt="n"
			,cex=cex
			#,lwd=cex
			,pch=19
			,las=2,cex.axis=cex
			,cex.main=cex
			,... )

#	axis(1,las=2,cex.axis=cex)
#	axis(2,las=2,cex.axis=cex)

	abline(coef=c(0,1),lty=2)
	### add axis labels
	mtext(side=1,xlab,las=1, line=4.8, cex=cex, ...)
	mtext(side=2,ylab,las=3, line=4.8, cex=cex, ...)

	if( "formula" %in% dispElements){
		legend("bottomright"
				,paste("log10(Est. CPC)"," = ", signif(coef(fit)[1],2)," + ",signif(coef(fit)[2],2)," * log10(iBAQ)" ,sep="")
				,box.lwd=0
				,box.col="white"
				,cex=cex-0.4
				,...
		)
	}

	if( "lowess" %in% dispElements){
		lws <- lowess(y ~ x)
		lines(10^lws$x, 10^lws$y,col="red",...)
	}

	if( "stats" %in% dispElements){

		df <- data.frame(cpc =  x,signal = y)
		medianFoldError <- median(abs(getLoocvFoldError(df)[,1]),na.rm=T)
		linRc <- as.vector(unlist(epi.ccc(x,y)$rho.c)[1])

		legend("topleft"
				,legend=c(as.expression(bquote(R^2*"" == .(round(summary(fit)$r.squared,2))))
						#,paste("Lin's Rc  = ",round(linRc,2))
						,as.expression(bquote(R[c]*"" == .(round(linRc,2))))
						,paste("Median Fold Error = ",round(medianFoldError,2))
				)
				#,text.col=c(1,2)
				,box.lwd=0
				,box.col="white"
				,cex=cex
				,...
		)
	}
	### reset margins
	par(mar=c(5.1,4.1,4.1,2.1))
}

## require pairide calibration mix eset
## @export
#.plotTMTRatioVsRefRatio <- function(esetCalibMixPair,ylim=c(-2.8,2.8), xlim=c(-2.8,2.8),cex.lab=2, cex.axis=2, ...){
#
#	calMixDilutionTag <- round(log10(fData(esetCalibMixPair)$calMixDilution)+0.6)
#	plot(0,0
#			,xlab="TMT Ratio (log2)"
#			,ylab="Reference Ratio (log2)"
#			,xlim=xlim
#			,ylim=ylim
#			,type="n"
#			,cex.lab=cex.lab
##			,xaxt="n"
##			,yaxt="n"
#			,cex.axis=cex.axis
#			,...
#	)
#
#	# add diagonal ref line
#	abline(coef=c(0,1), h=0, v=0,lty=2, col="grey")
#
#	points(	getRatios(esetCalibMixPair)[,1]
#			,log2(fData(esetCalibMixPair)$refRatio)+(calMixDilutionTag/10 - 0.2)
#			,pch=19
#			,col=calMixDilutionTag+1)
#
##	text(	getRatios(esetCalibMixPair)[,1]
##			,log2(fData(esetCalibMixPair)$refRatio)+(calMixDilutionTag/10 - 0.2)
##			,gsub("^[sptr]{2}\\|","",fData(esetCalibMixPair)$proteinName)
##			,col=calMixDilutionTag+1
##			,pos=1:4
##			,cex=0.6
##		)
#
#
#	fit <- getRatioCorrectionFactorModel(esetCalibMixPair)
#	abline(fit, lwd=1.5)
#	legend("topleft",c(expression(paste("4 (fmol/", mu,"g)")),"20","100"),title="Cal. Mix Concentration",fill=2:4,cex=1.5,box.lwd=F, box.col=0)
#
#}

## Plot (boxplot) spectrum ratio distributions per protein and Claibration Mix Concentration
## @param eset
## @export
## @note  No note
## @details No details
## @references NA
## @examples print("No examples")
#plotCalibrationMixRatios <- function(eset, cex.axis=1, cex.lab=1.1,...){
## assuming tenplex mix
#	for(targetProtein in names(CALIBMIXRATIOS)){
#
#		c1 <- exprs(eset)[(fData(eset)$proteinName %in% targetProtein),1]
#		c2 <- exprs(eset)[(fData(eset)$proteinName %in% targetProtein),2]
#		c3 <- exprs(eset)[(fData(eset)$proteinName %in% targetProtein),3]
#		c4 <- exprs(eset)[(fData(eset)$proteinName %in% targetProtein),4]
#		c5 <- exprs(eset)[(fData(eset)$proteinName %in% targetProtein),5]
#		c6 <- exprs(eset)[(fData(eset)$proteinName %in% targetProtein),6]
#		c7 <- exprs(eset)[(fData(eset)$proteinName %in% targetProtein),7]
#		c8 <- exprs(eset)[(fData(eset)$proteinName %in% targetProtein),8]
#		c9 <- exprs(eset)[(fData(eset)$proteinName %in% targetProtein),9]
#		c10 <- exprs(eset)[(fData(eset)$proteinName %in% targetProtein),10]
#
#		r1 <- c2 / c1
#		r2 <- c4 / c3
#		r3 <- c6 / c5
#		r4 <- c8 / c7
#		r5 <- c10 / c9
#
#		nbDp <- length(r1)
#		xlab <- expression(paste("Calibration Mix Conc. ","(fmol/",mu,"g)",sep=""))
#
#		ylab <- "Fold Change"
#		names <- c(4,4,20,20,100)
#		ylim <- c(0.2,5)
#
#		if(nbDp > 9){
#			# per channel pair
#			boxplot(r3,r5,r1,r4,r2, ylim=ylim, main=paste(targetProtein,"\nN=",length(c1),sep="")
#					, col=c(2,2,3,3,4),ylab=ylab,names=names,xlab=xlab,cex.axis=cex.axis,cex.lab=cex.lab,log="y", ...)
#			# per condiiton pair
#			#boxplot(log2(r3), log2(c(r1,r4)),log2(c(r2,r5)), ylim=c(-2,2), main=main, col=2:4,ylab="log2(RATIO)",names=c(4,20,100),xlab="Mix Conc.")
#		}else{ ### fewer than X data points, the plot points
#
#
#			##Plot first set of data.
#			##Need to check for sensible ranges
#			##Use the jitter function to spread data out.
#			plot(jitter(rep(0,nbDp),amount=0.2),r3 ,
#					xlim=range(-0.5,4.5)
#					,ylim=ylim
#					,xaxt="n"
#					,frame.plot=TRUE
#					,col=2
#					,main=targetProtein
#					,xlab=xlab
#					,ylab=ylab
#					,log="y"
#					,cex.lab=cex.lab
#					,cex.axis=cex.axis
#					, pch=19
#					,...
#			)
#			points(jitter(rep(1,nbDp), amount=0.2), r5, col=2, pch=19)
#			points(jitter(rep(2,nbDp), amount=0.2), r1, col=3, pch=19)
#			points(jitter(rep(3,nbDp), amount=0.2), r4, col=3, pch=19)
#			points(jitter(rep(4,nbDp), amount=0.2), r2, col=4, pch=19)
#
#			##Add in the x-axis
#			axis(1, at=0:4,labels=names	,cex.axis=cex.axis)
#
#			medianR1 <- median(r1,na.rm=T)
#			medianR2 <- median(r2,na.rm=T)
#			medianR3 <- median(r3,na.rm=T)
#			medianR4 <- median(r4,na.rm=T)
#			medianR5 <- median(r5,na.rm=T)
#
#			##Add in the medians
#			segments(-0.25, medianR3, 0.25, medianR3,lwd=2)
#			segments(0.75, medianR5, 1.25, medianR5,lwd=2)
#			segments(1.75, medianR1, 2.25, medianR1,lwd=2)
#			segments(2.75, medianR4, 3.25, medianR4,lwd=2)
#			segments(3.75, medianR2, 4.25, medianR2,lwd=2)
#
#		}
#
#		refRatio <- fData(eset)[fData(eset)$proteinName %in% targetProtein,]$refRatio[1]
#		abline(h=refRatio,lwd=1.5)
#		abline(h=1,lty=2,lwd=1.5)
#
#	}
#}


#' @title ma-plot
#' @description ma-plot
#' @note  No note
#' @references NA
#' @examples print("No examples")
#' @export
#' @param eset ExpressionSet
#' @param sample selected condition
#' @param lwd default 2
#' @param pch default 1
#' @param cex.lab default 1.5
#' @param cex.axis default 1.5
#' @param col green transparent
#' @param ... see plot
maPlotSQ <- function(eset
		,sample=colnames(exprs(eset))[1]
		,cex.lab=1.5
		,cex.axis=1.5
		,lwd=2
		,pch=1
		,col= rgb(0,100,0,50,maxColorValue=255)
		,...){


	A <- apply(log2(exprs(eset)),1,mean)
	M <- log2(exprs(eset)[,sample]) - A

	plot(M ~ A,
			,pch=pch
			, xlab="Mean Log2 Intensity"
			#, ylab= paste("M - " ,sample, sep="" )
			, ylab= "log2 Ratio"
			, main=sample
			, col = col
			,cex.lab=cex.lab
			,cex.axis=cex.axis
			,...)
#abline(h=c(-1,0,1),lwd=lwd,lty=2)
	abline(h=c(0),lwd=1,lty=2)
	lines(lowess(M ~A),lwd=lwd,col="black")


}

# import motifStack
#' Plot sequence logo
#' @param motif list of target residue centered motfis
#' @param bgPeptides peptides used to calculate residue background frequency (default uniform )
#' @param main see plot
#' @param targetResidues default [STY]
#' @param ic.scale   logical. If TRUE, the height of each column is proportional to its information content. Otherwise, all columns have the same height.
#' @note  No note
#' @references NA
#' @examples print("No examples")
#' @export
plotLogo <- function(motif, bgPeptides="ACDEFGHIKLMNPQRSTVWY", main="", targetResidues=c("S","T","Y"), ic.scale=F, ...){

	# discard short motifs
	nc = nchar(motif)
	motif = motif[nc == max(nc,na.rm=T)]

	motifMatrix = do.call(rbind, strsplit(motif,""))
	motifMatrix = motifMatrix[motifMatrix[,round((nchar(motif[1])+1)/2)] %in% targetResidues,]

	aas = c("A","C","D","E","F","G","H","I","K","L","M","N","P","Q","R","S","T","V","W","Y")
	aaFreq =  data.frame()
	for(pos in 1:ncol(motifMatrix)){
		aaFreq = rbind(aaFreq,table(motifMatrix[,pos])[aas])
	}
	names(aaFreq) = aas
	aaFreq[is.na(aaFreq)] = 0
	aaFreq = t(aaFreq)

	bgFreq = table(as.vector(unlist(strsplit(as.character(bgPeptides),""))))
	bgFreq = bgFreq / sum(bgFreq)

	pfm=pcm2pfm(aaFreq)
	pfm=new("pfm", mat=pfm, name=main,
			color=colorset(alphabet="AA",colorScheme="chemistry")
	)
	slot(pfm,"background")[aas] = as.numeric(bgFreq)

	plot(pfm, ic.scale=ic.scale,...)

}

#' Plot Cystein Frequency
#' @param peptides vector
#' @param ... see plot
#' @export
#' @note  No note
#' @details Selecting for peptides of length 7-19
#' @references NA
#' @examples print("No examples")
cysteinFreqBarplot = function(peptides, ...){

	peptides = peptides %>% unique %>% as.character
	peptidesLength = nchar(peptides)
	peptides = subset(peptides ,peptidesLength > 6 & peptidesLength < 20  )
	cFreq = str_detect(peptides,"C") %>% sum / length(peptides)

	barplot(c(cFreq,0.1148966,0.1326964,0.2126175)*100, col =c("blue",rep("lightgrey",3)), names=c("this sample","e.coli","yeast","human"), ylab ="Cystein containig peptides (%)", las=2, ... )
}


#' Hierarchical clustering heat map, cluster by runs intensity, features by ratio and display log2 ratios to control median DEPRECATED
#' @param eset ExpressionSet
#' @param conditionColors data.frame of colors per condition
#' @param breaks default seq(-2,2,length=20)
#' @param dendogram see heatmap.2 gplots
#' @param legendPos see legend
#' @param ... see plot
#' @import gplots
#' @importFrom graphics legend
#' @importFrom stats dist hclust as.dist as.dendrogram
#' @export
#' @note  No note
#' @details No details
#' @references NA
#' @examples print("No examples")
hClustHeatMapOld <- function(eset
                          ,conditionColors =.getConditionColors(eset)
                          ,breaks=seq(-2,2,length=20)
                          ,dendogram = "column"
                          ,legendPos="left"
                          ,...
){

  # need more than own data feature
  if(nrow(eset) <= 1 ) return()

  # do not plot filtered
  eset <- eset[!fData(eset)$isFiltered,]

  ### we need at least two conditions
  if(ncol(eset) == 1){
    return(cat("INFO: Only one condition no hClustHeatMap \n"))
  }

  #d <- log2(exprs(eset))
  ### log2 ratios to median of control condition
  # @TODO adapt to accomaodated paired expDesign
  log2RatioPerMsRun <- log2(exprs(eset)) - log2(getSignalPerCondition(eset,method="median")[,.getControlCondition(eset)])

  #	feature.cor = cor(t(log2RatioPerMsRun), use="pairwise.complete.obs", method="pearson")
  #	feature.cor.dist = as.dist(1-feature.cor)
  #	feature.cor.dist[is.na(feature.cor.dist)] <- 0
  #	feature.tree = hclust(feature.cor.dist, method="ward.D2")

  feature.tree = hclust(dist(log2RatioPerMsRun,method = "euclidean"), method="ward.D2")

  ### !!!! clustering runs based on ratios is not a good idea as ratio corrrealtion of control runs is not expected to be high
  msrun.cor.pearson = cor(log2(exprs(eset)), use="pairwise.complete.obs", method="pearson")
  msrun.cor.pearson.dist = as.dist(1-msrun.cor.pearson)
  ### to avoid error when replicates of the same condition are identical, DOES THIS EVER HAPPEN?
  msrun.cor.pearson.dist[is.na(msrun.cor.pearson.dist)] <- 0
  msrun.tree = hclust(msrun.cor.pearson.dist, method="ward.D2")

  ### sample colors
  samplecolors =  as.vector(unlist(conditionColors[pData(eset)$condition,]))
  labRow <- rownames(log2RatioPerMsRun)

  ### do not display feature names if too many
  if(nrow(log2RatioPerMsRun) > 300){
    labRow <- rep("",(nrow(log2RatioPerMsRun)))
  }

  hm <- heatmap.2(as.matrix(log2RatioPerMsRun)
                  , col=colorRampPalette(c(colors()[142],"black",colors()[128]))
                  , scale="none"
                  , key=TRUE
                  , symkey=FALSE
                  , trace="none"
                  , cexRow=0.5
                  , cexCol = 0.7
                  ,ColSideColors=samplecolors
                  ,labRow = labRow
                  ,Rowv=as.dendrogram(feature.tree)
                  ,Colv=as.dendrogram(msrun.tree)
                  ,dendrogram=dendogram
                  ,density.info="density"
                  #,KeyValueName="Prob. Response"
                  ,breaks=breaks
                  ,na.rm=T
                  ,key.title=""
                  ,key.xlab = "log2 Ratio"
                  ,key.ylab=""
                  , ...
  )

  legend(legendPos,levels(pData(eset)$condition), fill=as.character(conditionColors[,1]), cex=0.7, box.col=0)

}

#' Hierarchical clustering heat map, cluster by pearson correlation of run intensity, features by mean centered intenstiy
#' @param eset ExpressionSet
#' @param conditionColors data.frame of colors per condition
#' @param breaks default seq(-2,2,length=20)
#' @param dendogram see heatmap.2 gplots
#' @param legendPos see legend
#' @param ... see plot
#' @import gplots
#' @importFrom graphics legend
#' @importFrom stats dist hclust as.dist as.dendrogram
#' @export
#' @note  No note
#' @details No details
#' @references NA
#' @examples print("No examples")
hClustHeatMap <- function(eset
                          ,conditionColors =.getConditionColors(eset)
                          ,breaks=seq(-2,2,length=20)
                          ,dendogram = "column"
                          ,legendPos="left"
                          ,...
){

  conditionColors =.getConditionColors(eset)
  dendogram = "column"
  legendPos="left"

  # need more than own data feature
  if(nrow(eset) <= 1 ) return()

  # do not plot filtered
  eset <- eset[!fData(eset)$isFiltered,]

  ### we need at least two conditions
  if(ncol(eset) == 1){
    return(cat("INFO: Only one condition no hClustHeatMap \n"))
  }

  # mean center
  expMatrixMC = log2(exprs(eset))
  expMatrixMC = (expMatrixMC - rowMeans(expMatrixMC))

  # cluster rows
  feature.tree = hclust(dist(expMatrixMC,method = "euclidean"), method="ward.D2")

  # cluster columns based on log column intensites
  msrun.cor.pearson = cor(log2(exprs(eset)), use="pairwise.complete.obs", method="pearson")
  msrun.cor.pearson.dist = as.dist(1-msrun.cor.pearson)
  ### to avoid error when replicates of the same condition are identical, DOES THIS EVER HAPPEN?
  msrun.cor.pearson.dist[is.na(msrun.cor.pearson.dist)] <- 0
  msrun.tree = hclust(msrun.cor.pearson.dist, method="ward.D2")

  ### sample colors
  samplecolors =  as.vector(unlist(conditionColors[pData(eset)$condition,]))
  labRow <- rownames(expMatrixMC)

  ### do not display feature names if too many
  if(nrow(expMatrixMC) > 300){
    labRow <- rep("",(nrow(expMatrixMC)))
  }

  hm <- heatmap.2(as.matrix(expMatrixMC)
                  , col=colorRampPalette(c(colors()[142],"black",colors()[128]))
                  , scale="none"
                  , key=TRUE
                  , symkey=FALSE
                  , trace="none"
                  , cexRow=0.5
                  , cexCol = 0.7
                  ,ColSideColors=samplecolors
                  ,labRow = labRow
                  ,Rowv=as.dendrogram(feature.tree)
                  ,Colv=as.dendrogram(msrun.tree)
                  ,dendrogram=dendogram
                  ,density.info="density"
                  #,KeyValueName="Prob. Response"
                  ,breaks=breaks
                  ,na.rm=T
                  ,key.title=""
                  ,key.xlab = "Mean centered log2 intensity"
                  ,key.ylab=""
                  , ...
  )

  legend(legendPos,levels(pData(eset)$condition), fill=as.character(conditionColors[,1]), cex=0.7, box.col=0)


}