R/split_cell.R

#' Split a cell
#'
#' Splits a cell in a lineage tree.
#'
#' There is no guarantee that \code{cell} will be indeed splitted into \code{Nsplit} cells.
#' After the split operation, the resulting cells are sorted by ascending distance
#' between their centroid and the centroid of their mother, if it exists.
#' Otherwise, the order of the resulting cells is random.
#' The original \code{cell} is replaced in both \code{LT} and \code{cell_list} by the first listed cell, which implies best matching.
#' For each one of the remaining resulting cells a single-node (connected) lineage tree (motherless branch) is automatically created.
#'
#' @section Prerequisites:
#' This function can be used by \emph{BaSCA} users \bold{only},
#' importing the data with \code{\link{import_basca}}.
#'
#' @param LT The lineage tree where the cell specified in \code{cell} belongs, an object of class \code{"igraph"}.
#' @param cell The label of the cell in the \code{LT} to be splitted, a character string.
#' It can be any non-root cell, as returned from \code{\link{get_cells}}.
#' @param Nsplit Number of cells for the \code{cell} to be splitted, an integer value \code{>=2}.
#' The default value is \code{2}.
#' @param cell_list A list containing all the cell instants of the movie.
#' @param col_list A list containing all the colony instants of the movie.
#' @param Ncols Number of colonies in the movie, a non-zero positive integer value.
#' @param pixelR The pixel ratio in units of length, a non-zero positive numeric value.
#'
#' @param matFolder A character string naming the absolute path of the directory where
#' the \code{.mat} file generated by \emph{BaSCA} is saved (excluding the last \code{"/"}).
#' The default value is the current working directory \code{getwd()}.
#' \cr\cr
#' NOTE: The components should be separated by \code{"/"} on Windows.
#' @param matFileName A character string naming the \code{.mat} file generated by \emph{BaSCA}
#' (including the suffix \code{".mat"}).
#' The filename is relative to the \code{matFolder}.
#' @param exeFolder A character string naming the absolute path of the installation folder of the MATLAB executable
#' (excluding the last \code{"/"}).
#' \cr\cr
#' NOTE: The components should be separated by \code{"/"} on Windows.
#' @param mcrFolder A character string naming the absolute path of the installation folder of the Matlab Compiler Runtime (MCR)
#' (excluding the last \code{"/"}).
#' \cr\cr
#' NOTE: The components should be separated by \code{"/"} on Windows.
#'
#' @param show A logical value (\code{TRUE} or \code{FALSE}) indicating whether \code{\link{view_cell}}
#' will be called for the mother of the \code{cell} (if it exists), the \code{cell} before the split operation
#' and the resulting cells after the split operation.
#' This capability is useful in order to see the result of the function.
#' The default value is \code{TRUE}.
#'
#' @return A named list with the following components:
#' \item{LT}{The updated LT with the \code{cell} replaced, an object of class \code{"igraph"}.}
#' \item{cell_list}{The updated cell_list with the \code{cell} replaced and the rest resulting cells added.}
#' \item{branches}{A list with the single-node (connected) lineage trees
#' for each one of the rest resulting cells
#' Each branch (element of the list) is an object of class \code{"igraph"}.}
#'
#' @seealso \code{\link{merge_cells}} for the reverse,
#' \code{\link{add_branch}} for connecting a motherless branch to a lineage tree.
#' @export
#' @import igraph
#' @importFrom stats dist

split_cell <- function(LT, cell,
                       Nsplit = 2,
                       cell_list, col_list, Ncols,
                       pixelR,
                       matFolder = getwd(), matFileName,
                       exeFolder, mcrFolder,
                       show = TRUE) {

  #exeFolder = "/usr/BaSCA/mat_correct"
  #mcrFolder = "/usr/local/MATLAB/MATLAB_Compiler_Runtime"

  possible_cells <- get_cells(tree = LT, treeT = "LT", type = "nr")
  if (!(cell %in% possible_cells)) {
    stop(paste("Selected cell", paste("\"", cell, "\"", sep = ""), "does not exist\n"))
  }

  if (Nsplit < 2) {
    stop("Nsplit must be >= 2\n")
  }

  ########################################################

  createSingleCellLT <- function(cell_list_ID) {

    nodesLT <- as.data.frame(cell_list_ID + Ncols + 1)
    colnames(nodesLT) <- "cellID"

    # cellName attribute
    nodesLT[, "cellName"] <- cell_list[[cell_list_ID]]$cellName

    # numeric attributes
    numeric_attrs <- names(cell_list[[1]])[sapply(cell_list[[1]], function(x) class(x) == "numeric" || class(x) == "integer")]
    for (attr in numeric_attrs) {
      if (attr == "colId") {
        nodesLT[, "colId"] <- as.character(cell_list[[cell_list_ID]]$colId)
      } else {
        nodesLT[, attr] <- unlist(unname(cell_list[[cell_list_ID]][attr]))
      }
    }

    # boolean attributes
    boolean_attrs <- names(cell_list[[1]])[sapply(cell_list[[1]], function(x) is.logical(x))]
    for (attr in boolean_attrs) {
      nodesLT[, attr] <- unlist(unname(cell_list[[cell_list_ID]][attr]))
    }

    edgesLT <- data.frame(from = numeric(), to = numeric(), stringsAsFactors = FALSE)
    LT <- graph_from_data_frame(d = edgesLT, directed = TRUE, vertices = nodesLT)

    return(LT)

  }

  ########################################################

  branches <- list()

  mother_cell <- get_cell_fam(tree = LT, treeT = "LT", cell = cell, type = "m")

  if (!is.null(mother_cell)) { # has mother
    if (show) {
      view_cell(LT = LT, cells = mother_cell,
                cell_list = cell_list, col_list = col_list, Ncols = Ncols)
    }
  }

  if (show) {
    view_cell(LT = LT, cells = cell,
              cell_list = cell_list, col_list = col_list, Ncols = Ncols)
  }

  cell_list_ID <- as.numeric(cell) - (Ncols + 1)

  myFrame <- cell_list[[cell_list_ID]]$frame
  myColony <- cell_list[[cell_list_ID]]$colony
  myColId <- cell_list[[cell_list_ID]]$colId
  correctColony <- V(LT)[cell]$colony # if LT is mainLT, then it is already correct

  splitCellProps <- callBascaMatlabExe(Nsplit = Nsplit,
                                       cell_list_IDs = cell_list_ID,
                                       cell_list = cell_list,
                                       col_list = col_list,
                                       mat_folder = matFolder,
                                       mat_fileName = matFileName,
                                       install_folder_exe = exeFolder,
                                       install_folder_MCR = mcrFolder)

  N_splitted <- dim(splitCellProps$cellProps)[3]
  if (N_splitted == 1) {
    stop("Cell could not be splitted\n")
  }

  ####### find order of splitted cells by best-fit with mother ########

  if (!is.null(mother_cell)) { # has mother

    mother_cell_centroid <- cell_list[[as.numeric(mother_cell) - (Ncols + 1)]]$centroid # (row,col) in colImage
    distances <-  NULL
    for (i_cell in 1:N_splitted) {
      cell_centroid <- matrix(colMeans(splitCellProps$cellProps[, , i_cell]$pixelList), nrow = 1, ncol = 2) # (col,row) in colImage
      cell_centroid <- cbind(cell_centroid[, 2], cell_centroid[, 1]) # (row,col) in colImage
      distances <- c(distances, dist(rbind(mother_cell_centroid, cell_centroid)))
    }

    right_order <- order(distances)

  } else {
    right_order <- c(1:N_splitted)
  }

  ######################################################################

  for (i_cell in right_order) {

    if (i_cell == right_order[1]) { # update cell_list, LT for cell

      cell_list[[cell_list_ID]] <- splitCellProps$cellProps[, , i_cell]
      cell_list[[cell_list_ID]] <- updateCellInBascaList(x = cell_list[[cell_list_ID]],
                                                         frame = myFrame,
                                                         colony = myColony,
                                                         colId = myColId,
                                                         col_list = col_list,
                                                         pixelRatio = pixelR)

      LT <- updateCellInBascaLT(LT = LT,
                                cell = cell,
                                cell_list = cell_list,
                                Ncols = Ncols,
                                colony = correctColony)

      if (show) {
        view_cell(LT = LT, cells = cell,
                  cell_list = cell_list, col_list = col_list, Ncols = Ncols)
      }

    } else { # add extra cell to cell_list, add single node LT to list with motherless branches

      curr_cell <- length(cell_list) + 1

      cell_list[[curr_cell]] <- splitCellProps$cellProps[, , i_cell]
      cell_list[[curr_cell]] <- updateCellInBascaList(x = cell_list[[curr_cell]],
                                                      frame = myFrame,
                                                      colony = myColony,
                                                      colId = myColId,
                                                      col_list = col_list,
                                                      pixelRatio = pixelR)

      branches[[length(branches) + 1]] <- createSingleCellLT(cell_list_ID = curr_cell)

      if (show) {
        view_cell(LT = branches[[length(branches)]], cells = as.character(curr_cell + Ncols + 1),
                  cell_list = cell_list, col_list = col_list, Ncols = Ncols)
      }

    }

  }

  cat("Checking updated cell list...\n")
  cell_list <- checkCellList(cell_list = cell_list, col_list = col_list)

  return(list(LT = LT, cell_list = cell_list, branches = branches))

}
vicstefanou/ViSCA documentation built on May 31, 2019, 10:50 p.m.