R/refinement.R

In rhondabacher/methylscaper: Visualization of Methylation Data

#' Refinement
#' 
#' Reorders a subset of the methylation data. 
#' 
#' @param orderObject An object of class \code{orderObject}, generated 
#'  with the \code{initialOrder} function.
#' @param refineStart The index of the first sample (row) used 
#'  in the refinement.
#' @param refineEnd The index of the last sample (row) used in the refinement.
#' @param Method The seriation method used to perform the refinement.
#' 
#' @return The refinement reorders the cells/molecules (rows) between the indicated
#' start and end positions. The function returns the new complete ordering 
#' with the refinement applied.
#' @importFrom seriation seriate get_order
#' @export
#'
#' @examples 
#'  
#' data(singlemolecule_example)
#' 
#' orderObj <- initialOrder(singlemolecule_example, Method = "PCA")
#' # reordering the first 50 cells/molecules (rows)
#' orderObj$order1 <- refineFunction(orderObj, 1, 50) 

refineFunction <- function(orderObject, refineStart, refineEnd, 
                                    Method="PCA") {
  
  toClust <- orderObject$toClust
  order1 <- orderObject$order1
  
  if (refineEnd > max(order1) | refineStart > max(order1)) {
     message(paste0("Refine parameters are out of bounds. The maximum has been set to ", refineEnd))
     if (refineEnd > max(order1)) refineEnd <- max(order1)
     if (refineStart > max(order1)) refineStart <- max(order1)  
  }
  if (refineEnd < min(order1) | refineStart < min(order1)) {
     message(paste0("Refine parameters are out of bounds. The maximum has been set to ", refineEnd))
     if (refineEnd < min(order1)) refineEnd <- min(order1)
     if (refineStart < min(order1)) refineStart <- min(order1)  
  }
  
 
  
  toRefine_order <- order1[refineStart:refineEnd]
  
  
  toRefine_clust <- toClust[toRefine_order,]
  
  if (Method=="PCA") {
    col_centered <- apply(toRefine_clust, 2, function(x) x - mean(x))
    try1 <- svd(col_centered, nu = 1, nv = 0)
    order_new <- order(try1$u[,1])
  } else { # Methods available for refining are: ARSA, HC_complete, 
      ## HC_average, HC_ward.
    if (is.null(orderObject$distMat)) distMat <- dist(toRefine_clust,
                                                    method = "euclidean")
    else {
        distMat <- as.dist(as.matrix(orderObject$distMat)[toRefine_order,toRefine_order])
    }
    order_new <- seriate(distMat, method=Method, verbose=FALSE)
    order_new <- get_order(order_new)
  }
  
  
  # New order:
  order_new <- order1[seq(refineStart,refineEnd)][order_new]
  order_final <- order1
  order_final[refineStart:refineEnd] <- order_new
  
  return(order_final)
}

#' Force reversal of a subset of the ordering
#' 
#' This reverses a subset of the ordering, as determined by the user. 
#' By default, the entire ordering is reversed.
#' 
#' @param orderObject An object of class \code{orderObject}, 
#'  generated with the \code{initialOrder} function.
#' @param reverseStart The first index to be included in the reversal.
#' @param reverseEnd The last index to be included in the reversal.
#' 
#' @return The new complete ordering, with the reversal applied.
#' @export
#' @examples 
#'  
#' data(singlemolecule_example)
#' 
#' orderObj <- initialOrder(singlemolecule_example, Method = "PCA")
#' # reorder first 50 cells/molecules (rows)
#' orderObj$order1 <- refineFunction(orderObj, 1, 50) 
#' orderObj$order1 <- forceReverse(orderObj, 1, 50) 

forceReverse <- function(orderObject, reverseStart = 1, 
                        reverseEnd = length(orderObject$order1))
{
  order1 <- orderObject$order1
  order1[reverseStart:reverseEnd] <- rev(order1[reverseStart:reverseEnd])
  return(order1)
}