R/ht.R

In dyerlab/gstudio: Tools Related to the Spatial Analysis of Genetic Marker Data

#' Estimate expected heterozygosity
#'  
#' Returns the general expected total heterozygosity parameter
#' @param x A \code{data.frame} object with \code{locus} objects 
#' @param stratum  The name of the column representing the stratum variable (default=Population)
#' @export
#' @author Rodney J. Dyer \email{rjdyer@@vcu.edu}
#' @examples
#' loci1 <- c( locus( c("A","A") ), locus( c("A","A") ), locus( c("A","B")))
#' loci2 <- c( locus( c("A","A") ), locus( c("A","B")), locus(c("B","B")))
#' df <- data.frame( Population=c("One","One","One","Two","Two","Two"), Locus=c(loci1,loci2) )
#' Ht( df )
Ht <- function( x, stratum="Population" ) { 
  
  if( missing(x) || !(stratum %in% names(x))) 
    stop("You need to pass both a data.frame and the name of the stratum to this function")
  
  locus_names <- column_class(x,class="locus")
  if( length(locus_names)==0)
    stop("Cannot estimate expected total heterozygosity if there are no loci...")
  
  ret <- data.frame( Locus=locus_names, Ht=0 )

  ho <- Hos( x, stratum=stratum )
  hs <- Hes( x, stratum=stratum )
  freqs <- frequencies( x, stratum=stratum)
  cts <- genotype_counts(x, stratum)
  K <- nrow(cts)
  
  for( locus in locus_names ) {
    nbar <-  harmonic_mean(cts[[locus]])
    xki <- freqs[ freqs$Locus==locus, ]
    xibar <- unlist(by( xki$Frequency, xki$Allele, function(x) sum(x/K )))
    ht <- 1 - sum(xibar^2) + hs$Hes[hs$Locus==locus]/(K*nbar) - ho$Hos[ho$Locus==locus]/(K*nbar*2)
    ret$Ht[ ret$Locus == locus ] <- ht
  }
  
  return( ret )
}