R/3_buildMST.R

Defines functions QueryStarPlot PlotGroups TestOutliers NewData FlowSOMSubset PlotCenters PlotOverview2D PlotClusters2D PlotNode PlotStars plotStarQuery plotStarLegend shiftFunction legendContinuous mystar PlotPies PlotLabels PlotNumbers PlotSD PlotVariable PlotMarker PlotBackgroundLegend computeBackgroundColor UpdateNodeSize BuildMST

Documented in BuildMST computeBackgroundColor FlowSOMSubset NewData PlotCenters PlotClusters2D PlotGroups PlotLabels PlotMarker PlotNode PlotNumbers PlotOverview2D PlotPies PlotSD plotStarLegend PlotStars PlotVariable QueryStarPlot TestOutliers UpdateNodeSize

#' 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)
#' 
#' @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 <- tsne(fsom$map$codes)   
        #library(RDRToolbox)
        #fsom$MST$l2 <- Isomap(fsom$map$codes,dims=2,k=3)[[1]]
    }
    
    UpdateNodeSize(fsom)
}

####################
## UpdateNodeSize ##
####################
#' Update nodesize of FlowSOM object
#' 
#' Add size property to the graph based on cellcount for each node
#' 
#' @param fsom     FlowSOM object, as generated by \code{\link{BuildMST}}
#' @param count    Absolute cell count of the sample
#' @param reset    Logical. If \code{TRUE}, all nodes get the same size
#' @param transform Transformation function. Use e.g. square root to let counts
#'                 correspond with area of node instead of radius
#' @param maxNodeSize Maximum node size after rescaling. Default: 15
#' @param shift    Shift of the counts, defaults to 0
#' @param scale    Scaling of the counts, defaults to the maximum of the value
#'                 minus the shift. With shift and scale set as default, the 
#'                 largest node will be maxNodesize and an empty node will have
#'                 size 0
#' 
#' @return Updated FlowSOM object
#' @seealso \code{\link{BuildMST}}
#' 
#' @examples 
#' # Read from file, build self-organizing map and minimal spanning tree
#' 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))
#' flowSOM.res <- BuildMST(flowSOM.res)
#'
#' # Give all nodes same size
#' flowSOM.res <- UpdateNodeSize(flowSOM.res,reset=TRUE)
#' PlotStars(flowSOM.res)
#'
#' # Node sizes relative to amount of cells assigned to the node
#' flowSOM.res <- UpdateNodeSize(flowSOM.res)
#' PlotStars(flowSOM.res)
#' 
#' @export
UpdateNodeSize <- function(fsom, count = NULL, reset=FALSE, transform=sqrt,
                            maxNodeSize = 15,
                            shift = 0, scale = NULL){
    if(reset){
        fsom$MST$size <- rep(maxNodeSize, nrow(fsom$map$grid))
    } else {
        t <- rep(0, fsom$map$nNodes)
        names(t) <- as.character(seq_len(fsom$map$nNodes))
        t_tmp <- table(fsom$map$mapping[, 1]) / nrow(fsom$map$mapping)
        t[names(t_tmp)] <- t_tmp
        
        if(!is.null(count)){
            t <- (t/sum(t)) * count
        }
        
        if(!is.null(transform)){
            t <- transform(t)
        }
        
        if(is.null(shift)) shift <- min(t)
        if(is.null(scale)) scale <- max(t - shift)
        rescaled <- maxNodeSize * (t - shift)/scale
        rescaled[rescaled == 0] <- 0.0001
        fsom$MST$size <- rescaled   
    }
    fsom
}

############################
## computeBackgroundColor ##
############################
#' Internal function for computing background nodes
#' @param backgroundValues  Values to be used for background coloring, either
#'                          numerical values or something that can be made into
#'                          a factor (e.g. a clustering)
#' @param backgroundColor   Colorpalette to be used for the background coloring
#'                          . Can be either a function or an array specifying
#'                          colors
#' @param backgroundLim     Only used when backgroundValues are numerical. 
#'                          Defaults to min and max of the backgroundValues.
#' @param backgroundBreaks  Breaks to pass on to \code{\link{cut}}, to split
#'                          numerical background values. If NULL, the length of
#'                          backgroundColor will be used (default 100).
computeBackgroundColor <- function(backgroundValues, 
                                    backgroundColor,
                                    backgroundLim = NULL,
                                    backgroundBreaks = NULL){
    # Choose background colour
    backgroundList <- list()
    backgroundColors <- NULL
    if(!is.null(backgroundValues)){
        if(is.numeric(backgroundValues)){
            backgroundList <- as.list(seq_along(backgroundValues))
            
            if(class(backgroundColor)=="function" & 
               !is.null(backgroundBreaks) & 
               length(backgroundBreaks)>1)
            {
                backgroundColor <- backgroundColor(length(backgroundBreaks))
            } else if (class(backgroundColor)=="function"){
                backgroundColor <- backgroundColor(100)
                backgroundBreaks <- length(backgroundColor)
            } else if (is.null(backgroundBreaks)){
                backgroundBreaks <- length(backgroundColor)
            }
            
            if(length(backgroundLim) > 0){
                ids <- cut(c(backgroundLim,backgroundValues),
                           backgroundBreaks
                )[-c(seq_along(backgroundLim))]
            } else {
                ids <- cut(backgroundValues,
                           backgroundBreaks)
            }
            backgroundValues <- ids
#             backgroundColors <- backgroundColor[ids]    
        } else {
            if(! is.factor(backgroundValues)){
                backgroundValues <- as.factor(backgroundValues)
            }
            
            backgroundList <- as.list(seq_along(backgroundValues))
            
            if(class(backgroundColor)=="function"){
                backgroundColor <- backgroundColor(
                    length(levels(backgroundValues)))
            }
            
            if(length(backgroundColor) < length(levels(backgroundValues))){
                stop("You specified less backgroundcolours than groups.")
            }
            
        }    
    }
    backgroundColors <- backgroundColor[backgroundValues]    

    list(values=backgroundValues,
         col=backgroundColor,
         groups=backgroundList)
}

##########################
## PlotBackgroundLegend ##
##########################
PlotBackgroundLegend <- function(backgroundValues, background, 
                                 main="Background"){
    graphics::plot.new()
    if(is.numeric(backgroundValues)) {
        legendContinuous(background$col,
                         as.numeric(gsub(".*,","",
                                         gsub("].*","",
                                              levels(background$values)))))
    } else {
        graphics::legend("center", legend=levels(background$values),
               fill=background$col, 
               cex=0.7, 
               ncol =  ceiling(length(levels(background$values)) / 10),
               bty="n",
               title=main)       
    }
}

################
## PlotMarker ##
################
#' Plot marker values
#' 
#' Plot FlowSOM grid or tree, coloured by node values for a specific marker
#' 
#' @param fsom         FlowSOM object, as generated by \code{\link{BuildMST}}
#' @param marker       Name or index of marker to plot
#' @param view     Preferred view, options: "MST" (default), "grid" or 
#'                  "tSNE" (if this option was selected while building the MST)
#' @param main         Title of the plot
#' @param colorPalette Color palette to use
#' 
#' @param backgroundValues  Values to be used for background coloring, either
#'                          numerical values or something that can be made into
#'                          a factor (e.g. a clustering)
#' @param backgroundColor   Colorpalette to be used for the background coloring
#'                          . Can be either a function or an array specifying
#'                          colors
#' @param backgroundLim     Only used when backgroundValues are numerical. 
#'                          Defaults to min and max of the backgroundValues.
#' @param backgroundBreaks  Breaks to pass on to \code{\link{cut}}, to split
#'                          numerical background values. If NULL, the length of
#'                          backgroundColor will be used (default 100).
#'                          
#' @return Nothing is returned. A plot is drawn in which each node 
#'         is coloured depending on its median value for the given marker
#' @references This visualization technique resembles SPADE results. 
#' M. Linderman, P. Qiu, E. Simonds and Z. Bjornson (). spade: SPADE -- An 
#' analysis and visualization tool for Flow Cytometry. R package version 1.12.2.
#' http://cytospade.org
#' @seealso \code{\link{PlotStars}},\code{\link{PlotPies}},
#'          \code{\link{PlotCenters}},\code{\link{BuildMST}}
#' 
#' @examples 
#' # Read from file, build self-organizing map and minimal spanning tree
#' 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))
#' flowSOM.res <- BuildMST(flowSOM.res)
#' 
#' # Plot one marker
#' PlotMarker(flowSOM.res,"FSC-A")
#' 
#' @export
PlotMarker <- function(fsom, marker=NULL, view="MST",main=NULL, 
                        colorPalette=grDevices::colorRampPalette(c("#00007F",
                            "blue","#007FFF", "cyan", "#7FFF7F", "yellow",
                            "#FF7F00","red", "#7F0000")),
                        backgroundValues = NULL,
                        backgroundColor = function(n){
                            grDevices::rainbow(n,alpha=0.3)},
                        backgroundBreaks = NULL,
                        backgroundLim = NULL){
    switch(view,
           MST  = { 
             layout <- fsom$MST$l 
             lty <- 1
           }, grid = { 
               layout <- as.matrix(fsom$map$grid)
               lty <- 0
           }, tSNE = { 
             layout <- fsom$MST$l2
             lty <- 0
           }, stop("The view should be MST, grid or tSNE. tSNE will only work
                   if you specified this when building the MST.")
    )
    
    # Choose background colour
    if(!is.null(backgroundValues)){
        background <- computeBackgroundColor(backgroundValues,backgroundColor,
                                             backgroundLim, backgroundBreaks)
    } else {
        background <- NULL
    }
    
    oldpar <- graphics::par(no.readonly = TRUE)
    
    if(is.null(main)) main <- marker
    
    if(is.null(marker)){
        igraph::plot.igraph(fsom$MST$graph, layout=layout, 
            vertex.size=fsom$MST$size, 
            vertex.label=NA,edge.lty=lty)
    }else{
        f <- fsom
        igraph::V(f$MST$graph)$color <- colorPalette(100)[as.numeric(cut(
            fsom$map$medianValues[, marker], breaks = 100))] 
        igraph::plot.igraph(f$MST$graph, 
            layout=layout, 
            vertex.size=fsom$MST$size, 
            vertex.label=NA, 
            main=main,
            edge.lty=lty,
            mark.groups=background$groups, 
            mark.col=background$col[background$values], 
            mark.border=background$col[background$values])
        
        graphics::par(fig=c(0,0.2,0,1),mar=c(0,0,0,0),new=TRUE)
        legendContinuous(colorPalette(100),
                   fsom$map$medianValues[, marker])
    }
    graphics::par(oldpar)
}

##################
## PlotVariable ##
##################
#' Plot a variable for all nodes
#' 
#' Plot FlowSOM grid or tree, coloured by node values given in variable
#' 
#' @param fsom         FlowSOM object, as generated by \code{\link{BuildMST}}
#' @param variable     Vector containing a value for each node
#' @param view         Preferred view, options: "MST", "grid" or "tSNE" (if this
#'                     option was selected while building the MST)
#' @param main         Title of the plot
#' @param colorPalette Color palette to use
#' @param symmetric    Plot colours symmetric around zero
#' @param lim          Variable limits
#' 
#' @param backgroundValues  Values to be used for background coloring, either
#'                          numerical values or something that can be made into
#'                          a factor (e.g. a clustering)
#' @param backgroundColor   Colorpalette to be used for the background coloring
#'                          . Can be either a function or an array specifying
#'                          colors
#' @param backgroundLim     Only used when backgroundValues are numerical. 
#'                          Defaults to min and max of the backgroundValues.
#' @param backgroundBreaks  Breaks to pass on to \code{\link{cut}}, to split
#'                          numerical background values. If NULL, the length of
#'                          backgroundColor will be used (default 100).
#' 
#' @return Nothing is returned. A plot is drawn in which each node 
#'         is coloured depending on its value for the given variable
#' @seealso \code{\link{PlotMarker}},\code{\link{PlotStars}},
#'          \code{\link{PlotPies}},\code{\link{PlotCenters}},
#'          \code{\link{BuildMST}}
#' 
#' @examples 
#' # Read from file, build self-organizing map and minimal spanning tree
#' 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))
#' flowSOM.res <- BuildMST(flowSOM.res)
#' 
#' # Plot some random values
#' rand <- runif(flowSOM.res$map$nNodes)
#' PlotVariable(flowSOM.res,rand)
#' 
#' @export
PlotVariable <- function(fsom, variable, view="MST",main=NULL, 
                        colorPalette=grDevices::colorRampPalette(c("#00007F",
                        "blue", "#007FFF", "cyan", "#7FFF7F", "yellow",
                        "#FF7F00","red", "#7F0000")),symmetric=FALSE,
                        lim=NULL,
                        backgroundValues = NULL,
                        backgroundColor = function(n){
                            grDevices::rainbow(n,alpha=0.3)},
                        backgroundLim = NULL,
                        backgroundBreaks = NULL){
    switch(view,
            MST  = { layout <- fsom$MST$l 
            lty <- 1},
            grid = { layout <- as.matrix(fsom$map$grid)
            lty <- 0},
            tSNE = { layout <- fsom$MST$l2
            lty <- 0}, 
           stop("The view should be MST, grid or tSNE. tSNE will only work
                   if you specified this when building the MST.")
    )
    
    # Choose background colour
    if(!is.null(backgroundValues)){
        background <- computeBackgroundColor(backgroundValues,backgroundColor,
                                             backgroundLim, backgroundBreaks)
    } else {
        background <- NULL
    }
    
    oldpar <- graphics::par(no.readonly = TRUE)
    graphics::par(mar=c(1,1,1,1))
    graphics::layout(matrix(c(1,2), 1, 2, byrow = TRUE), 
            widths=c(1,2), heights=c(1))
    graphics::plot.new()
    
    if(is.null(lim)){
        if(symmetric){
            lim <- c(-max(abs(variable)),max(abs(variable)))
        } else {
            lim <- c(min(variable),max(variable))
        }
    }
    
    legendContinuous(colorPalette(100),lim)
    
    
    f <- fsom

    igraph::V(f$MST$graph)$color <- colorPalette(100)[as.numeric(cut(
        c(lim,variable), 
        breaks = 100))[-c(1,2)]] 
    
    igraph::plot.igraph(f$MST$graph, 
        layout=layout, 
        vertex.size=fsom$MST$size, 
        vertex.label=NA, 
        main=main,
        edge.lty=lty,
        mark.groups=background$groups, 
        mark.col=background$col[background$values], 
        mark.border=background$col[background$values])
    
    
    graphics::layout(1)
    graphics::par(oldpar)
}





#' Plot SD
#' 
#' --- Function in development, use with caution ---
#' Plot FlowSOM grid or tree, coloured by standard deviaton
#' 
#' From suggestion in email:
#' I am currently considering a way to summarize for each node all the SD as one value. 
#' After computing the SD matrix (nrow = # nodes, ncol = # markers), I compute the median value per column, 
#' then divide the SD matrix by it, and finally take the maximum ratio of each line (aka node). 
#' Doing so I got a unique dispersion score per node. 
#' 
#' @param fsom         FlowSOM object, as generated by \code{\link{BuildMST}}
#' @param marker       If a marker is given, the sd for this marker is shown.
#'                     Otherwise, the maximum ratio is used.
#' @param view         Preferred view, options: "MST", "grid" or "tSNE" (if this
#'                     option was selected while building the MST)
#' @param main         Title of the plot
#' @param colorPalette Color palette to use
#' @param symmetric    Plot colours symmetric around zero
#' @param lim          Variable limits
#' 
#' @param backgroundValues  Values to be used for background coloring, either
#'                          numerical values or something that can be made into
#'                          a factor (e.g. a clustering)
#' @param backgroundColor   Colorpalette to be used for the background coloring
#'                          . Can be either a function or an array specifying
#'                          colors
#' @param backgroundLim     Only used when backgroundValues are numerical. 
#'                          Defaults to min and max of the backgroundValues.
#' @param backgroundBreaks  Breaks to pass on to \code{\link{cut}}, to split
#'                          numerical background values. If NULL, the length of
#'                          backgroundColor will be used (default 100).
#' 
#' @return Nothing is returned. A plot is drawn in which each node 
#'         is coloured depending on its standard deviation
#' @seealso \code{\link{PlotMarker}},\code{\link{PlotStars}},
#'          \code{\link{PlotPies}},\code{\link{PlotCenters}},
#'          \code{\link{BuildMST}}
#' 
#' @examples 
#' # Read from file, build self-organizing map and minimal spanning tree
#' 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))
#' flowSOM.res <- BuildMST(flowSOM.res)
#' 
#' PlotSD(flowSOM.res)
#' 
#' @export
PlotSD <- function(fsom, 
                   marker= NULL,
                   view="MST", main=NULL, 
                   colorPalette=grDevices::colorRampPalette(c(
                     "#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow",
                     "#FF7F00","red", "#7F0000")),symmetric=FALSE,
                   lim=NULL,
                   backgroundValues = NULL,
                   backgroundColor = function(n){
                     grDevices::rainbow(n,alpha=0.3)},
                   backgroundLim = NULL,
                   backgroundBreaks = NULL){
  switch(view,
         MST  = { layout <- fsom$MST$l 
         lty <- 1},
         grid = { layout <- as.matrix(fsom$map$grid)
         lty <- 0},
         tSNE = { layout <- fsom$MST$l2
         lty <- 0}, 
         stop("The view should be MST, grid or tSNE. tSNE will only work
                   if you specified this when building the MST.")
  )
  
  # Choose background colour
  if(!is.null(backgroundValues)){
    background <- computeBackgroundColor(backgroundValues,backgroundColor,
                                         backgroundLim, backgroundBreaks)
  } else {
    background <- NULL
  }
  
  stdevs <- apply(fsom$data[, fsom$map$colsUsed], 
                  2, 
                  function(x){
                    tapply(x, fsom$map$mapping[, 1], stats::sd)
                  })
  stdev_medians <- apply(stdevs, 2, stats::median)
  if(is.null(marker)){
    variable <- apply(stdevs, 1, function(x){max(x / stdev_medians)})
  } else{
    variable <- stdevs[,marker] / stdev_medians[marker]
  }
  oldpar <- graphics::par(no.readonly = TRUE)
  graphics::par(mar=c(1,1,1,1))
  graphics::layout(matrix(c(1,2), 1, 2, byrow = TRUE), 
                   widths=c(1,2), heights=c(1))
  graphics::plot.new()
  
  if(is.null(lim)){
    if(symmetric){
      lim <- c(-max(abs(variable)),max(abs(variable)))
    } else {
      lim <- c(min(variable),max(variable))
    }
  }
  
  legendContinuous(colorPalette(100),lim)
  
  
  f <- fsom
  
  igraph::V(f$MST$graph)$color <- colorPalette(100)[as.numeric(cut(
    c(lim,variable), 
    breaks = 100))[-c(1,2)]] 
  
  igraph::plot.igraph(f$MST$graph, 
                      layout=layout, 
                      vertex.size=fsom$MST$size, 
                      vertex.label=NA, 
                      main=main,
                      edge.lty=lty,
                      mark.groups=background$groups, 
                      mark.col=background$col[background$values], 
                      mark.border=background$col[background$values])
  
  
  graphics::layout(1)
  graphics::par(oldpar)
}


                                                                                                                      
#################
## PlotNumbers ##
#################
#' Plot the index of each node
#' 
#' Plot FlowSOM grid or tree, with in each node a number indicating its index
#' 
#' @param fsom         FlowSOM object, as generated by \code{\link{BuildMST}}
#' @param view     Preferred view, options: "MST", "grid" or "tSNE" (if this
#'                 option was selected while building the MST)
#' @param main         Title of the plot
#' @param nodeSize   Nodesize. The plot might be easier to read if this is a 
#'                   constant number, e.g. 10 or 15
#' @param fontSize  Fontsize, passed to label.cex
#' @param backgroundValues  Values to be used for background coloring, either
#'                          numerical values or something that can be made into
#'                          a factor (e.g. a clustering)
#' @param backgroundColor   Colorpalette to be used for the background coloring
#'                          . Can be either a function or an array specifying
#'                          colors
#' @param backgroundLim     Only used when backgroundValues are numerical. 
#'                          Defaults to min and max of the backgroundValues.
#' @param backgroundBreaks  Breaks to pass on to \code{\link{cut}}, to split
#'                          numerical background values. If NULL, the length of
#'                          backgroundColor will be used (default 100).
#'                          
#' @return Nothing is returned. A plot is drawn in which each node 
#'         is assigned a number
#' @seealso \code{\link{PlotMarker}},\code{\link{PlotStars}},
#'          \code{\link{PlotPies}},\code{\link{PlotCenters}},
#'          \code{\link{BuildMST}}
#' 
#' @examples 
#' # Read from file, build self-organizing map and minimal spanning tree
#' 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))
#' flowSOM.res <- BuildMST(flowSOM.res)
#' 
#' # Plot the node IDs
#' PlotNumbers(flowSOM.res)
#' 
#' # Adapt node size for easier readability
#' PlotNumbers(flowSOM.res, nodeSize=14)
#' 
#' @export
PlotNumbers <- function(fsom, view="MST",main=NULL,nodeSize=fsom$MST$size,
                        fontSize = 1,
                        backgroundValues = NULL,
                        backgroundColor = function(n){
                            grDevices::rainbow(n,alpha=0.3)},
                        backgroundLim = NULL,
                        backgroundBreaks = NULL){
    switch(view,
            MST  = { layout <- fsom$MST$l 
            lty <- 1},
            grid = { layout <- as.matrix(fsom$map$grid)
            lty <- 0},
            tSNE = { layout <- fsom$MST$l2
            lty <- 0},
           stop("The view should be MST, grid or tSNE. tSNE will only work
                   if you specified this when building the MST.")
    )
    
    # Choose background colour
    if(!is.null(backgroundValues)){
        background <- computeBackgroundColor(backgroundValues,backgroundColor,
                                             backgroundLim, backgroundBreaks)
    } else {
        background <- NULL
    }
    
    igraph::plot.igraph(fsom$MST$graph, 
        layout=layout, 
        vertex.size=nodeSize, 
        vertex.label=seq_len(nrow(fsom$map$codes)),
        vertex.label.cex = fontSize,
        edge.lty=lty,
        mark.groups=background$groups,
        mark.col=background$col[background$values],
        mark.border=background$col[background$values],
        main=main)

}

#' Plot a label in each node
#' 
#' Plot FlowSOM grid or tree, with in each node a label. Especially useful to
#' show metacluster numbers
#' 
#' @param fsom         FlowSOM object, as generated by \code{\link{BuildMST}}
#' @param labels    A label for every node
#' @param view     Preferred view, options: "MST", "grid" or "tSNE" (if this
#'                 option was selected while building the MST)
#' @param main         Title of the plot
#' @param nodeSize   Nodesize. The plot might be easier to read if this is a 
#'                   constant number, e.g. 10 or 15
#' @param fontSize  Fontsize, passed to label.cex
#' @param backgroundValues  Values to be used for background coloring, either
#'                          numerical values or something that can be made into
#'                          a factor (e.g. a clustering)
#' @param backgroundColor   Colorpalette to be used for the background coloring
#'                          . Can be either a function or an array specifying
#'                          colors
#' @param backgroundLim     Only used when backgroundValues are numerical. 
#'                          Defaults to min and max of the backgroundValues.
#' @param backgroundBreaks  Breaks to pass on to \code{\link{cut}}, to split
#'                          numerical background values. If NULL, the length of
#'                          backgroundColor will be used (default 100).
#'                          
#' @return Nothing is returned. A plot is drawn in which each node 
#'         is assigned a label
#' @seealso \code{\link{PlotNumbers}}
#' 
#' @examples 
#' # Read from file, build self-organizing map and minimal spanning tree
#' fileName <- system.file("extdata", "68983.fcs", package="FlowSOM")
#' ff <- flowCore::read.FCS(fileName)
#' ff <- flowCore::compensate(ff, ff@description$SPILL)
#' ff <- flowCore::transform(ff, flowCore::estimateLogicle(ff,
#'                                                flowCore::colnames(ff)[8:18]))
#' flowSOM.res <- FlowSOM(ff,
#'                        scale=TRUE,
#'                        colsToUse=c(9,12,14:18),
#'                        nClus = 10,
#'                        seed = 1)
#' 
#' # Plot the node IDs
#' PlotLabels( flowSOM.res$FlowSOM, flowSOM.res$metaclustering, nodeSize=15)
#' 
#' @export
PlotLabels <- function(fsom, 
                       labels,
                       view="MST",
                       main=NULL,
                       nodeSize=fsom$MST$size,
                       fontSize = 1,
                       backgroundValues = NULL,
                       backgroundColor = function(n){
                         grDevices::rainbow(n,alpha=0.3)},
                       backgroundLim = NULL,
                       backgroundBreaks = NULL){
  switch(view,
         MST  = { layout <- fsom$MST$l 
         lty <- 1},
         grid = { layout <- as.matrix(fsom$map$grid)
         lty <- 0},
         tSNE = { layout <- fsom$MST$l2
         lty <- 0}, 
         stop("The view should be MST, grid or tSNE. tSNE will only work
                   if you specified this when building the MST.")
  )
  
  # Choose background colour
  if(!is.null(backgroundValues)){
    background <- computeBackgroundColor(backgroundValues,backgroundColor,
                                         backgroundLim, backgroundBreaks)
  } else {
    background <- NULL
  }
  
  igraph::plot.igraph(fsom$MST$graph, 
                      layout=layout, 
                      vertex.size=nodeSize, 
                      vertex.label=labels,
                      vertex.label.cex = fontSize,
                      edge.lty=lty,
                      mark.groups=background$groups,
                      mark.col=background$col[background$values],
                      mark.border=background$col[background$values],
                      main=main)

}

##############
## PlotPies ##
##############
#' Plot comparison with other clustering
#' 
#' Plot FlowSOM grid or tree, with pies indicating another clustering
#'    or manual gating result
#'
#' @param fsom      FlowSOM object, as generated by \code{\link{BuildMST}}
#' @param cellTypes Array of factors indicating the celltypes
#' @param view     Preferred view, options: "MST", "grid" or "tSNE" (if this
#'                 option was selected while building the MST)
#' @param colorPalette      Colorpalette to be used for the markers
#' @param backgroundValues  Values to be used for background coloring, either
#'                          numerical values or something that can be made into
#'                          a factor (e.g. a clustering)
#' @param backgroundColor   Colorpalette to be used for the background coloring
#'                          . Can be either a function or an array specifying
#'                          colors
#' @param backgroundLim     Only used when backgroundValues are numerical. 
#'                          Defaults to min and max of the backgroundValues.
#' @param backgroundBreaks  Breaks to pass on to \code{\link{cut}}, to split
#'                          numerical background values. If NULL, the length of
#'                          backgroundColor will be used (default 100).
#' @param legend   Logicle, if T add a legend
#' @param main     Title of the plot
#' 
#' @return Nothing is returned. A plot is drawn in which each node is 
#' represented by a pie chart indicating the percentage of cells present 
#' of each cell type. At the end, the layout is set to 1 figure again.
#' 
#' @seealso \code{\link{PlotStars}},\code{\link{PlotMarker}},
#' \code{\link{PlotCenters}},\code{\link{BuildMST}}
#'
#' @examples 
#'#' # Identify the files
#' fcs_file <- system.file("extdata", "68983.fcs", package = "FlowSOM")
#' wsp_file <- system.file("extdata", "gating.wsp", package = "FlowSOM")
#' 
#' # Specify the cell types of interest for assigning one label per cell
#' cell_types <- c("B cells",
#'                 "gd T cells", "CD4 T cells", "CD8 T cells",
#'                 "NK cells","NK T cells")
#'
#' # Parse the FlowJo workspace   
#' library(flowWorkspace)             
#' gatingResult <- GetFlowJoLabels(fcs_file, wsp_file,
#'                                 cell_types = cell_types)
#'
#' # Check the number of cells assigned to each gate
#' colSums(gatingResult$matrix)
#' 
#' # Build a FlowSOM tree
#' flowSOM.res <- FlowSOM(fcs_file, 
#'                        compensate = TRUE, 
#'                        transform = TRUE,
#'                        toTransform = 8:18, 
#'                        colsToUse = c(9,12,14:18),
#'                        nClus = 10,
#'                        seed = 1)
#'    
#'  # Plot pies indicating the percentage of cell types present in the nodes
#'  PlotPies(flowSOM.res$FlowSOM,
#'           gatingResult$manual,
#'           backgroundValues = flowSOM.res$metaclustering)
#'
#' @export
PlotPies <- function(fsom, 
                    cellTypes, 
                    view="MST", #"grid","tSNE"
                    colorPalette=grDevices::colorRampPalette(c("white",
                        "#00007F","blue", "#007FFF", "cyan", "#7FFF7F", 
                        "yellow","#FF7F00", "red")),
                    backgroundValues = NULL,
                    backgroundColor = function(n){
                        grDevices::rainbow(n,alpha=0.3)},
                    backgroundLim = NULL,
                    backgroundBreaks = NULL,
                    # colorRampPalette(c("#FFFFFF","#FF000077"),alpha=TRUE)
                    legend=TRUE,
                    main=""){

    # Choose pie values
    if(!is.factor(cellTypes)){
        cellTypes <- as.factor(cellTypes)
    }
    t <- table(factor(fsom$map$mapping[,1], 
                        levels=seq_along(fsom$MST$size)),
                factor(cellTypes,levels=c(levels(cellTypes),"empty")))
    t[rowSums(t)==0, "empty"] <- 1
    data <- unlist(apply(t, 1, list), recursive = FALSE)
    colors <- list(c(colorPalette(length(levels(cellTypes))),"#000000"))
    colors <- rep(colors,length(fsom$MST$size))
    
    # Choose layout type
    switch(view,
            MST  = { layout <- fsom$MST$l 
            lty <- 1},
            grid = { layout <- as.matrix(fsom$map$grid)
            lty <- 0},
            tSNE = { layout <- fsom$MST$l2
            lty <- 0}, 
           stop("The view should be MST, grid or tSNE. tSNE will only work
                   if you specified this when building the MST.")
    )
    
    # Choose background colour
    if(!is.null(backgroundValues)){
        background <- computeBackgroundColor(backgroundValues,backgroundColor,
                                         backgroundLim, backgroundBreaks)
    } else {
        background <- NULL
    }
    
    # Save plot window settings and minimize margin
    oldpar <- graphics::par(no.readonly = TRUE)
    graphics::par(mar=c(1,1,1,1))
    
    # Add legend
    if(legend){
        if(!is.null(backgroundValues)){
            # Make plot with folowing layout
            # 1 3
            # 2 3
            graphics::layout(matrix(c(1,3,2,3), 2, 2, byrow = TRUE), 
                    widths=c(1,2), heights=c(1))
        } else {
            graphics::layout(matrix(c(1,2), 1, 2, byrow = TRUE), 
                   widths=c(1,2), heights=c(1))
        }
        
        graphics::plot.new()
        legend("center",legend=levels(cellTypes), 
                fill=colors[[1]], 
                cex=0.7, ncol=1, bty="n")
        
        if(!is.null(backgroundValues)){
            PlotBackgroundLegend(backgroundValues,background)
        }
    }
    
    # Plot the actual graph
    igraph::plot.igraph(fsom$MST$g, 
        vertex.shape="pie", 
        vertex.label=NA, 
        vertex.size=fsom$MST$size, 
        vertex.pie=data,
        vertex.pie.color=colors,
        layout=layout, 
        edge.lty=lty, 
        mark.groups=background$groups, 
        mark.col=background$col[background$values], 
        mark.border=background$col[background$values],
        main=main
    )
    
    # Reset plot window
    graphics::par(oldpar)
    graphics::layout(1)
}

############
## mystar ##
############
# Internal use only:
# Add a new vertex shape to iGraph to make star charts
mystar <- function(coords, v=NULL, params) {
    vertex.color <- params("vertex", "color")
    if (length(vertex.color) != 1 && !is.null(v)) {
        vertex.color <- vertex.color[v]
    }
    vertex.size    <- 1/200 * params("vertex", "size")
    if (length(vertex.size) != 1 && !is.null(v)) {
        vertex.size <- vertex.size[v]
    }
    data <- params("vertex", "data")
    cP <- params("vertex","cP")
    scale <- params("vertex","scale")
    bg <- params("vertex","bg")
    graphics::symbols(coords[, 1], coords[, 2], circles = vertex.size, 
                      inches = FALSE, bg = bg, bty='n', add=TRUE) 
    graphics::stars(data, locations = coords, labels = NULL,scale=scale, 
            len = vertex.size, col.segments = cP, 
            draw.segments = TRUE, mar = c(0, 0, 0, 0), add=TRUE, 
            inches=FALSE)
    
}

######################
## legendContinuous ##
######################
# Internal use only, color gradient legend
legendContinuous <- function(col, lev){
    n <- length(col)
    bx <- graphics::par("usr")
    
    box.cx <- c((bx[2]+bx[1])/2 - (bx[2] - bx[1]) / 50,
                (bx[2]+bx[1])/2 + (bx[2] - bx[1]) / 50)
    box.cy <- c(bx[3], bx[3])
    box.sy <- (bx[4] - bx[3]) / n
    
    xx <- rep(box.cx, each = 2)
    
    graphics::par(xpd = TRUE)
    yval <- seq(min(lev), max(lev),by=(max(lev)-min(lev))/n)

    for(i in 1:n){
        
        yy <- c(box.cy[1] + (box.sy * (i - 1)),
                box.cy[1] + (box.sy * (i)),
                box.cy[1] + (box.sy * (i)),
                box.cy[1] + (box.sy * (i - 1)))
        graphics::polygon(xx, yy, col = col[i], border = col[i])
        if((i==1) | (i%%ceiling(n/5)) == 0){
            graphics::text(x = box.cx[2]+ (bx[2] - bx[1]) / 50,
                y = box.cy[1]+(box.sy * (i - 0.5)),
                labels = format(yval[i], digits=3,nsmall=3),
                adj = c(0,0.5))
        }
    }
    yy <- c(box.cy[1],
            box.cy[1] + (box.sy * (n)),
            box.cy[1] + (box.sy * (n)),
            box.cy[1])
    graphics::polygon(xx, yy, col = "#00000000", border = "#000000FF")
}

####################
## shiftFunction ##
####################
shiftFunction <- function(x,n){
    c(x[(n+1):length(x)],x[1:n])
}

####################
## plotStarLegend ##
####################
#' Plot legend for star plot
#' 
#' Plot a single star chart, annotated with labels
#'  
#' @param labels Names to show in the legend
#' @param colors Corresponding colors
#' @param main   Title of the legend
#'
#' @return Nothing is returned. A plot is drawn with 1 star chart, which is 
#' filled completely and annotated with the given labels.
#' @seealso \code{\link{PlotStars}}
plotStarLegend <- function(labels, colors=grDevices::rainbow(length(labels)),
                            main=""){
    graphics::plot(1, type="n", xlab="", ylab="", 
        xlim=c(-10, 10), ylim=c(-3, 3),asp=1,
        bty="n",xaxt="n",yaxt="n",main=main)
    
    graphics::stars(matrix(c(1:(2*length(labels))),nrow=2),col.segments=colors,
        locations=c(0,0),draw.segments = TRUE,add=TRUE,
        inches=FALSE)
    n <- length(labels)
    angle <- 2*pi / n
    angles <- seq(angle/2,2*pi,by=angle)
    
    left <- (angles > (pi/2) & angles < (3*pi/2))
    x <- c(2,-2)[left+1]
    y_tmp <- c(seq(-2,2,by= 4/(sum(!left)+1))[-c(1,sum(!left)+2)],
                seq(2,-2,by=-4/(sum(left)+1))[-c(1,sum(left)+2)])
    y <- shiftFunction(y_tmp,max((cummax(y_tmp)<0)*seq_along(y_tmp)))
    
    for(i in seq_along(labels)){
        graphics::text(x= x[i], 
            y= y[i],
            labels=labels[i],
            adj = c(as.numeric(left)[i],0.5),
            cex = 0.5)
        
        graphics::lines(x=c(x[i]+c(-0.2,0.2)[left[i]+1],
                c(1.5,-1.5)[left[i]+1],
                cos(angles[i])),
            y=c(y[i],
                y[i],
                sin(angles[i])),
            col=colors[i],
            lwd=2)    
    }
}

###################
## plotStarQuery ##
###################
plotStarQuery <- function(labels,values, 
                            colors=grDevices::rainbow(length(labels)),main=""){
    #browser()
    graphics::plot(1, type="n", xlab="", ylab="", 
        xlim=c(-10, 10), ylim=c(-3, 3),asp=1,
        bty="n",xaxt="n",yaxt="n",main=main)
    #graphics::stars(matrix(c(1:(2*length(labels))),nrow=2),col.segments=colors,
    #      locations=c(0,0),draw.segments = TRUE,add=TRUE,
    #      inches=FALSE)
    
    graphics::stars(matrix(c(rep(0,length(labels)),values),nrow=2,byrow=TRUE),
            col.segments=colors,
            locations=c(0,0),draw.segments = TRUE,add=TRUE,
            inches=FALSE)
    n <- length(labels)
    angle <- 2*pi / n
    angles <- seq(angle/2,2*pi,by=angle)
    
    left <- (angles > (pi/2) & angles < (3*pi/2))
    x <- c(2,-2)[left+1]
    y_tmp <- c(seq(-2,2,by= 4/(sum(!left)+1))[-c(1,sum(!left)+2)],
                seq(2,-2,by=-4/(sum(left)+1))[-c(1,sum(left)+2)])
    y <- shiftFunction(y_tmp,max((cummax(y_tmp)<0)*seq_along(y_tmp)))
    
    for(i in seq_along(labels)){
        graphics::text(x= x[i], 
            y= y[i],
            labels=labels[i],
            adj = c(as.numeric(left)[i],0.5),
            cex = 0.5)
        
        graphics::lines(x=c(x[i]+c(-0.2,0.2)[left[i]+1],
                c(1.5,-1.5)[left[i]+1],
                cos(angles[i])),
            y=c(y[i],
                y[i],
                sin(angles[i])),
            col=colors[i],
            lwd=2)    
    }
}

###############
## PlotStars ##
###############
#' Plot star charts
#' 
#' Plot FlowSOM grid or tree, where each node is represented by 
#' a star chart indicating median marker values
#' 
#' @param fsom     FlowSOM object, as generated by \code{\link{BuildMST}}
#' @param markers  Array of markers to use. Default: the markers used to build 
#'                 the tree
#' @param view     Preferred view, options: "MST", "grid" or "tSNE" (if this
#'                 option was selected while building the MST)
#' @param colorPalette      Colorpalette to be used for the markers
#' @param starBg Background color inside the star circle. Default is "white".
#'               Can also be put to  "transparent" (as was the case for older 
#'               versions).
#' @param backgroundValues  Values to be used for background coloring, either
#'                          numerical values or something that can be made into
#'                          a factor (e.g. a clustering)
#' @param backgroundColor   Colorpalette to be used for the background coloring
#'                          . Can be either a function or an array specifying
#'                          colors
#' @param backgroundLim     Only used when backgroundValues are numerical. 
#'                          Defaults to min and max of the backgroundValues.
#' @param backgroundBreaks  Breaks to pass on to \code{\link{cut}}, to split
#'                          numerical background values. If NULL, the length of
#'                          backgroundColor will be used (default 100).
#' @param backgroundSize Size of the background circles. Default 15.
#' @param thresholds    Optional. Array containing a number for each of the 
#'                      markers to be used as the split between high/low. 
#'                      If provided, the percentage of positive cells is
#'                      indicated instead of the MFI
#' @param legend   Logical, if TRUE add a legend
#' @param query    Show a low/high profile for certain markers in the legend.
#'                 See also \code{\link{QueryStarPlot}}
#' @param range    If "all" (default), range is computed on all markers passed,
#'                 if "one", range is computed on every marker separately. The
#'                 height of the pie pieces will be computed relative to this 
#'                 range.
#' @param main     Title of the plot
#' 
#' @return Nothing is returned. A plot is drawn in which each node is 
#' represented by a star chart indicating the median fluorescence intensities.
#' Resets the layout back to 1 plot at the end.
#' 
#' @seealso \code{\link{PlotPies}},\code{\link{PlotMarker}},
#' \code{\link{PlotCenters}}, \code{\link{BuildMST}}
#' 
#' @examples
#'    # Read from file, build self-organizing map and minimal spanning tree
#'    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))
#'    flowSOM.res <- BuildMST(flowSOM.res)
#'    
#'    # Plot stars indicating the MFI of the cells present in the nodes
#'    PlotStars(flowSOM.res)
#'
#' @export
PlotStars <- function(fsom, 
                      markers=fsom$map$colsUsed, 
                      view="MST", #"grid","tSNE"
                      colorPalette=grDevices::colorRampPalette(
                        c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", 
                          "yellow", "#FF7F00", "red", "#7F0000")),
                      starBg = "white",
                      backgroundValues = NULL,
                      backgroundColor = function(n){
                        grDevices::rainbow(n, alpha=0.3)},
                      backgroundLim = NULL,
                      backgroundBreaks = NULL,
                      backgroundSize = NULL,
                      thresholds=NULL,
                      legend=TRUE,
                      query=NULL,
                      range = "all",
                      main=""){
    # Add star chart option to iGraph
    add.vertex.shape("star", clip=igraph.shape.noclip, plot=mystar, 
                    parameters=list(vertex.data=NULL,vertex.cP = colorPalette,
                                    vertex.scale=FALSE, vertex.bg = starBg))
    
    if (is.null(thresholds)) {
        # Use MFIs
        data <- fsom$map$medianValues[, markers, drop=FALSE]
        if (range == "all") {
          min_data <- min(data, na.rm = TRUE)
          max_data <- max(data, na.rm = TRUE)
          data <- (data - min_data) / (max_data - min_data)
        } else if (range == "one"){
          data <- apply(data, 2, function(x){
            min_x <- min(x, na.rm = TRUE)
            max_x <- max(x, na.rm = TRUE)
            (x - min_x)/ (max_x - min_x) 
          })
        }
    } else {
        # scale thresholds same as data
        if(fsom$transform){
            warning("Thresholds should be given in the transformed space")
        }
        if(!is.null(fsom$scaled.center)){
          thresholds <- scale(t(thresholds), 
                              center = fsom$scaled.center[markers],
                              scale = fsom$scaled.scale[markers])
        }
        # Use pctgs of cells above threshold as star plot values
        data <-
            t(sapply(seq_len(fsom$map$nNodes), function(i) {
                res = NULL
                for(m in seq_along(markers)){
                    res = c(res, 
                            sum(subset(fsom$data, 
                               fsom$map$mapping[,1] == i)[,
                                                  markers[m]] > thresholds[m])/
                            sum(fsom$map$mapping[,1] == i))
                    
                }
                res
            }))
    }
    
    # Choose layout type
    switch(view,
        MST  = { layout <- fsom$MST$l 
                    lty <- 1},
        grid = { layout <- as.matrix(fsom$map$grid)
                    lty <- 0},
        tSNE = { layout <- fsom$MST$l2
                    lty <- 0}, 
        stop("The view should be MST, grid or tSNE. tSNE will only work
                   if you specified this when building the MST.")
    )
    
    # Choose background colour
    if (!is.null(backgroundValues)) {
        background <- computeBackgroundColor(backgroundValues,backgroundColor,
                                             backgroundLim, backgroundBreaks)
        if (is.null(backgroundSize)) { 
          backgroundSize <- fsom$MST$size
          backgroundSize[backgroundSize == 0] <- 3
        }
    } else {
        background <- NULL
    }
    
    # Save plot window settings and minimize margin
    oldpar <- graphics::par(no.readonly = TRUE)
    graphics::par(mar=c(1,1,1,1))
    
    # Add legend
    if(legend){
        if(!is.null(backgroundValues)){
            # Make plot with folowing layout
            # 1 3
            # 2 3
            graphics::layout(matrix(c(1,3,2,3), 2, 2, byrow = TRUE), 
                    widths=c(1,2), heights=c(1))
        } else {
            graphics::layout(matrix(c(1,2), 1, 2, byrow = TRUE), 
                   widths=c(1,2), heights=c(1))
        }
         
       if(is.null(query)){
            plotStarLegend(fsom$prettyColnames[markers], 
                            colorPalette(ncol(data)))
        } else {
            plotStarQuery(fsom$prettyColnames[markers],
                            values=query == "high",
                            colorPalette(ncol(data)))
        }
        
        if(!is.null(backgroundValues)){
            PlotBackgroundLegend(backgroundValues,background)
        }
    }
    
    # Plot the actual graph
    igraph::plot.igraph(fsom$MST$g, 
                        vertex.shape = "star", 
                        vertex.label = NA, 
                        vertex.size = fsom$MST$size, 
                        vertex.data = data,
                        vertex.cP = colorPalette(ncol(data)),
                        vertex.scale = FALSE,
                        layout = layout, 
                        edge.lty = lty,  
                        mark.groups = background$groups, 
                        mark.col = background$col[background$values], 
                        mark.border = background$col[background$values],
                        mark.expand	= backgroundSize,
                        main=main
    )
    # Reset plot window
    graphics::par(oldpar)
    graphics::layout(1)
}




##############
## PlotNode ##
##############
#' Plot star chart
#' 
#' Plot a star chart indicating median marker values of a single node
#' 
#' @param fsom     FlowSOM object, as generated by \code{\link{BuildMST}} or
#'                 the first element of the list returned by 
#'                 \code{\link{FlowSOM}}
#' @param id       Id of the node to plot (check PlotNumbers to get the ids)
#' @param markers  Array of markers to use. Default: the markers used to build 
#'                 the tree
#' @param colorPalette      Colorpalette to be used for the markers
#' @param main     Title of the plot
#' 
#' @return Nothing is returned. A plot is drawn in which the node is 
#' represented by a star chart indicating the median fluorescence intensities.
#' 
#' @seealso \code{\link{PlotStars}},\code{\link{PlotNumbers}},
#'  \code{\link{FlowSOM}}
#' 
#' @examples
#'    # Read from file, build self-organizing map and minimal spanning tree
#'    fileName <- system.file("extdata", "68983.fcs", package="FlowSOM")
#'    flowSOM.res <- FlowSOM(fileName, compensate=TRUE,transform=TRUE,
#'                             scale=TRUE,colsToUse=c(9,12,14:18),nClus=10)
#'    
#'    # Plot stars indicating the MFI of the cells present in the nodes
#'    PlotNode(flowSOM.res$FlowSOM,1)
#'
#' @export
PlotNode <- function(fsom, id, 
                     markers=fsom$map$colsUsed, 
                     colorPalette=grDevices::colorRampPalette(
                         c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", 
                           "yellow", "#FF7F00", "red", "#7F0000")),
                     main=paste0("Cluster ",id)){
    
    graphics::plot(1, type="n", xlab="", ylab="", 
                   xlim=c(-10, 10), ylim=c(-3, 3),asp=1,
                   bty="n",xaxt="n",yaxt="n",main=main)
    
    labels <- fsom$prettyColnames[markers]
    n <- length(markers)
    
    if(is.function(colorPalette)){colorPalette <- colorPalette(n)}
    data <- rbind(apply(fsom$map$medianValues[, markers,drop=FALSE],2,min, 
                        na.rm = TRUE),
                  fsom$map$medianValues[id, markers,drop=FALSE],
                  apply(fsom$map$medianValues[, markers,drop=FALSE],2,max, 
                        na.rm = TRUE))
    coords <- matrix(c(100,0,100,100,0,100),nrow=3)
    
    print(data)
    graphics::stars(data,col.segments=colorPalette,
                    locations=coords,
                    draw.segments = TRUE,add=TRUE,
                    inches=FALSE)
    
    graphics::symbols(coords[, 1], coords[, 2], circles = c(1,1,1), 
                      inches = FALSE, bg = "transparent", bty='n', add=TRUE) 
    
    angle <- 2*pi / n
    angles <- seq(angle/2,2*pi,by=angle)
    
    left <- (angles > (pi/2) & angles < (3*pi/2))
    x <- c(2,-2)[left+1]
    y_tmp <- c(seq(-2,2,by= 4/(sum(!left)+1))[-c(1,sum(!left)+2)],
               seq(2,-2,by=-4/(sum(left)+1))[-c(1,sum(left)+2)])
    y <- shiftFunction(y_tmp,max((cummax(y_tmp)<0)*seq_along(y_tmp)))
    
    for(i in seq_along(labels)){
        graphics::text(x= x[i], 
                       y= y[i],
                       labels=labels[i],
                       adj = c(as.numeric(left)[i],0.5),
                       cex = 0.5)
        
        graphics::lines(x=c(x[i]+c(-0.2,0.2)[left[i]+1],
                            c(1.5,-1.5)[left[i]+1],
                            cos(angles[i])),
                        y=c(y[i],
                            y[i],
                            sin(angles[i])),
                        col=colorPalette[i],
                        lwd=2)    
    }
}
# #################
# ## PlotStarsSD ##
# #################
# #' Plot standard deviation star charts
# #' 
# #' Plot FlowSOM grid or tree, where each node is represented by 
# #' a star chart indicating standard deviation of the marker values
# #' 
# #' @param fsom     FlowSOM object, as generated by \code{\link{BuildMST}}
# #' @param markers  Array of markers to use. Default: the markers used to build 
# #'                 the tree
# #' @param view     Preferred view, options: "MST", "grid" or "tSNE" (if this
# #'                 option was selected while building the MST)
# #' @param colorPalette      Colorpalette to be used for the markers
# #' @param backgroundValues  Values to be used for background coloring, either
# #'                          numerical values or something that can be made 
# #'                          into a factor (e.g. a clustering)
# #' @param backgroundColor   Colorpalette to be used for the background 
# #'                          coloring. Can be either a function or an array 
# #'                          specifying colors
# #' @param backgroundLim     Only used when backgroundValues are numerical. 
# #'                          Defaults to min and max of the backgroundValues.
# #' @param legend   Logical, if TRUE add a legend
# #' @param query    Show a low/high profile for certain markers in the legend.
# #'                 See also \code{\link{QueryStarPlot}}
# #' @param main     Title of the plot
# #' 
# #' @return Nothing is returned. A plot is drawn in which each node is 
# #' represented by a star chart indicating the standard deviation of the
# #' fluorescence intensities. Resets the layout back to 1 plot at the end.
# #' 
# #' @seealso \code{\link{PlotStars}}
# #' 
# #' @examples
# #'    # Read from file, build self-organizing map and minimal spanning tree
# #'    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))
# #'    flowSOM.res <- BuildMST(flowSOM.res)
# #'    
# #'    # Plot stars indicating the MFI of the cells present in the nodes
# #'    PlotStars(flowSOM.res)
# #'    # Plot stars indicating the standard deviations of the MFIs
# #'    PlotStarsSD(flowSOM.res)
# #'
# #' @export
# PlotStarsSD <- function(fsom, 
#                         markers=fsom$map$colsUsed, 
#                         view="MST", #"grid","tSNE"
#                         colorPalette=grDevices::colorRampPalette(
#                             c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", 
#                             "yellow", "#FF7F00", "red", "#7F0000")),
#                         backgroundValues = NULL,
#                         backgroundColor = function(n){grDevices::rainbow(n,
#                                                                 alpha=0.3)},
#                         backgroundLim = NULL,
#                         #colorRampPalette(c("#FFFFFF","#FF000077"),alpha=TRUE)
#                         legend=TRUE,
#                         query=NULL,
#                         main=""){
#     # Add star chart option to iGraph
#     add.vertex.shape("star", clip=igraph.shape.noclip, plot=mystar, 
#                     parameters=list(vertex.data=NULL,
#                                     vertex.cP = colorPalette,
#                                     vertex.scale = TRUE))
#     
#     # Choose star values
#     data <- fsom$map$sdValues[, markers,drop=FALSE]
#     for(i in seq_along(markers)){
#         data[,i] <- data[,i]/diff(range(
#             fsom$map$medianValues[,markers[i]] , na.rm = TRUE))
#     }
#     # Choose layout type
#     switch(view,
#             MST  = { layout <- fsom$MST$l 
#                     lty <- 1},
#             grid = { layout <- as.matrix(fsom$map$grid)
#                     lty <- 0},
#             tSNE = { layout <- fsom$MST$l2
#                     lty <- 0}
#     )
#     
#     # Choose background colour
#     backgroundList <- list()
#     backgroundColors <- NULL
#     if(!is.null(backgroundValues)){
#         if(is.numeric(backgroundValues)){
#             backgroundList <- as.list(seq_along(backgroundValues))
#             
#             if(class(backgroundColor)=="function"){
#                 backgroundColor <- backgroundColor(100)
#             }
#             
#             if(length(backgroundLim) > 0){
#                 ids <- cut(c(backgroundLim,backgroundValues),
#                                 length(backgroundColor)
#                             )[-c(seq_along(backgroundLim))]
#             } else {
#                 ids <- cut(backgroundValues,
#                             length(backgroundColor))
#                 
#             }
#             backgroundColors <- backgroundColor[ids]    
#         } else {
#             if(! is.factor(backgroundValues)){
#                 backgroundValues <- as.factor(backgroundValues)
#             }
#             
#             backgroundList <- as.list(seq_along(backgroundValues))
#             
#             if(class(backgroundColor)=="function"){
#                 backgroundColor <- backgroundColor(
#                     length(levels(backgroundValues)))
#             }
#             
#             if(length(backgroundColor) < length(levels(backgroundValues))){
#                 stop("You specified less backgroundcolours than groups.")
#             }
#             
#             backgroundColors <- backgroundColor[backgroundValues]    
#             
#         }    
#     }
#     
#     # Save plot window settings and minimize margin
#     oldpar <- par(no.readonly = TRUE)
#     par(mar=c(1,1,1,1))
#     
#     # Add legend
#     if(legend){
#         if(!is.null(backgroundValues)){
#             # Make plot with folowing layout
#             # 1 3
#             # 2 3
#             graphics::layout(matrix(c(1,3,2,3), 2, 2, byrow = TRUE), 
#                     widths=c(1,2), heights=c(1))
#             if(is.null(query)){
#                 plotStarLegend(fsom$prettyColnames[markers], 
#                                 colorPalette(ncol(data)))
#             } else {
#                 plotStarQuery(fsom$prettyColnames[markers],
#                                 values=query == "high",
#                                 colorPalette(ncol(data)))
#             }
#             graphics::plot.new()
#             if(is.factor(backgroundValues)){
#                 legend("center", legend=levels(backgroundValues),
#                         fill=backgroundColor, 
#                         cex=0.7, ncol=1, bty="n",title="Background")       
#             } else {
#                 legendContinuous(backgroundColor,c(backgroundLim,
#                                             backgroundValues))
#             }
#         } else {
#             graphics::layout(matrix(c(1,2), 1, 2, byrow = TRUE), 
#                     widths=c(1,2), heights=c(1))
#             #graphics::layout(matrix(c(1,2), 2, 1, byrow = TRUE), 
#             #       widths=c(1), heights=c(1,2))
#             plotStarLegend(fsom$prettyColnames[markers], 
#                             colorPalette(ncol(data)))
#         }
#     }
#     
#     # Plot the actual graph
#     igraph::plot.igraph(fsom$MST$g, 
#         vertex.shape="star", 
#         vertex.label=NA, 
#         vertex.size=fsom$MST$size, 
#         vertex.data=data,
#         vertex.cP=colorPalette(ncol(data)),
#         vertex.scale=FALSE,
#         layout=layout, 
#         edge.lty=lty, 
#         mark.groups=backgroundList, 
#         mark.col=backgroundColors, 
#         mark.border=backgroundColors,
#         main=main
#     )
#     
#     # Reset plot window
#     par(oldpar)
#     graphics::slayout(1)
# }

####################
## PlotClusters2D ##
####################
#' Plot nodes on scatter plot
#' 
#' Plot a 2D scatter plot. All cells of fsom$data are plotted in
#' black, and those of the selected nodes are plotted in red.
#' The nodes in the grid are indexed starting from the left bottom,
#' first going right, then up. E.g. In a 10x10 grid, the node at
#' top left will have index 91.
#'  
#' @param fsom    FlowSOM object, as generated by \code{\link{BuildMST}}
#' @param marker1 Marker to plot on the x-axis
#' @param marker2 Marker to plot on the y-axis
#' @param nodes   Nodes of which the cells should be plotted in red
#' @param main    Title of the plot
#' @param col     Colors for all the cells in the selected nodes (ordered array)
#' @param maxBgPoints Maximum number of background points to plot
#' @param pchBackground Character to use for background cells
#' @param pchCluster Character to use for cells in cluster
#' @param xlab    Label for the x axis
#' @param ylab    Label for the y axis
#' @param xlim    Limits for the x axis
#' @param ylim    Limits for the y axis
#' @param ...     Other parameters to pass on to plot
#'
#' @return Nothing is returned. A plot is drawn in which all cells are plotted
#'         in black and the cells of the selected nodes in red.
#' @seealso \code{\link{PlotNumbers}},\code{\link{PlotCenters}},
#'          \code{\link{BuildMST}}
#' 
#' @examples
#'    # Read from file, build self-organizing map and minimal spanning tree
#'    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))
#'   flowSOM.res <- BuildMST(flowSOM.res)
#'    
#'    # Plot cells
#'    PlotClusters2D(flowSOM.res,1,2,91)
#'
#' @export
PlotClusters2D <- function(fsom, marker1, marker2, nodes, 
                            col = "#FF0000",maxBgPoints=10000,
                            pchBackground=".",
                            pchCluster=".",
                            main="",
                            xlab=fsom$prettyColnames[marker1],
                            ylab=fsom$prettyColnames[marker2],
                            xlim=c(min(fsom$data[,marker1]),
                                   max(fsom$data[,marker1])),
                            ylim=c(min(fsom$data[,marker2]),
                                   max(fsom$data[,marker2])),
                            ...){
    if(!is.null(maxBgPoints)){
        background <- sample(seq_len(nrow(fsom$data)),
                            min(maxBgPoints,nrow(fsom$data)))
    } else{
        background <- seq_len(nrow(fsom$data))
    }
    graphics::plot(fsom$data[background, c(marker1, marker2)], 
        pch=pchBackground, col="#000000AA",
        main=main,xlab=xlab,ylab=ylab,
        xlim=xlim,ylim=ylim,
        ...)
    graphics::points(fsom$data[fsom$map$mapping[,1] %in% nodes,
            c(marker1, marker2)], 
            pch=pchCluster, 
            col=col)
    graphics::points(fsom$map$medianValues[nodes,marker1],
            fsom$map$medianValues[nodes,marker2],
            pch="x",col="blue")
    #cat(nodes,": \n",table(m[fsom$map$mapping[,1] %in% nodes]),"\n")
}

####################
## PlotOverview2D ##
####################
#' Plot metaclusters on scatter plots
#' 
#' Write multiple 2D scatter plots to a png file. 
#' All cells of fsom$data are plotted in black, and those of the selected 
#' metaclusters are plotted in color.
#'  
#' @param fsom         FlowSOM object, as generated by \code{\link{FlowSOM}}.
#'                     If using a FlowSOM object as generated by 
#'                     \code{\link{BuildMST}}, it needs to be wrapped in a list,
#'                     list(FlowSOM = fsom, metaclustering = metaclustering).
#' @param markerlist   List in which each element is a pair of marker names
#' @param metaclusters Metaclusters of interest
#' @param colors       Named vector with color value for each metacluster. 
#'                     If NULL (default) colorbrewer "paired" is interpolated
#' @param ff           flowFrame to use as reference for the marker names
#' @param ...     Other parameters to pass on to PlotClusters2D
#'
#' @return Nothing is returned, but a plot is drawn for every markerpair and
#'         every metacluster. The individual cells are colored, and the 
#'         center of each FlowSOM cluster is indicated with a blue cross.
#' @seealso \code{\link{PlotClusters2D}}
#' 
#' @examples
#'    # Read from file, build self-organizing map and minimal spanning tree
#'    fileName <- system.file("extdata", "68983.fcs", package="FlowSOM")
#'    flowSOM.res <- FlowSOM(fileName, 
#'                           compensate=TRUE, transform=TRUE, scale=TRUE,
#'                           colsToUse=c(9,12,14:18),
#'                           nClus = 10,
#'                           seed = 1)
#'                           
#'    # Plot cells
#'    markers_of_interest = list(c("FSC-A", "SSC-A"),
#'                               c("CD3", "CD19"),
#'                               c("TCRb", "TCRyd"),
#'                               c("CD4", "CD8"))
#'    metaclusters_of_interest = 1:10
#'    
#'    # Recommended to write to png
#'    
#'    png("Markeroverview.png",
#'        width = 500 * length(markers_of_interest),
#'        height = 500 * length(metaclusters_of_interest))
#'    PlotOverview2D(flowSOM.res,
#'                   markerlist = markers_of_interest,
#'                   metaclusters = metaclusters_of_interest,
#'                   pchCluster = 19,
#'                   ff = flowCore::read.FCS(fileName))
#'    dev.off()
#'
#' @export
PlotOverview2D <- function(fsom, 
                           markerlist, 
                           metaclusters,
                           colors = NULL,
                           ff,
                           ...){
  graphics::layout(matrix(seq_len(length(markerlist) * length(metaclusters)), 
                nrow = length(metaclusters)))
  if(is.null(colors)){
    colors <- RColorBrewer::brewer.pal(12,
                                       "Paired")
    colors <- rep(colors, length.out = length(metaclusters))
    names(colors) <- as.character(metaclusters)
  }
  
  for(marker_pair in markerlist){
    for(metacluster in metaclusters){
      PlotClusters2D(fsom$FlowSOM,
                     marker1 = get_channels(ff, marker_pair[1]),
                     marker2 = get_channels(ff, marker_pair[2]),
                     nodes = which(fsom$metaclustering == metacluster),
                     main = paste0("Metacluster ", metacluster),
                     col = colors[as.character(metacluster)],
                     ...)
    }
  }
  
  graphics::layout(1)
}

#################
## PlotCenters ##
#################
#' Plot cluster centers on a 2D plot
#' 
#' Plot FlowSOM nodes on a 2D scatter plot of the data
#' 
#' @param fsom    FlowSOM object, as generated by \code{\link{BuildMST}}
#' @param marker1 Marker to show on the x-axis
#' @param marker2 Marker to show on the y-axis
#' @param MST     Type of visualization, if 1 plot tree, else plot grid
#' @return Nothing is returned. A 2D scatter plot is drawn on which the nodes 
#'         of the grid are indicated
#'         
#' @seealso \code{\link{PlotStars}},\code{\link{PlotPies}},
#'          \code{\link{PlotMarker}},\code{\link{BuildMST}}
#'          
#' @examples
#'    # Read from file, build self-organizing map and minimal spanning tree
#'    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))
#'    flowSOM.res <- BuildMST(flowSOM.res)
#'    
#'    # Plot centers
#'    PlotCenters(flowSOM.res,"FSC-A","SSC-A")
#'    PlotCenters(flowSOM.res,2,5)
#'
#' @export
PlotCenters <- function(fsom, marker1, marker2, MST=TRUE){
    
    graphics::plot(fsom$data[, c(marker1, marker2)], pch=".", col="#000000AA")
    d <- fsom$map$medianValues[, c(marker1, marker2)]
    graphics::points(d,col="red")
    
    if(MST==1){
        g <- fsom$MST$graph
        e <- get.edges(g, E(g))
    } else {
        e <- which(as.matrix(
            stats::dist(fsom$map$grid,method = "manhattan"))==1,
            arr.ind=TRUE)
    }
    
    for(i in seq_len(nrow(e))){
        graphics::lines(d[c(e[i, 1], e[i, 2]), 1], 
                        d[c(e[i, 1], e[i, 2]), 2], 
                        col="red")
    }
}

###################
## 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,]
#'    ff1@@description$FIL <- "File1"
#'    ff2 <- flowCore::read.FCS(fileName)[1001:2000,]
#'    ff2@@description$FIL <- "File2"
#'    
#'    flowSOM.res <- FlowSOM(flowCore::flowSet(c(ff1,ff2)), compensate=TRUE,
#'                           transform=TRUE, scale=TRUE,
#'                           colsToUse=c(9,12,14:18), maxMeta=10)
#'    fSOM <- flowSOM.res[[1]]
#'    
#'    # see $metadata for subsets:
#'    fSOM$metaData
#'    
#'    # Use only the second file, without changing the map
#'    fSOM2 <- FlowSOMSubset(fSOM,
#'                           (fSOM$metaData[[2]][1]):(fSOM$metaData[[2]][2]))
#'
#' @export
FlowSOMSubset <- function(fsom,ids){
    fsom_tmp <- fsom
    fsom_tmp$data <- fsom$data[ids,]
    fsom_tmp$map$mapping <- fsom$map$mapping[ids,]
    fsom_tmp <- UpdateDerivedValues(fsom_tmp)
    fsom_tmp <- UpdateNodeSize(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 mad_allowed   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 + mad_allowed * 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 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.
#'        
#' @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,ff@@description$SPILL)
#'    ff <- flowCore::transform(ff,
#'              flowCore::transformList(colnames(ff@@description$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,
                    mad_allowed = 4,
                    compensate = NULL, 
                    spillover = NULL,
                    transform = NULL, 
                    toTransform = NULL, 
                    transformFunction = NULL,
                    scale = NULL, 
                    scaled.center = NULL, 
                    scaled.scale = NULL){
  
  if (class(fsom) == "list" & !is.null(fsom$FlowSOM)) {
    fsom_o <- fsom
    fsom <- fsom$FlowSOM 
  }
  if (class(fsom) != "FlowSOM") {
    stop("fsom should be a FlowSOM object.")
  }
  
  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(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, 
                          scale = scale, scaled.center = scaled.center,
                          scaled.scale = scaled.scale)
    
    fsom_new$map <- fsom$map
    fsom_new$MST <- fsom$MST
    
    fsom_new$map$mapping <- MapDataToCodes(fsom$map$codes,fsom_new$data)
    fsom_new <- UpdateDerivedValues(fsom_new)
    fsom_new <- UpdateNodeSize(fsom_new)
    
    test_outliers <- TestOutliers(fsom_new,
                                  mad_allowed = mad_allowed,
                                  fsom_reference = 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.")
    }
   
    if(exists("fsom_o")){
      return(list(FlowSOM = fsom_new,
                  metaclustering = fsom_o$metaclustering))
    } else {
      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{mad_allowed} times the median absolute
#' deviation of the distances. 
#'
#' @param fsom  FlowSOM object
#' @param mad_allowed Number of median absolute deviations allowed. Default = 4.
#' @param fsom_reference FlowSOM object to use as reference. If NULL (default),
#'                       the original fsom object is used.
#' @param plot Should a plot be generated showing the distribution of the
#'             distances. Default is FALSE.
#' @param img_file If plot is TRUE, the output will be written to this file.
#'                 Default is "testOutliers.pdf"
#'        
#' @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)
#'
#' @export
TestOutliers <- function(fsom, 
                         mad_allowed = 4,
                         fsom_reference = NULL,
                         plot = FALSE,
                         img_file = "testOutliers.pdf"){
  
  if (class(fsom) == "list" & !is.null(fsom$FlowSOM)) {
    fsom <- fsom$FlowSOM 
  }
  if (class(fsom) != "FlowSOM") {
    stop("fsom should be a FlowSOM object.")
  }
  
  if(is.null(fsom_reference)){
    fsom_reference <- fsom
  } else {
    if (class(fsom_reference) == "list" & !is.null(fsom_reference$FlowSOM)) {
      fsom_reference <- fsom_reference$FlowSOM 
    }
    if (class(fsom_reference) != "FlowSOM") {
      stop("fsom should be a FlowSOM object.")
    }
  }
  
  distances_median <- sapply(seq_len(fsom_reference$map$nNodes),
                             function(x){
                               ids <- which(GetClusters(fsom_reference) == x)
                               if(length(ids) > 0){
                                 m <- stats::median(
                                   fsom_reference$map$mapping[ids, 2])
                               } else {
                                 m <- 0
                               }
                               return(m)
                             })
  
  distances_mad <- sapply(seq_len(fsom_reference$map$nNodes),
                             function(x){
                               ids <- which(GetClusters(fsom_reference) == x)
                               if(length(ids) > 0){
                                 m <- stats::mad(
                                   fsom_reference$map$mapping[ids, 2])
                               } else {
                                 m <- 0
                               }
                               return(m)
                             })
  
  thresholds <- distances_median + mad_allowed * 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(plot){
    grDevices::pdf(img_file, width = 20, height = 20)
    xdim <- fsom$map$xdim
    ydim <- fsom$map$ydim
    graphics::layout(matrix(1:(xdim*ydim), nrow = xdim))
    for(i in 1:(xdim*ydim)){
      ids <- which(GetClusters(fsom) == i)
      values <- fsom$map$mapping[ids, 2]
        if(length(values) > 1){
          nOutliers <- sum(values > thresholds[i])
          graphics::hist(values, main = paste0(i," (", nOutliers,")"), 
                         xlab = "")
          graphics::abline(v = distances_median[i], col = "black", lwd = 2)
          graphics::abline(v = thresholds[i], col = "red", lwd = 2)
        } else {
          graphics::plot.new()
        }
    }
    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)
}

################
## PlotGroups ##
################
#' Plot differences between groups
#' 
#' Plot FlowSOM trees, where each node is represented by 
#' a star chart indicating mean marker values, the size of the node is relative
#' to the mean percentage of cells present in each
#'
#' @param fsom    FlowSOM object, as generated by \code{\link{BuildMST}} or
#'        the first list item of \code{\link{FlowSOM}}
#' @param groups  groups result as generated by \code{\link{CountGroups}}
#' @param tresh   Relative difference in groups before the node is coloured
#' @param p_tresh Threshold on p-value from wilcox-test before the node is
#'        coloured. If this is not NULL, tresh will be ignored.
#' @param heatmap   If TRUE, the scores are plotted in a gradient instead
#'                  of only the selection that passes the threshold
#' @param ... Other parameters to pass to \code{\link{PlotStars}}
#' @return A vector containing the labels assigned to the nodes for
#'         all groups except the first
#' @seealso \code{\link{PlotStars}},\code{\link{CountGroups}}
#' @examples
#'    ## Use the wrapper function to build a flowSOM object (saved in fsom[[1]])
#'    ## and a metaclustering (saved in fsom[[2]])
#'    # fsom <- FlowSOM(ff,compensate = FALSE, transform = FALSE,scale = TRUE,
#'    #              colsToUse = colsToUse, nClus = 10, silent = FALSE,
#'    #             xdim=7, ydim=7)
#'    
#'    ## Make a list with for each group a list of files
#'    ## The reference group should be the first
#'    # groups <- list("C"=file.path(workDir,grep("C",files,value = TRUE)),
#'    #             "D"=file.path(workDir,grep("D",files,value=TRUE)))
#'    
#'    ## Compute the percentages for all groups
#'    # groups_res <- CountGroups(fsom[[1]],groups)
#'   
#'    ## Plot the groups. For all groups except the first, differences with the
#'    ## first group are indicated
#'    # annotation <- PlotGroups(fsom[[1]],groups_res)
#'    
#' @export
PlotGroups <- function(fsom, groups,
                       tresh = NULL, p_tresh = 0.05, 
                       heatmap = FALSE, ...){
    
  groupnames <- rownames(groups$means)
  annotation <- list()
  
  # Compare all groups with the first one
  for(group in groupnames[-1]){
    if (!is.null(p_tresh) & is.null(tresh)) {
      # Comparing the percentages of the individual samples,
      # using a wilcox test
      values <- groups$pctgs
      p_v <- rep(NA, ncol(values))
      for(i in seq_len(ncol(values))){
        test <- stats::wilcox.test(
          values[groups$groups %in% groupnames[1], i],
          values[groups$groups %in% group, i],
          exact = FALSE)
        p_v[i] <- test$p.value
      }
      diff <- groups$means[group,] - groups$means[groupnames[1],]
      adj_p <- stats::p.adjust(p_v, "BH")
      
      if(heatmap){ 
        score <- log10(adj_p) * (-1)^(diff < 0)
        annotation[[group]] <- score
      } else{
        score <- 1 + (adj_p < p_tresh) + (diff > 0 & adj_p < p_tresh)
        annotation[[group]] <- factor(c("--",
                                        groupnames[1],
                                        group)[score],
                                      levels=c("--",groupnames[1],group))
      }
      
    } else if (is.null(p_tresh) & !is.null(tresh)) {
      # Using fold change of the mean per group
      diff <- groups$means[group,] - groups$means[1,]
      values <- (apply(groups$means[c(groupnames[1],group),],2,max) / 
                   apply(groups$means[c(groupnames[1],group),],2,min))
      
      if (heatmap) {
        annotation[[group]] <- values * (-1)^(diff < 0)
      } else {
        annotation[[group]] <- as.factor(c("--",groupnames[1],group)[
          1 + (values > tresh) + (diff > 0 & values > tresh)])
      }
    } else {
      stop("Please use only tresh or p_tresh, not both.")
    }
  }
  
  # Plot the first group
  fsom <- UpdateNodeSize(fsom, reset=TRUE)
  fsom$MST$size <- fsom$MST$size * groups$means_norm[[groupnames[1]]]
  PlotStars(fsom,
            main = groupnames[1],
            ...)
  
  if(!heatmap){
    backgroundColors <- grDevices::colorRampPalette(
      c("#FFFFFF22","#00FFFF55","#FF000055"),alpha=TRUE)
  } else {
    backgroundColors <- grDevices::colorRampPalette(c("#FFFFFF","#FF000077"),
                                               alpha=TRUE)(49)[
                                                 order(annotation[[group]])]
  }
  # Plot the others with background annotation
  for(group in groupnames[-1]){
    fsom <- UpdateNodeSize(fsom,reset=TRUE)
    fsom$MST$size <- fsom$MST$size * groups$means_norm[[group]]
    PlotStars(fsom, backgroundValues = annotation[[group]],
              main = group,
              backgroundColor = backgroundColors,
              ...)
  }
  
  annotation
}

###################
## QueryStarPlot ##
###################
#' Query a certain cell type
#' 
#' Identify nodes in the tree which resemble a certain profile of "high"
#' or "low" marker expressions.
#'
#' @param fsom    FlowSOM object, as generated by \code{\link{BuildMST}} or
#'        the first list item of \code{\link{FlowSOM}}
#' @param query  Array containing "high" or "low" for the specified column
#'               names of the FlowSOM data
#' @param plot   If true, a plot with a gradient of scores for the nodes is 
#'               shown
#' @param color  Color to use for nodes with a high score in the plot
#' @param debug  If TRUE, some extra output will be printed
#' @param ...    Other parameters to pass to \code{\link{PlotStars}}
#' @return A list, containing the ids of the selected nodes, the individual
#'         scores for all nodes and the scores for each marker for each node
#' @examples
#'    file <- system.file("extdata", "68983.fcs", package="FlowSOM")
#'    # Use the wrapper function to build a flowSOM object (saved in fsom[[1]])
#'    # and a metaclustering (saved in fsom[[2]])
#'    fsom <- FlowSOM(file,compensate = TRUE, transform = TRUE,scale = TRUE,
#'                   colsToUse = c(9,12,14:18), nClus = 10, silent = FALSE,
#'                   xdim=7, ydim=7)
#'    query <- c("PE-Cy7-A" = "high", #CD3
#'               "APC-Cy7-A" = "high", #TCRb
#'               "Pacific Blue-A" = "high") #CD8
#'    query_res <- QueryStarPlot(UpdateNodeSize(fsom[[1]],reset=TRUE), query)
#'    
#'    cellTypes <- factor(rep("Unknown",49),levels=c("Unknown","CD8 T cells"))
#'    cellTypes[query_res$selected] <- "CD8 T cells"
#'    PlotStars(fsom[[1]],
#'                  backgroundValues=cellTypes,
#'                  backgroundColor=c("#FFFFFF00","#ca0020aa"))
#'    
#' @export
QueryStarPlot <- function(fsom, query, plot=TRUE,
                            color="#ca0020",debug=FALSE,...){
    
    scores <- matrix(NA,
                    ncol=length(query),
                    nrow=nrow(fsom$map$medianValues),
                    dimnames=list(NULL,names(query)))
    
    for(marker in names(query)){
        data <- fsom$map$medianValues[,marker]
        if(query[marker]=="high"){
            scores[,marker] = (max(data,na.rm = TRUE) - data)^2
        } else if(query[marker]=="low"){
            scores[,marker] = (data - min(data,na.rm = TRUE))^2
        } else {
            stop("Only high or low marker expressions are 
                supported at the moment")
        }
    }
    
    if(debug){
        message("Scores")
        print(scores)
    }
    
    # Normalize between 0 and 1 and make highest score best
    scores <- apply(scores,2,function(x){1-((x-min(x,na.rm=TRUE))/
                                    (max(x,na.rm=TRUE)-min(x,na.rm=TRUE)))})
    if(debug){
        message("Normalized")
        print(scores)
    }
    
    nodeScores <- apply(scores,1,mean)
    nodeScores[is.na(nodeScores)] <- 0
    if(debug){
        message("nodeScores")
        print(nodeScores)
    }
    
    if(plot){
        PlotStars(fsom, 
                markers = names(query),
                backgroundValues = nodeScores,
                backgroundColor = grDevices::colorRampPalette(
                    c("#ffffff00","#ffffff00","#ffffff00",color)),
                query = query,
                ...)
    }

    cutoff <- max(nodeScores) *0.95

    scores_cutoff <- nodeScores > cutoff

    return(list("selected"=which(scores_cutoff),
                "nodeScores"=nodeScores,
                "fullScores"=scores))
}

Try the FlowSOM package in your browser

Any scripts or data that you put into this service are public.

FlowSOM documentation built on Nov. 8, 2020, 6:40 p.m.