R/Phenot.getProgenyStatusLookup.R

In aleczwart/PolyPatEx: Paternity Exclusion in Autopolyploid Species

##  Function Phenot.getProgenyStatusLookup.R
##
## For a specified ploidy, generate the progeny-mother relationship
## lookup table.
##
## The \sQuote{progeny status lookup table} describes whether a
## particular non-trivial combination of mother and progeny
## phenotypes are compatible, in the sense that at least one maternal
## gamete can be generated from the maternal phenotype which is
## compatible with at least one progeny genotype that can be
## generated from the progeny phenotype.
##
## \code{Phenot.getProgenyStatusLookup} generates the progeny status
## lookup table for the non-trivial combinations of maternal and
## progeny allele counts generated by
## \code{\link{getPossMPPhenotCounts}}.
##
## Status codes produced by \code{Phenot.getProgenyStatusLookup} are:
##
## \describe{
##
## \item{\code{MP.noMatch}}{Mother and progeny's allele sets fail to match
##                   this mother cannot generate a gamete compatible
##                   with the alleles in the progeny's phenotype.}
##
## \item{\code{MAO}}{Mother and progeny are compatible and the progeny's
##            phenotype contains only alleles that are present in the
##            mother phenotype (\sQuote{Mother Alleles Only}).}
##
## \item{\code{NMA}}{Mother and progeny are compatible and the progeny
##            contains alleles that are not present in the mother
##            (\sQuote{Non-Mother Alleles}).}
##
## }
##
## @title Get the progeny-mother relationship lookup table
## @param p the ploidy, an even number greater than or equal to 4.
## @return The progeny status table, as a data frame.
##
## The data frame contains the following columns:
##
## \describe{
##
## \item{\code{nM}}{the number of alleles in the maternal phenotype.}
##
## \item{\code{nP}}{the number of alleles in the progeny's
## phenotype.}
##
## \item{\code{nMP}}{the number of alleles that appear in both
## maternal and progeny phenotypes.}
##
## \item{\code{MPStatus}}{Mother-progeny comparison status.  Either
## \sQuote{MP.noMatch},\sQuote{MAO}, or \sQuote{NMA} - see the
## \sQuote{details} section for details.}
##
## }
##
## \sQuote{Trivial} cases where nM, nP, or nMP are zero, are excluded
## from the table.
##
## @author Alexander Zwart (alec.zwart at csiro.au)
## @examples
## \dontrun{
##
## Phenot.getProgenyStatusLookup(p=4)
##
## }
##
Phenot.getProgenyStatusLookup <- function(p) {
  ##Phenot prefix to avoid confusion between dataType cases
  aa <- array(NA_character_,
              dim=c(p,p,p),
              dimnames=list(
              M=paste("nM",1:p,sep=""),
              P=paste("nP",1:p,sep=""),
              MP=paste("nMP",1:p,sep="")
              ))
  ## Get all possible genotype combos
  apMPg <- Phenot.getAllPossMPGenots(p)
  ## Loop over all possible genotype combos
  for (i in 1:length(apMPg)) {
    ll <- apMPg[[i]]  ## For brevity of code
    matchFound <- FALSE
    for (iM in 1:length(ll$MGenots)) {
      for (iP in 1:length(ll$PGenots)) {
        ## Compare genotypes - can the mum provide a valid gamete?
        UMGenot <- make.unique(ll$MGenots[[iM]])
        UPGenot <- make.unique(ll$PGenots[[iP]])
        if (sum(UPGenot %in% UMGenot) >= p/2) {
          matchFound <- TRUE
          break  ##Break out of iP loop
        }
      }
      if (matchFound) break  ##Break out of iM loop
    }
    if (matchFound) {
      if (ll$nMP < ll$nP) { ## Non-maternal alleles in progeny
        aa[ll$nM,ll$nP,ll$nMP] <- "NMA"
      } else {  ## nMP = nP, only maternal alleles in progeny
        aa[ll$nM,ll$nP,ll$nMP] <- "MAO"
      }
    } else {  ##match not found
      aa[ll$nM,ll$nP,ll$nMP] <- "MP.noMatch"
    }
  }
  ##Inefficient, but safe!
  bb <- data.frame(nM = rep(1:p,p*p),
                   nP = rep(rep(1:p,each=p),p),
                   nMP = rep(1:p,each=p*p),
                   MPStatus = NA)
  for (i in 1:p^3) {
    bb$MPStatus[i] <- aa[bb$nM[i],bb$nP[i],bb$nMP[i]]
  }
  ## Strip out the NA rows - these correspond to nMP > min(nM,nP)
  bb <- bb[!is.na(bb$MPStatus),]
  ## Add in the cases nMP = 0:
  bb <- rbind(bb,
              data.frame(nM = rep(1:p,p),
                         nP = rep(1:p,each=p),
                         nMP = rep(0,p*p),
                         MPStatus = rep("MP.noMatch",p*p)))
  rownames(bb) <- 1:(dim(bb)[1])
  return(bb)
}