R/findLocalMax.R

In chengvt/EEM: Read and Preprocess Fluorescence Excitation-Emission Matrix (EEM) Data

#' Find local maximum peaks
#' 
#' Find local maximum peaks in EEM data
#' 
#' @param data EEM data generated by \code{\link{readEEM}} function, unfolded EEM data 
#' generated by \code{\link{unfold}} function or  a vector of numeric values which have names in 
#' the format of EX...EM...
#' @param n sample number. The number should not exceed \code{length(EEM)}.
#' @param threshold threshold value in between 0 and 1. Lower the value to cover low peaks.
#' @param showprint logical value whether to print out the results or not
#' @param ... (optional) further arguments passed to other methods
#' 
#' @return return a character vector of peak names. If showprint = TRUE, it will also 
#' print a dataframe of indicating the value of local maximum peaks.
#' 
#' @examples
#' data(applejuice)
#' findLocalMax(applejuice, 1) 
#' 
#' applejuice_uf <- unfold(applejuice)
#' findLocalMax(applejuice_uf, 1) 
#' 
#' @export
#' 
#' @importFrom reshape2 melt
#' @importFrom sp point.in.polygon
#' @importFrom grDevices contourLines
#' 
findLocalMax <- function(data, ...) UseMethod("findLocalMax", data)

#' @describeIn findLocalMax for EEM data created by \code{\link{readEEM}} function
#' @export
findLocalMax.EEM <- function(data, n, threshold = 0.7, showprint = TRUE, ...){
    
        # make sure that n is given
        if (nargs() < 2) stop('Require at least 2 arguments')
    
        # get information from EEM
        data <- data[[n]]
        
        # calculate for regional max
        peak_names <- findLocalMax.internal(data, threshold = threshold, showprint = showprint)
        return(peak_names)
}

#' @describeIn findLocalMax for unfolded EEM data created by \code{\link{unfold}} function
#' @export
findLocalMax.matrix <- function(data, n, threshold = 0.7, showprint = TRUE, ...){
    
    # make sure that n is given
    if (nargs() < 2) stop('Require at least 2 arguments')
    
    # get information from data
    data <- fold(data[n,])[[1]]
    
    # calculate for regional max
    peak_names <- findLocalMax.internal(data, threshold = threshold, showprint = showprint)
    return(peak_names)
}

#' @describeIn findLocalMax for a vector of numeric values which have names in 
#' the format of EX...EM...
#' @export
findLocalMax.numeric <- function(data, threshold = 0.7, showprint = TRUE, ...){

    # convert data to matrix form
    name <- names(data)
    EX <- getEX(name)
    EM <- getEM(name)
    data <- data.frame(ex = as.numeric(EX), em = as.numeric(EM), value = data)
    data.casted <- acast(data, em ~ ex, value.var = "value")
    
    # calculate for regional max
    peak_names <- findLocalMax.internal(data, threshold = threshold, showprint = showprint)
    return(peak_names)
}

#' @export
findLocalMax.internal <- function(data, threshold = 0.7, showprint = TRUE, ...){
    
    # data = a matrix with columns being 
    # excitation wavelength and rows being emission wavelength
    
    if (threshold == 1) stop("Threshold value has to be lower than 1.")
    
    # get contourLines
    x <- as.numeric(colnames(data)) # EX
    y <- as.numeric(rownames(data)) # EM
    z <- t(as.matrix(data))
    cLine <- contourLines(x, y, z, nlevels = 100) 
    
    # melt z
    z_melted <- melt(z)
    names(z_melted)[1:2] <- c("EX", "EM")
    
    # retrieve high contours
    level <- sapply(cLine, function(x) x$level)
    dif <- diff(range(level))
    cutoff <- dif * threshold
    value <- min(level[level > cutoff])
    cLine_selected <- cLine[level == value]
    
    # find max value for each contour
    local_max <- data.frame(EX = numeric(), EM = numeric(), value = numeric(), 
                            stringsAsFactors = FALSE)
    I <- length(cLine_selected)
    for (i in 1:I) {
        pol <- cLine_selected[i]
        pol.x <- as.numeric(sapply(pol, function(x) x$x))
        pol.y <- as.numeric(sapply(pol, function(x) x$y))
        index <- point.in.polygon(z_melted$EX, z_melted$EM, pol.x, pol.y) > 0
        if (sum(index) == 0){ 
            # when there is no real point within polygon, give NA value
            local_max[i,] <- NA
        } else {
            MAX <- max(z_melted$value[index])
            MAX.index <- which(z_melted$value == MAX)
            local_max[i,] <- z_melted[MAX.index,]
        }
    }
    row.names(local_max) <- NULL
    if (showprint) print(local_max)
    peak_names <- paste0("EX", local_max$EX, "EM", local_max$EM) 
    return(peak_names)
}