Nothing
R/VolcanoPlotW.R
In wrGraph: Graphics in the Context of Analyzing High-Throughput Data
Defines functions
.levIndex
.colorByPvalue
.sampNoDeMArrayLM
VolcanoPlotW
Documented in
.colorByPvalue
.levIndex
.sampNoDeMArrayLM
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).
#'
#' 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.
#'
#' @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 ProjNa (character) custom title
#' @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 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 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, ProjNa=NULL, FCthrs=NULL, FdrList=NULL, FdrThrs=NULL, FdrType=NULL,
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, 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
plotByFdr <- TRUE #
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
splNa <- annot <- ptType <- colPass <- ptBg <- grpMeans <- pcol <- FDRvalue <- NULL # initialize
multiComp <- TRUE # initialize
if(debug) silent <- FALSE
if(length(Mvalue) <1) message(fxNa,"Nothing to do, 'Mvalue' seems to be empty !") else {
## data seem valid to make MAplot
if(length(cexTxLab) <0) cexTxLab <- 0.7
if(inherits(Mvalue, "MArrayLM")) {
## try working based on MArrayLM-object (Mvalue)
## 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, setting to 1") }
if(debug) {message(fxNa,"VPW2")}
if(length(pValue) <1) {
## no spearate pValues provided : extract from MArrayLM-object (Mvalue)
pcol <- wrMisc::naOmit(match(c("p.value","pvalue","pval","p"), tolower(names(Mvalue))))
if(length(pcol) >0) pValue <- Mvalue[[pcol[1]]] else stop("can't find suitable element for p-Values from MArrayLM-object")
if(length(dim(pValue)) >0) if(colnames(pValue)[1] =="(Intercept)" & ncol(pValue) >1) {
## extract 2nd col if result from wrMisc::moderTest2grp()
pNa <- rownames(pValue)
pValue <- as.numeric(pValue[,2])
names(pValue) <- pNa
multiComp <- FALSE
} else {
## select corresponding of multiple comparisons
if(useComp > ncol(pValue)) { useComp <- 1
if(!silent) message(fxNa,"Argument 'useComp' for pValues invalid or too high; reset to 1") }
names(useComp) <- colnames(pValue)[useComp]
pNa <- rownames(pValue)
pValue <- as.numeric(pValue[,useComp])
if(length(pNa) >0) names(pValue) <- pNa }
} ## .. otherwise spearate pValue was provided
if(debug) {message(fxNa,"VPW3")}
## look for M-values (need to create if not available - using useComp checked when extracting pValue)
Melem <- wrMisc::naOmit(match(c("mvalues","mvalue","m"), tolower(names(Mvalue)))) # which list-element
Fcol <- wrMisc::naOmit(match(if(length(FdrType) <1) c("fdr","bh","lfdr","by","p.value") else tolower(FdrType), tolower(names(Mvalue)))) # needed for matching means to pair-wise names and for extracting FDR values
if(debug) {message(fxNa,"Fcol : ",wrMisc::pasteC(Fcol)," VPW4")}
if(length(Fcol) <1) stop("Can't find elment suitable for statistical values", if(length(FdrType)==1) c(" to '",FdrType,"'")) else Fcol <- Fcol[1]
if(!silent & length(FdrType) <1) message(fxNa,"Using element '",names(Mvalue)[Fcol],"' as FDR-values for plot")
if("lfdr" %in% names(Mvalue)[Fcol]) { # switch from directly plotting FDR-values to uncorrected p-values
plotByFdr <- FALSE }
## look for group-means & identify column association to current question/pairwise comparison
if("means" %in% names(Mvalue)) {
## identify sammple-groups to comparsison(s) - needed lateron
pairwCol <- wrMisc::sampNoDeMArrayLM(Mvalue, useComp, lstMeans="means", lstP=Fcol, callFrom=fxNa, silent=silent)
grpMeans <- cbind(mean1=Mvalue$means[,pairwCol[1]], mean2=Mvalue$means[,pairwCol[2]])
## are all group-means needed (for exporting) ??
} else {grpMeans <- NULL; warning(fxNa,"Could not find suitable field '$means' in '",namesIn[1],"'") }
if(length(Melem) >0) { ## M-values are available
Mvalue$Mval <- Mvalue[[Melem]]
if(length(dim(Mvalue$Mval)) >0) if(ncol(Mvalue$Mval) >1) {
if(useComp[1] > ncol(Mvalue$Mval)) { if(!silent) message(fxNa," 'useComp' is too high, set to 1"); useComp <- 1 }
Mvalue$Mval <- Mvalue$Mval[,useComp]}
} else { # need to construct M-values based on means
if("means" %in% names(Mvalue)) {
## construct Mvalue based on means (only one/current pairwise comparison needed)
Mvalue$Mval <- grpMeans[,2] - grpMeans[,1]
Melem <- which(names(Mvalue) =="Mval") # update
} else stop("Can't construct M-values since suitable field '$means' missing in '",namesIn[1],"' !")
}
if(debug) {message(fxNa,"VPW4b")}
## now one can check if 'pValue' & Mvalue match
chPM <- length(pValue) >0 && length(as.numeric(pValue)) == length(as.numeric(Mvalue[[Melem]]))
if(!chPM) {
if(length(pcol) >0) { # pVal avilable in Mvalue
if(length(as.numeric(Mvalue[[pcol[1]]][,useComp])) ==length(as.numeric(Mvalue[[Melem]]))) { # pVal from Mvalue seems to fit => use
pValue <- as.numeric(Mvalue[[pcol[1]]][,useComp])
} else stop("'pValue' & 'Mvalue' don't match")
} else stop("'pValue' & 'Mvalue' don't match (no field in MArrayLM-object available)")
}
if(debug) {message(fxNa,"Fcol : ",wrMisc::pasteC(Fcol)," VPW5"); VPW5 <- list(FdrList=FdrList,Mvalue=Mvalue,Fcol=Fcol,Melem=Melem,FdrType=FdrType,plotByFdr=plotByFdr,pValue=pValue,useComp=useComp,pairwCol=pairwCol,grpMeans=grpMeans)}
## extract FDR from MArrayLM-object
if(length(FdrList) <1) {
## no explicit pValue, try to extract from MArrayLM-object (Mvalue)
if(length(Fcol) >0) { FDRvalue <- Mvalue[[Fcol[1]]]
## extract 2nd col if result from wrMisc::moderTest2grp()
if(length(dim(FDRvalue)) >0) if(colnames(FDRvalue)[1] =="(Intercept)" && ncol(FDRvalue) >1) {
pNa <- rownames(FDRvalue)
FDRvalue <- as.numeric(FDRvalue[,2])
names(FDRvalue) <- pNa
multiComp <- FALSE } #else FDRvalue <- as.numeric(FDRvalue[,useComp])
} else {
FDRvalue <- stats::p.adjust(pValue)
if(!silent) message(fxNa,"No FDR data found, generating BH-FDR")}
## need to find corresponding of multiple comparisons
if(length(dim(FDRvalue)) >0) {
pNa <- rownames(FDRvalue)
if(ncol(FDRvalue) >1) { if(useComp > ncol(FDRvalue)) { useComp <- 1
if(!silent) message(fxNa," argument 'useComp' for FDRvalues invalid or too high; reset to 1") }
FDRvalue <- as.numeric(FDRvalue[,useComp])
if(length(pNa) >0) names(FDRvalue) <- pNa }}
}
if(debug) {message(fxNa,"VPW6")}
## recuperate filtering - if present, but only when no custom filtering provided
if(length(filtFin) <1 || isFALSE(filtFin)) {
Fcol <- wrMisc::naOmit(match(c("filtfin","filter","filt","finfilt"), tolower(names(Mvalue))))
filtFin <- if(length(Fcol) >0) Mvalue[[Fcol[1]]] else rep(TRUE,length(pValue))
if(length(dim(filtFin)) >1) filtFin <- filtFin[,useComp]
}
## recuperate $annot if present and use for symbol
if("annot" %in% names(Mvalue)) {
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
} }
if(length(pch)==1) pch <- rep(as.integer(pch), length(pValue))
if(debug) {message(fxNa,"VPW7")}
## recuperate M values (& dismiss rest of MArrayLM-object)
Mvalue <- Mvalue$Mval
if(length(dim(Mvalue)) >1) { MNa <- rownames(Mvalue)
Mvalue <- as.numeric(Mvalue)
if(length(MNa) >0) names(Mvalue) <- MNa }
## additional check for length
chpM <- length(Mvalue)==length(pValue)
if(!chpM & !silent) message(fxNa,"Trouble ahead ? p- and M- values have different length !! (M=",length(Mvalue)," vs p=",length(pValue),")")
## done with extracing MArrayLM-object
if(!silent) message(fxNa,"Successfully extracted ",length(Mvalue)," Mvalues and ",length(pValue)," pValues", if(length(annot) >0) c(" plus anotation"))
} else {
## thus argument 'Mvalue' is not 'MArrayLM'-object
## ... 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")}
## need to introduce -log10 to pValue
chNA <- is.na(pValue)
if(all(chNA)) stop(fxNa," All p-values are NA, nothing to draw !")
if(any(pValue <0, na.rm=TRUE)) warning(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,plotByFdr=plotByFdr,pValue=pValue,useComp=useComp,pairwCol=pairwCol,grpMeans=grpMeans, filtFin=filtFin,annotColumn=annotColumn,annot=annot, pch=pch,sortLeg=sortLeg,batchFig=batchFig)}
## 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(merg))
## 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")}
## 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")}
## 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]})
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 <- paste(c(if(!batchFig) c(ProjNa, if(!is.null(ProjNa)) ": ","Volcano-Plot"),
if(!is.null(compNa)) 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")}
## 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")}
## 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)]
}
## 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
## adjust fill color for open symbols
chPch <- pch %in% c(21:25)
if(any(chPch)) { 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,"VPW14")} }
if(length(Mvalue) >0) {
## rather directly plot FDR
if(!"FDR" %in% colnames(merg)) plotByFdr <- FALSE
pl1 <- try(graphics::plot(Mvalue, if(plotByFdr) -1*log10(merg[,"FDR"]) else merg[,"pValue"], pch=pch, cex=useCex, main=tit1,
ylab=if(plotByFdr) "- log10 FDR" 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(length(Mvalue) >0) {
sTxt <- if(length(subTxt) ==1) subTxt else { if(multiComp) paste0(if(length(names(useComp)) >0) names(useComp) else c("useComp=",useComp),"; ",collapse="")}
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)
graphics::abline(v=c(-1,1)*(log2(FCthrs) + diff(graphics::par("usr")[1:2])/500), 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,plotByFdr=plotByFdr,limM=limM,limp=limp,ptBg=ptBg,cexLa=cexLa)}
if(sum(passFdr, na.rm=TRUE) >0) {
if(plotByFdr) {
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 {
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=-1.4,cex=0.7)
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 -1*log10(merg[useLi[useL2], if(plotByFdr) "FDR" else "pValue"]),
names(tmP)[useL2], cex=cexTxLab, col=NbestCol, adj=0) }
}
}
if(debug) {message(fxNa,"VPW17")}
## 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 }
} }
}
#' Locate sample index from index or name of pair-wise comparisons
#'
#' This function helps locating sample index from index or name of pair-wise comparisons
#'
#' @param MArrayObj (MArray type objct) main input
#' @param useComp (matrix) types of pair-wise comparisons to be performed
#' @param groupSep (character) separator used with pair-wise grouping
#' @param lstMeans (character) type of summarization, default is 'means'
#' @param lstP (character) type of multiple testing correction data to choose from \code{MArrayObj}
#' @param silent (logical) suppress messages
#' @param debug (logical) supplemental messages for debugging
#' @param callFrom (character) allow easier tracking of messages produced
#' @return This function returns a integer vector of indexes
#' @seealso (for PCA) \code{\link{plotPCAw}})
#' @examples
#' aa <- 1:5
#' @export
.sampNoDeMArrayLM <- function(MArrayObj, useComp, groupSep="-", lstMeans="means", lstP="BH", silent=FALSE, debug=FALSE, callFrom=NULL) {
## old version,not used any more
## locate sample index from index or name of pair-wise comparisons in list or MArrayLM-object
fxNa <- wrMisc::.composeCallName(callFrom, newNa=".sampNoDeMArrayLM")
if(!isTRUE(silent)) silent <- FALSE
if(isTRUE(debug)) silent <- FALSE else debug <- FALSE
errMsg <- c("Argument 'MArrayObj' is ","empty","doesn't contain the list-element needed ('",lstMeans,"') !")
if(length(MArrayObj) <1) stop(errMsg[1:2])
if(length(MArrayObj[[lstMeans]]) <1) stop(errMsg[-2])
if(length(colnames(MArrayObj[[lstMeans]])) <1) stop(" Problem with 'MArrayObj$lstMeans' (does not contain matrix of means)")
if(ncol(MArrayObj[[lstMeans]]) ==2) { # only 2 mean-values, no other choice, don't need to try matching anything
if(!identical(as.character(useComp),"1") && !silent) message(fxNa,"Only 2 columns of mean-values available, can't interpret properly 'useComp=",useComp,"'")
out <- 1:2
} else {
if(length(lstP) <0) stop(" 'lstP' is empty !")
if(length(colnames(MArrayObj[[lstP]])) <1) stop(" Problem with 'MArrayObj' (does not contain matrix of p-values)")
## convert/locate names to index
if(is.character(useComp) && length(grep("[[:alpha:]]",useComp)) >0) useComp <- wrMisc::naOmit(match(useComp, MArrayObj[[lstP]] ))
if(length(useComp) <1) stop("Argument 'useComp' is empty or can't locate in comparison-names")
if(length(useComp) >1) { useComp <- useComp[1]; message(fxNa,"Using only 1st instance of argument 'useComp'")}
if(useComp > ncol(MArrayObj[[lstP]])) {useComp <- 1; message(fxNa,"Argument 'useComp' is too high (max ",ncol(MArrayObj[[lstP]]),"), resetting to 1")}
## main
out <- if(ncol(MArrayObj[[lstP]])==2 && colnames(MArrayObj[[lstP]])[1] =="(Intercept)") 1:2 else {
wrMisc::matchSampToPairw(grpNa=colnames(MArrayObj[[lstMeans]]), pairwNa=colnames(MArrayObj[[lstP]])[useComp], sep=groupSep,silent=silent,callFrom=fxNa)}
if(length(useComp)==1) out <- as.integer(out)}
out }
#' 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 incormporated to wtMisc::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 }
#' Transform levels into index
#'
#' This function transforms levels into index. This function has been depreciated, please use wrMisc::levIndex() insted.
#'
#' @param dat (numeric) initial levels (main input)
#' @param asSortedLevNa (logical)
#' @return This function retruns a color-gradient based on p-Values
#' @seealso (for PCA) \code{\link{plotPCAw}})
#' @examples
#' library(wrMisc)
#' @export
.levIndex <- function(dat, asSortedLevNa=FALSE) {
## transform levels into index; should get integrated to wrMisc
## depreciated, please use wrMisc::levIndex() instead
fxNa <- ".levIndex"
.Deprecated("levIndex()", package="wrMisc", msg="The function .levIndex() has been depreciated, please use wrMisc::levIndex() instead")
if(asSortedLevNa) out <- as.integer(as.factor(dat)) else {
out <- dat
levU <- wrMisc::naOmit(unique(out)) # levels in orig order (non-alpahbetical)
for(i in 1:length(levU)) out[which(out==levU[i])] <- i
out <- as.integer(out)}
out }
Try the wrGraph package in your browser
Any scripts or data that you put into this service are public.
wrGraph documentation built on Oct. 20, 2023, 5:08 p.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 .levIndex .colorByPvalue .sampNoDeMArrayLM VolcanoPlotW
Documented in .colorByPvalue .levIndex .sampNoDeMArrayLM 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).
#'
#' 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.
#'
#' @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 ProjNa (character) custom title
#' @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 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 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, ProjNa=NULL, FCthrs=NULL, FdrList=NULL, FdrThrs=NULL, FdrType=NULL,
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, 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
plotByFdr <- TRUE #
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
splNa <- annot <- ptType <- colPass <- ptBg <- grpMeans <- pcol <- FDRvalue <- NULL # initialize
multiComp <- TRUE # initialize
if(debug) silent <- FALSE
if(length(Mvalue) <1) message(fxNa,"Nothing to do, 'Mvalue' seems to be empty !") else {
## data seem valid to make MAplot
if(length(cexTxLab) <0) cexTxLab <- 0.7
if(inherits(Mvalue, "MArrayLM")) {
## try working based on MArrayLM-object (Mvalue)
## 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, setting to 1") }
if(debug) {message(fxNa,"VPW2")}
if(length(pValue) <1) {
## no spearate pValues provided : extract from MArrayLM-object (Mvalue)
pcol <- wrMisc::naOmit(match(c("p.value","pvalue","pval","p"), tolower(names(Mvalue))))
if(length(pcol) >0) pValue <- Mvalue[[pcol[1]]] else stop("can't find suitable element for p-Values from MArrayLM-object")
if(length(dim(pValue)) >0) if(colnames(pValue)[1] =="(Intercept)" & ncol(pValue) >1) {
## extract 2nd col if result from wrMisc::moderTest2grp()
pNa <- rownames(pValue)
pValue <- as.numeric(pValue[,2])
names(pValue) <- pNa
multiComp <- FALSE
} else {
## select corresponding of multiple comparisons
if(useComp > ncol(pValue)) { useComp <- 1
if(!silent) message(fxNa,"Argument 'useComp' for pValues invalid or too high; reset to 1") }
names(useComp) <- colnames(pValue)[useComp]
pNa <- rownames(pValue)
pValue <- as.numeric(pValue[,useComp])
if(length(pNa) >0) names(pValue) <- pNa }
} ## .. otherwise spearate pValue was provided
if(debug) {message(fxNa,"VPW3")}
## look for M-values (need to create if not available - using useComp checked when extracting pValue)
Melem <- wrMisc::naOmit(match(c("mvalues","mvalue","m"), tolower(names(Mvalue)))) # which list-element
Fcol <- wrMisc::naOmit(match(if(length(FdrType) <1) c("fdr","bh","lfdr","by","p.value") else tolower(FdrType), tolower(names(Mvalue)))) # needed for matching means to pair-wise names and for extracting FDR values
if(debug) {message(fxNa,"Fcol : ",wrMisc::pasteC(Fcol)," VPW4")}
if(length(Fcol) <1) stop("Can't find elment suitable for statistical values", if(length(FdrType)==1) c(" to '",FdrType,"'")) else Fcol <- Fcol[1]
if(!silent & length(FdrType) <1) message(fxNa,"Using element '",names(Mvalue)[Fcol],"' as FDR-values for plot")
if("lfdr" %in% names(Mvalue)[Fcol]) { # switch from directly plotting FDR-values to uncorrected p-values
plotByFdr <- FALSE }
## look for group-means & identify column association to current question/pairwise comparison
if("means" %in% names(Mvalue)) {
## identify sammple-groups to comparsison(s) - needed lateron
pairwCol <- wrMisc::sampNoDeMArrayLM(Mvalue, useComp, lstMeans="means", lstP=Fcol, callFrom=fxNa, silent=silent)
grpMeans <- cbind(mean1=Mvalue$means[,pairwCol[1]], mean2=Mvalue$means[,pairwCol[2]])
## are all group-means needed (for exporting) ??
} else {grpMeans <- NULL; warning(fxNa,"Could not find suitable field '$means' in '",namesIn[1],"'") }
if(length(Melem) >0) { ## M-values are available
Mvalue$Mval <- Mvalue[[Melem]]
if(length(dim(Mvalue$Mval)) >0) if(ncol(Mvalue$Mval) >1) {
if(useComp[1] > ncol(Mvalue$Mval)) { if(!silent) message(fxNa," 'useComp' is too high, set to 1"); useComp <- 1 }
Mvalue$Mval <- Mvalue$Mval[,useComp]}
} else { # need to construct M-values based on means
if("means" %in% names(Mvalue)) {
## construct Mvalue based on means (only one/current pairwise comparison needed)
Mvalue$Mval <- grpMeans[,2] - grpMeans[,1]
Melem <- which(names(Mvalue) =="Mval") # update
} else stop("Can't construct M-values since suitable field '$means' missing in '",namesIn[1],"' !")
}
if(debug) {message(fxNa,"VPW4b")}
## now one can check if 'pValue' & Mvalue match
chPM <- length(pValue) >0 && length(as.numeric(pValue)) == length(as.numeric(Mvalue[[Melem]]))
if(!chPM) {
if(length(pcol) >0) { # pVal avilable in Mvalue
if(length(as.numeric(Mvalue[[pcol[1]]][,useComp])) ==length(as.numeric(Mvalue[[Melem]]))) { # pVal from Mvalue seems to fit => use
pValue <- as.numeric(Mvalue[[pcol[1]]][,useComp])
} else stop("'pValue' & 'Mvalue' don't match")
} else stop("'pValue' & 'Mvalue' don't match (no field in MArrayLM-object available)")
}
if(debug) {message(fxNa,"Fcol : ",wrMisc::pasteC(Fcol)," VPW5"); VPW5 <- list(FdrList=FdrList,Mvalue=Mvalue,Fcol=Fcol,Melem=Melem,FdrType=FdrType,plotByFdr=plotByFdr,pValue=pValue,useComp=useComp,pairwCol=pairwCol,grpMeans=grpMeans)}
## extract FDR from MArrayLM-object
if(length(FdrList) <1) {
## no explicit pValue, try to extract from MArrayLM-object (Mvalue)
if(length(Fcol) >0) { FDRvalue <- Mvalue[[Fcol[1]]]
## extract 2nd col if result from wrMisc::moderTest2grp()
if(length(dim(FDRvalue)) >0) if(colnames(FDRvalue)[1] =="(Intercept)" && ncol(FDRvalue) >1) {
pNa <- rownames(FDRvalue)
FDRvalue <- as.numeric(FDRvalue[,2])
names(FDRvalue) <- pNa
multiComp <- FALSE } #else FDRvalue <- as.numeric(FDRvalue[,useComp])
} else {
FDRvalue <- stats::p.adjust(pValue)
if(!silent) message(fxNa,"No FDR data found, generating BH-FDR")}
## need to find corresponding of multiple comparisons
if(length(dim(FDRvalue)) >0) {
pNa <- rownames(FDRvalue)
if(ncol(FDRvalue) >1) { if(useComp > ncol(FDRvalue)) { useComp <- 1
if(!silent) message(fxNa," argument 'useComp' for FDRvalues invalid or too high; reset to 1") }
FDRvalue <- as.numeric(FDRvalue[,useComp])
if(length(pNa) >0) names(FDRvalue) <- pNa }}
}
if(debug) {message(fxNa,"VPW6")}
## recuperate filtering - if present, but only when no custom filtering provided
if(length(filtFin) <1 || isFALSE(filtFin)) {
Fcol <- wrMisc::naOmit(match(c("filtfin","filter","filt","finfilt"), tolower(names(Mvalue))))
filtFin <- if(length(Fcol) >0) Mvalue[[Fcol[1]]] else rep(TRUE,length(pValue))
if(length(dim(filtFin)) >1) filtFin <- filtFin[,useComp]
}
## recuperate $annot if present and use for symbol
if("annot" %in% names(Mvalue)) {
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
} }
if(length(pch)==1) pch <- rep(as.integer(pch), length(pValue))
if(debug) {message(fxNa,"VPW7")}
## recuperate M values (& dismiss rest of MArrayLM-object)
Mvalue <- Mvalue$Mval
if(length(dim(Mvalue)) >1) { MNa <- rownames(Mvalue)
Mvalue <- as.numeric(Mvalue)
if(length(MNa) >0) names(Mvalue) <- MNa }
## additional check for length
chpM <- length(Mvalue)==length(pValue)
if(!chpM & !silent) message(fxNa,"Trouble ahead ? p- and M- values have different length !! (M=",length(Mvalue)," vs p=",length(pValue),")")
## done with extracing MArrayLM-object
if(!silent) message(fxNa,"Successfully extracted ",length(Mvalue)," Mvalues and ",length(pValue)," pValues", if(length(annot) >0) c(" plus anotation"))
} else {
## thus argument 'Mvalue' is not 'MArrayLM'-object
## ... 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")}
## need to introduce -log10 to pValue
chNA <- is.na(pValue)
if(all(chNA)) stop(fxNa," All p-values are NA, nothing to draw !")
if(any(pValue <0, na.rm=TRUE)) warning(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,plotByFdr=plotByFdr,pValue=pValue,useComp=useComp,pairwCol=pairwCol,grpMeans=grpMeans, filtFin=filtFin,annotColumn=annotColumn,annot=annot, pch=pch,sortLeg=sortLeg,batchFig=batchFig)}
## 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(merg))
## 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")}
## 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")}
## 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]})
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 <- paste(c(if(!batchFig) c(ProjNa, if(!is.null(ProjNa)) ": ","Volcano-Plot"),
if(!is.null(compNa)) 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")}
## 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")}
## 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)]
}
## 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
## adjust fill color for open symbols
chPch <- pch %in% c(21:25)
if(any(chPch)) { 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,"VPW14")} }
if(length(Mvalue) >0) {
## rather directly plot FDR
if(!"FDR" %in% colnames(merg)) plotByFdr <- FALSE
pl1 <- try(graphics::plot(Mvalue, if(plotByFdr) -1*log10(merg[,"FDR"]) else merg[,"pValue"], pch=pch, cex=useCex, main=tit1,
ylab=if(plotByFdr) "- log10 FDR" 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(length(Mvalue) >0) {
sTxt <- if(length(subTxt) ==1) subTxt else { if(multiComp) paste0(if(length(names(useComp)) >0) names(useComp) else c("useComp=",useComp),"; ",collapse="")}
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)
graphics::abline(v=c(-1,1)*(log2(FCthrs) + diff(graphics::par("usr")[1:2])/500), 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,plotByFdr=plotByFdr,limM=limM,limp=limp,ptBg=ptBg,cexLa=cexLa)}
if(sum(passFdr, na.rm=TRUE) >0) {
if(plotByFdr) {
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 {
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=-1.4,cex=0.7)
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 -1*log10(merg[useLi[useL2], if(plotByFdr) "FDR" else "pValue"]),
names(tmP)[useL2], cex=cexTxLab, col=NbestCol, adj=0) }
}
}
if(debug) {message(fxNa,"VPW17")}
## 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 }
} }
}
#' Locate sample index from index or name of pair-wise comparisons
#'
#' This function helps locating sample index from index or name of pair-wise comparisons
#'
#' @param MArrayObj (MArray type objct) main input
#' @param useComp (matrix) types of pair-wise comparisons to be performed
#' @param groupSep (character) separator used with pair-wise grouping
#' @param lstMeans (character) type of summarization, default is 'means'
#' @param lstP (character) type of multiple testing correction data to choose from \code{MArrayObj}
#' @param silent (logical) suppress messages
#' @param debug (logical) supplemental messages for debugging
#' @param callFrom (character) allow easier tracking of messages produced
#' @return This function returns a integer vector of indexes
#' @seealso (for PCA) \code{\link{plotPCAw}})
#' @examples
#' aa <- 1:5
#' @export
.sampNoDeMArrayLM <- function(MArrayObj, useComp, groupSep="-", lstMeans="means", lstP="BH", silent=FALSE, debug=FALSE, callFrom=NULL) {
## old version,not used any more
## locate sample index from index or name of pair-wise comparisons in list or MArrayLM-object
fxNa <- wrMisc::.composeCallName(callFrom, newNa=".sampNoDeMArrayLM")
if(!isTRUE(silent)) silent <- FALSE
if(isTRUE(debug)) silent <- FALSE else debug <- FALSE
errMsg <- c("Argument 'MArrayObj' is ","empty","doesn't contain the list-element needed ('",lstMeans,"') !")
if(length(MArrayObj) <1) stop(errMsg[1:2])
if(length(MArrayObj[[lstMeans]]) <1) stop(errMsg[-2])
if(length(colnames(MArrayObj[[lstMeans]])) <1) stop(" Problem with 'MArrayObj$lstMeans' (does not contain matrix of means)")
if(ncol(MArrayObj[[lstMeans]]) ==2) { # only 2 mean-values, no other choice, don't need to try matching anything
if(!identical(as.character(useComp),"1") && !silent) message(fxNa,"Only 2 columns of mean-values available, can't interpret properly 'useComp=",useComp,"'")
out <- 1:2
} else {
if(length(lstP) <0) stop(" 'lstP' is empty !")
if(length(colnames(MArrayObj[[lstP]])) <1) stop(" Problem with 'MArrayObj' (does not contain matrix of p-values)")
## convert/locate names to index
if(is.character(useComp) && length(grep("[[:alpha:]]",useComp)) >0) useComp <- wrMisc::naOmit(match(useComp, MArrayObj[[lstP]] ))
if(length(useComp) <1) stop("Argument 'useComp' is empty or can't locate in comparison-names")
if(length(useComp) >1) { useComp <- useComp[1]; message(fxNa,"Using only 1st instance of argument 'useComp'")}
if(useComp > ncol(MArrayObj[[lstP]])) {useComp <- 1; message(fxNa,"Argument 'useComp' is too high (max ",ncol(MArrayObj[[lstP]]),"), resetting to 1")}
## main
out <- if(ncol(MArrayObj[[lstP]])==2 && colnames(MArrayObj[[lstP]])[1] =="(Intercept)") 1:2 else {
wrMisc::matchSampToPairw(grpNa=colnames(MArrayObj[[lstMeans]]), pairwNa=colnames(MArrayObj[[lstP]])[useComp], sep=groupSep,silent=silent,callFrom=fxNa)}
if(length(useComp)==1) out <- as.integer(out)}
out }
#' 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 incormporated to wtMisc::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 }
#' Transform levels into index
#'
#' This function transforms levels into index. This function has been depreciated, please use wrMisc::levIndex() insted.
#'
#' @param dat (numeric) initial levels (main input)
#' @param asSortedLevNa (logical)
#' @return This function retruns a color-gradient based on p-Values
#' @seealso (for PCA) \code{\link{plotPCAw}})
#' @examples
#' library(wrMisc)
#' @export
.levIndex <- function(dat, asSortedLevNa=FALSE) {
## transform levels into index; should get integrated to wrMisc
## depreciated, please use wrMisc::levIndex() instead
fxNa <- ".levIndex"
.Deprecated("levIndex()", package="wrMisc", msg="The function .levIndex() has been depreciated, please use wrMisc::levIndex() instead")
if(asSortedLevNa) out <- as.integer(as.factor(dat)) else {
out <- dat
levU <- wrMisc::naOmit(unique(out)) # levels in orig order (non-alpahbetical)
for(i in 1:length(levU)) out[which(out==levU[i])] <- i
out <- as.integer(out)}
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.