R/kinematic_analysis.R

In impres-lab/leafkin: Leaf kinematic analysis in R

#' Perform the kinematic analysis.
#'
#' Combines calculated leaf elongation rates, fitted cell sizes and meristem
#' sizes to perform a kinematic analysis, resulting in data about cell division
#' and cell elongation.
#'
#' This function takes in the data generated by \link{calculate_LER},
#' \link{get_all_fitted_cell_lengths} and a data.frame or tibble containing the
#' meristem size of each plant involved in the microscopy analysis. The first
#' column of the meristem size file should contain plant_id's and should be
#' named exactly plant_id. The second column should contain meristem sizes in
#' micrometre and should be named exactly mer_length_um. It are these column
#' names that the function will look for and use to select the data.
#'
#' @param LER_means contains a tidy tibble with the previoulsy calculated mean leaf elongation rates for each plant.
#' @param tidy_cell_lengths contains a tidy tibble with the previously fitted cell lengths for each plant.
#' @param meristem_size_micrometre contains a tidy tibble with the meristem sizes in micrometre for each plant.
#'
#' @author Jonas Bertels
#'
#' @import tidyr
#' @import tibble
#' @import dplyr
#' @import tools
#'
#' @return A tidy tibble, containing all calculated kinematic parameters.
#' @export
#'


kinematic_analysis <- function(LER_means,
                               tidy_cell_lengths,
                               meristem_size_micrometre) {

  # check whether the user has provided a data.frame or tbl (tibble) for meristem sizes
  if(sum(class(meristem_size_micrometre) %in% c("data.frame", "tbl", "tbl_df")) == 0) { # If file does not exist
    stop("\n ##### \n Meristem data is not provided as data.frame or tbl (tibble).\n #####")
  }

  # Check column headers
  if(names(meristem_size_micrometre)[1] != "plant_id") { # Check whether the first column is named plant_id
    stop("\n ##### \n The first column of the meristem size file should contain plant_id's and should be named exactly plant_id.\n #####")
  }

  if(names(meristem_size_micrometre)[2] != "mer_length_um") { # Check whether the second column is named mer_length_um
    stop("\n ##### \n The second column of the meristem size file should contain meristem lengths in micrometer and should be named exactly mer_length_um. \n #####")
  }


  # Check whether mer_length_um is all nummeric
  if(!is.numeric(meristem_size_micrometre$mer_length_um)) {
    stop("\n ##### \n Please check the entered meristem sizes. Not all of them seem to be numeric.\n #####")
  }
  # check whether all plant_id's are unique
  if(length(meristem_size_micrometre$plant_id) != length(unique(meristem_size_micrometre$plant_id))){
    stop("\n ##### \n There are duplicate plant_id's in meristem size file. Please use unique plant_id's for each plant.\n #####")
  }
  # check whether the meristem plant_id's match the cell length plant_id's
  if(!setequal(meristem_size_micrometre$plant_id, unique(tidy_cell_lengths$plant_id))){
    stop("\n ##### \n The plant_id's of the meristem size file do not match with the plant_id's of the cell length measurements. The set of plants for the meristem length measurements and cell length measurements must be identical.\n #####")
  }

  # Check whether LER_means has the correct format
  if(!setequal(colnames(LER_means), c("plant_id", "mean_plant_LER"))) {
    stop("\n ##### \n Please check the LER_mean columns. They should contain the following columns: 'plant_id' and 'mean_plant_LER'. \n #####")
  }
  # Check whether tidy_cell_lengths has the correct format
  if(!setequal(colnames(tidy_cell_lengths), c("plant_id","position","cell_length"))) {
    stop("\n ##### \n Please check the tidy_cell_lengths columns. They should contain the following columns: 'plant_id', 'position' and 'cell_length'. \n #####")
  }
  # Check whether all tidy_cell_lengths plant_id's are also in the LER_means plant_id's
  if(!all(unique(tidy_cell_lengths$plant_id) %in% unique(LER_means$plant_id))) {
    stop("\n ##### \n Not all plant_id's in tidy_cell_lengths are represented in the plant_id's of LER_means. Both measurements are needed to perform a kinematic analysis. Please make sure plant_id's are exactly the same as well.  \n #####")
  }


  # Initiate tibble where all the results will be stored
  no_of_plants <- length(levels(as.factor(tidy_cell_lengths$plant_id)))
  all_results <- tibble(
    plant_id = levels(as.factor(tidy_cell_lengths$plant_id)),
    LER = vector(mode = "numeric", length = length(no_of_plants)),
    meristem_size = vector(mode = "numeric", length = length(no_of_plants)),
    length_of_growth_zone = vector(mode = "numeric", length = length(no_of_plants)),
    length_of_elongation_zone = vector(mode = "numeric", length = length(no_of_plants)),
    mature_cell_length = vector(mode = "numeric", length = length(no_of_plants)),
    cell_production_rate = vector(mode = "numeric", length = length(no_of_plants)),
    number_of_cells_in_the_meristem = vector(mode = "numeric", length = length(no_of_plants)),
    number_of_cells_in_the_elongation_zone = vector(mode = "numeric", length = length(no_of_plants)),
    number_of_cells_in_the_growth_zone = vector(mode = "numeric", length = length(no_of_plants)),
    cell_division_rate = vector(mode = "numeric", length = length(no_of_plants)),
    cell_cycle_duration = vector(mode = "numeric", length = length(no_of_plants)),
    time_in_elongation_zone = vector(mode = "numeric", length = length(no_of_plants)),
    time_in_division_zone = vector(mode = "numeric", length = length(no_of_plants)),
    length_of_cells_leaving_meristem = vector(mode = "numeric", length = length(no_of_plants)),
    relative_cell_elongation_rate = vector(mode = "numeric", length = length(no_of_plants))
  )

  # Loop over all plants of the microscopy study to calculate all parameters

  for (plant in levels(as.factor(all_results$plant_id))) {

    # Filter cell length profile
    filtered_plant <- tidy_cell_lengths %>%
      filter(plant_id == plant)

    # Find the row number of this plant in the all_results tibble
    row_number_of_plant <- which(all_results$plant_id == plant)

    # Find the matching LER and enter it in all_results tibble
    all_results$LER[row_number_of_plant] <- LER_means$mean_plant_LER[which(LER_means$plant_id == plant)]

    # Find the matching meristem size and enter it in all_results tibble
    all_results$meristem_size[row_number_of_plant] <- meristem_size_micrometre$mer_length_um[which(meristem_size_micrometre$plant_id == plant)]

    # Calculate cummulative cell number and return it as a vector
    result_cumm_cell_no <- cumm_cell_no(filtered_plant$cell_length, filtered_plant$position)
    result_cumm_cell_no

    # Calculate 95 perc value
    result_perc_95 <- perc_95(size_vect <- filtered_plant$cell_length)
    result_perc_95

    # Find growth zone size
    result_growth_zone <- growth_zone(pos_vect = filtered_plant$position,
                                      size_vect = filtered_plant$cell_length,
                                      perc95_vect = result_perc_95)
    all_results$length_of_growth_zone[row_number_of_plant] <- result_growth_zone

    # Calculate elongation zone size
    all_results$length_of_elongation_zone[row_number_of_plant] <- all_results$length_of_growth_zone[row_number_of_plant] - (all_results$meristem_size[row_number_of_plant]/1000)

    # Mature cell length
    result_mat_cell_length <- mat_cell_length(growth_zone_size = result_growth_zone,
                                              pos_vect = filtered_plant$position,
                                              cell_lengths = filtered_plant$cell_length)
    all_results$mature_cell_length[row_number_of_plant] <- result_mat_cell_length

    # Cell production rate
    result_cell_prod_rate <- cell_prod_rate(LER = all_results$LER[row_number_of_plant],
                                            mat_cell_length = result_mat_cell_length)
    all_results$cell_production_rate[row_number_of_plant] <- result_cell_prod_rate

    # number of cells in the meristem
    result_mer_cell_no <- value_at_end_of_mer(mer_size = all_results$meristem_size[row_number_of_plant],
                                              value = result_cumm_cell_no,
                                              pos_vect = filtered_plant$position)
    all_results$number_of_cells_in_the_meristem[row_number_of_plant] <- result_mer_cell_no

    # length of cells leaving the meristem
    result_mer_cell_length <- value_at_end_of_mer(mer_size = all_results$meristem_size[row_number_of_plant],
                                                  value = filtered_plant$cell_length,
                                                  pos_vect = filtered_plant$position)
    all_results$length_of_cells_leaving_meristem[row_number_of_plant] <- result_mer_cell_length

    # number of cells in the growth zone
    result_growth_zone_cell_no <- growth_zone_cell_no(pos_vect = filtered_plant$position,
                                                      cum_cell_no = result_cumm_cell_no,
                                                      growth_zone_size = result_growth_zone)
    all_results$number_of_cells_in_the_growth_zone[row_number_of_plant] <- result_growth_zone_cell_no

    # number of cells in the elongation zone
    result_elong_zone_cell_no <- elong_zone_cell_no(growth_zone_cell_no = result_growth_zone_cell_no,
                                                    mer_cell_no = result_mer_cell_no)
    all_results$number_of_cells_in_the_elongation_zone[row_number_of_plant] <- result_elong_zone_cell_no

    # Average cell division rate
    result_avg_cell_div_rate <- avg_cell_div_rate(cell_prod_rate = result_cell_prod_rate,
                                                  mer_cell_no = result_mer_cell_no)
    all_results$cell_division_rate[row_number_of_plant] <- result_avg_cell_div_rate

    # Cell cycle duration
    result_cell_cycle_duration <- cell_cycle_duration(avg_cell_div_rate = result_avg_cell_div_rate)
    all_results$cell_cycle_duration[row_number_of_plant] <- result_cell_cycle_duration

    # Time in elongation zone
    result_time_in_elo_zone <- time_in_elo_zone(elong_zone_cell_no = result_elong_zone_cell_no,
                                                cell_prod_rate = result_cell_prod_rate)
    all_results$time_in_elongation_zone[row_number_of_plant] <- result_time_in_elo_zone

    # Time in division zone
    result_time_in_div_zone <-  time_in_div_zone(mer_cell_no = result_mer_cell_no,
                                                 cell_cycle_duration = result_cell_cycle_duration)
    all_results$time_in_division_zone[row_number_of_plant] <- result_time_in_div_zone

    # Average cell expansion rate
    result_avg_cell_exp_rate <- avg_cell_exp_rate(mat_cell_length = result_mat_cell_length,
                                                  mer_cell_length = result_mer_cell_length,
                                                  time_in_elo_zone = result_time_in_elo_zone)
    all_results$relative_cell_elongation_rate[row_number_of_plant] <- result_avg_cell_exp_rate
  }
  message("\n ##### \n Analysis finished. All calculated kinematics parameters were returned by the function.\n #####")
  return(all_results)

}