R/Sequoia_F90wrappers.R

In sequoia: Pedigree Inference from SNPs

#=====================================================================
#=====================================================================

#' @title Fortran Wrapper for Pedigree Reconstruction
#'
#' @description Call main Fortran part of sequoia, and convert its output to a
#'   list with dataframes.
#'
#' @param ParSib either "par" to call parentage assignment, or "sib" to call the
#'   rest of the algorithm.
#' @param FortPARAM a named vector with parameter values, as generated by
#'   \code{ParamToSpecs}.
#' @param GenoM matrix with genotype data, size nInd x nSnp.
#' @param LhIN  life history data: ID - sex - birth year.
#' @param AgePriors matrix with agepriors, size `FortPARAM["nAgeClasses"]` by 8.
#' @param Parents  matrix with rownumbers of assigned parents, size nInd by 2.
#' @param mtDif matrix indicating whether individuals have definitely a
#'   different mitochondrial haplotype (1), or (possibly) the same (0). Size
#'   nInd x nInd.
#' @param DumPfx dummy prefixes
#' @param quiet suppress messages.
#'
#' @return A list with
#' \item{PedigreePar or Pedigree}{the pedigree}
#' \item{DummyIDs}{Info on dummies (not included if parentage-only)}
#' \item{TotLikParents or TotLikSib}{Total log-likelihood per iteration}
#' \item{AgePriorExtra}{Ageprior including columns for grandparental and
#'   avuncular relationships}
#' \item{LifeHistPar or LifeHistSib}{Includes sex and birthyear estimate
#'   inferred from the pedigree for individuals with initially unknown sex
#'   and/or birthyear}.
#'
#' For a detailed description of the output see \code{\link{sequoia}}.
#'
#' @useDynLib sequoia, .registration = TRUE
#'
#' @keywords internal

SeqParSib <- function(ParSib,
                      FortPARAM,
                      GenoM,
                      LhIN,
                      AgePriors,
                      Parents,
                      mtDif,
                      DumPfx,
                      quiet)
{
  on.exit(.Fortran(deallocall), add=TRUE) #, PACKAGE = "sequoia"))

  Ng <- nrow(GenoM)
  gID <- rownames(GenoM)
  LHF <- orderLH(LhIN, gID)
  PedN <- PedToNum(Parents, gID, DoDummies = "no")
  MaxMaxAgePO <- 100  # if changed, change in Fortran too!

  SpecsIntMkPed <- c(switch(ParSib,
                            par = 1,
                            sib = 2),
                     FortPARAM$SpecsIntMkPed)

  TMP <- .Fortran(makeped,
                  # IN:
                  ng = as.integer(Ng),
                  specsintglb = as.integer(FortPARAM$SpecsInt),
                  specsintmkped = as.integer(SpecsIntMkPed),
                  specsdbl = as.double(FortPARAM$SpecsDbl),
                  errv = as.double(FortPARAM$ErrM),
                  genofr = as.integer(GenoM),
                  sexrf = as.integer(LHF$Sex),
                  byrf = as.integer(c(LHF$BirthYear, LHF$BY.min, LHF$BY.max)),
                  lyrf = as.integer(LHF$Year.last),
                  aprf = as.double(AgePriors),
                  mtdif_rf = as.integer(mtDif),
                  # IN/OUT:
                  parentsrf = as.integer(PedN$PedPar),
                  # OUT:
                  lrrf = double(3*Ng),
                  ohrf = integer(3*Ng),
                  nd = integer(2),
                  dumparrf = integer(2*Ng),
                  dumlrrf = double(3*Ng),
                  dumbyrf = integer(3*Ng),
                  totll = double(42),
                  apout = double((3*MaxMaxAgePO)*5*3))
  #                  PACKAGE = "sequoia")

  TMP$lrrf[abs(TMP$lrrf - 999) < 0.1] <- NA
  TMP$dumlrrf[abs(TMP$dumlrrf - 999) < 0.1] <- NA
  TMP$lrrf <- round(TMP$lrrf, 2)
  TMP$dumlrrf <- round(TMP$dumlrrf, 2)
  TMP$ohrf[TMP$ohrf < 0] <- NA

  #=========================
  # Pedigree ----

  PedColNames <- c("id", "dam", "sire", "LLRdam", "LLRsire", "LLRpair",
                   "OHdam", "OHsire", "MEpair", "Sex")

  # parents of genotyped individuals
  NumPed <- data.frame(id = seq_along(gID),
                       VtoM(TMP$parentsrf),
                       VtoM(TMP$lrrf, nc=3),
                       VtoM(TMP$ohrf, nc=3),
                       Sex = 3,   # only relevant for dummies
                       stringsAsFactors=FALSE)
  names(NumPed) <- PedColNames

  # parents of dummies
  if (grepl("sib", ParSib) && any(TMP$nd>0)) {
    NgOdd <- Ng%%2==1
    DumPed <- data.frame(id = -c(seq_len(TMP$nd[1]), seq_len(TMP$nd[2])),
                         VtoM(TMP$dumparrf, sum(TMP$nd), 2, NgOdd),
                         VtoM(TMP$dumlrrf, sum(TMP$nd), 3, NgOdd),
                         OHdam = NA, OHsire = NA, MEpair = NA,
                         Sex = rep(1:2, TMP$nd),
                         stringsAsFactors=FALSE)
    names(DumPed) <- PedColNames

    if (any(NumPed$dam < -1e6)) {   # hermaphrodite dummies
      IsHermDumDam <- with(DumPed, id %in% (NumPed$dam + 1e6) & Sex==1)
      DumPed$id[IsHermDumDam] <- DumPed$id[IsHermDumDam] - 1e6
      IsHermDumSire <- with(DumPed, (!id %in% NumPed$sire) & Sex==2)
      DumPed <- DumPed[!IsHermDumSire, ]
    }
    NumPed <- rbind(NumPed, DumPed)
  }

  # numeric to character IDs
  Pedigree <- NumPed
  Pedigree[,"id"] <- NumToID(Pedigree[,"id"], k = Pedigree[,"Sex"], gID, DumPfx)
  for (k in 1:2) Pedigree[, k+1] <- NumToID(Pedigree[, k+1], k, gID, DumPfx)

  #=========================
  # all info for each dummy
  if (grepl("sib", ParSib) && any(TMP$nd>0)) {
    NumOff <- list("mat" = table(Pedigree$dam[NumPed$dam < 0]),
                   "pat" = table(Pedigree$sire[NumPed$sire < 0]))
    MaxOff <- max(unlist(NumOff))

    OffIDs <- c(plyr::dlply(Pedigree, "dam", function(df) df$id),
                plyr::dlply(Pedigree, "sire", function(df) df$id))
    OffIDs <- OffIDs[c(names(NumOff[[1]]), names(NumOff[[2]]))]
    # includes dummy offspring

    DummyIDs <- plyr::join(data.frame(id = c(names(NumOff[["mat"]]),
                                             names(NumOff[["pat"]])),
                                      stringsAsFactors=FALSE),
                           Pedigree[, c("id", "dam", "sire")],
                           by="id")
    DummyIDs <- cbind(DummyIDs,
                      Sex = rep(1:2, times=TMP$nd),
                      VtoM(TMP$dumbyrf, sum(TMP$nd),3, NgOdd),
                      unlist(NumOff),
                      plyr::laply(OffIDs, .fun = function(x) x[1:MaxOff], .drop=FALSE),
                      row.names=NULL)
    names(DummyIDs)[5:ncol(DummyIDs)] <- c("BY.est", "BY.lo", "BY.hi", "NumOff",
                                           paste0("O", 1:MaxOff))
  } else  DummyIDs <- NULL

  Pedigree <- Pedigree[, colnames(Pedigree)!="Sex"]


  #=========================
  # returned ageprior w columns for GP + avuncular

  if (grepl("par", ParSib) | grepl("sib", ParSib)) {
    APM <- VtoM(TMP$apout, nc=15)
    colnames(APM) <- c("M", "P", "FS", "MS", "PS",
                       "MGM", "MGF", "MFA", "MMA", "MPA",
                       "PGM", "PGF", "PFA", "PMA", "PPA")
    MaxAgePO <- ifelse(any(APM[,"FS"]==0),
                       mean(which(APM[,"FS"]>0)) -1,  # sibs are always symmetrical around 0
                       min(which(APM[,"M"]>0)) -2)  # may be flat, tail may be cut off
    rownames(APM) <- c(1:nrow(APM)) -MaxAgePO -1
    MaxRow <- min(max(which(apply(APM, 1, function(x) any(x>0)))) +1, nrow(APM))
    APM <- round(APM[1:MaxRow, ], 3)
  }


  #=========================
  # update lifehist w. inferred sex & estimated birth year

  LhOUT <- data.frame(LHF,
                      Sexx = TMP$sexrf,
                      BY.est = TMP$byrf[1:Ng],
                      BY.lo = TMP$byrf[1:Ng + Ng],
                      BY.hi = TMP$byrf[1:Ng + 2*Ng],
                      stringsAsFactors = FALSE)
  LhOUT$BY.est[LhOUT$BY.est < 0] <- NA
  names(LhOUT)[names(LhOUT)=="ID"] <- "id"

  #=========================
  # output

  if (!quiet) {
    nAss <- apply(Pedigree[,c("dam", "sire")], 2, function(x) sum(!is.na(x)))
    cli::cli_alert_success(cli::col_green(
      c("assigned {nAss[1]} dams and {nAss[2]} sires to {nrow(unique(GenoM))}",
        '{ifelse(ParSib=="sib", paste(" +", sum(TMP$nd)), "")} individuals',
        '{ifelse(ParSib=="sib", " (real + dummy)", "")} \n')))
  }

  rownames(Pedigree) <- 1:nrow(Pedigree)

  TotLik <- TMP$totll[seq_len(sum(TMP$totll!=0))]
  names(TotLik) <- seq_along(TotLik) -1

  if (grepl("par", ParSib)) {
    OUT <- list(PedigreePar = Pedigree,
                TotLikPar = TotLik,
                AgePriorExtra = APM,
                LifeHistPar = LhOUT)

  } else if (grepl("sib", ParSib)) {
    OUT <- list(Pedigree = Pedigree,
                DummyIDs = DummyIDs,
                TotLikSib = TotLik,
                AgePriorExtra = APM,
                LifeHistSib = LhOUT)
  }
  return(OUT[!sapply(OUT, is.null)])
}


#============================================================================
#============================================================================