mismatchplot.R 3.69 KB
Newer Older
Eric Coissac committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
#'@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)
#'               
29 30 31 32
#' @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.
#' 
Eric Coissac committed
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
#' @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,
55 56 57
                        legend.cex=0.7,inset=c(0.02,0.02),
                        view.legend=TRUE,
                        maxerror=NA) {
Eric Coissac committed
58 59 60
  
  maxforward_error = max(ecopcr$forward_mismatch)
  maxreverse_error = max(ecopcr$reverse_mismatch)
61 62
  if (is.na(maxerror))
    maxerror=max(maxforward_error,maxreverse_error)
Eric Coissac committed
63 64 65 66 67 68 69 70 71 72 73 74

  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
  
75
  if (length(levels(group)) > 1 & view.legend)
Eric Coissac committed
76 77 78 79 80 81 82 83 84 85 86
    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)
  
  
87 88
  data = aggregate(group,by=list(forward=factor(ecopcr$forward_mismatch,levels=0:maxerror),
                                 reverse=factor(ecopcr$reverse_mismatch,levels=0:maxerror)),
Eric Coissac committed
89 90 91 92 93 94 95 96 97 98
                   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]
99 100
  x = as.integer(data[,1]) - 1
  y = as.integer(data[,2]) - 1
Eric Coissac committed
101 102 103 104 105 106 107 108 109
  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))
  
  
110
  if ((length(levels(group))) > 1 & view.legend)
Eric Coissac committed
111 112 113
    legend('topright',legend=legend,fill=col, cex=legend.cex, inset=inset)
  
}