R/reportGV.R

In rmsharp/nprcmanager: Genetic Tools for Colony Management

#' Generates a genetic value report for a provided pedigree.
#'
## Copyright(c) 2017-2020 R. Mark Sharp
## This file is part of nprcgenekeepr
#' This is the main function for the Genetic Value Analysis.
#'
#' @return A dataframe with the genetic value report. Animals are ranked
#' in order of descending value.
#'
#' @examples
#' \donttest{
#' library(nprcgenekeepr)
#' examplePedigree <- nprcgenekeepr::examplePedigree
#' breederPed <- qcStudbook(examplePedigree, minParentAge = 2,
#'                          reportChanges = FALSE,
#'                          reportErrors = FALSE)
#' focalAnimals <- breederPed$id[!(is.na(breederPed$sire) &
#'                                   is.na(breederPed$dam)) &
#'                                 is.na(breederPed$exit)]
#' ped <- setPopulation(ped = breederPed, ids = focalAnimals)
#' trimmedPed <- trimPedigree(focalAnimals, breederPed)
#' probands <- ped$id[ped$population]
#' ped <- trimPedigree(probands, ped, removeUninformative = FALSE,
#'                     addBackParents = FALSE)
#' geneticValue <- reportGV(ped, guIter = 50, # should be >= 1000
#'                          guThresh = 3,
#'                          byID = TRUE,
#'                          updateProgress = NULL)
#' trimmedGeneticValue <- reportGV(trimmedPed, guIter = 50, # should be >= 1000
#'                                 guThresh = 3,
#'                                 byID = TRUE,
#'                                 updateProgress = NULL)
#' rpt <- trimmedGeneticValue[["report"]]
#' kmat <- trimmedGeneticValue[["kinship"]]
#' f <- trimmedGeneticValue[["total"]]
#' mf <- trimmedGeneticValue[["maleFounders"]]
#' ff <- trimmedGeneticValue[["femaleFounders"]]
#' nmf <- trimmedGeneticValue[["nMaleFounders"]]
#' nff <- trimmedGeneticValue[["nFemaleFounders"]]
#' fe <- trimmedGeneticValue[["fe"]]
#' fg <- trimmedGeneticValue[["fg"]]
#' }
#'
#' @param ped The pedigree information in data.frame format
#' @param guIter Integer indicating the number of iterations for the gene-drop
#'  analysis. Default is 5000 iterations
#' @param guThresh Integer indicating the threshold number of animals for
#' defining a unique allele. Default considers an allele "unique"
#' if it is found in only 1 animal.
#' @param pop Character vector with animal IDs to consider as the population of
#' interest. The default is NULL.
#' @param byID Logical variable of length 1 that is passed through to
#' eventually be used by \code{alleleFreq()}, which calculates the count of each
#'  allele in the provided vector. If \code{byID} is TRUE and ids are provided,
#'  the function will only count the unique alleles for an individual
#'   (homozygous alleles will be counted as 1).
#' @param updateProgress Function or NULL. If this function is defined, it
#' will be called during each iteration to update a
#' \code{shiny::Progress} object.
#' @export
reportGV <- function(ped, guIter = 5000, guThresh = 1, pop = NULL,
                     byID = TRUE, updateProgress = NULL) {
  # Generates a genetic value report for a provided pedigree

  ## If user has limited the population of interest by defining 'pop',
  ## that information is incorporated via the 'population' column.
  ped$population <- getGVPopulation(ped, pop)

  # Get the list of animals in the population to consider
  probands <- as.character(ped$id[ped$population])

  ## Extract genotype data if available otherwise NULL is returned.
  genotype <- getGVGenotype(ped)

  # Generate the kinship matrix and filter down to the animals of interest
  kmat <- filterKinMatrix(probands, kinship(ped$id, ped$sire, ped$dam,
                                            ped$gen))

  # Calculate the mean kinship, and convert to z-scores
  indivMeanKin <- meanKinship(kmat)
  indivMeanKin <- indivMeanKin[probands] # making sure the order is correct
  zScores <- scale(indivMeanKin)

  # Perform the gene drop simulation
  alleles <- geneDrop(ids = ped$id, sires = ped$sire, dams = ped$dam,
                      gen = ped$gen, genotype = genotype, n = guIter,
                      updateProgress = updateProgress)

  if (!is.null(updateProgress)) {
    updateProgress(detail = "Calculating Genome Uniqueness", value = 1,
                   reset = TRUE)
  }

  # Calculate genome uniqueness and order the rows of the returned data.frame
  gu <- calcGU(alleles, threshold = guThresh, byID = byID, pop = probands)
  gu <- gu[probands, , drop = FALSE]

  if (!is.null(updateProgress)) {
    updateProgress(detail = "Calculating Numbers of Offspring", value = 1,
                   reset = TRUE)
  }

  # Get a data.frame of offspring counts for the probands
  offspring <- offspringCounts(probands, ped, considerPop = TRUE)

  include.cols <- intersect(getIncludeColumns(), names(ped))

  # Subsetting out the needed demographic information from the pedigree
  rownames(ped) <- ped$id
  demographics <- ped[probands, include.cols]

  if (!is.null(updateProgress)) {
    updateProgress(detail = "Calculating Founder Equivalents", value = 1,
                   reset = TRUE)
  }

  # Calculating founder equivalents and founder genome equivalents
  feFg <- calcFEFG(ped, alleles)

  # Calculating known founders
  founders <- ped[is.na(ped$sire) & is.na(ped$dam), ]
  males <- founders[(founders$sex == "M") & !grepl("^U", founders$id,
                                                   ignore.case = TRUE), ]
  females <- founders[(founders$sex == "F") & !grepl("^U", founders$id,
                                                     ignore.case = TRUE), ]

  finalData <- cbind(demographics, indivMeanKin, zScores, gu, offspring)
  finalData <- list(report = orderReport(finalData, ped),
                    kinship = kmat,
                    gu = gu,
                    fe = feFg$FE,
                    fg = feFg$FG,
                    maleFounders = males,
                    femaleFounders = females,
                    nMaleFounders = nrow(males),
                    nFemaleFounders = nrow(females),
                    total = (nrow(males) + nrow(females)))
  class(finalData) <- append(class(finalData),"nprcgenekeeprGV")

  return(finalData)
}