Nothing
R/VolcanoPlotW.R
In wrGraph: Graphics in the Context of Analyzing High-Throughput Data
Defines functions
.colorByPvalue
VolcanoPlotW
Documented in
.colorByPvalue
VolcanoPlotW
#' Volcano-Plot (Statistical Test Outcome versus Relative Change)
#'
#' This type of plot is very common in high-throughput biology, see \href{https://en.wikipedia.org/wiki/Volcano_plot_(statistics)}{Volcano-plot}.
#' Basically, this plot allows comparing the outcome of a statistical test to the relative change (ie log fold-change, M-value).
#'
#' @details
#' In high-throughput biology data are typically already transformed to log2 and thus, the 'M'-values (obtained by subtrating two group means) represent a relative change.
#' Output from statistical testing by \code{\link[wrMisc]{moderTest2grp}} or \code{\link[wrMisc]{moderTestXgrp}} can be directly read to produce Volcano plots for diagnostic reasons.
#' Please note, that plotting a very high number of points (eg >10000) in transparency may take several seconds.
#'
#' Using the argument \code{plotAsFdr} FDR value will be used for plaotting instead of p-values. Since selection of 'significant' hits is usually done based on FDR-values, it makes sense to plot them directly.
#' However, FDR values are not fully linear to initial p-values and plateau-phenomena may occur frequently.
#' Thus, resulting plots may appear less spread, so finally drawing p-values (but cutting based on FDR-thersholds) many times results in nicer plots.
#'
#'
#' @param Mvalue (MArrayLM-object, numeric or matrix) data to plot; M-values are typically calculated as difference of log2-abundance values and 'pValue' the mean of log2-abundance values;
#' M-values and p-values may be given as 2 columsn of a matrix, in this case the argument \code{pValue} should remain NULL.
#' One may also furnish MArrayLM-objects created bpackage wrProteo or limma.
#' @param pValue (numeric, list or data.frame) if \code{NULL} it is assumed that 2nd column of 'Mvalue' contains the p-values to be used
#' @param useComp (integer, length=1) choice of which of multiple comparisons to present in \code{Mvalue} (if generated using \code{moderTestXgrp()})
#' @param filtFin (matrix or logical) The data may get filtered before plotting: If \code{FALSE} no filtering will get applied; if matrix of \code{TRUE}/\code{FALSE} it will be used as optional custom filter, otherwise (if \code{Mvalue} if an \code{MArrayLM}-object eg from limma) a default filtering based on the \code{filtFin} element will be applied
#' @param tit (character) custom title (has priority over \code{ProjNa})
#' @param ProjNa (character) will be added to automatic title (if custom title \code{tit} was not specified)
#' @param FCthrs (numeric) Fold-Change threshold (display as line) give as Fold-change and NOT log2(FC), default at 1.5, set to \code{NA} for omitting
#' @param FdrList (numeric) FDR data or name of list-element
#' @param FdrThrs (numeric) FDR threshold (display as line), default at 0.05, set to \code{NA} for omitting
#' @param FdrType (character) FDR-type to extract if \code{Mvalue} is 'MArrayLM'-object (eg produced by from \code{moderTest2grp} etc);
#' if \code{NULL} it will search for suitable fields/values in this order : 'FDR','BH','lfdr' and 'BY'
#' @param plotAsFdr (logical) rather plot using FDR-values (if available) instead of p-Values
#' @param subTxt (character) custom sub-title
#' @param grayIncrem (logical) if \code{TRUE}, display overlay of points (not exceeding thresholds) as increased shades of gray
#' @param col (character) custom color(s) for points of plot (see also \code{\link[graphics]{par}})
#' @param pch (integer) type of symbol(s) to plot (default=16) (see also \code{\link[graphics]{par}})
#' @param compNa (character) names of groups compared
#' @param batchFig (logical) if \code{TRUE} figure title and axes legends will be kept shorter for display on fewer splace
#' @param cexMa (numeric) font-size of title, as expansion factor (see also \code{cex} in \code{\link[graphics]{par}})
#' @param cexLa (numeric) size of axis-labels, as expansion factor (see also \code{cex} in \code{\link[graphics]{par}})
#' @param limM (numeric, length=2) range of axis M-values
#' @param limp (numeric, length=2) range of axis FDR / p-values
#' @param annotColumn (character) column names of annotation to be extracted (only if \code{Mvalue} is \code{MArrayLM}-object containing matrix $annot).
#' The first entry (typically 'SpecType') is used for different symbols in figure, the second (typically 'GeneName') is used as prefered text for annotating the best points (if \code{namesNBest} allows to do so.)
#' @param annColor (character or integer) colors for specific groups of annotation (only if \code{Mvalue} is \code{MArrayLM}-object containing matrix $annot)
#' @param expFCarrow (logical, character or numeric) optional adding arrow for expected fold-change; if \code{TRUE} the expected ratio will be extracted from numeric concentration-indications from sample-names
#' if \code{numeric} an arrow will be drawn (M-value as 1st position, color of 2nd position of vector).
#' @param cexPt (numeric) size of points, as expansion factor (see also \code{cex} in \code{\link[graphics]{par}})
#' @param cexSub (numeric) size of subtitle, as expansion factor (see also \code{cex} in \code{\link[graphics]{par}})
#' @param cexTxLab (numeric) size of text-labels for points, as expansion factor (see also \code{cex} in \code{\link[graphics]{par}})
#' @param namesNBest (integer or character) for display of labels to points in figure: if 'pass','passThr' or 'signif' all points passing thresholds; if numeric (length=1) this number of best points will get labels
#' if the initial object \code{Mvalue} contains a list-element called 'annot' the second of the column specified in argument \code{annotColumn} will be used as text
#' @param NbestCol (character or integer) colors for text-labels of best points, also used for arrow
#' @param colBySpecType (logical) incase arument \code{Mvalue} is MArrayLM-object it is possible to use different color-codes for points passing thresholds based on Mvalue$annot[,"SpecType"]; use this with multi-species benchmark tests
#' @param sortLeg (character) sorting of 'SpecType' annotation either ascending ('ascend') or descending ('descend'), no sorting if \code{NULL}
#' @param NaSpecTypeAsContam (logical) consider lines/proteins with \code{NA} in Mvalue$annot[,"SpecType"] as contaminants (if a 'SpecType' for contaminants already exits)
#' @param useMar (numeric, length=4) custom margings (see also \code{\link[graphics]{par}})
#' @param returnData (logical) optional returning data.frame with (ID, Mvalue, pValue, FDRvalue, passFilt)
#' @param silent (logical) suppress messages
#' @param callFrom (character) allow easier tracking of messages produced
#' @param debug (logical) additional messages for debugging
#' @return This function simply plots an MA-plot (to the current graphical device), if \code{returnData=TRUE} an optional data.frame with (ID, Mvalue, pValue, FDRvalue, passFilt) can be returned
#' @seealso (for PCA) \code{\link{plotPCAw}})
#' @examples
#' library(wrMisc)
#' set.seed(2005); mat <- matrix(round(runif(900),2), ncol=9)
#' rownames(mat) <- paste0(rep(letters[1:25], each=4), rep(letters[2:26],4))
#' mat[1:50,4:6] <- mat[1:50,4:6] + rep(c(-1,1)*0.1,25)
#' mat[3:7,4:9] <- mat[3:7,4:9] + 0.7
#' mat[11:15,1:6] <- mat[11:15,1:6] - 0.7
#' ## assume 2 groups with 3 samples each
#' gr3 <- gl(3, 3, labels=c("C","A","B"))
#' tRes2 <- moderTest2grp(mat[,1:6], gl(2,3), addResults = c("FDR","means"))
#' # Note: due to the small number of lines only FDR chosen to calculate
#' VolcanoPlotW(tRes2)
#' ## Add names of points passing custom filters
#' VolcanoPlotW(tRes2, FCth=1.3, FdrThrs=0.2, namesNBest="passThr")
#'
#' ## assume 3 groups with 3 samples each
#' tRes <- moderTestXgrp(mat, gr3, addResults = c("FDR","means"))
#' # Note: due to the small number of lines only FDR chosen to calculate
#' VolcanoPlotW(tRes)
#' VolcanoPlotW(tRes, FCth=1.3, FdrThrs=0.2)
#' VolcanoPlotW(tRes, FCth=1.3, FdrThrs=0.2, useComp=2)
#'
#' @export
VolcanoPlotW <- function(Mvalue, pValue=NULL, useComp=1, filtFin=NULL, tit=NULL, ProjNa=NULL, FCthrs=NULL, FdrList=NULL, FdrThrs=NULL, FdrType=NULL, plotAsFdr=FALSE,
subTxt=NULL, grayIncrem=TRUE, col=NULL, pch=16, compNa=NULL, batchFig=FALSE, cexMa=1.8, cexLa=1.1, limM=NULL, limp=NULL,
annotColumn=c("SpecType","GeneName","EntryName","Accession","Species","Contam"), annColor=NULL, expFCarrow=FALSE,cexPt=NULL, cexSub=NULL,
cexTxLab=0.7, namesNBest=NULL, NbestCol=1, colBySpecType=FALSE, sortLeg="descend", NaSpecTypeAsContam=TRUE, useMar=c(6.2,4,4,2), returnData=FALSE, callFrom=NULL, silent=FALSE, debug=FALSE) {
## MA plot
## optional arguments for explicit title in batch-mode
fxNa <- wrMisc::.composeCallName(callFrom, newNa="VolcanoPlotW")
opar <- graphics::par(no.readonly=TRUE)
opar2 <- opar[-which(names(opar) %in% c("fig","fin","font","mfcol","mfg","mfrow","oma","omd","omi"))] #
on.exit(graphics::par(opar2)) # progression ok
## during function changes in $mar,$cex.main,$cex.lab,$las
namesIn <- c(deparse(substitute(Mvalue)), deparse(substitute(pValue)), deparse(substitute(filtFin)))
basRGB <- c(0.3,0.3,0.3) # grey
fcRGB <- c(1,0,0) # red for points passing FC filt line
multiComp <- TRUE # initialize
inversedComp <- singleCompSetup <- FALSE # inversedComp: if useComp is given as text it may be inverse of p-values given; singleCompSetup : cases sof single comparison BUT 2 cols of p.values
splNa <- annot <- ptType <- colPass <- ptBg <- grpMeans <- pcol <- pwComb <- pwNames <- FDRvalue <- FdrList <- FDRty <- useCompNa <- NULL # initialize
.extrSep2 <- function(comb , indiv) { ## extract separator by striping indiv elements (iniv) from combined pairwise names (comb) .. move to wrMisc ??
fx2 <- function(x, y) sub(paste0(wrMisc::protectSpecChar(x),"$"),"", sub(paste0("^",wrMisc::protectSpecChar(x)),"", y)) # remove x from head and tail (eg rm C from end)
for(i in indiv) comb <- fx2(i, comb)
unique(comb[which(nchar(comb) >0)]) # will/might return NULL when sep==""
}
if(length(pch) <1) pch <- 16
if(isTRUE(debug)) silent <- FALSE else debug <- FALSE
if(!isTRUE(silent)) silent <- FALSE
if(debug) message(fxNa,"Length Mvalue ",length(Mvalue)," ; length pValue ",length(pValue)," ; useComp ",useComp)
if(identical(col,"FDR")) {FDR4color <- TRUE; col <- NULL} else FDR4color <- FALSE
if(length(Mvalue) <1) message("Nothing to do, 'Mvalue' seems to be empty !") else {
## data seem valid to make MAplot
if(length(cexTxLab) <0) cexTxLab <- 0.7
if("MArrayLM" %in% class(Mvalue) || "list" %in% class(Mvalue)) {
if(debug) {message(fxNa,"VPW0"); VPW0 <- list(Mvalue=Mvalue,pValue=pValue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,fxNa=fxNa )}
## try working based on MArrayLM-object (Mvalue)
if(debug) message(fxNa," '",namesIn[1],"' is list or MArrayLM-object ")
## initial check of useComp
if(length(useComp) >1) { useComp <- wrMisc::naOmit(useComp)[1]
if(!silent) message(fxNa,"Argument 'useComp' should be integer of length=1; using only 1st entry") }
if(length(useComp) <1) { useComp <- 1
if(!silent) message(fxNa,"Argument 'useComp' invalid with this dataset, setting to 1") }
## address multiple questions (via useComp): need to check 'useComp', thus need to locate FDRvalues or pValues for group-names ..
## extract data: find suitable columns
pcol <- wrMisc::naOmit(match(c("p.value","pvalue","pval","p"), tolower(names(Mvalue))))
if(length(pcol) >0) names(pcol) <- c("p.value","pvalue","pval","p")[pcol]
## NEED WAY TO CUSTOM CHOOSE WHICH FDR TO USE !!
FDRcol <- wrMisc::naOmit(match(c("fdr","bh","lfdr","by","bonferroni"), tolower(names(Mvalue))))
if(length(pcol) >0) names(pcol) <- names(Mvalue)[pcol]
if(length(FDRcol) >0) names(FDRcol) <- names(Mvalue)[FDRcol]
#if(length(FDRcol) >0) names(FDRcol) <- c("FDR","BH","lfdr","BY","bonferroni")[FDRcol]
if(debug) {message(fxNa,"VPW1"); VPW1 <- list(Mvalue=Mvalue,pValue=pValue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol )}
if(length(pcol) ==0) stop(fxNa,"Cannot find suitable element for p-values") else pcol <- pcol[1]
if(length(FDRcol) ==0) stop(fxNa,"Cannot find suitable element for FDR-values") else FDRcol <- FDRcol[1]
if("means" %in% names(Mvalue)) grpMeans <- Mvalue$means else stop(fxNa,"Need element 'means' (group-means) in object 'Mvalue'")
## need to know which of pairwise comparisons get addressed: need to understand setup -either from character useComp or index & p-values
if(debug) {message(fxNa,"VPW1b"); VPW1b <- list(Mvalue=Mvalue,pValue=pValue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol,grpMeans=grpMeans )}
if("design" %in% names(Mvalue)) {
if(length(Mvalue$design)==0 || length(dim(Mvalue$design)) !=2 || nrow(Mvalue$design) ==0) stop(fxNa,"Unable to locate element with names of pairwise groups ! (ie $design)")
## note $contrasts NOT present in result from wrMisc::moderTest2grp() !!
if(colnames(Mvalue$design)[1] =="(Intercept)" && ncol(Mvalue$design) ==2) {
## single comparison case (colnames of Mvalue$design are fixed); determine optimal sep and create pairwise name(s)
singleCompSetup <- TRUE
#Mvalue$design <- Mvalue$design[,2, drop=FALSE]
grp <- colnames(Mvalue$means) #
if(utils::packageVersion("wrMisc") < "2.0.0") {
pwSep <- " "
} else {
pwSep <- wrMisc::getPWseparator(grp=colnames(Mvalue$design), includeGrp=FALSE, silent=silent, debug=debug, callFrom=fxNa)
colnames(Mvalue$design) <- grp[1:ncol(Mvalue$design)]
}
pwNames <- if(length(Mvalue$setup$pwNames)==1) Mvalue$setup$pwNames else {if(length(grp)==2) as.matrix(grp) else utils::combn(grp, 2)}
allCompNa <- if(length(Mvalue$setup$allCompNa) >0) Mvalue$setup$allCompNa else {if(length(grp)==2) paste(grp, collapse=pwSep) else paste(pwNames[1,], pwNames[2,], sep=pwSep)}
for(i in c(pcol,FDRcol)) {
if(colnames(Mvalue[[i]])[1] =="(Intercept)") { Mvalue[[i]] <- Mvalue[[i]][,2, drop=FALSE] # remove '(Intercept)' column
colnames( Mvalue[[i]]) <- allCompNa[1]
}
}
} else {
## regular multi-comparison case
grp <- unique(colnames(Mvalue$design))
if(utils::packageVersion("wrMisc") < "2.0.0") {
pwSep <- " "
} else {
pwSep <- if(length(Mvalue$setup$sep)==1) Mvalue$setup$sep else wrMisc::getPWseparator(colnames(Mvalue[[c(pcol,FDRcol)[1]]]), grp=colnames(Mvalue$design), includeGrp=FALSE, silent=silent, debug=debug, callFrom=fxNa)
}
allCompNa <- if(length(Mvalue$setup$allCompNa) >0) Mvalue$setup$allCompNa else colnames(Mvalue[[c(pcol,FDRcol)[1] ]])
pwGrInd <- t(wrMisc::matchSampToPairw(grpNa=grp, pairwNa=allCompNa, sep=pwSep, silent=FALSE, debug=FALSE, callFrom=fxNa))
pwNames <- matrix(grp[pwGrInd], ncol=length(grp))
}
## created in this block : pwSep, pwNames
## investigate pairwise names from testing (& get separator) - need to find out which element of Mvalue
if(length(allCompNa) ==0 ) stop(fxNa,"Unable to locate element with names of pairwise groups !")
if(debug) {message(fxNa,"VPW1c"); VPW1c <- list(Mvalue=Mvalue,pValue=pValue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol,grpMeans=grpMeans,grp=grp,pwSep=pwSep,pwNames=pwNames,allCompNa=allCompNa,useComp=useComp,useCompNa=useCompNa )}
if(is.numeric(useComp)) { ## simple matching on allCompNa
if(singleCompSetup && useComp != 1) {useComp <- 1; names(useComp) <- paste(pwNames[,1], sep=pwSep)
if(!silent) message(fxNa,"Adjusting useComp to single comparison setup") }
names(useComp) <- allCompNa[useComp]
useCompNa2 <- pwNames[,useComp]
if(debug) message(fxNa," (useComp numeric) useCompNa2 ",paste(useCompNa2,collapse=" & "))
} else {
## check if name of useComp matches any of (default) pairwise
chComp <- allCompNa %in% names(useComp)
if(any(chComp, na.rm=TRUE)) {
useCompNa2 <- pwNames[,which(chComp)]
useCompNa <- useComp
useComp <- which(chComp)
names(useComp) <- useCompNa
}
if(any(chComp, na.rm=TRUE)) {
useCompNa2 <- pwNames[,which(chComp)]
useCompNa <- useComp
useComp <- which(chComp)
names(useComp) <- useCompNa
} else {
## possible reversed written - split (text-)useComp
useCompIn2 <- try(wrMisc::matchSampToPairw(grpNa=grp, pairwNa=useComp, sep=pwSep, silent=FALSE, debug=FALSE, callFrom=fxNa), silent=TRUE)
if(inherits(useCompIn2, "try-error")) stop(fxNa,"Unable to understand argument useComp ('",useComp,"') - can't match to any of the groups of replicates")
useCompNa2 <- grp[useCompIn2]
chComp <- which(colSums(matrix(pwNames %in% useCompNa2, ncol=length(grp)))==2)
if(length(chComp) !=1) stop(fxNa,"Unable to understand argument useComp ('",useComp,"') - can't match to any of pairwise groups")
useCompNa <- useComp
useComp <- chComp
names(useComp) <- useCompNa
inversedComp <- TRUE
}
if(debug) message(fxNa," (useComp char) useCompNa2 ",paste(useCompNa2,collapse=" & "))
}
## created in this block : update useComp, useCompNa2 (separate names)
if(debug) {message(fxNa,"v"); VPW1d <- list(Mvalue=Mvalue,pValue=pValue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol,grpMeans=grpMeans,pwSep=pwSep,annotColumn=annotColumn,useCompNa2=useCompNa2 )}
Mvalue$Mval <- grpMeans[,useCompNa2[1]] - grpMeans[,useCompNa2[2]]
if(inversedComp) Mvalue$Mval <- -1 * Mvalue$Mval
pValue <- if(length(pcol)==1) Mvalue[[pcol]][,useComp] else NULL
if(debug) {message(fxNa,"VPW2"); VPW2 <- list(Mvalue=Mvalue,pValue=pValue,FDRvalue=FDRvalue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol,grpMeans=grpMeans,pwSep=pwSep,useCompNa2=useCompNa2,annotColumn=annotColumn )}
} else stop(fxNa,"Argument 'Mvalue' is missing key element $design")
## recuperate FDRvalue
if(length(FDRcol)==1 && names(FDRcol) %in% names(Mvalue)) {
FDRvalue <- Mvalue[[FDRcol]][,useComp]
}
## recuperate $annot if present and use for symbol
if("annot" %in% names(Mvalue) && isTRUE(colBySpecType)) {
useAnnCol <- match(annotColumn, colnames(Mvalue$annot))
if(!is.na(useAnnCol[1])) { # annotation (for multiple groups) exists
ptType <- Mvalue$annot[,useAnnCol[1]] # SpecType
chNA <- is.na(ptType)
## associate NAs from 'SpecType' in ptType with conta ?
if(NaSpecTypeAsContam) {
chConta <- tolower(ptType) %in% c("contaminant","contam","conta","cont")
if(any(chConta)) ptType[which(is.na(ptType))] <- unique(ptType[which(chConta)])[1]}
if(any(is.na(ptType))) ptType[which(chNA)] <- "NA"
if(length(pch) < length(pValue) && length(unique(wrMisc::naOmit(ptType))) >1) {
if(length(pch) >1 && !silent) message(fxNa," (invalid pch) using default 'pch' oriented by $annot and starting from 15")
pch <- 14 + as.integer(as.factor(ptType))
}
useAnnCol <- wrMisc::naOmit(useAnnCol) }
annot <- Mvalue$annot[,useAnnCol]
if(annotColumn[1] %in% colnames(annot)) annot[,annotColumn[1]] <- ptType # update with NAs trasformed to "NA"
if(debug) {message(fxNa,"VPW2b"); VPW2b <- list(Mvalue=Mvalue,pValue=pValue,FDRvalue=FDRvalue,annot=annot,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol,grpMeans=grpMeans,pwSep=pwSep,useCompNa2=useCompNa2,annotColumn=annotColumn,useCompNa=useCompNa,pch=pch )}
}
## filtering
if("filter" %in% names(Mvalue) && ((length(filtFin)==1 && isTRUE(filtFin)) || length(filtFin)==0)) {
filtFin <- Mvalue$filter[,useComp]
#filtFin <- Mvaue$filter[,which(colnames(Mvalue$filter)==useCompNa)]
} else { ## check custom filtFin
chFi <- is.logical(filtFin) && length(filtFin)==length(pValue)
if(!chFi) { filtFin <- NULL
if(!silent) message(fxNa,"Argument 'filtFin' seems invalid (should be logical vector of length of number of elements for plotting )")
} else if(debug) message(fxNa,"Using custom filtFin")
}
if(length(filtFin) >0) {
if(debug) message(fxNa,"Filtering : Setting ",sum(!filtFin & !is.na(pValue))," additional values as NA")
pValue[which(!filtFin)] <- NA # better to set NA so that annotation, filtering etc stays same
if(length(FDRvalue) ==length(pValue)) FDRvalue[which(!filtFin)] <- NA
## Mvalue : if only 2 groups $Mval should be OK
if("Mval" %in% names(Mvalue)) Mvalue$Mval[which(!filtFin)] <- NA
}
if(debug) {message(fxNa,"VPW2c"); VPW2c <- list(Mvalue=Mvalue,pValue=pValue,FDRvalue=FDRvalue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol,grpMeans=grpMeans,pwSep=pwSep,useCompNa2=useCompNa2,annotColumn=annotColumn )}
Mvalue <- Mvalue$Mval # DISMISS rest of object !!! .....
} else {
## ... case of explicit pValue argument
if(length(pValue) <1) stop("Argument 'pValue' is required (if 'Mvalue' not 'MArrayLM'-type object) !")
if(length(dim(pValue)) >1) if(ncol(pValue) >1) {
if(!silent) message(fxNa," Note, ",namesIn[2]," has ",ncol(pValue)," columns, using last column")
pNa <- rownames(pValue)
pValue <- as.numeric(pValue[,ncol(pValue)] )
names(pValue) <- pNa }
FDRvalue <- if(length(FdrList) <1) NULL else FdrList
}
if(debug) {message(fxNa,"VPW8"); VPW8 <- list(FdrList=FdrList,FDRvalue=FDRvalue, Mvalue=Mvalue,pValue=pValue,FdrType=FdrType,plotAsFdr=plotAsFdr,pValue=pValue,useComp=useComp,grpMeans=grpMeans)}
## need to introduce -log10 to pValue
chNA <- is.na(pValue)
if(all(chNA)) stop(fxNa,"All p-values are NA, nothing to draw !")
if(debug && any(pValue <0, na.rm=TRUE)) message(fxNa,"Some p-values are negative, this should not be ! Maybe log values were given by error ?")
pValue <- -log10(pValue)
## check for (same) order, adjust Mvalue & pValue according to names
chNa <- list(MNa=if(length(dim(Mvalue)) >1) rownames(Mvalue) else names(Mvalue),
pNa=if(length(dim(pValue)) >1) rownames(pValue) else names(pValue))
nIni <- c(M=length(Mvalue), p=length(pValue))
if(length(chNa$MNa) >0 && length(chNa$pNa) >0) { # ie both have names, so one can match names
if(!all(chNa$MNa==chNa$pNa, na.rm=TRUE)) {
matchNa <- wrMisc::naOmit(match(chNa$MNa, chNa$pNa))
if(length(matchNa) <1) stop("Both 'Mvalue' and 'pValue' have names, but none of them match !!")
if(!all(matchNa, 1:length(pValue), na.rm=TRUE)) {
pValue <- pValue[matchNa]
Mvalue <- wrMisc::naOmit(Mvalue[match(names(pValue), names(Mvalue))]) }
}
} else {
if(length(Mvalue) != length(pValue)) stop("p- and M- values have different length, but no names to match !! (M=",length(Mvalue)," vs p=",length(pValue),")")
}
if(length(grpMeans) <1) grpMeans <- matrix(rep(NA,2*length(Mvalue)), ncol=2, dimnames=list(names(Mvalue),c("mean1","mean2")))
if(debug) {message(fxNa,"VPW9"); VPW9 <- list(FdrList=FdrList,Mvalue=Mvalue,FdrType=FdrType,FDRvalue=FDRvalue,plotAsFdr=plotAsFdr,pValue=pValue,useComp=useComp,grpMeans=grpMeans, filtFin=filtFin,annotColumn=annotColumn,annot=annot, pch=pch,sortLeg=sortLeg,batchFig=batchFig)} # pairwCol=pairwCol,
## prepare/integrate FILTERING
if(length(filtFin) >0) {
## if filtFin is matrix use each line with min 1 instance of TRUE,
if(length(dim(filtFin)) >1) filtFin <- as.logical(as.matrix(filtFin)[,useComp]) # use rows with >= 1 TRUE
if(length(names(filtFin)) >0) {
matchNa <- wrMisc::naOmit(match(rownames(pValue), names(filtFin)))
if(length(matchNa)==length(pValue)) filtFin <- as.logical(filtFin[matchNa])
}
} else filtFin <- rep(TRUE, nrow(grpMeans))
if(debug) {message(fxNa,"VPW9b"); VPW9b <- list(FdrList=FdrList,Mvalue=Mvalue,FdrType=FdrType,FDRvalue=FDRvalue,plotAsFdr=plotAsFdr,pValue=pValue,useComp=useComp,grpMeans=grpMeans, filtFin=filtFin,annotColumn=annotColumn,annot=annot, pch=pch,sortLeg=sortLeg,batchFig=batchFig)} # pairwCol=pairwCol,
## start creating merged data for plot (& export)
merg <- if(length(annot) >0) data.frame(ID=NA, grpMeans, Mvalue=Mvalue, pValue=pValue,
FDR=if(length(FDRvalue) >0) as.numeric(FDRvalue) else rep(NA,length(pValue)),
filtFin=filtFin, annot, pch=pch, stringsAsFactors=FALSE) else {
data.frame(ID=NA, grpMeans, Mvalue=Mvalue, pValue=pValue, FDR=if(length(FDRvalue) >0) as.numeric(FDRvalue) else rep(NA,length(pValue)),
filtFin=filtFin, pch=pch, stringsAsFactors=FALSE) }
if(length(names(Mvalue)) >0) merg[,1] <- names(Mvalue) else {if(length(names(pValue)) >0) merg[,1] <- names(pValue)}
## replace NA in 'SpecType' by 'NA'
if(annotColumn[1] %in% colnames(merg)) { chNa <- is.na(merg[,annotColumn[1]]) # replace NAs in col "SpecType" by "NA"
if(any(chNa)) merg[which(chNa),annotColumn[1]] <- "NA"
} else { merg <- cbind(merg, rep(1,nrow(merg))) # add colum for 'SpecType'
colnames(merg)[ncol(merg)] <- annotColumn[1] }
if(debug) {message(fxNa,"VPW10"); VPW10 <- list(merg=merg,FdrList=FdrList,Mvalue=Mvalue,FdrType=FdrType,FDRvalue=FDRvalue,plotAsFdr=plotAsFdr,pValue=pValue,useComp=useComp,grpMeans=grpMeans, filtFin=filtFin,annotColumn=annotColumn,annot=annot, pch=pch,sortLeg=sortLeg,batchFig=batchFig)}
#check#plot(pValue ~Mvalue, data=VPW10$merg) # ok
## adjust col & pch
if(!any(c(1,length(Mvalue)) %in% length(pch))) {
if(!silent) message(fxNa,"Argument 'pch' should be either length=1 or correspond to length of data, reset to default=16")
pch <- 16 }
if(length(col) >1 && length(col) <length(Mvalue)) {
if(!silent) message(fxNa,"Argument 'col' should be either length=1 or correspond to length of data, reset to default=NULL")
col <- NULL }
if(debug) message(fxNa," ++ DONE extracting columns : ",wrMisc::pasteC(colnames(merg),quo="'"))
## apply filtering
msg <- " data provided in 'Mvalue' and 'pValue' "
if(debug) {message(fxNa,"VPW10")}
if(!silent && nrow(merg) < round(length(Mvalue)/10)) message(fxNa," .. note : less than 10% of",msg," were matched") else {
if(!silent && nrow(merg) < length(Mvalue)/2) message(fxNa," .. NOTE : less than 50% of",msg," were matched !!")}
if(debug) message(fxNa,msg," were matched to ",nrow(merg)," common entries")
## apply filtering (keep all lines where at least one condition passes)
if(length(filtFin) >0 && any(isFALSE(filtFin), na.rm=TRUE)) { # use filtering provided
if(sum(filtFin) >0 && sum(filtFin) < nrow(merg)) {
whFilt <- which(merg$filtFin)
if(length(pch) >1) pch <- pch[whFilt]
if(length(col) >1) col <- col[whFilt]
merg <- merg[whFilt,]
if(!silent) message(fxNa,"Filtered (based on 'filtFin') from ",length(filtFin)," to ",nrow(merg)," lines")
}
} else filtFin <- rep(TRUE, nrow(merg))
if(debug) {message(fxNa,"VPW11"); VPW11 <- list(merg=merg,FdrList=FdrList,Mvalue=Mvalue,FdrType=FdrType,FDRvalue=FDRvalue,plotAsFdr=plotAsFdr,pValue=pValue,useComp=useComp,grpMeans=grpMeans, filtFin=filtFin,annotColumn=annotColumn,annot=annot, pch=pch,sortLeg=sortLeg,batchFig=batchFig)}
#check#plot(pValue ~Mvalue, data=VPW11$merg) # ok
## sort merg, so that legend always gets constructed the same order, ascending ('ascend') or descending ('descend')
sortLeg <- if(identical(sortLeg,"ascend")) FALSE else {if(identical(sortLeg,"descend")) TRUE else NULL}
if(length(sortLeg) ==1 && annotColumn[1] %in% colnames(merg)) merg <- merg[order(merg[,annotColumn[1]], decreasing=sortLeg),]
## update ..
nIDco <- sum(c("ID","nredID","uniqID") %in% colnames(merg)) # number of heading columns in 'merg'
Mvalue <- as.numeric(if("Mvalue" %in% colnames(merg)) merg[,"Mvalue"] else merg[,nIDco+1])
pValue <- as.numeric(if("pValue" %in% colnames(merg)) merg[,"pValue"] else {
if(length(dim(Mvalue)) >0) merg[,ncol(Mvalue) +nIDco +1] else merg[,nIDco+2]})
FDRvalue <- as.numeric(if("FDR" %in% colnames(merg)) merg[,"FDR"] else merg[,grep("(fdr)|(bh)", tolower(colnames(merg)))[1]])
if("Lfdr" %in% colnames(merg)) FdrList <- merg[,"Lfdr"] else {
if("lfdr" %in% colnames(merg)) FdrList <- merg[,"lfdr"]}
pch <- merg[,"pch"] # update
ptType <- if(annotColumn[1] %in% colnames(merg)) merg[,annotColumn[1]] else rep(1,nrow(merg)) # update "SpecType"
if(debug) {message(fxNa,"VPW12"); VPW12 <- list(merg=merg,Mvalue=Mvalue,filtFin=filtFin)}
## prepare for plotting
if(is.null(cexSub)) cexSub <- cexLa +0.05
xLab <- "M-value (log2 fold-change)"
tit1 <- if(length(tit)==1 && nchar(tit) >0) tit else paste(c(if(!batchFig) c(ProjNa, if(!is.null(ProjNa)) ": ","Volcano-Plot"),
if(!is.null(compNa)) {if(length(compNa)==1) compNa else c(compNa[1]," vs ",compNa[2])}), collapse=" ") # but what title if batchFig=NULL & compNa=NULL -> only "Volcano-plot"
if(length(FCthrs) <1) FCthrs <- 1.5
if(length(FdrThrs) <1) FdrThrs <- 0.05
if(debug) {message(fxNa,"VPW12b"); VPW12b <- list(merg=merg,FdrList=FdrList,Mvalue=Mvalue,FdrType=FdrType,FDRvalue=FDRvalue,plotAsFdr=plotAsFdr,pValue=pValue,useComp=useComp,grpMeans=grpMeans, filtFin=filtFin,annotColumn=annotColumn,annot=annot, pch=pch,sortLeg=sortLeg,batchFig=batchFig,FCthrs=FCthrs,FdrThrs=FdrThrs)}
## count no of passing
passFC <- if(length(FCthrs) ==1 && !any(is.na(FCthrs))) abs(merg[,"Mvalue"]) >= log2(FCthrs) else merg[,"filtFin"] ## convert FCthrs to log2
passFdr <- if(length(FdrThrs) ==1 && !any(is.na(FdrThrs)) && "FDR" %in% colnames(merg)) {merg[,"FDR"] <= FdrThrs} else merg[,"filtFin"]
passAll <- merg[,"filtFin"] & passFC & passFdr
if(debug) {message(fxNa,"VPW12c")}
chNA <- is.na(passAll) # passFdr may contain NAs
if(any(chNA)) passAll[which(chNA)] <- FALSE
if(debug) message(fxNa," ",sum(passFC,na.rm=TRUE)," passing FCthrs ; ",sum(passFdr,na.rm=TRUE)," passing FdrThrs ; combined ",sum(passAll,na.rm=TRUE))
if(debug) {message(fxNa,"VPW13"); VPW13 <- list(merg=merg,Mvalue=Mvalue,filtFin=filtFin,col=col,passAll=passAll,cexPt=cexPt,annotColumn=annotColumn,annColor=annColor)}
## color for points passing filter
if(length(col) >0) if(length(col) != nrow(merg)) { col <- NULL
if(!silent) message(fxNa," invalid entry for 'col', should be of length=",nrow(merg),", resetting to default")}
if(length(col) <1) {
alph <- sort(c(0.14, round(0.6/log10(length(Mvalue)),2), 0.8))[2] # alph <- round(12/sqrt(nrow(eBayesLst$pValue)),2)
alph2 <- sort(c(round(7/(5 +sum(passAll)^0.7),2), alph,0.9))[2] # for points passing thresholds
useCol <- if(grayIncrem) grDevices::rgb(0.35,0.35,0.35,alph) else grDevices::rgb(0.7,0.7,0.7) # basic color
useCex <- if(length(cexPt) >0) cexPt else max(round(0.8 + 2/(1 +sum(filtFin, na.rm=TRUE))^0.28,2), 1.1)
chCol <- unique(merg[, annotColumn[1]]) # check how many different colors may be needed
chNaC <- is.na(chCol)
if(any(chNaC)) chCol[which(chNaC)] <- "NA"
if(length(annColor) >0) {colPass <- annColor} else if(length(chCol) >4) {
colPass <- cbind(red=c(141,72,90,171, 220,253,244,255), green=c(129,153,194,221, 216,174,109,0), blue=c(194,203,185,164, 83,97,67,0))
colPass <- grDevices::rgb(red=colPass[,1], green=colPass[,2], blue=colPass[,3], alph2, maxColorValue=255)
if(length(chCol) >8) { colPass <- c(colPass, rep(colPass[8], length(chCol) -8))
if(!silent) message(fxNa," > 8 different groups found, using 8th color after 7th group")}
} else colPass <- grDevices::rgb(c(0.95,0.2,0,0.75), c(0.15,0.2,0.9,0.35), c(0.15,0.95,0,0.8), alph2) # red, blue, green, purple (luminosity adjusted)
useCol <- rep(useCol[1], nrow(merg)) # fuse basic gray to colors for different types
## integrate names of annColor as order of colPass
if(length(names(annColor)) >0) {
uniTy <- unique(merg[which(passAll), annotColumn[1]])
colPass <- colPass[match(names(annColor), uniTy)]
}
if(debug) {message(fxNa,"VPW13a"); VPW13a <- list(useCol=useCol,colPass=colPass, merg=merg,Mvalue=Mvalue,filtFin=filtFin,col=col,passAll=passAll,cexPt=cexPt,annotColumn=annotColumn,annColor=annColor)}
## assign color for those passing
useCol[which(passAll)] <- if(any(passAll, na.rm=TRUE)) colPass[if(length(unique(merg[which(passAll),annotColumn[1]])) >1) wrMisc::levIndex(merg[which(passAll),annotColumn[1]]) else rep(1,sum(passAll))] # assign colors for those passing
} else useCol <- col
if(debug) {message(fxNa,"VPW13b"); VPW13b <- list(merg=merg,Mvalue=Mvalue,filtFin=filtFin,col=col,passAll=passAll,cexPt=cexPt,annotColumn=annotColumn,annColor=annColor,pch=pch,useCol=useCol)}
#}} #ok
## adjust fill color for open symbols
chPch <- pch %in% c(21:25)
if(any(chPch, na.rm=TRUE)) { ptBg <- useCol
ptBg[which(chPch)] <- useCol[which(chPch)] # background color for filled symbols
useCol[which(chPch)] <- 1 # contour as black
}
## main graphics
if(length(Mvalue) >0) {
pl1 <- try(graphics::par(mar=if(length(useMar)==4) useMar else c(6.5,4,4,2), cex.main=cexMa, las=1), silent=TRUE)
if(inherits(pl1, "try-error")) {Mvalue <- NULL; message("UNABLE TO SET PLOT MARGINS !! check plotting device ...")}
if(debug) {message(fxNa,"VPW14a"); VPW14a <- list(pl1=pl1,merg=merg,Mvalue=Mvalue,filtFin=filtFin,useMar=useMar,FCthrs=FCthrs,FdrThrs=FdrThrs,FdrType=FdrType,plotAsFdr=plotAsFdr)}
}
if(length(Mvalue) >0) {
## rather directly plot FDR
if(!"FDR" %in% colnames(merg)) plotAsFdr <- FALSE
#check# plot(pValue ~Mvalue, data=VPW14a$merg)
pl1 <- try(graphics::plot(Mvalue, if(plotAsFdr) log10(merg[,"FDR"]) -min(log10(merg[,"FDR"]), na.rm=TRUE) else merg[,"pValue"], pch=pch, cex=useCex, main=tit1,
ylab=if(plotAsFdr) paste("- log10",if(is.null(FdrType)) "(probability)" else FdrType) else "- log10 p-value (uncorrected)", col=useCol, xlab=xLab, cex.lab=cexLa, xlim=limM,ylim=limp, pt.bg=ptBg), silent=TRUE)
if(inherits(pl1, "try-error")) {Mvalue <- NULL; message("UNABLE TO PLOT !! check plotting device ...")} }
if(debug) {message(fxNa,"VPW14b"); VPW14b <- list(merg=merg,Mvalue=Mvalue,filtFin=filtFin,passAll=passAll,useMar=useMar,FCthrs=FCthrs,FdrThrs=FdrThrs,useCol=useCol,multiComp=multiComp,useComp=useComp,useCompNa=useCompNa,subTxt=subTxt,plotAsFdr=plotAsFdr )}
if(length(Mvalue) >0) { # useCompNa
sTxt <- if(length(subTxt) ==1) subTxt else { if(multiComp) if(length(names(useComp)) >0) paste0(names(useComp),"; ",collapse="") else paste0("comp=",useComp,"; ")}
sTxt <- paste0(sTxt,"n=",length(Mvalue),
if(!all(is.na(c(FCthrs,FdrThrs)))) paste(";",sum(passAll, na.rm=TRUE),"(color) points passing",
if(!is.na(FCthrs)) paste0("(FCthr=", as.character(FCthrs),", ") else " (", paste0("FdrThrs=",as.character(FdrThrs),")")))
graphics::mtext(sTxt,cex=0.75,line=0.2)
if(!all(is.na(c(FCthrs,FdrThrs)))) {
if(debug) message(fxNa," n=",length(Mvalue)," FCthrs=",as.character(FCthrs)," filt.ini=", sum(filtFin, na.rm=TRUE),
" passAll=",sum(passAll,na.rm=TRUE)," ; range Mva ",wrMisc::pasteC(signif(range(Mvalue,na.rm=TRUE),3))," ; alph=",alph," useCex=",useCex," alph2=",alph2)
## previous versions : adjusted abline by adding diff(graphics::par("usr")[1:2])/1000)
graphics::abline(v=c(-1,1)* log2(FCthrs) , col=grDevices::rgb(0.87,0.72,0.72), lty=2) }
if(debug) {message(fxNa,"VPW15"); VPW15 <- list(merg=merg,Mvalue=Mvalue,filtFin=filtFin,FCthrs=FCthrs,FdrThrs=FdrThrs,passFdr=passFdr,annotColumn=annotColumn,xLab=xLab,tit1=tit1,col=col,useCol=useCol,assAll=passAll,grayIncrem=grayIncrem,cexPt=cexPt,annColor=annColor,plotAsFdr=plotAsFdr,limM=limM,limp=limp,ptBg=ptBg,cexLa=cexLa,passAll=passAll)}
if(sum(passFdr, na.rm=TRUE) >0) {
if(plotAsFdr) {
if(any(passAll, na.rm=TRUE)) graphics::abline(h=-1*log10(max(merg[passAll,"FDR"], na.rm=TRUE)) -diff(graphics::par("usr")[3:4])/400, col=grDevices::rgb(0.87,0.72,0.72), lty=2)
} else {
if(plotAsFdr) graphics::mtext("Note, that FDR and p-value may not correlate perfectly, thus points may appear at good p-value but finally don't get retained",line=5,cex=0.65, side=1)
pRa <- range(merg[which(passFdr),"pValue"], na.rm=TRUE)
graphics::abline(h=pRa[1] +diff(graphics::par("usr")[3:4])/400, col=grDevices::rgb(0.87,0.72,0.72), lty=2) }}
## add names to best points
if(length(namesNBest) >0) {
if(any(sapply( c("pass","passThr","passFDR","signif"), identical, namesNBest))) namesNBest <- sum(passAll)
if(!is.integer(namesNBest)) namesNBest <- try(as.integer(namesNBest), silent=TRUE)
if(inherits(namesNBest, "try-error")) { namesNBest <- NULL
message(fxNa,"Unable to understand argument 'namesNBest', must be integer or 'pass','passThr','passFDR' or 'signif' (for display of labels to points in figure)") }
}
if(debug) {message(fxNa,"VPW16")}
if(length(namesNBest) >0) {
if(namesNBest >0 && any(passAll)) {
useLi <- if(any(!passAll)) which(passAll) else 1:nrow(merg)
tmP <- as.numeric(merg[useLi,"pValue"])
names(tmP) <- rownames(merg)[useLi]
## look for more informative names to display
if(length(annot) >0) {
proNa <- annot[match(names(tmP), rownames(annot)), annotColumn[2]] # normally 'Description'
chNa <- is.na(proNa)
if(!all(chNa)) names(tmP)[which(!chNa)] <- proNa[which(!chNa)]
}
useL2 <- order(tmP, decreasing=TRUE)[1:min(namesNBest,sum(passAll))]
xOffs <- signif(diff(graphics::par("usr")[1:2])/170,3)
yOffs <- signif(diff(graphics::par("usr")[3:4])/90,3)
noNa <- if(is.null(names(tmP[useL2]))) 1:length(tmP) else which(is.na(names(tmP)[useL2]))
if(length(noNa) >0 && all(annotColumn %in% colnames(merg))) names(tmP)[useL2[noNa]] <- merg[useLi[useL2[noNa]], wrMisc::naOmit(match(annotColumn[-1], colnames(merg)))[1]]
if(length(NbestCol) <1) NbestCol <- 1
if(is.null(names(tmP[useL2]))) {if(!silent) message(fxNa,"No names available for displaying names of best in plot")
} else { displTx <- wrMisc::naOmit(match(annotColumn[2:4], colnames(merg)))
displTx <- if(length(displTx) >0) cbind(merg[useLi[useL2], displTx], rownames(merg)[useLi[useL2]]) else as.matrix(useLi[useL2])
chNa <- rowSums(!is.na(displTx)) >0
if(any(chNa)) displTx[which(chNa),1] <- "unknown"
displTx <- apply(displTx,1, function(x) wrMisc::naOmit(x)[1])
graphics::text(Mvalue[useLi[useL2]] +xOffs, yOffs +if(plotAsFdr) log10(merg[useLi[useL2],"FDR"]) -min(log10(merg[,"FDR"]), na.rm=TRUE) else merg[useLi[useL2],"pValue"],
names(tmP)[useL2], cex=cexTxLab, col=NbestCol, adj=0) }
}
}
if(debug) {message(fxNa,"VPW17"); VPW17 <- list(pch=pch,passAll=passAll,merg=merg,annot=annot,annotColumn=annotColumn)}
## legend (if multiple symbols)
pch[which(is.na(pch))] <- -2
ch1 <- unique(pch)
if(length(ch1) >1 && sum(passAll) >0) {
legInd <- which(!duplicated(merg[which(passAll), annotColumn[1]], fromLast=FALSE))
if(length(legInd) <1 && !silent) message(fxNa,"Trouble ahead : Can't find non-duplicated ",annotColumn[1]," for ",sum(passAll)," points passing thresholds ! (ie as 'legInd')")
legPch <- pch[which(passAll)[legInd]]
legCol <- useCol[which(passAll)[legInd]]
legBg <- ptBg[which(passAll)[legInd]]
if(alph2 <1) {legCol <- substr(legCol,1,7); legBg <- substr(legBg,1,7)} # reset to no transparency
legLab <- merg[which(passAll)[legInd], annotColumn[1]]
chNa <- is.na(legLab)
if(any(chNa)) legLab[chNa] <- "NA"
legOr <- if(length(legLab) >1) order(legLab) else 1 # not used so far
legLoc <- checkForLegLoc(cbind(Mvalue, pValue), sampleGrp=legLab, showLegend=FALSE)
legCex <- stats::median(c(useCex, cexTxLab, 1.2), na.rm=TRUE)
if(length(legLab) >0) { chLeg <- try(graphics::legend(legLoc$loc, legend=legLab, col=legCol, text.col=1, pch=legPch,
if(length(ptBg) >0) pt.bg=ptBg, cex=legCex, pt.cex=1.2*legCex, xjust=0.5, yjust=0.5), silent=TRUE) # as points
if(inherits(chLeg, "try-error") && !silent) message(fxNa,"Note: Failed to add legend .. ",chLeg) }
} else legCol <- NULL
if(debug) {message(fxNa,"VPW18")}
## arrow for expected ratio
## how to extract in smartest way ??
drawArrow <- if(length(expFCarrow) >0) {isTRUE(as.logical(expFCarrow[1])) || grepl("^[[:digit:]]+$", expFCarrow[1])} else FALSE
if(drawArrow) {
regStr <-"[[:space:]]*[[:alpha:]]+[[:punct:]]*[[:alpha:]]*"
if(isTRUE(expFCarrow)) {
## automatic extraction of FC
if(debug) message(fxNa," .. ",names(useComp)," -> ",paste(unlist(strsplit(names(useComp), "-")), collapse=" "),"\n")
expM <- sub(paste0("^",regStr),"", sub(paste0(regStr,"$"), "", unlist(strsplit(names(useComp), "-")))) # assume '-' separator from pairwise comparison
arrCol <- 1
chN2 <- try(as.numeric(expM), silent=TRUE)
if(inherits(chN2, "try-error")) { drawArrow <- FALSE #expM <- arrCol <- NA
} else {
expM <- log2(chN2[2] / chN2[1]); arrCol <- 1 # transform automatic extracted to ratio # expFCarrow <- c(expM, 1)
}
msg <- "Argument 'expFCarrow' was set to TRUE; "
if(!silent) message(fxNa, msg, if(drawArrow) "Extract concentration values of group-names for calculating ratios" else "Failed automatic extraction of FC-values (possibly no numeric values in setup)")
} else {
expM <- try(as.numeric(expFCarrow[1]), silent=TRUE)
if(inherits(expM, "try-error")) {
if(!silent) message(fxNa,"Argument 'expFCarrow' was given, but no numeric content in 1st place found, can't draw FC-arrow") #understand 1st value of'expFCarrow' (should be numeric-like)
drawArrow <- FALSE
} else arrCol <- if(length(expFCarrow) >1) expFCarrow[2] else 1
}
if(length(expFCarrow) >1) {arrCol <- expFCarrow[2] # 2nd value of expFCarrow for color
} else arrCol <- 1 # autom : use 3rd if possible
}
if(drawArrow) {
if(is.finite(expM)) {
figCo <- graphics::par("usr") # c(x1, x2, y1, y2)
figRa <- diff(range(figCo[1:2]))*0.1
if(expM > figCo[2] +figRa || expM < figCo[1] -figRa) {
if(!silent) message(fxNa,"Can't draw arrow, ",round(expM,2)," is too far outside the plotting frame")
expM <- NA }
}
if(is.finite(expM)) {
arr <- c(0.019,0.14) # start- and end-points of arrow (as relative to entire plot)
graphics::arrows(expM, figCo[3] + arr[1]*(figCo[4] -figCo[3]), expM, figCo[3] + arr[2]*(figCo[4] -figCo[3]),
col=arrCol, lwd=1, length=0.1)
graphics::mtext(paste("expect at",signif(expM,3)), at=expM, side=1, adj=0.5, col=arrCol, cex=cexLa*0.7, line=-0.9)
} else { if(!silent) message(fxNa,"Unable to draw arrow for expexted M-value)") }
}
if(debug) {message(fxNa,"VPW19")}
tmp <- try(graphics::par(mar=opar$mar, cex.main=opar$cex.main, las=opar$las), silent=TRUE) # resetto previous
## export data used for plotting
if(returnData) {
merg <- merg[,-1*c(1,ncol(merg))] # remove col 'ID' 'redundant' & 'pch'
annCo <- wrMisc::naOmit(match(annotColumn, colnames(merg)))
if(length(annCo) >0) cbind(merg[,annCo], merg[,-annCo]) else merg }
} }
}
#' Colors based on p-Values
#'
#' This function helps defining color-gradient based on p-Values.
#' This fuction requires package RColorBrewer being installed
#'
#' @param x (numeric) p-values (main input)
#' @param br (numeric) custom breaks (used with cut)
#' @param col custom colors (must be of length(br) -1)
#' @param asIndex (logical) custom breaks (used with cut)
#' @param silent (logical) suppress messages
#' @param callFrom (character) allow easier tracking of messages produced
#' @param debug (logical) supplemental messages for debugging
#' @return This function retruns a color-gradient based on p-Values
#' @seealso (for PCA) \code{\link{plotPCAw}})
#' @examples
#' .colorByPvalue((1:10)/10)
#' @export
.colorByPvalue <- function(x, br=NULL, col=NULL, asIndex=FALSE, silent=FALSE, debug=FALSE, callFrom=NULL) {
## this function should ultimately get incorporated to wrMisc::colorAccording2 as option "pValue" or "FDR"
## colors significant (***,**,*) as red-tones, 0.05-0.075 as pale orange, others as blue-tones
## 'br' .. custom breaks (used with cut)
## 'col' .. custom colors (must be of length(br) -1)
## 'asIndex' .. optional returning of index to \code{col} instead of final color designation
fxNa <- ".colorByPvalue"
if(!isTRUE(silent)) silent <- FALSE
if(isTRUE(debug)) silent <- FALSE else debug <- FALSE
if(length(br) >0) br <- sort(unique(wrMisc::naOmit(br)))
if(length(col) != length(br) -1 && length(br) >0) { br <- col <- NULL
if(!silent) message(fxNa, "Entries for 'col' and 'br' don't match, ignoring") }
if(length(col) <1) {
hasDep <- requireNamespace("RColorBrewer", quietly=TRUE)
col <- if(hasDep) grDevices::colorRampPalette(RColorBrewer::brewer.pal(n=7, name="RdYlBu"))( 22 )[c(2:4,7,17:20)] else {
grDevices::heat.colors(n=length(unique(x)) +2)
if(!silent) message(fxNa,"Please install first from CRAN the package 'RColorBrewer'")
}
if(length(col) != length(br)) br <- sort(c(1e-20,10^(-3:0), 0.05,0.075,0.2,0.3))
}
out <- if(asIndex) as.numeric(cut(x, br)) else col[as.numeric(cut(x, br))]
if(asIndex) out <- as.numeric(cut(x, br)) else {
out <- col[as.numeric(cut(x, br))]
chNa <- is.na(x)
if(any(chNa)) out[which(chNa)] <- grDevices::grey(0.4) }
out }
Try the wrGraph package in your browser
Any scripts or data that you put into this service are public.
wrGraph documentation built on March 17, 2026, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .colorByPvalue VolcanoPlotW
Documented in .colorByPvalue VolcanoPlotW
#' Volcano-Plot (Statistical Test Outcome versus Relative Change)
#'
#' This type of plot is very common in high-throughput biology, see \href{https://en.wikipedia.org/wiki/Volcano_plot_(statistics)}{Volcano-plot}.
#' Basically, this plot allows comparing the outcome of a statistical test to the relative change (ie log fold-change, M-value).
#'
#' @details
#' In high-throughput biology data are typically already transformed to log2 and thus, the 'M'-values (obtained by subtrating two group means) represent a relative change.
#' Output from statistical testing by \code{\link[wrMisc]{moderTest2grp}} or \code{\link[wrMisc]{moderTestXgrp}} can be directly read to produce Volcano plots for diagnostic reasons.
#' Please note, that plotting a very high number of points (eg >10000) in transparency may take several seconds.
#'
#' Using the argument \code{plotAsFdr} FDR value will be used for plaotting instead of p-values. Since selection of 'significant' hits is usually done based on FDR-values, it makes sense to plot them directly.
#' However, FDR values are not fully linear to initial p-values and plateau-phenomena may occur frequently.
#' Thus, resulting plots may appear less spread, so finally drawing p-values (but cutting based on FDR-thersholds) many times results in nicer plots.
#'
#'
#' @param Mvalue (MArrayLM-object, numeric or matrix) data to plot; M-values are typically calculated as difference of log2-abundance values and 'pValue' the mean of log2-abundance values;
#' M-values and p-values may be given as 2 columsn of a matrix, in this case the argument \code{pValue} should remain NULL.
#' One may also furnish MArrayLM-objects created bpackage wrProteo or limma.
#' @param pValue (numeric, list or data.frame) if \code{NULL} it is assumed that 2nd column of 'Mvalue' contains the p-values to be used
#' @param useComp (integer, length=1) choice of which of multiple comparisons to present in \code{Mvalue} (if generated using \code{moderTestXgrp()})
#' @param filtFin (matrix or logical) The data may get filtered before plotting: If \code{FALSE} no filtering will get applied; if matrix of \code{TRUE}/\code{FALSE} it will be used as optional custom filter, otherwise (if \code{Mvalue} if an \code{MArrayLM}-object eg from limma) a default filtering based on the \code{filtFin} element will be applied
#' @param tit (character) custom title (has priority over \code{ProjNa})
#' @param ProjNa (character) will be added to automatic title (if custom title \code{tit} was not specified)
#' @param FCthrs (numeric) Fold-Change threshold (display as line) give as Fold-change and NOT log2(FC), default at 1.5, set to \code{NA} for omitting
#' @param FdrList (numeric) FDR data or name of list-element
#' @param FdrThrs (numeric) FDR threshold (display as line), default at 0.05, set to \code{NA} for omitting
#' @param FdrType (character) FDR-type to extract if \code{Mvalue} is 'MArrayLM'-object (eg produced by from \code{moderTest2grp} etc);
#' if \code{NULL} it will search for suitable fields/values in this order : 'FDR','BH','lfdr' and 'BY'
#' @param plotAsFdr (logical) rather plot using FDR-values (if available) instead of p-Values
#' @param subTxt (character) custom sub-title
#' @param grayIncrem (logical) if \code{TRUE}, display overlay of points (not exceeding thresholds) as increased shades of gray
#' @param col (character) custom color(s) for points of plot (see also \code{\link[graphics]{par}})
#' @param pch (integer) type of symbol(s) to plot (default=16) (see also \code{\link[graphics]{par}})
#' @param compNa (character) names of groups compared
#' @param batchFig (logical) if \code{TRUE} figure title and axes legends will be kept shorter for display on fewer splace
#' @param cexMa (numeric) font-size of title, as expansion factor (see also \code{cex} in \code{\link[graphics]{par}})
#' @param cexLa (numeric) size of axis-labels, as expansion factor (see also \code{cex} in \code{\link[graphics]{par}})
#' @param limM (numeric, length=2) range of axis M-values
#' @param limp (numeric, length=2) range of axis FDR / p-values
#' @param annotColumn (character) column names of annotation to be extracted (only if \code{Mvalue} is \code{MArrayLM}-object containing matrix $annot).
#' The first entry (typically 'SpecType') is used for different symbols in figure, the second (typically 'GeneName') is used as prefered text for annotating the best points (if \code{namesNBest} allows to do so.)
#' @param annColor (character or integer) colors for specific groups of annotation (only if \code{Mvalue} is \code{MArrayLM}-object containing matrix $annot)
#' @param expFCarrow (logical, character or numeric) optional adding arrow for expected fold-change; if \code{TRUE} the expected ratio will be extracted from numeric concentration-indications from sample-names
#' if \code{numeric} an arrow will be drawn (M-value as 1st position, color of 2nd position of vector).
#' @param cexPt (numeric) size of points, as expansion factor (see also \code{cex} in \code{\link[graphics]{par}})
#' @param cexSub (numeric) size of subtitle, as expansion factor (see also \code{cex} in \code{\link[graphics]{par}})
#' @param cexTxLab (numeric) size of text-labels for points, as expansion factor (see also \code{cex} in \code{\link[graphics]{par}})
#' @param namesNBest (integer or character) for display of labels to points in figure: if 'pass','passThr' or 'signif' all points passing thresholds; if numeric (length=1) this number of best points will get labels
#' if the initial object \code{Mvalue} contains a list-element called 'annot' the second of the column specified in argument \code{annotColumn} will be used as text
#' @param NbestCol (character or integer) colors for text-labels of best points, also used for arrow
#' @param colBySpecType (logical) incase arument \code{Mvalue} is MArrayLM-object it is possible to use different color-codes for points passing thresholds based on Mvalue$annot[,"SpecType"]; use this with multi-species benchmark tests
#' @param sortLeg (character) sorting of 'SpecType' annotation either ascending ('ascend') or descending ('descend'), no sorting if \code{NULL}
#' @param NaSpecTypeAsContam (logical) consider lines/proteins with \code{NA} in Mvalue$annot[,"SpecType"] as contaminants (if a 'SpecType' for contaminants already exits)
#' @param useMar (numeric, length=4) custom margings (see also \code{\link[graphics]{par}})
#' @param returnData (logical) optional returning data.frame with (ID, Mvalue, pValue, FDRvalue, passFilt)
#' @param silent (logical) suppress messages
#' @param callFrom (character) allow easier tracking of messages produced
#' @param debug (logical) additional messages for debugging
#' @return This function simply plots an MA-plot (to the current graphical device), if \code{returnData=TRUE} an optional data.frame with (ID, Mvalue, pValue, FDRvalue, passFilt) can be returned
#' @seealso (for PCA) \code{\link{plotPCAw}})
#' @examples
#' library(wrMisc)
#' set.seed(2005); mat <- matrix(round(runif(900),2), ncol=9)
#' rownames(mat) <- paste0(rep(letters[1:25], each=4), rep(letters[2:26],4))
#' mat[1:50,4:6] <- mat[1:50,4:6] + rep(c(-1,1)*0.1,25)
#' mat[3:7,4:9] <- mat[3:7,4:9] + 0.7
#' mat[11:15,1:6] <- mat[11:15,1:6] - 0.7
#' ## assume 2 groups with 3 samples each
#' gr3 <- gl(3, 3, labels=c("C","A","B"))
#' tRes2 <- moderTest2grp(mat[,1:6], gl(2,3), addResults = c("FDR","means"))
#' # Note: due to the small number of lines only FDR chosen to calculate
#' VolcanoPlotW(tRes2)
#' ## Add names of points passing custom filters
#' VolcanoPlotW(tRes2, FCth=1.3, FdrThrs=0.2, namesNBest="passThr")
#'
#' ## assume 3 groups with 3 samples each
#' tRes <- moderTestXgrp(mat, gr3, addResults = c("FDR","means"))
#' # Note: due to the small number of lines only FDR chosen to calculate
#' VolcanoPlotW(tRes)
#' VolcanoPlotW(tRes, FCth=1.3, FdrThrs=0.2)
#' VolcanoPlotW(tRes, FCth=1.3, FdrThrs=0.2, useComp=2)
#'
#' @export
VolcanoPlotW <- function(Mvalue, pValue=NULL, useComp=1, filtFin=NULL, tit=NULL, ProjNa=NULL, FCthrs=NULL, FdrList=NULL, FdrThrs=NULL, FdrType=NULL, plotAsFdr=FALSE,
subTxt=NULL, grayIncrem=TRUE, col=NULL, pch=16, compNa=NULL, batchFig=FALSE, cexMa=1.8, cexLa=1.1, limM=NULL, limp=NULL,
annotColumn=c("SpecType","GeneName","EntryName","Accession","Species","Contam"), annColor=NULL, expFCarrow=FALSE,cexPt=NULL, cexSub=NULL,
cexTxLab=0.7, namesNBest=NULL, NbestCol=1, colBySpecType=FALSE, sortLeg="descend", NaSpecTypeAsContam=TRUE, useMar=c(6.2,4,4,2), returnData=FALSE, callFrom=NULL, silent=FALSE, debug=FALSE) {
## MA plot
## optional arguments for explicit title in batch-mode
fxNa <- wrMisc::.composeCallName(callFrom, newNa="VolcanoPlotW")
opar <- graphics::par(no.readonly=TRUE)
opar2 <- opar[-which(names(opar) %in% c("fig","fin","font","mfcol","mfg","mfrow","oma","omd","omi"))] #
on.exit(graphics::par(opar2)) # progression ok
## during function changes in $mar,$cex.main,$cex.lab,$las
namesIn <- c(deparse(substitute(Mvalue)), deparse(substitute(pValue)), deparse(substitute(filtFin)))
basRGB <- c(0.3,0.3,0.3) # grey
fcRGB <- c(1,0,0) # red for points passing FC filt line
multiComp <- TRUE # initialize
inversedComp <- singleCompSetup <- FALSE # inversedComp: if useComp is given as text it may be inverse of p-values given; singleCompSetup : cases sof single comparison BUT 2 cols of p.values
splNa <- annot <- ptType <- colPass <- ptBg <- grpMeans <- pcol <- pwComb <- pwNames <- FDRvalue <- FdrList <- FDRty <- useCompNa <- NULL # initialize
.extrSep2 <- function(comb , indiv) { ## extract separator by striping indiv elements (iniv) from combined pairwise names (comb) .. move to wrMisc ??
fx2 <- function(x, y) sub(paste0(wrMisc::protectSpecChar(x),"$"),"", sub(paste0("^",wrMisc::protectSpecChar(x)),"", y)) # remove x from head and tail (eg rm C from end)
for(i in indiv) comb <- fx2(i, comb)
unique(comb[which(nchar(comb) >0)]) # will/might return NULL when sep==""
}
if(length(pch) <1) pch <- 16
if(isTRUE(debug)) silent <- FALSE else debug <- FALSE
if(!isTRUE(silent)) silent <- FALSE
if(debug) message(fxNa,"Length Mvalue ",length(Mvalue)," ; length pValue ",length(pValue)," ; useComp ",useComp)
if(identical(col,"FDR")) {FDR4color <- TRUE; col <- NULL} else FDR4color <- FALSE
if(length(Mvalue) <1) message("Nothing to do, 'Mvalue' seems to be empty !") else {
## data seem valid to make MAplot
if(length(cexTxLab) <0) cexTxLab <- 0.7
if("MArrayLM" %in% class(Mvalue) || "list" %in% class(Mvalue)) {
if(debug) {message(fxNa,"VPW0"); VPW0 <- list(Mvalue=Mvalue,pValue=pValue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,fxNa=fxNa )}
## try working based on MArrayLM-object (Mvalue)
if(debug) message(fxNa," '",namesIn[1],"' is list or MArrayLM-object ")
## initial check of useComp
if(length(useComp) >1) { useComp <- wrMisc::naOmit(useComp)[1]
if(!silent) message(fxNa,"Argument 'useComp' should be integer of length=1; using only 1st entry") }
if(length(useComp) <1) { useComp <- 1
if(!silent) message(fxNa,"Argument 'useComp' invalid with this dataset, setting to 1") }
## address multiple questions (via useComp): need to check 'useComp', thus need to locate FDRvalues or pValues for group-names ..
## extract data: find suitable columns
pcol <- wrMisc::naOmit(match(c("p.value","pvalue","pval","p"), tolower(names(Mvalue))))
if(length(pcol) >0) names(pcol) <- c("p.value","pvalue","pval","p")[pcol]
## NEED WAY TO CUSTOM CHOOSE WHICH FDR TO USE !!
FDRcol <- wrMisc::naOmit(match(c("fdr","bh","lfdr","by","bonferroni"), tolower(names(Mvalue))))
if(length(pcol) >0) names(pcol) <- names(Mvalue)[pcol]
if(length(FDRcol) >0) names(FDRcol) <- names(Mvalue)[FDRcol]
#if(length(FDRcol) >0) names(FDRcol) <- c("FDR","BH","lfdr","BY","bonferroni")[FDRcol]
if(debug) {message(fxNa,"VPW1"); VPW1 <- list(Mvalue=Mvalue,pValue=pValue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol )}
if(length(pcol) ==0) stop(fxNa,"Cannot find suitable element for p-values") else pcol <- pcol[1]
if(length(FDRcol) ==0) stop(fxNa,"Cannot find suitable element for FDR-values") else FDRcol <- FDRcol[1]
if("means" %in% names(Mvalue)) grpMeans <- Mvalue$means else stop(fxNa,"Need element 'means' (group-means) in object 'Mvalue'")
## need to know which of pairwise comparisons get addressed: need to understand setup -either from character useComp or index & p-values
if(debug) {message(fxNa,"VPW1b"); VPW1b <- list(Mvalue=Mvalue,pValue=pValue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol,grpMeans=grpMeans )}
if("design" %in% names(Mvalue)) {
if(length(Mvalue$design)==0 || length(dim(Mvalue$design)) !=2 || nrow(Mvalue$design) ==0) stop(fxNa,"Unable to locate element with names of pairwise groups ! (ie $design)")
## note $contrasts NOT present in result from wrMisc::moderTest2grp() !!
if(colnames(Mvalue$design)[1] =="(Intercept)" && ncol(Mvalue$design) ==2) {
## single comparison case (colnames of Mvalue$design are fixed); determine optimal sep and create pairwise name(s)
singleCompSetup <- TRUE
#Mvalue$design <- Mvalue$design[,2, drop=FALSE]
grp <- colnames(Mvalue$means) #
if(utils::packageVersion("wrMisc") < "2.0.0") {
pwSep <- " "
} else {
pwSep <- wrMisc::getPWseparator(grp=colnames(Mvalue$design), includeGrp=FALSE, silent=silent, debug=debug, callFrom=fxNa)
colnames(Mvalue$design) <- grp[1:ncol(Mvalue$design)]
}
pwNames <- if(length(Mvalue$setup$pwNames)==1) Mvalue$setup$pwNames else {if(length(grp)==2) as.matrix(grp) else utils::combn(grp, 2)}
allCompNa <- if(length(Mvalue$setup$allCompNa) >0) Mvalue$setup$allCompNa else {if(length(grp)==2) paste(grp, collapse=pwSep) else paste(pwNames[1,], pwNames[2,], sep=pwSep)}
for(i in c(pcol,FDRcol)) {
if(colnames(Mvalue[[i]])[1] =="(Intercept)") { Mvalue[[i]] <- Mvalue[[i]][,2, drop=FALSE] # remove '(Intercept)' column
colnames( Mvalue[[i]]) <- allCompNa[1]
}
}
} else {
## regular multi-comparison case
grp <- unique(colnames(Mvalue$design))
if(utils::packageVersion("wrMisc") < "2.0.0") {
pwSep <- " "
} else {
pwSep <- if(length(Mvalue$setup$sep)==1) Mvalue$setup$sep else wrMisc::getPWseparator(colnames(Mvalue[[c(pcol,FDRcol)[1]]]), grp=colnames(Mvalue$design), includeGrp=FALSE, silent=silent, debug=debug, callFrom=fxNa)
}
allCompNa <- if(length(Mvalue$setup$allCompNa) >0) Mvalue$setup$allCompNa else colnames(Mvalue[[c(pcol,FDRcol)[1] ]])
pwGrInd <- t(wrMisc::matchSampToPairw(grpNa=grp, pairwNa=allCompNa, sep=pwSep, silent=FALSE, debug=FALSE, callFrom=fxNa))
pwNames <- matrix(grp[pwGrInd], ncol=length(grp))
}
## created in this block : pwSep, pwNames
## investigate pairwise names from testing (& get separator) - need to find out which element of Mvalue
if(length(allCompNa) ==0 ) stop(fxNa,"Unable to locate element with names of pairwise groups !")
if(debug) {message(fxNa,"VPW1c"); VPW1c <- list(Mvalue=Mvalue,pValue=pValue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol,grpMeans=grpMeans,grp=grp,pwSep=pwSep,pwNames=pwNames,allCompNa=allCompNa,useComp=useComp,useCompNa=useCompNa )}
if(is.numeric(useComp)) { ## simple matching on allCompNa
if(singleCompSetup && useComp != 1) {useComp <- 1; names(useComp) <- paste(pwNames[,1], sep=pwSep)
if(!silent) message(fxNa,"Adjusting useComp to single comparison setup") }
names(useComp) <- allCompNa[useComp]
useCompNa2 <- pwNames[,useComp]
if(debug) message(fxNa," (useComp numeric) useCompNa2 ",paste(useCompNa2,collapse=" & "))
} else {
## check if name of useComp matches any of (default) pairwise
chComp <- allCompNa %in% names(useComp)
if(any(chComp, na.rm=TRUE)) {
useCompNa2 <- pwNames[,which(chComp)]
useCompNa <- useComp
useComp <- which(chComp)
names(useComp) <- useCompNa
}
if(any(chComp, na.rm=TRUE)) {
useCompNa2 <- pwNames[,which(chComp)]
useCompNa <- useComp
useComp <- which(chComp)
names(useComp) <- useCompNa
} else {
## possible reversed written - split (text-)useComp
useCompIn2 <- try(wrMisc::matchSampToPairw(grpNa=grp, pairwNa=useComp, sep=pwSep, silent=FALSE, debug=FALSE, callFrom=fxNa), silent=TRUE)
if(inherits(useCompIn2, "try-error")) stop(fxNa,"Unable to understand argument useComp ('",useComp,"') - can't match to any of the groups of replicates")
useCompNa2 <- grp[useCompIn2]
chComp <- which(colSums(matrix(pwNames %in% useCompNa2, ncol=length(grp)))==2)
if(length(chComp) !=1) stop(fxNa,"Unable to understand argument useComp ('",useComp,"') - can't match to any of pairwise groups")
useCompNa <- useComp
useComp <- chComp
names(useComp) <- useCompNa
inversedComp <- TRUE
}
if(debug) message(fxNa," (useComp char) useCompNa2 ",paste(useCompNa2,collapse=" & "))
}
## created in this block : update useComp, useCompNa2 (separate names)
if(debug) {message(fxNa,"v"); VPW1d <- list(Mvalue=Mvalue,pValue=pValue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol,grpMeans=grpMeans,pwSep=pwSep,annotColumn=annotColumn,useCompNa2=useCompNa2 )}
Mvalue$Mval <- grpMeans[,useCompNa2[1]] - grpMeans[,useCompNa2[2]]
if(inversedComp) Mvalue$Mval <- -1 * Mvalue$Mval
pValue <- if(length(pcol)==1) Mvalue[[pcol]][,useComp] else NULL
if(debug) {message(fxNa,"VPW2"); VPW2 <- list(Mvalue=Mvalue,pValue=pValue,FDRvalue=FDRvalue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol,grpMeans=grpMeans,pwSep=pwSep,useCompNa2=useCompNa2,annotColumn=annotColumn )}
} else stop(fxNa,"Argument 'Mvalue' is missing key element $design")
## recuperate FDRvalue
if(length(FDRcol)==1 && names(FDRcol) %in% names(Mvalue)) {
FDRvalue <- Mvalue[[FDRcol]][,useComp]
}
## recuperate $annot if present and use for symbol
if("annot" %in% names(Mvalue) && isTRUE(colBySpecType)) {
useAnnCol <- match(annotColumn, colnames(Mvalue$annot))
if(!is.na(useAnnCol[1])) { # annotation (for multiple groups) exists
ptType <- Mvalue$annot[,useAnnCol[1]] # SpecType
chNA <- is.na(ptType)
## associate NAs from 'SpecType' in ptType with conta ?
if(NaSpecTypeAsContam) {
chConta <- tolower(ptType) %in% c("contaminant","contam","conta","cont")
if(any(chConta)) ptType[which(is.na(ptType))] <- unique(ptType[which(chConta)])[1]}
if(any(is.na(ptType))) ptType[which(chNA)] <- "NA"
if(length(pch) < length(pValue) && length(unique(wrMisc::naOmit(ptType))) >1) {
if(length(pch) >1 && !silent) message(fxNa," (invalid pch) using default 'pch' oriented by $annot and starting from 15")
pch <- 14 + as.integer(as.factor(ptType))
}
useAnnCol <- wrMisc::naOmit(useAnnCol) }
annot <- Mvalue$annot[,useAnnCol]
if(annotColumn[1] %in% colnames(annot)) annot[,annotColumn[1]] <- ptType # update with NAs trasformed to "NA"
if(debug) {message(fxNa,"VPW2b"); VPW2b <- list(Mvalue=Mvalue,pValue=pValue,FDRvalue=FDRvalue,annot=annot,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol,grpMeans=grpMeans,pwSep=pwSep,useCompNa2=useCompNa2,annotColumn=annotColumn,useCompNa=useCompNa,pch=pch )}
}
## filtering
if("filter" %in% names(Mvalue) && ((length(filtFin)==1 && isTRUE(filtFin)) || length(filtFin)==0)) {
filtFin <- Mvalue$filter[,useComp]
#filtFin <- Mvaue$filter[,which(colnames(Mvalue$filter)==useCompNa)]
} else { ## check custom filtFin
chFi <- is.logical(filtFin) && length(filtFin)==length(pValue)
if(!chFi) { filtFin <- NULL
if(!silent) message(fxNa,"Argument 'filtFin' seems invalid (should be logical vector of length of number of elements for plotting )")
} else if(debug) message(fxNa,"Using custom filtFin")
}
if(length(filtFin) >0) {
if(debug) message(fxNa,"Filtering : Setting ",sum(!filtFin & !is.na(pValue))," additional values as NA")
pValue[which(!filtFin)] <- NA # better to set NA so that annotation, filtering etc stays same
if(length(FDRvalue) ==length(pValue)) FDRvalue[which(!filtFin)] <- NA
## Mvalue : if only 2 groups $Mval should be OK
if("Mval" %in% names(Mvalue)) Mvalue$Mval[which(!filtFin)] <- NA
}
if(debug) {message(fxNa,"VPW2c"); VPW2c <- list(Mvalue=Mvalue,pValue=pValue,FDRvalue=FDRvalue,useComp=useComp,filtFin=filtFin,ProjNa=ProjNa,FCthrs=FCthrs,pcol=pcol,FDRcol=FDRcol,grpMeans=grpMeans,pwSep=pwSep,useCompNa2=useCompNa2,annotColumn=annotColumn )}
Mvalue <- Mvalue$Mval # DISMISS rest of object !!! .....
} else {
## ... case of explicit pValue argument
if(length(pValue) <1) stop("Argument 'pValue' is required (if 'Mvalue' not 'MArrayLM'-type object) !")
if(length(dim(pValue)) >1) if(ncol(pValue) >1) {
if(!silent) message(fxNa," Note, ",namesIn[2]," has ",ncol(pValue)," columns, using last column")
pNa <- rownames(pValue)
pValue <- as.numeric(pValue[,ncol(pValue)] )
names(pValue) <- pNa }
FDRvalue <- if(length(FdrList) <1) NULL else FdrList
}
if(debug) {message(fxNa,"VPW8"); VPW8 <- list(FdrList=FdrList,FDRvalue=FDRvalue, Mvalue=Mvalue,pValue=pValue,FdrType=FdrType,plotAsFdr=plotAsFdr,pValue=pValue,useComp=useComp,grpMeans=grpMeans)}
## need to introduce -log10 to pValue
chNA <- is.na(pValue)
if(all(chNA)) stop(fxNa,"All p-values are NA, nothing to draw !")
if(debug && any(pValue <0, na.rm=TRUE)) message(fxNa,"Some p-values are negative, this should not be ! Maybe log values were given by error ?")
pValue <- -log10(pValue)
## check for (same) order, adjust Mvalue & pValue according to names
chNa <- list(MNa=if(length(dim(Mvalue)) >1) rownames(Mvalue) else names(Mvalue),
pNa=if(length(dim(pValue)) >1) rownames(pValue) else names(pValue))
nIni <- c(M=length(Mvalue), p=length(pValue))
if(length(chNa$MNa) >0 && length(chNa$pNa) >0) { # ie both have names, so one can match names
if(!all(chNa$MNa==chNa$pNa, na.rm=TRUE)) {
matchNa <- wrMisc::naOmit(match(chNa$MNa, chNa$pNa))
if(length(matchNa) <1) stop("Both 'Mvalue' and 'pValue' have names, but none of them match !!")
if(!all(matchNa, 1:length(pValue), na.rm=TRUE)) {
pValue <- pValue[matchNa]
Mvalue <- wrMisc::naOmit(Mvalue[match(names(pValue), names(Mvalue))]) }
}
} else {
if(length(Mvalue) != length(pValue)) stop("p- and M- values have different length, but no names to match !! (M=",length(Mvalue)," vs p=",length(pValue),")")
}
if(length(grpMeans) <1) grpMeans <- matrix(rep(NA,2*length(Mvalue)), ncol=2, dimnames=list(names(Mvalue),c("mean1","mean2")))
if(debug) {message(fxNa,"VPW9"); VPW9 <- list(FdrList=FdrList,Mvalue=Mvalue,FdrType=FdrType,FDRvalue=FDRvalue,plotAsFdr=plotAsFdr,pValue=pValue,useComp=useComp,grpMeans=grpMeans, filtFin=filtFin,annotColumn=annotColumn,annot=annot, pch=pch,sortLeg=sortLeg,batchFig=batchFig)} # pairwCol=pairwCol,
## prepare/integrate FILTERING
if(length(filtFin) >0) {
## if filtFin is matrix use each line with min 1 instance of TRUE,
if(length(dim(filtFin)) >1) filtFin <- as.logical(as.matrix(filtFin)[,useComp]) # use rows with >= 1 TRUE
if(length(names(filtFin)) >0) {
matchNa <- wrMisc::naOmit(match(rownames(pValue), names(filtFin)))
if(length(matchNa)==length(pValue)) filtFin <- as.logical(filtFin[matchNa])
}
} else filtFin <- rep(TRUE, nrow(grpMeans))
if(debug) {message(fxNa,"VPW9b"); VPW9b <- list(FdrList=FdrList,Mvalue=Mvalue,FdrType=FdrType,FDRvalue=FDRvalue,plotAsFdr=plotAsFdr,pValue=pValue,useComp=useComp,grpMeans=grpMeans, filtFin=filtFin,annotColumn=annotColumn,annot=annot, pch=pch,sortLeg=sortLeg,batchFig=batchFig)} # pairwCol=pairwCol,
## start creating merged data for plot (& export)
merg <- if(length(annot) >0) data.frame(ID=NA, grpMeans, Mvalue=Mvalue, pValue=pValue,
FDR=if(length(FDRvalue) >0) as.numeric(FDRvalue) else rep(NA,length(pValue)),
filtFin=filtFin, annot, pch=pch, stringsAsFactors=FALSE) else {
data.frame(ID=NA, grpMeans, Mvalue=Mvalue, pValue=pValue, FDR=if(length(FDRvalue) >0) as.numeric(FDRvalue) else rep(NA,length(pValue)),
filtFin=filtFin, pch=pch, stringsAsFactors=FALSE) }
if(length(names(Mvalue)) >0) merg[,1] <- names(Mvalue) else {if(length(names(pValue)) >0) merg[,1] <- names(pValue)}
## replace NA in 'SpecType' by 'NA'
if(annotColumn[1] %in% colnames(merg)) { chNa <- is.na(merg[,annotColumn[1]]) # replace NAs in col "SpecType" by "NA"
if(any(chNa)) merg[which(chNa),annotColumn[1]] <- "NA"
} else { merg <- cbind(merg, rep(1,nrow(merg))) # add colum for 'SpecType'
colnames(merg)[ncol(merg)] <- annotColumn[1] }
if(debug) {message(fxNa,"VPW10"); VPW10 <- list(merg=merg,FdrList=FdrList,Mvalue=Mvalue,FdrType=FdrType,FDRvalue=FDRvalue,plotAsFdr=plotAsFdr,pValue=pValue,useComp=useComp,grpMeans=grpMeans, filtFin=filtFin,annotColumn=annotColumn,annot=annot, pch=pch,sortLeg=sortLeg,batchFig=batchFig)}
#check#plot(pValue ~Mvalue, data=VPW10$merg) # ok
## adjust col & pch
if(!any(c(1,length(Mvalue)) %in% length(pch))) {
if(!silent) message(fxNa,"Argument 'pch' should be either length=1 or correspond to length of data, reset to default=16")
pch <- 16 }
if(length(col) >1 && length(col) <length(Mvalue)) {
if(!silent) message(fxNa,"Argument 'col' should be either length=1 or correspond to length of data, reset to default=NULL")
col <- NULL }
if(debug) message(fxNa," ++ DONE extracting columns : ",wrMisc::pasteC(colnames(merg),quo="'"))
## apply filtering
msg <- " data provided in 'Mvalue' and 'pValue' "
if(debug) {message(fxNa,"VPW10")}
if(!silent && nrow(merg) < round(length(Mvalue)/10)) message(fxNa," .. note : less than 10% of",msg," were matched") else {
if(!silent && nrow(merg) < length(Mvalue)/2) message(fxNa," .. NOTE : less than 50% of",msg," were matched !!")}
if(debug) message(fxNa,msg," were matched to ",nrow(merg)," common entries")
## apply filtering (keep all lines where at least one condition passes)
if(length(filtFin) >0 && any(isFALSE(filtFin), na.rm=TRUE)) { # use filtering provided
if(sum(filtFin) >0 && sum(filtFin) < nrow(merg)) {
whFilt <- which(merg$filtFin)
if(length(pch) >1) pch <- pch[whFilt]
if(length(col) >1) col <- col[whFilt]
merg <- merg[whFilt,]
if(!silent) message(fxNa,"Filtered (based on 'filtFin') from ",length(filtFin)," to ",nrow(merg)," lines")
}
} else filtFin <- rep(TRUE, nrow(merg))
if(debug) {message(fxNa,"VPW11"); VPW11 <- list(merg=merg,FdrList=FdrList,Mvalue=Mvalue,FdrType=FdrType,FDRvalue=FDRvalue,plotAsFdr=plotAsFdr,pValue=pValue,useComp=useComp,grpMeans=grpMeans, filtFin=filtFin,annotColumn=annotColumn,annot=annot, pch=pch,sortLeg=sortLeg,batchFig=batchFig)}
#check#plot(pValue ~Mvalue, data=VPW11$merg) # ok
## sort merg, so that legend always gets constructed the same order, ascending ('ascend') or descending ('descend')
sortLeg <- if(identical(sortLeg,"ascend")) FALSE else {if(identical(sortLeg,"descend")) TRUE else NULL}
if(length(sortLeg) ==1 && annotColumn[1] %in% colnames(merg)) merg <- merg[order(merg[,annotColumn[1]], decreasing=sortLeg),]
## update ..
nIDco <- sum(c("ID","nredID","uniqID") %in% colnames(merg)) # number of heading columns in 'merg'
Mvalue <- as.numeric(if("Mvalue" %in% colnames(merg)) merg[,"Mvalue"] else merg[,nIDco+1])
pValue <- as.numeric(if("pValue" %in% colnames(merg)) merg[,"pValue"] else {
if(length(dim(Mvalue)) >0) merg[,ncol(Mvalue) +nIDco +1] else merg[,nIDco+2]})
FDRvalue <- as.numeric(if("FDR" %in% colnames(merg)) merg[,"FDR"] else merg[,grep("(fdr)|(bh)", tolower(colnames(merg)))[1]])
if("Lfdr" %in% colnames(merg)) FdrList <- merg[,"Lfdr"] else {
if("lfdr" %in% colnames(merg)) FdrList <- merg[,"lfdr"]}
pch <- merg[,"pch"] # update
ptType <- if(annotColumn[1] %in% colnames(merg)) merg[,annotColumn[1]] else rep(1,nrow(merg)) # update "SpecType"
if(debug) {message(fxNa,"VPW12"); VPW12 <- list(merg=merg,Mvalue=Mvalue,filtFin=filtFin)}
## prepare for plotting
if(is.null(cexSub)) cexSub <- cexLa +0.05
xLab <- "M-value (log2 fold-change)"
tit1 <- if(length(tit)==1 && nchar(tit) >0) tit else paste(c(if(!batchFig) c(ProjNa, if(!is.null(ProjNa)) ": ","Volcano-Plot"),
if(!is.null(compNa)) {if(length(compNa)==1) compNa else c(compNa[1]," vs ",compNa[2])}), collapse=" ") # but what title if batchFig=NULL & compNa=NULL -> only "Volcano-plot"
if(length(FCthrs) <1) FCthrs <- 1.5
if(length(FdrThrs) <1) FdrThrs <- 0.05
if(debug) {message(fxNa,"VPW12b"); VPW12b <- list(merg=merg,FdrList=FdrList,Mvalue=Mvalue,FdrType=FdrType,FDRvalue=FDRvalue,plotAsFdr=plotAsFdr,pValue=pValue,useComp=useComp,grpMeans=grpMeans, filtFin=filtFin,annotColumn=annotColumn,annot=annot, pch=pch,sortLeg=sortLeg,batchFig=batchFig,FCthrs=FCthrs,FdrThrs=FdrThrs)}
## count no of passing
passFC <- if(length(FCthrs) ==1 && !any(is.na(FCthrs))) abs(merg[,"Mvalue"]) >= log2(FCthrs) else merg[,"filtFin"] ## convert FCthrs to log2
passFdr <- if(length(FdrThrs) ==1 && !any(is.na(FdrThrs)) && "FDR" %in% colnames(merg)) {merg[,"FDR"] <= FdrThrs} else merg[,"filtFin"]
passAll <- merg[,"filtFin"] & passFC & passFdr
if(debug) {message(fxNa,"VPW12c")}
chNA <- is.na(passAll) # passFdr may contain NAs
if(any(chNA)) passAll[which(chNA)] <- FALSE
if(debug) message(fxNa," ",sum(passFC,na.rm=TRUE)," passing FCthrs ; ",sum(passFdr,na.rm=TRUE)," passing FdrThrs ; combined ",sum(passAll,na.rm=TRUE))
if(debug) {message(fxNa,"VPW13"); VPW13 <- list(merg=merg,Mvalue=Mvalue,filtFin=filtFin,col=col,passAll=passAll,cexPt=cexPt,annotColumn=annotColumn,annColor=annColor)}
## color for points passing filter
if(length(col) >0) if(length(col) != nrow(merg)) { col <- NULL
if(!silent) message(fxNa," invalid entry for 'col', should be of length=",nrow(merg),", resetting to default")}
if(length(col) <1) {
alph <- sort(c(0.14, round(0.6/log10(length(Mvalue)),2), 0.8))[2] # alph <- round(12/sqrt(nrow(eBayesLst$pValue)),2)
alph2 <- sort(c(round(7/(5 +sum(passAll)^0.7),2), alph,0.9))[2] # for points passing thresholds
useCol <- if(grayIncrem) grDevices::rgb(0.35,0.35,0.35,alph) else grDevices::rgb(0.7,0.7,0.7) # basic color
useCex <- if(length(cexPt) >0) cexPt else max(round(0.8 + 2/(1 +sum(filtFin, na.rm=TRUE))^0.28,2), 1.1)
chCol <- unique(merg[, annotColumn[1]]) # check how many different colors may be needed
chNaC <- is.na(chCol)
if(any(chNaC)) chCol[which(chNaC)] <- "NA"
if(length(annColor) >0) {colPass <- annColor} else if(length(chCol) >4) {
colPass <- cbind(red=c(141,72,90,171, 220,253,244,255), green=c(129,153,194,221, 216,174,109,0), blue=c(194,203,185,164, 83,97,67,0))
colPass <- grDevices::rgb(red=colPass[,1], green=colPass[,2], blue=colPass[,3], alph2, maxColorValue=255)
if(length(chCol) >8) { colPass <- c(colPass, rep(colPass[8], length(chCol) -8))
if(!silent) message(fxNa," > 8 different groups found, using 8th color after 7th group")}
} else colPass <- grDevices::rgb(c(0.95,0.2,0,0.75), c(0.15,0.2,0.9,0.35), c(0.15,0.95,0,0.8), alph2) # red, blue, green, purple (luminosity adjusted)
useCol <- rep(useCol[1], nrow(merg)) # fuse basic gray to colors for different types
## integrate names of annColor as order of colPass
if(length(names(annColor)) >0) {
uniTy <- unique(merg[which(passAll), annotColumn[1]])
colPass <- colPass[match(names(annColor), uniTy)]
}
if(debug) {message(fxNa,"VPW13a"); VPW13a <- list(useCol=useCol,colPass=colPass, merg=merg,Mvalue=Mvalue,filtFin=filtFin,col=col,passAll=passAll,cexPt=cexPt,annotColumn=annotColumn,annColor=annColor)}
## assign color for those passing
useCol[which(passAll)] <- if(any(passAll, na.rm=TRUE)) colPass[if(length(unique(merg[which(passAll),annotColumn[1]])) >1) wrMisc::levIndex(merg[which(passAll),annotColumn[1]]) else rep(1,sum(passAll))] # assign colors for those passing
} else useCol <- col
if(debug) {message(fxNa,"VPW13b"); VPW13b <- list(merg=merg,Mvalue=Mvalue,filtFin=filtFin,col=col,passAll=passAll,cexPt=cexPt,annotColumn=annotColumn,annColor=annColor,pch=pch,useCol=useCol)}
#}} #ok
## adjust fill color for open symbols
chPch <- pch %in% c(21:25)
if(any(chPch, na.rm=TRUE)) { ptBg <- useCol
ptBg[which(chPch)] <- useCol[which(chPch)] # background color for filled symbols
useCol[which(chPch)] <- 1 # contour as black
}
## main graphics
if(length(Mvalue) >0) {
pl1 <- try(graphics::par(mar=if(length(useMar)==4) useMar else c(6.5,4,4,2), cex.main=cexMa, las=1), silent=TRUE)
if(inherits(pl1, "try-error")) {Mvalue <- NULL; message("UNABLE TO SET PLOT MARGINS !! check plotting device ...")}
if(debug) {message(fxNa,"VPW14a"); VPW14a <- list(pl1=pl1,merg=merg,Mvalue=Mvalue,filtFin=filtFin,useMar=useMar,FCthrs=FCthrs,FdrThrs=FdrThrs,FdrType=FdrType,plotAsFdr=plotAsFdr)}
}
if(length(Mvalue) >0) {
## rather directly plot FDR
if(!"FDR" %in% colnames(merg)) plotAsFdr <- FALSE
#check# plot(pValue ~Mvalue, data=VPW14a$merg)
pl1 <- try(graphics::plot(Mvalue, if(plotAsFdr) log10(merg[,"FDR"]) -min(log10(merg[,"FDR"]), na.rm=TRUE) else merg[,"pValue"], pch=pch, cex=useCex, main=tit1,
ylab=if(plotAsFdr) paste("- log10",if(is.null(FdrType)) "(probability)" else FdrType) else "- log10 p-value (uncorrected)", col=useCol, xlab=xLab, cex.lab=cexLa, xlim=limM,ylim=limp, pt.bg=ptBg), silent=TRUE)
if(inherits(pl1, "try-error")) {Mvalue <- NULL; message("UNABLE TO PLOT !! check plotting device ...")} }
if(debug) {message(fxNa,"VPW14b"); VPW14b <- list(merg=merg,Mvalue=Mvalue,filtFin=filtFin,passAll=passAll,useMar=useMar,FCthrs=FCthrs,FdrThrs=FdrThrs,useCol=useCol,multiComp=multiComp,useComp=useComp,useCompNa=useCompNa,subTxt=subTxt,plotAsFdr=plotAsFdr )}
if(length(Mvalue) >0) { # useCompNa
sTxt <- if(length(subTxt) ==1) subTxt else { if(multiComp) if(length(names(useComp)) >0) paste0(names(useComp),"; ",collapse="") else paste0("comp=",useComp,"; ")}
sTxt <- paste0(sTxt,"n=",length(Mvalue),
if(!all(is.na(c(FCthrs,FdrThrs)))) paste(";",sum(passAll, na.rm=TRUE),"(color) points passing",
if(!is.na(FCthrs)) paste0("(FCthr=", as.character(FCthrs),", ") else " (", paste0("FdrThrs=",as.character(FdrThrs),")")))
graphics::mtext(sTxt,cex=0.75,line=0.2)
if(!all(is.na(c(FCthrs,FdrThrs)))) {
if(debug) message(fxNa," n=",length(Mvalue)," FCthrs=",as.character(FCthrs)," filt.ini=", sum(filtFin, na.rm=TRUE),
" passAll=",sum(passAll,na.rm=TRUE)," ; range Mva ",wrMisc::pasteC(signif(range(Mvalue,na.rm=TRUE),3))," ; alph=",alph," useCex=",useCex," alph2=",alph2)
## previous versions : adjusted abline by adding diff(graphics::par("usr")[1:2])/1000)
graphics::abline(v=c(-1,1)* log2(FCthrs) , col=grDevices::rgb(0.87,0.72,0.72), lty=2) }
if(debug) {message(fxNa,"VPW15"); VPW15 <- list(merg=merg,Mvalue=Mvalue,filtFin=filtFin,FCthrs=FCthrs,FdrThrs=FdrThrs,passFdr=passFdr,annotColumn=annotColumn,xLab=xLab,tit1=tit1,col=col,useCol=useCol,assAll=passAll,grayIncrem=grayIncrem,cexPt=cexPt,annColor=annColor,plotAsFdr=plotAsFdr,limM=limM,limp=limp,ptBg=ptBg,cexLa=cexLa,passAll=passAll)}
if(sum(passFdr, na.rm=TRUE) >0) {
if(plotAsFdr) {
if(any(passAll, na.rm=TRUE)) graphics::abline(h=-1*log10(max(merg[passAll,"FDR"], na.rm=TRUE)) -diff(graphics::par("usr")[3:4])/400, col=grDevices::rgb(0.87,0.72,0.72), lty=2)
} else {
if(plotAsFdr) graphics::mtext("Note, that FDR and p-value may not correlate perfectly, thus points may appear at good p-value but finally don't get retained",line=5,cex=0.65, side=1)
pRa <- range(merg[which(passFdr),"pValue"], na.rm=TRUE)
graphics::abline(h=pRa[1] +diff(graphics::par("usr")[3:4])/400, col=grDevices::rgb(0.87,0.72,0.72), lty=2) }}
## add names to best points
if(length(namesNBest) >0) {
if(any(sapply( c("pass","passThr","passFDR","signif"), identical, namesNBest))) namesNBest <- sum(passAll)
if(!is.integer(namesNBest)) namesNBest <- try(as.integer(namesNBest), silent=TRUE)
if(inherits(namesNBest, "try-error")) { namesNBest <- NULL
message(fxNa,"Unable to understand argument 'namesNBest', must be integer or 'pass','passThr','passFDR' or 'signif' (for display of labels to points in figure)") }
}
if(debug) {message(fxNa,"VPW16")}
if(length(namesNBest) >0) {
if(namesNBest >0 && any(passAll)) {
useLi <- if(any(!passAll)) which(passAll) else 1:nrow(merg)
tmP <- as.numeric(merg[useLi,"pValue"])
names(tmP) <- rownames(merg)[useLi]
## look for more informative names to display
if(length(annot) >0) {
proNa <- annot[match(names(tmP), rownames(annot)), annotColumn[2]] # normally 'Description'
chNa <- is.na(proNa)
if(!all(chNa)) names(tmP)[which(!chNa)] <- proNa[which(!chNa)]
}
useL2 <- order(tmP, decreasing=TRUE)[1:min(namesNBest,sum(passAll))]
xOffs <- signif(diff(graphics::par("usr")[1:2])/170,3)
yOffs <- signif(diff(graphics::par("usr")[3:4])/90,3)
noNa <- if(is.null(names(tmP[useL2]))) 1:length(tmP) else which(is.na(names(tmP)[useL2]))
if(length(noNa) >0 && all(annotColumn %in% colnames(merg))) names(tmP)[useL2[noNa]] <- merg[useLi[useL2[noNa]], wrMisc::naOmit(match(annotColumn[-1], colnames(merg)))[1]]
if(length(NbestCol) <1) NbestCol <- 1
if(is.null(names(tmP[useL2]))) {if(!silent) message(fxNa,"No names available for displaying names of best in plot")
} else { displTx <- wrMisc::naOmit(match(annotColumn[2:4], colnames(merg)))
displTx <- if(length(displTx) >0) cbind(merg[useLi[useL2], displTx], rownames(merg)[useLi[useL2]]) else as.matrix(useLi[useL2])
chNa <- rowSums(!is.na(displTx)) >0
if(any(chNa)) displTx[which(chNa),1] <- "unknown"
displTx <- apply(displTx,1, function(x) wrMisc::naOmit(x)[1])
graphics::text(Mvalue[useLi[useL2]] +xOffs, yOffs +if(plotAsFdr) log10(merg[useLi[useL2],"FDR"]) -min(log10(merg[,"FDR"]), na.rm=TRUE) else merg[useLi[useL2],"pValue"],
names(tmP)[useL2], cex=cexTxLab, col=NbestCol, adj=0) }
}
}
if(debug) {message(fxNa,"VPW17"); VPW17 <- list(pch=pch,passAll=passAll,merg=merg,annot=annot,annotColumn=annotColumn)}
## legend (if multiple symbols)
pch[which(is.na(pch))] <- -2
ch1 <- unique(pch)
if(length(ch1) >1 && sum(passAll) >0) {
legInd <- which(!duplicated(merg[which(passAll), annotColumn[1]], fromLast=FALSE))
if(length(legInd) <1 && !silent) message(fxNa,"Trouble ahead : Can't find non-duplicated ",annotColumn[1]," for ",sum(passAll)," points passing thresholds ! (ie as 'legInd')")
legPch <- pch[which(passAll)[legInd]]
legCol <- useCol[which(passAll)[legInd]]
legBg <- ptBg[which(passAll)[legInd]]
if(alph2 <1) {legCol <- substr(legCol,1,7); legBg <- substr(legBg,1,7)} # reset to no transparency
legLab <- merg[which(passAll)[legInd], annotColumn[1]]
chNa <- is.na(legLab)
if(any(chNa)) legLab[chNa] <- "NA"
legOr <- if(length(legLab) >1) order(legLab) else 1 # not used so far
legLoc <- checkForLegLoc(cbind(Mvalue, pValue), sampleGrp=legLab, showLegend=FALSE)
legCex <- stats::median(c(useCex, cexTxLab, 1.2), na.rm=TRUE)
if(length(legLab) >0) { chLeg <- try(graphics::legend(legLoc$loc, legend=legLab, col=legCol, text.col=1, pch=legPch,
if(length(ptBg) >0) pt.bg=ptBg, cex=legCex, pt.cex=1.2*legCex, xjust=0.5, yjust=0.5), silent=TRUE) # as points
if(inherits(chLeg, "try-error") && !silent) message(fxNa,"Note: Failed to add legend .. ",chLeg) }
} else legCol <- NULL
if(debug) {message(fxNa,"VPW18")}
## arrow for expected ratio
## how to extract in smartest way ??
drawArrow <- if(length(expFCarrow) >0) {isTRUE(as.logical(expFCarrow[1])) || grepl("^[[:digit:]]+$", expFCarrow[1])} else FALSE
if(drawArrow) {
regStr <-"[[:space:]]*[[:alpha:]]+[[:punct:]]*[[:alpha:]]*"
if(isTRUE(expFCarrow)) {
## automatic extraction of FC
if(debug) message(fxNa," .. ",names(useComp)," -> ",paste(unlist(strsplit(names(useComp), "-")), collapse=" "),"\n")
expM <- sub(paste0("^",regStr),"", sub(paste0(regStr,"$"), "", unlist(strsplit(names(useComp), "-")))) # assume '-' separator from pairwise comparison
arrCol <- 1
chN2 <- try(as.numeric(expM), silent=TRUE)
if(inherits(chN2, "try-error")) { drawArrow <- FALSE #expM <- arrCol <- NA
} else {
expM <- log2(chN2[2] / chN2[1]); arrCol <- 1 # transform automatic extracted to ratio # expFCarrow <- c(expM, 1)
}
msg <- "Argument 'expFCarrow' was set to TRUE; "
if(!silent) message(fxNa, msg, if(drawArrow) "Extract concentration values of group-names for calculating ratios" else "Failed automatic extraction of FC-values (possibly no numeric values in setup)")
} else {
expM <- try(as.numeric(expFCarrow[1]), silent=TRUE)
if(inherits(expM, "try-error")) {
if(!silent) message(fxNa,"Argument 'expFCarrow' was given, but no numeric content in 1st place found, can't draw FC-arrow") #understand 1st value of'expFCarrow' (should be numeric-like)
drawArrow <- FALSE
} else arrCol <- if(length(expFCarrow) >1) expFCarrow[2] else 1
}
if(length(expFCarrow) >1) {arrCol <- expFCarrow[2] # 2nd value of expFCarrow for color
} else arrCol <- 1 # autom : use 3rd if possible
}
if(drawArrow) {
if(is.finite(expM)) {
figCo <- graphics::par("usr") # c(x1, x2, y1, y2)
figRa <- diff(range(figCo[1:2]))*0.1
if(expM > figCo[2] +figRa || expM < figCo[1] -figRa) {
if(!silent) message(fxNa,"Can't draw arrow, ",round(expM,2)," is too far outside the plotting frame")
expM <- NA }
}
if(is.finite(expM)) {
arr <- c(0.019,0.14) # start- and end-points of arrow (as relative to entire plot)
graphics::arrows(expM, figCo[3] + arr[1]*(figCo[4] -figCo[3]), expM, figCo[3] + arr[2]*(figCo[4] -figCo[3]),
col=arrCol, lwd=1, length=0.1)
graphics::mtext(paste("expect at",signif(expM,3)), at=expM, side=1, adj=0.5, col=arrCol, cex=cexLa*0.7, line=-0.9)
} else { if(!silent) message(fxNa,"Unable to draw arrow for expexted M-value)") }
}
if(debug) {message(fxNa,"VPW19")}
tmp <- try(graphics::par(mar=opar$mar, cex.main=opar$cex.main, las=opar$las), silent=TRUE) # resetto previous
## export data used for plotting
if(returnData) {
merg <- merg[,-1*c(1,ncol(merg))] # remove col 'ID' 'redundant' & 'pch'
annCo <- wrMisc::naOmit(match(annotColumn, colnames(merg)))
if(length(annCo) >0) cbind(merg[,annCo], merg[,-annCo]) else merg }
} }
}
#' Colors based on p-Values
#'
#' This function helps defining color-gradient based on p-Values.
#' This fuction requires package RColorBrewer being installed
#'
#' @param x (numeric) p-values (main input)
#' @param br (numeric) custom breaks (used with cut)
#' @param col custom colors (must be of length(br) -1)
#' @param asIndex (logical) custom breaks (used with cut)
#' @param silent (logical) suppress messages
#' @param callFrom (character) allow easier tracking of messages produced
#' @param debug (logical) supplemental messages for debugging
#' @return This function retruns a color-gradient based on p-Values
#' @seealso (for PCA) \code{\link{plotPCAw}})
#' @examples
#' .colorByPvalue((1:10)/10)
#' @export
.colorByPvalue <- function(x, br=NULL, col=NULL, asIndex=FALSE, silent=FALSE, debug=FALSE, callFrom=NULL) {
## this function should ultimately get incorporated to wrMisc::colorAccording2 as option "pValue" or "FDR"
## colors significant (***,**,*) as red-tones, 0.05-0.075 as pale orange, others as blue-tones
## 'br' .. custom breaks (used with cut)
## 'col' .. custom colors (must be of length(br) -1)
## 'asIndex' .. optional returning of index to \code{col} instead of final color designation
fxNa <- ".colorByPvalue"
if(!isTRUE(silent)) silent <- FALSE
if(isTRUE(debug)) silent <- FALSE else debug <- FALSE
if(length(br) >0) br <- sort(unique(wrMisc::naOmit(br)))
if(length(col) != length(br) -1 && length(br) >0) { br <- col <- NULL
if(!silent) message(fxNa, "Entries for 'col' and 'br' don't match, ignoring") }
if(length(col) <1) {
hasDep <- requireNamespace("RColorBrewer", quietly=TRUE)
col <- if(hasDep) grDevices::colorRampPalette(RColorBrewer::brewer.pal(n=7, name="RdYlBu"))( 22 )[c(2:4,7,17:20)] else {
grDevices::heat.colors(n=length(unique(x)) +2)
if(!silent) message(fxNa,"Please install first from CRAN the package 'RColorBrewer'")
}
if(length(col) != length(br)) br <- sort(c(1e-20,10^(-3:0), 0.05,0.075,0.2,0.3))
}
out <- if(asIndex) as.numeric(cut(x, br)) else col[as.numeric(cut(x, br))]
if(asIndex) out <- as.numeric(cut(x, br)) else {
out <- col[as.numeric(cut(x, br))]
chNa <- is.na(x)
if(any(chNa)) out[which(chNa)] <- grDevices::grey(0.4) }
out }
Try the wrGraph package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.