R/plotAllele.R

In sequinstandards/RAnaquin: Statistical analysis of sequins

.plotAllele <- function(x, y, label, legs, cols, xlab, ylab, title)
{
    data <- data.frame(ExpFreq=x, ObsFreq=y, Label=label)
    
    # False positives don't have expected allele frequency
    data$ExpFreq <- as.numeric.factor(data$ExpFreq)
    
    # False negatives don't have measured allele frequency
    data$ObsFreq <- as.numeric.factor(data$ObsFreq)
    
    # Only detected sequins
    data <- data[!is.na(data$ObsFreq),]
    
    if (nrow(data) == 0)
    {
        return(emptyPlot(xlab, ylab, title))
    }
    
    data <- data[,c('ExpFreq', 'ObsFreq', 'Label')]
    
    x1 <- min(data$ExpFreq[data$ExpFreq > 0], na.rm=TRUE)
    x2 <- max(data$ExpFreq[data$ExpFreq > 0], na.rm=TRUE)
    y1 <- min(data$ObsFreq[data$ObsFreq > 0], na.rm=TRUE)
    y2 <- max(data$ObsFreq[data$ObsFreq > 0], na.rm=TRUE)
    
    if (nrow(data[data$ObsFreq == 0,]))
    {
        data[data$ObsFreq == 0,]$ObsFreq <- y1
    }
    
    data[data$ExpFreq != 0,]$ExpFreq <- log2(data[data$ExpFreq != 0,]$ExpFreq)
    data$ObsFreq <- log2(data$ObsFreq)
    
    breaks <- c(0.0001, 0.001, 0.01, 0.1, 0.5)
    xlims  <- log2(c(x1, x2))
    ylims  <- log2(c(y1, y2)) #log2(c(min(breaks[breaks >= y1]), max(breaks[breaks >= y1])))
    
    p <- ggplot(data=data, aes_string(x='data$ExpFreq', y='data$ObsFreq')) +
           xlab(xlab) +
           ylab(ylab) +
           ggtitle(title) +
           labs(colour='') +
           geom_point(aes_string(colour='Label'), alpha=0.3) +
           scale_x_continuous(breaks=c(log2(breaks)), labels=breaks, limits=xlims) +
           scale_y_continuous(breaks=c(log2(breaks)), labels=breaks, limits=ylims) +
           theme_bw() +
           theme(legend.position="none", plot.title = element_text(hjust = 0.5))
    
    minX <- min(data$ExpFreq)
    maxY <- max(data$ObsFreq)
    
    if (!is.null(legs)) { p <- p + scale_color_manual(labels=legs, values=cols) }
    
    if ("TP" %in% data$Label || "FP" %in% data$Label)
    {
        p <- p + annotation_raster(AnaquinSupImage1, ymin=maxY - 2.5, ymax=maxY, xmin=minX, xmax=minX + 4)
    }
    
    p <- p + geom_smooth(data=data[data$Label == 'TP',], aes(x=ExpFreq, y=ObsFreq),
                         colour='black', linetype='dashed', size=0.5,
                         method='lm',
                         formula=y~x)
    
    suppressWarnings(print(.transformPlot(p)))
}

plotAllele <- function(x, y, label, legs, cols, xl="Expected Allele Frequency", yl="Measured Allele Frequency", title="Allele Frequency Ladder")
{
    tryCatch({ .plotAllele(x, y, label, legs, cols, xl, yl, title) }, error=function(x) { print(x); emptyPlot(xl, yl, title) })
}