R/GetReplStats_func_20200914.R

In MHCtools: Analysis of MHC Data in Non-Model Species

#' GetReplStats function
#'
#' \code{\link{GetReplStats}} uses the output files produced by the ReplMatch()
#' function to calculate statistics on the agreement between replicated samples
#' in the sequencing experiment.
#'
#' If you publish data or results produced with MHCtools, please cite both of
#' the following references:
#' Roved, J. (2022). MHCtools: Analysis of MHC data in non-model species. Cran.
#' Roved, J. (2024). MHCtools 1.5: Analysis of MHC sequencing data in R. In S.
#' Boegel (Ed.), HLA Typing: Methods and Protocols (2nd ed., pp. 275–295).
#' Humana Press. https://doi.org/10.1007/978-1-0716-3874-3_18
#'
#' @param filepath is a user defined path to the folder where the output files
#'   from the ReplMatch() function have been saved.
#' @return  A list containing the number of replicate sets with zero incongruent
#'   sequences, the proportion of replicate sets with zero incongruent
#'   sequences, the mean of the mean proportion of incongruent sequences across
#'   all replicate sets, and the repeatability of the sequencing experiment.
#' @seealso \code{\link{ReplMatch}}; \code{\link{GetReplTable}}
#' @examples
#' filepath <- system.file("extdata/ReplMatchOut/", package="MHCtools")
#' GetReplStats(filepath)
#' @export

GetReplStats <- function(filepath) {

  # Get the file names of the .Rds output generated by the ReplMatch function
  file_names <- dir(filepath)

  mean_props <- numeric()

  # Extract the mean proportion of incongruent sequences for each replicate set
  for(i in 1:length(file_names)) {

    mean_props[i] <- readRDS(file.path(filepath, file_names[i]))$Mean_prop_incongr_seqs

  }

  # Calculate the number of replicate sets with zero incongruent sequences
  N_zeros <- length(which(mean_props==0))

  # Calculate the proportion of replicate sets with zero incongruent sequences
  Pr_zeros <- length(which(mean_props==0))/length(mean_props)

  # Calculate the mean of the mean proportions of incongruent sequences for all replicate sets
  MoM_props <- mean(mean_props)

  # Calculate the repeatability
  Rep <- 1-mean(mean_props)

  list(Number_repl_sets_w_zero_inc=c(N_zeros), Prop_repl_sets_w_zero_inc=c(Pr_zeros), Mean_of_mean_props=c(MoM_props), Repeatability=c(Rep))

}