3_buildMST.R

#' BuildMST
#' 
#' Build Minimal Spanning Tree
#' 
#' Add minimal spanning tree description to the FlowSOM object
#' 
#' @param fsom   FlowSOM object, as generated by \code{\link{BuildSOM}}
#' @param silent If \code{TRUE}, no progress updates will be printed
#' @param tSNE   If \code{TRUE}, an alternative tSNE layout is computed as well
#' 
#' @return FlowSOM object containing MST description
#' 
#' @seealso \code{\link{BuildSOM}}, \code{\link{PlotStars}}
#' 
#' @examples 
#' # Read from file, build self-organizing map
#' fileName <- system.file("extdata", "68983.fcs", package="FlowSOM")
#' flowSOM.res <- ReadInput(fileName, compensate=TRUE, transform = TRUE,
#'                          scale = TRUE)
#' flowSOM.res <- BuildSOM(flowSOM.res, colsToUse = c(9, 12, 14:18))
#' 
#' # Build the Minimal Spanning Tree
#' flowSOM.res <- BuildMST(flowSOM.res)
#' 
#' @importFrom Rtsne Rtsne
#' 
#' @export
BuildMST <- function(fsom, silent = FALSE, tSNE = FALSE){
  
  fsom$MST <- list()
  if(!silent) message("Building MST\n")
  
  adjacency <- as.matrix(stats::dist(fsom$map$codes, method = "euclidean"))
  fullGraph <- igraph::graph.adjacency(adjacency,
                                       mode = "undirected",
                                       weighted = TRUE)
  fsom$MST$graph <- igraph::minimum.spanning.tree(fullGraph)
  ws <- igraph::edge.attributes(fsom$MST$graph)$weight
  #normalize edge weights to match the grid size in coords (see below)
  ws <- ws / mean(ws)
  igraph::edge.attributes(fsom$MST$graph)$weight <- ws
  fsom$MST$l <- igraph::layout.kamada.kawai(
    coords = as.matrix(fsom$map$grid),
    fsom$MST$graph)
  
  
  if(tSNE){
    fsom$MST$l2 <- Rtsne::Rtsne(fsom$map$codes)$Y   
    #library(RDRToolbox)
    #fsom$MST$l2 <- Isomap(fsom$map$codes, dims = 2, k = 3)[[1]]
  }
  
  return(fsom)
}

#' FlowSOMSubset 
#' 
#' FlowSOM subset
#' 
#' Take a subset from a FlowSOM object
#' 
#' @param fsom FlowSOM object, as generated by \code{\link{BuildMST}}
#' @param ids  Array containing the ids to keep
#' 
#' @return FlowSOM object containg updated data and medianvalues, 
#'    but with the same grid
#' @seealso \code{\link{BuildMST}}
#' 
#' @examples
#'    # Read two files (Artificially, as we just split 1 file in 2 subsets)
#'    fileName <- system.file("extdata", "68983.fcs", package = "FlowSOM")
#'    ff1 <- flowCore::read.FCS(fileName)[1:1000, ]
#'    flowCore::keyword(ff1)[["FIL"]] <- "File1"
#'    ff2 <- flowCore::read.FCS(fileName)[1001:2000, ]
#'    flowCore::keyword(ff2)[["FIL"]] <- "File2"
#'    
#'    flowSOM.res <- FlowSOM(flowCore::flowSet(c(ff1, ff2)), compensate = TRUE,
#'                           transform = TRUE, scale = TRUE,
#'                           colsToUse = c(9, 12, 14:18), maxMeta = 10)
#'    
#'    # see $metadata for subsets:
#'    flowSOM.res$metaData
#'    
#'    # Use only the second file, without changing the map
#'    fSOM2 <- FlowSOMSubset(flowSOM.res,
#'                           (flowSOM.res$metaData[[2]][1]):
#'                            (flowSOM.res$metaData[[2]][2]))
#'
#' @export
FlowSOMSubset <- function(fsom, ids){
  fsom_tmp <- fsom
  fsom_tmp$data <- fsom$data[ids, , drop = FALSE]
  fsom_tmp$map$mapping <- fsom$map$mapping[ids, , drop = FALSE]
  fsom_tmp <- UpdateDerivedValues(fsom_tmp)
  return(fsom_tmp)
}

#' NewData 
#' 
#' Map new data to a FlowSOM grid
#'
#' New data is mapped to an existing FlowSOM object. The input is similar to the
#' readInput function.
#' A new FlowSOM object is created, with the same grid, but a new
#' mapping, node sizes and mean values. The same preprocessing steps
#' (compensation, tranformation and scaling) will happen to this file as was 
#' specified in the original FlowSOM call. The scaling parameters from the 
#' original grid will be used.
#'
#' @param fsom          FlowSOM object
#' @param input         A flowFrame, a flowSet or an array of paths to files 
#'                      or directories   
#' @param madAllowed   A warning is generated if the distance of the new
#'                      data points to their closest cluster center is too
#'                      big. This is computed based on the typical distance
#'                      of the points from the original dataset assigned to
#'                      that cluster, the threshold being set to
#'                      median + madAllowed * MAD. Default is 4.
#' @param compensate    logical, does the data need to be compensated. If NULL,
#'                      the same value as in the original FlowSOM call will be 
#'                      used.
#' @param spillover     spillover matrix to compensate with. If NULL,
#'                      the same value as in the original FlowSOM call will be 
#'                      used.
#' @param transform     logical, does the data need to be transformed. If NULL,
#'                      the same value as in the original FlowSOM call will be 
#'                      used.
#' @param toTransform   column names or indices that need to be transformed.
#'                      If NULL,  the same value as in the original FlowSOM 
#'                      call will be used.
#' @param transformFunction  If NULL, the same value as in the original FlowSOM 
#'                      call will be used.
#' @param transformList  If NULL, the same value as in the original FlowSOM 
#'                      call will be used.
#' @param scale         Logical, does the data needs to be rescaled. If NULL, 
#'                      the same value as in the original FlowSOM call will be 
#'                      used.
#' @param scaled.center See \code{\link{scale}}. If NULL, the same value as in 
#'                      the original FlowSOM call will be used.
#' @param scaled.scale  See \code{\link{scale}}. If NULL, the same value as in 
#'                      the original FlowSOM call will be used.
#' @param silent        Logical. If \code{TRUE}, print progress messages. 
#'                      Default = \code{FALSE}.
#'        
#' @return A new FlowSOM object
#' @seealso \code{\link{FlowSOMSubset}} if you want to get a subset of the
#'          current data instead of a new dataset
#' @examples 
#'  # Build FlowSom result
#'  fileName <- system.file("extdata", "68983.fcs", package = "FlowSOM")
#'  ff <- flowCore::read.FCS(fileName)
#'  ff <- flowCore::compensate(ff, flowCore::keyword(ff)[["SPILL"]])
#'  ff <- flowCore::transform(ff,
#'          flowCore::transformList(colnames(flowCore::keyword(ff)[["SPILL"]]),
#'                                 flowCore::logicleTransform()))
#'    flowSOM.res <- FlowSOM(ff[1:1000, ], 
#'                           scale = TRUE, 
#'                           colsToUse = c(9, 12, 14:18),
#'                           nClus = 10)
#'    
#'    # Map new data
#'    fSOM2 <- NewData(flowSOM.res, ff[1001:2000, ])
#'
#' @export
NewData <- function(fsom, 
                    input,
                    madAllowed = 4,
                    compensate = NULL, 
                    spillover = NULL,
                    transform = NULL, 
                    toTransform = NULL, 
                    transformFunction = NULL,
                    transformList = NULL,
                    scale = NULL, 
                    scaled.center = NULL, 
                    scaled.scale = NULL,
                    silent = FALSE){
  fsom <- UpdateFlowSOM(fsom)
  if(is.null(compensate)){
    compensate <- fsom$compensate
  }
  if(is.null(spillover)){
    spillover <- fsom$spillover
  }
  if(is.null(transform)){
    transform <- fsom$transform
  }
  if(is.null(toTransform)){
    toTransform <- fsom$toTransform
  }
  if(is.null(transformFunction)){
    transformFunction <- fsom$transformFunction
  }
  if(is.null(transformList)){
    transformList <- fsom$transformList
  }
  if(is.null(scale)){
    scale <- fsom$scale
  }
  if(is.null(scaled.center)){
    scaled.center <- fsom$scaled.center
  }
  if(is.null(scaled.scale)){
    scaled.scale <- fsom$scaled.scale
  }
  
  fsom_new <- ReadInput(input, 
                        compensate = compensate, spillover = spillover,
                        transform = transform, toTransform = toTransform,
                        transformFunction = transformFunction, 
                        transformList = transformList,
                        scale = scale, scaled.center = scaled.center,
                        scaled.scale = scaled.scale, silent = silent)
  
  fsom_new$map <- fsom$map
  fsom_new$MST <- fsom$MST
  
  fsom_new$metaclustering <- fsom$metaclustering
  fsom_new$map$mapping <- MapDataToCodes(fsom$map$codes, fsom_new$data)
  fsom_new <- UpdateDerivedValues(fsom_new)
  
  test_outliers <- TestOutliers(fsom_new,
                                madAllowed = madAllowed,
                                fsomReference = fsom)
  max_outliers <- max(test_outliers$Number_of_outliers) 
  n_outliers <- sum(test_outliers$Number_of_outliers) 
  if(max_outliers > 100){
    warning(n_outliers, 
            " cells (",
            round(n_outliers / nrow(fsom_new$data) * 100, 2),
            "%) seem far from their cluster centers.")
  }
  
  fsom_new$outliers <- test_outliers
  return(fsom_new)
}

#' TestOutliers 
#' 
#' Test if any cells are too far from their cluster centers
#'
#' For every cluster, the distance from the cells to the cluster centers is
#' used to label cells which deviate too far as outliers. The threshold is
#' chosen as the median distance + \code{madAllowed} times the median absolute
#' deviation of the distances. 
#'
#' @param fsom  FlowSOM object
#' @param madAllowed Number of median absolute deviations allowed. Default = 4.
#' @param fsomReference FlowSOM object to use as reference. If NULL (default),
#'                       the original fsom object is used.
#' @param plotFile If \code{NULL} (default), no plot will be created. If a 
#'                 filepath is given for a pdf, the plot will be written in the
#'                 corresponding file
#'        
#' @return A new FlowSOM object
#' @seealso \code{\link{FlowSOMSubset}} if you want to get a subset of the
#'          current data instead of a new dataset
#' @examples 
#'  # Build FlowSom result
#'  fileName <- system.file("extdata", "68983.fcs", package = "FlowSOM")
#'  ff <- flowCore::read.FCS(fileName)
#'  flowSOM.res <- FlowSOM(ff,
#'                         compensate = TRUE, transform = TRUE, scale = TRUE,
#'                         colsToUse = c(9, 12, 14:18),
#'                         maxMeta = 10)
#'    
#'  # Map new data
#'  outlier_report <- TestOutliers(flowSOM.res)
#'  
#' @import     ggplot2
#' @importFrom grDevices pdf dev.off
#' @importFrom ggpubr ggarrange
#'
#' @export
TestOutliers <- function(fsom, 
                         madAllowed = 4,
                         fsomReference = NULL,
                         plotFile = NULL){
  
  fsom <- UpdateFlowSOM(fsom)
  
  if(is.null(fsomReference)){
    fsomReference <- fsom
  } else {
    fsomReference <- UpdateFlowSOM(fsomReference)
  }
  
  distances_median <- sapply(seq_len(fsomReference$map$nNodes),
                             function(x){
                               ids <- which(GetClusters(fsomReference) == x)
                               if(length(ids) > 0){
                                 m <- stats::median(
                                   fsomReference$map$mapping[ids, 2])
                               } else {
                                 m <- 0
                               }
                               return(m)
                             })
  
  distances_mad <- sapply(seq_len(fsomReference$map$nNodes),
                          function(x){
                            ids <- which(GetClusters(fsomReference) == x)
                            if(length(ids) > 0){
                              m <- stats::mad(
                                fsomReference$map$mapping[ids, 2])
                            } else {
                              m <- 0
                            }
                            return(m)
                          })
  
  thresholds <- distances_median + madAllowed * distances_mad
  
  max_distances_new <- sapply(seq_len(fsom$map$nNodes),
                              function(x){
                                ids <- which(GetClusters(fsom) == x)
                                if(length(ids) > 0){
                                  m <- max(fsom$map$mapping[ids, 2])
                                } else {
                                  m <- 0
                                }
                                return(m)
                              })
  
  outliers <- sapply(seq_len(fsom$map$nNodes),
                     function(x){
                       ids <- which(GetClusters(fsom) == x)
                       distances <- fsom$map$mapping[ids, 2]
                       return(sum(distances > thresholds[x]))
                     })
  if(!is.null(plotFile)){
    xdim <- fsom$map$xdim
    ydim <- fsom$map$ydim
    graphics::layout(matrix(1:(xdim*ydim), nrow = xdim))
    plotList <- lapply(seq_len(xdim * ydim), function(i){
      ids <- which(GetClusters(fsom) == i)
      values <- fsom$map$mapping[ids, 2]
      if(length(values) > 1){
        nOutliers <- sum(values > thresholds[i])
        p <- suppressMessages(ggplot2::ggplot() +
                                ggplot2::geom_histogram(ggplot2::aes(values), fill = "grey90",
                                                        col = "black", size = 0.2) +
                                ggplot2::geom_vline(ggplot2::aes(xintercept = distances_median[i]),
                                                    col = "black") +
                                ggplot2::geom_vline(ggplot2::aes(xintercept = thresholds[i]), 
                                                    col = "red") +
                                ggplot2::ggtitle(paste0(i, " (", nOutliers, ")")) +
                                ggplot2::theme_minimal() +
                                ggplot2::theme(axis.title.x = ggplot2::element_blank()) +
                                ggplot2::ylab("Frequency"))
        return(p)
      }
    })
    p <- suppressMessages(ggpubr::ggarrange(plotlist = plotList))
    grDevices::pdf(plotFile, width = 20, height = 20)
    print(p)
    grDevices::dev.off()
  }
  result <- data.frame(
    "Median_distance" = distances_median, 
    "Median_absolute_deviation" = distances_mad, 
    "Threshold" = thresholds, 
    "Number_of_outliers" = outliers,
    "Maximum_outlier_distance" = max_distances_new)[outliers > 0, ]
  
  result <- result[order(outliers[outliers > 0], decreasing = TRUE), ]
  
  return(result)
}

#' NClusters
#' 
#' Extracts the number of clusters from a FlowSOM object
#' 
#' @param fsom FlowSOM object
#' 
#' @return The number of clusters
#' 
#' @examples 
#'  # Build FlowSom result
#'  fileName <- system.file("extdata", "68983.fcs", package = "FlowSOM")
#'  ff <- flowCore::read.FCS(fileName)
#'  flowSOM.res <- FlowSOM(ff,
#'                         compensate = TRUE, transform = TRUE, scale = TRUE,
#'                         colsToUse = c(9, 12, 14:18),
#'                         maxMeta = 10)
#'  NClusters(flowSOM.res)
#' 
#' @export
NClusters <- function(fsom){
  fsom <- UpdateFlowSOM(fsom)
  return(fsom$map$nNodes)
}

#' NMetaclusters
#' 
#' Extracts the number of metaclusters from a FlowSOM object
#' 
#' @param fsom FlowSOM object
#' 
#' @return The number of metaclusters
#' 
#' @examples 
#'  # Build FlowSom result
#'  fileName <- system.file("extdata", "68983.fcs", package = "FlowSOM")
#'  ff <- flowCore::read.FCS(fileName)
#'  flowSOM.res <- FlowSOM(ff,
#'                         compensate = TRUE, transform = TRUE, scale = TRUE,
#'                         colsToUse = c(9, 12, 14:18),
#'                         maxMeta = 10)
#'  NMetaclusters(flowSOM.res)
#' @export
NMetaclusters <- function(fsom){
  fsom <- UpdateFlowSOM(fsom)
  nMetaclusters <- fsom$map$nMetaclusters
  return(nMetaclusters)
}