Files
ROBIBarcodes/R/mismatchplot.R

113 lines
3.7 KiB
R

#'@include ROBIBarcodes.R
#'@include logo.R
NULL
#' Draw a scatter plot of the reverse mismatches as a function of forward mismatches.
#'
#' The \code{mismatchplot} function draws a scatter plot of the number of mismatches
#' observed in an ecoPCR result for the reverse primer as a function of the mismatches
#' for the reverse primer. Each point for a pair (forward_mismatch,reverse_mismatch) is
#' drawn as a circle having a surface proportional to the aboundance of this pair in the
#' results. If a grouping factor is specified, then the circle is replaced by a pie chart.
#'
#' @param ecopcr an ecoPCR result data.frame as returned by the \code{\link{read.ecopcr.result}}
#' function.
#'
#' @param group a factor decribing classes amongst the amplicon described in the ecoPCR
#' result
#'
#' @param col a vector describing the colored used for the bubble or the pie charts
#'
#' @param legend a character vector describing the legend for each modality of the
#' grouping factor. By default the factor levels are used for the legend
#'
#' @param legend.cex the expension factor for the legend text
#'
#' @param inset the distance to the margin of the legend box (see the \code{\link{legend}}
#' documentation)
#'
#' @param view.legend if set to \code{FALSE} the legend corresponding to the groups is not dispayed.
#'
#' @param maxerror allows for specifying the maximum of errors to display on the graph.
#'
#' @examples
#'
#' # Load the ROBITools library
#' library(ROBITools)
#'
#' # Load the default taxonomy
#' taxo = default.taxonomy()
#'
#' # Load the sample ecoPCR data file
#' data(GH.ecopcr)
#'
#' # Computes classes associated to each taxid
#' orders = as.factor(taxonatrank(taxo,GH.ecopcr$taxid,'order',name=T))
#'
#' # Plot the graph
#' mismatchplot(GH.ecopcr,group=orders)
#'
#' @seealso \code{\link{read.ecopcr.result}}
#' @author Eric Coissac
#' @export
mismatchplot = function(ecopcr,group=NULL,
col=NULL,legend=NULL,
legend.cex=0.7,inset=c(0.02,0.02),
view.legend=TRUE,
maxerror=NA) {
maxforward_error = max(ecopcr$forward_mismatch)
maxreverse_error = max(ecopcr$reverse_mismatch)
if (is.na(maxerror))
maxerror=max(maxforward_error,maxreverse_error)
if (is.null(group))
group=factor(rep("all",dim(ecopcr)[1]))
else
group=as.factor(group)
if (is.null(legend))
legend = levels(group)
actualheight= maxerror + 1
actualwidth = maxerror + 1
if (length(levels(group)) > 1 & view.legend)
actualwidth = actualwidth + 2
whitepaper(actualwidth,actualheight,xmin=-0.5,ymin=-0.5,asp=1)
axis(1,at=0:maxerror,
labels=0:maxerror)
axis(2,at=0:maxerror,
labels=0:maxerror)
data = aggregate(group,by=list(forward=factor(ecopcr$forward_mismatch,levels=0:maxerror),
reverse=factor(ecopcr$reverse_mismatch,levels=0:maxerror)),
table)
data <- data[rowSums(data[,c(-1,-2),drop=FALSE])>0, , drop=FALSE]
if (is.null(col))
col <- c("white", "lightblue", "mistyrose", "lightcyan",
"lavender", "cornsilk")
value=data[,c(-1,-2),drop=FALSE]
x = as.integer(data[,1]) - 1
y = as.integer(data[,2]) - 1
diam = sqrt(rowSums(value))
radius = diam / max(diam) / 2
hide = mapply(pie.xy,x,y,
data=lapply(1:(dim(value)[1]),function(y) value[y,]),
radius=radius,
label="",MoreArgs=list(col=col))
if ((length(levels(group))) > 1 & view.legend)
legend('topright',legend=legend,fill=col, cex=legend.cex, inset=inset)
}