R/plotLinear.R

In sequinstandards/RAnaquin: Statistical analysis of sequins

.plotLinear <- function(seqs, x, y, title=NULL, xl=NULL, yl=NULL,
                        std=NULL,
                        showCol=TRUE,
                        showSD=TRUE,
                        showLOQ=TRUE,
                        showStats=FALSE,
                        xBreaks=NULL,
                        yBreaks=NULL,
                        errors=NULL,
                        shouldLog=TRUE,
                        showLinear=TRUE,
                        showAxis=FALSE)
{
    stopifnot(is.null(errors) | errors == 'SD' | errors == 'Range')
    
    stopifnot(!is.null(seqs))
    stopifnot(!is.null(x))
    stopifnot(!is.null(y))

    data <- data.frame(row.names=seqs, x=x, y=y)
    data <- data[data$x != "NA" & data$y != "NA" & !is.na(data$x) & !is.na(data$y),]

    if (nrow(data) == 0) { stop("No sequin reads detected, please check your library and ensure you are using the correct reference file.") }

    data$x <- as.numeric.factor(data$x)
    data$y <- as.numeric.factor(data$y)

    if (shouldLog) { data$x <- log2(data$x); data$y <- log2(data$y) }

    if (!is.null(std)) { data$sd <- std }
    
    if (ncol(data) == 2)
    {
        data <- data[!is.na(data$y),]
        data <- data[!is.infinite(data$y),]
    }
    
    # For mutliple measurments
    data$y <- data[,c(2:ncol(data))]
    
    data <- data[!is.na(data$x),]
    data$grp <- as.factor(abs(data$x))
    
    #stopifnot(length(data$x) > 0)
    stopifnot(length(data$x) == length((data$y)) ||
              length(data$x) == nrow((data$y)))
    
    isMultiDF <- is(data$y, 'data.frame')
    
    # Should we show standard deviation?
    isMultiSD <- sum(data$sd) > 0 & showSD
    
    isMulti <- isMultiDF | isMultiSD
    
    data$ymax <- NULL
    data$ymin <- NULL
    
    # All data points (including all replicates)
    data.all <- NULL
    
    if (isMultiDF)
    {
        for (i in 1:ncol(data$y))
        {
            data.all <- rbind(data.all, data.frame(x=data$x, y=(data$y)[,i]))
        }
        
        if (is.null(data$sd))
        {
            data$sd <- apply(data$y, 1, sd, na.rm=TRUE)
        }
        
        data$min <- do.call(pmin, c(as.data.frame(data$y), na.rm=TRUE))
        data$max <- do.call(pmax, c(as.data.frame(data$y), na.rm=TRUE))
        data$y   <- rowMeans(data$y, na.rm=TRUE)
    }
    else
    {
        data.all <- data.frame(x=data$x, y=data$y)
    }
    
    if (isMulti)
    {
        #
        # Quote from: "https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2064100"
        #
        #   About two thirds of the data points will lie within the region of
        #   mean ± 1 SD, and ∼95% of the data points will be within 2 SD of
        #   the mean.
        #
        # We want to establish 95% confidence interval.
        #
        if (is.null(errors) || errors == 'SD')
        {
            data$ymax <- data$y + 2*data$sd
            data$ymin <- data$y - 2*data$sd
        }
        
        #
        # Range error encompass the lowest and highest values
        #
        else if (errors == 'Range')
        {
            data$ymax <- data$max
            data$ymin <- data$min
        }
        
        data <- data[!is.na(data$ymax),]
        data <- data[!is.na(data$ymin),]    
    }
    else
    {
        data$sd <- NULL
    }
    
    data <- data[!is.na(data$y),]
    
    p <- ggplot(data=data, aes_string(x='data$x', y='data$y', label='row.names(data)')) +
                    xlab(xl) +
                    ylab(yl) +
                    ggtitle(title) +
                    labs(colour='Ratio') +
                    geom_point(aes_string(colour='grp'), alpha=0.5, size=0.7) +
                    guides(colour=FALSE) +
                    theme_bw() +
                    theme(plot.title=element_text(hjust=0.5))

    if (!showCol) { p <- p + geom_point(alpha=0.5, size=0.7) }
    
    if (showLinear)
    {
       p <- p + geom_smooth(method='lm',
                           formula=y~x,
                          linetype='11',
                             color='black',
                              size=0.7)
    }
        
    if (showSD & !is.null(data$sd))
    {
        p <- p + geom_segment(aes_string(x='data$x',
                                         y='data$ymax',
                                         xend='data$x',
                                         yend='data$ymin'),
                              data=data,
                              size=0.2,
                              alpha=0.5)
    }
    
    y_off <- ifelse(max(data$y) - min(data$y) <= 10, 0.7, 0.7)
    
    if (showAxis)
    {
        p <- p + geom_vline(xintercept=c(0), linetype='solid', size=0.1)
        p <- p + geom_hline(yintercept=c(0), linetype='solid', size=0.1)
    }

    if (!is.null(xBreaks)) { p <- p + scale_x_continuous(breaks=xBreaks) }
    if (!is.null(yBreaks)) { p <- p + scale_y_continuous(breaks=yBreaks) }
    
    overall <- if(nrow(data) > 0) { .lm2str(data) } else { NULL }
    above   <- NULL
    
    LOQ <- NULL
    
    if (showLOQ)
    {
        tryCatch( {
            LOQ <- estimateLOQ(data.all$x, data.all$y)
        }, error = function(x) {})
        
        if (!is.null(LOQ))
        {
            # Print out the regression above LOQ
            above <- .m2str(LOQ$model$rModel)
            
            # Assume the break-point is on the log2-scale. Convert it back.
            label <- 2^LOQ$breaks$k
            
            p <- p + geom_vline(xintercept=c(LOQ$breaks$k),
                                linetype='33',
                                size=0.6)
            p <- p + geom_label(aes_string(x='max(LOQ$breaks$k)', y='min(y)'),
                                label=paste('LOQ:', signif(label, 3)),
                                colour='black',
                                show.legend=FALSE,
                                hjust=0.1,
                                vjust=0.7)                
        }
    }
    
    if (nrow(data) > 0 && showStats)
    {
        size <- ifelse(!is.null(.env("annotate.text")), as.numeric(.env("annotate.text")), 4)

        if (showLOQ & !is.null(LOQ))
        {
            overall <- paste(c('bold(Overall): ', overall), collapse='')
            above   <- paste(c('bold(Above)~bold(LOQ):', above), collapse='')
            label   <- paste(c('atop(', overall, ',', above, ')'), collapse='')
            
            p <- p + annotate("text",
                              label=label,
                              x=min(data$x),
                              y=max(data$y),
                              size=size,
                              colour='grey24',
                              parse=TRUE,
                              hjust=0,
                              vjust=1)
        }
        else
        {
            p <- p + annotate("text",
                              label=overall,
                              x=min(data$x),
                              y=max(data$y),
                              size=size,
                              colour='grey24',
                              parse=TRUE,
                              hjust=0,
                              vjust=1)
        }
    }
    
    if (!is.null(.env("report")))
    {
        print(lm(data$y ~ data$x)$coefficients[[2]])
        print(cor(data$x,data$y)[[1]]^2)
    }
    
    suppressWarnings(print(.transformPlot(p)))
    LOQ
}

plotLinear <- function(seqs, x, y, title=NULL, xl=NULL, yl=NULL,
                       std=NULL,
                       showCol=TRUE,
                       showSD=TRUE,
                       showLOQ=TRUE,
                       showStats=FALSE,
                       xBreaks=NULL,
                       yBreaks=NULL,
                       errors=NULL,
                       shouldLog=TRUE,
                       showLinear=TRUE,
                       showAxis=FALSE)
{
    tryCatch({
        .plotLinear(seqs, x, y, title, xl, yl,
                    std,
                    showCol,
                    showSD,
                    showLOQ,
                    showStats,
                    xBreaks,
                    yBreaks,
                    errors,
                    shouldLog,
                    showLinear,
                    showAxis) },
    error=function(x) {  print(x); emptyPlot(xl, yl, title) })
}