R/ped_stats_summary.R

In pedtricks: Visualize, Summarize and Simulate Data from Pedigrees

#' Post-processes output from ped_stats
#'
#' Generates a manageable summary of pedigree-wide statistics
#' reported by ped_stats, either for a single pedigree or
#' for a comparison between pedigrees
#'
#' @param object An object of class ped_stats generated by \code{ped_stats}
#' @param ... extra arguments
#'
#' @return
#'   Returns a table of numbers of records, maternities,
#' paternities, pairwise sibship relationships, numbers
#' of different classes of grand-parental relationships,
#' pedigree depth, number of founders, mean sibship sizes,
#' simple statistics of numbers of inbred and non-inbred
#' individuals, and proportions of pairwise relationship
#' coefficients equal to or greater than several
#' thresholds.
#'
#' @export

summary.ped_stats <- function(object, ...) {
  sumData <- array(NA, dim = 22)

  names(sumData) <- c("records ", "maternities", "paternities", "full sibs", "maternal sibs", "maternal half sibs", "paternal sibs", "paternal half sibs", "maternal grandmothers", "maternal grandfathers", "paternal grandmothers", "paternal grandfathers", "maximum pedigree depth", "founders", "mean maternal sibsip size", "mean paternal sibsip size", "non-zero F", "F > 0.125", "mean pairwise relatedness", "pairwise relatedness>=0.125", "pairwise relatedness>=0.25", "pairwise relatedness>=0.5")

  sumData[1] <- object$totalSampleSize
  sumData[2] <- object$totalMaternities
  sumData[3] <- object$totalPaternities
  sumData[4] <- object$totalFullSibs
  sumData[5] <- object$totalMaternalSibs
  sumData[6] <- object$totalMaternalSibs - object$totalFullSibs
  sumData[7] <- object$totalPaternalSibs
  sumData[8] <- object$totalPaternalSibs - object$totalFullSibs
  sumData[9] <- object$totalMaternalGrandmothers
  sumData[10] <- object$totalMaternalGrandfathers
  sumData[11] <- object$totalPaternalGrandmothers
  sumData[12] <- object$totalPaternalGrandfathers
  sumData[13] <- max(as.numeric(names(object$pedigreeDepth)))
  sumData[14] <- object$pedigreeDepth[1]
  sumData[15] <- mean(object$maternalSibships[, 2])
  sumData[16] <- mean(object$paternalSibships[, 2])
  sumData[17] <- sum(object$inbreedingCoefficients != 0 + 0)
  sumData[18] <- sum(object$inbreedingCoefficients > 0.125 + 0)
  rc <- subset(object$relatednessCategories, is.na(object$relatednessCategories) == FALSE)
  sumData[19] <- weighted.mean(as.numeric(names(rc)), rc, na.rm = TRUE)
  sumData[20] <- sum(subset(rc, as.numeric(names(rc)) >= 0.125)) / sum(rc)
  sumData[21] <- sum(subset(rc, as.numeric(names(rc)) >= 0.25)) / sum(rc)
  sumData[22] <- sum(subset(rc, as.numeric(names(rc)) >= 0.5)) / sum(rc)

  data.frame(sumData)
}

#' @rdname pedantics-deprecated
#' @section \code{pedStatSummary}: the function has been simplified but only the functionality are still available via \code{ped_stats} and its summary and plot methods
#' @export
pedStatSummary <- function() {
  .Deprecated(summary.ped_stats,
    msg = "this function from pedantics is deprecated and not working anymore. Please use 'ped_stats()' with its summary instead",
  )
}

#' Plot output from ped_stats
#'
#' Generates a manageable summary of pedigree-wide statistics
#' reported by ped_stats, either for a single pedigree or
#' for a comparison between pedigrees
#'
#' @param x An object of class ped_stats generated by \code{ped_stats}
#' @param lowMem If TRUE, then stats based on calculation of A are not performed.
#' @param grContrast If TRUE, then uglier shades of red and blue are used to denote male and female statistics in graphical reports, but these colours provide better contrast in greyscale.
#' @param ... extra arguments
#'
#' @return
#'   Returns a table of numbers of records, maternities,
#' paternities, pairwise sibship relationships, numbers
#' of different classes of grand-parental relationships,
#' pedigree depth, number of founders, mean sibship sizes,
#' simple statistics of numbers of inbred and non-inbred
#' individuals, and proportions of pairwise relationship
#' coefficients equal to or greater than several
#' thresholds.
#' @keywords plot
#' @export
#'
plot.ped_stats <- function(x, lowMem = FALSE, grContrast = FALSE, ...) {
  col1 <- "red"
  col2 <- "blue"
  if (grContrast == TRUE) {
    col1 <- colors()[117]
    col2 <- colors()[109]
  }
  cohorts <- x$cohorts

  if (is.null(cohorts) == FALSE & lowMem == FALSE) {
    cohortRelatedness <- as.data.frame(as.matrix(x$meanRelatednessAmongCohorts))
    cohortTakeOneRelatedness <- array(dim = length(cohortRelatedness[1, ]))
    for (i in 1:(length(cohortTakeOneRelatedness) - 1)) {
      cohortTakeOneRelatedness[i + 1] <- cohortRelatedness[i, i + 1]
    }
    cohortTakeTwoRelatedness <- array(dim = length(cohortRelatedness[1, ]))
    for (i in 1:(length(cohortTakeOneRelatedness) - 2)) {
      cohortTakeTwoRelatedness[i + 2] <- cohortRelatedness[i, i + 2]
    }
    cohortTakeThreeRelatedness <- array(dim = length(cohortRelatedness[1, ]))
    for (i in 1:(length(cohortTakeOneRelatedness) - 3)) {
      cohortTakeThreeRelatedness[i + 3] <- cohortRelatedness[i, i + 3]
    }
    cohortTakeFourRelatedness <- array(dim = length(cohortRelatedness[1, ]))
    for (i in 1:(length(cohortTakeOneRelatedness) - 4)) {
      cohortTakeFourRelatedness[i + 4] <- cohortRelatedness[i, i + 4]
    }
    
    opar <- par(no.readonly = TRUE)
    on.exit(par(opar))

    par(oma = c(5, 1, 1, 1))

    plot(as.numeric(names(cohortRelatedness)), cohortTakeOneRelatedness, type = "l", xlab = "Cohort", ylab = "Pairwise mean cohort relatedness")
    lines(as.numeric(names(cohortRelatedness)), cohortTakeTwoRelatedness, lty = "dashed")
    lines(as.numeric(names(cohortRelatedness)), cohortTakeThreeRelatedness, col = "gray")
    lines(as.numeric(names(cohortRelatedness)), cohortTakeFourRelatedness, col = "gray", lty = "dashed")
    mtext("Mean relatedness between individuals born 1", side = 1, line = 5)
    mtext("(black, solid), 2 (black, dashed), 3 (gray,", side = 1, line = 6)
    mtext("solid) and 4 (gray, dashed) cohorts apart.", side = 1, line = 7)
    inp <- readline(prompt = "Press <s> to save current plot or press <Enter> to continue...")
    if (inp == "s") {
      s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
      savePlot(paste(s, ".jpeg", sep = ""), type = "jpeg")
      readline(prompt = "File saved.  Press <Enter> to continue...")
    }
  }

  #  inbredSubset<-as.data.frame(cbind(x$cohorts,x$inbreedingCoefficients))
  inbredSubset <- subset(x$inbreedingCoefficients, x$inbreedingCoefficients > 0)
  proportionInbred <- length(inbredSubset) / x$totalSampleSize
  hist(as.numeric(as.character(inbredSubset)), xlab = "Inbreeding coefficient", ylab = "Count", main = "")
  mtext("Distribution of inbreeding coefficients, as ", side = 1, line = 5)
  mtext("estimated from the pedigree, among the", side = 1, line = 6)
  mtext("ndividuals in the pedigree with F>0.", side = 1, line = 7)
  inp <- readline(prompt = "Press <s> to save current plot or press <Enter> to continue...")
  if (inp == "s") {
    s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
    savePlot(paste(s, ".jpeg", sep = ""), type = "jpeg")
    readline(prompt = "File saved.  Press <Enter> to continue...")
  }


  if (is.null(cohorts) == FALSE) {
    plot(as.numeric(as.character(x$cohorts)), as.numeric(as.character(x$inbreedingCoefficients)), xlab = "Inbreeding coefficient", ylab = "Count (note: zeros excluded)")
    mtext("Distribution of inbreeding coefficients, as ", side = 1, line = 5)
    mtext("estimated from the pedigree, by cohort.", side = 1, line = 6)

    inp <- readline(prompt = "Press <s> to save current plot or press <Enter> to continue...")
    if (inp == "s") {
      s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
      savePlot(paste(s, ".jpeg", sep = ""), type = "jpeg")
      readline(prompt = "File saved.  Press <Enter> to continue...")
    }
  }

  plot(as.numeric(names(x$pedigreeDepth)), as.numeric(as.character(x$pedigreeDepth)), type = "l", xlab = "Pedigree depth", ylab = "Count")
  mtext("Distribution of pedigree depth.", side = 1, line = 5)
  inp <- readline(prompt = "Press <s> to save current plot or press <Enter> to continue...")
  if (inp == "s") {
    s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
    savePlot(paste(s, ".jpeg", sep = ""), type = "jpeg")
    readline(prompt = "File saved.  Press <Enter> to continue...")
  }

  if (is.null(cohorts) == FALSE) {
    plot(as.numeric(colnames(x$cumulativePedigreeDepth)),
      as.numeric(as.character(x$cumulativePedigreeDepth[1, ])),
      type = "l",
      ylim = c(0, max(x$cumulativePedigreeDepth)), xlab = "Pedigree depth", ylab = "Count"
    )
    for (i in 2:length(x$cumulativePedigreeDepth[, 1])) {
      lines(
        as.numeric(colnames(x$cumulativePedigreeDepth)),
        as.numeric(as.character(x$cumulativePedigreeDepth[i, ]))
      )
    }
    mtext("Cumulative distribution of pedigree depth.", side = 1, line = 5)
    mtext("Earliest cohorts are represented by the lowest lines.", side = 1, line = 6)
    inp <- readline(prompt = "Press <s> to save current plot or press <Enter> to continue...")
    if (inp == "s") {
      s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
      savePlot(paste(s, ".jpeg", sep = ""), type = "jpeg")
      readline(prompt = "File saved.  Press <Enter> to continue...")
    }
  }

  if (is.null(cohorts) == FALSE) {
    thickness <- 1
    if (grContrast == TRUE) thickness <- 1.7
    sibCohortData <- as.data.frame(cbind(
      x$fullSibsByCohort,
      x$maternalSibsByCohort,
      x$paternalSibsByCohort,
      x$maternalSibsByCohort
        - x$fullSibsByCohort,
      x$paternalSibsByCohort
        - x$fullSibsByCohort
    ))
    plot(as.numeric(rownames(x$fullSibsByCohort)), sibCohortData[, 1], type = "l", ylim = c(0, max(sibCohortData)), xlab = "Cohort", ylab = "Count", lwd = thickness)
    lines(as.numeric(rownames(x$fullSibsByCohort)), sibCohortData[, 2], lty = "dashed", col = col1, lwd = thickness)
    lines(as.numeric(rownames(x$fullSibsByCohort)), sibCohortData[, 4], col = col1, lwd = thickness)
    lines(as.numeric(rownames(x$fullSibsByCohort)), sibCohortData[, 3], lty = "dashed", col = col2, lwd = thickness)
    lines(as.numeric(rownames(x$fullSibsByCohort)), sibCohortData[, 5], col = col2, lwd = thickness)
    mtext("Known sib pairs by cohort.  Black line: full sibs, red: maternal sibs,", side = 1, line = 5)
    mtext("blue: paternal sibs, dashed: total maternal and paternal", side = 1, line = 6)
    mtext("bs, solid coloured: maternal and paternal half sibs.", side = 1, line = 7)
    inp <- readline(prompt = "Press <s> to save current plot or press <Enter> to continue...")
    if (inp == "s") {
      s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
      savePlot(paste(s, ".jpeg", sep = ""), type = "jpeg")
      readline(prompt = "File saved.  Press <Enter> to continue...")
    }

    if (inp == "e") {
      s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
      postscript(paste(s, ".eps", sep = ""), width = 8, height = 8, horizontal = FALSE)
      par(oma = c(5, 1, 1, 1))
      plot(as.numeric(rownames(x$fullSibsByCohort)), sibCohortData[, 1], type = "l", ylim = c(0, max(sibCohortData)), xlab = "Cohort", ylab = "Count", lwd = thickness)
      lines(as.numeric(rownames(x$fullSibsByCohort)), sibCohortData[, 2], lty = "dashed", col = col1, lwd = thickness)
      lines(as.numeric(rownames(x$fullSibsByCohort)), sibCohortData[, 4], col = col1, lwd = thickness)
      lines(as.numeric(rownames(x$fullSibsByCohort)), sibCohortData[, 3], lty = "dashed", col = col2, lwd = thickness)
      lines(as.numeric(rownames(x$fullSibsByCohort)), sibCohortData[, 5], col = col2, lwd = thickness)
      mtext("Known sib pairs by cohort.  Black line: full sibs, red: maternal sibs,", side = 1, line = 5)
      mtext("blue: paternal sibs, dashed: total maternal and paternal", side = 1, line = 6)
      mtext("bs, solid coloured: maternal and paternal half sibs.", side = 1, line = 7)
      dev.off()
      readline(prompt = "File saved.  Press <Enter> to continue...")
    }
  }

  if (is.null(cohorts) == FALSE) {
    thickness <- 1
    if (grContrast == TRUE) thickness <- 1.7
    ymax <- max(max(x$maternalGrandmothersByCohort), max(x$maternalGrandfathersByCohort), max(x$paternalGrandmothersByCohort), max(x$paternalGrandfathersByCohort))
    plot(as.numeric(rownames(x$fullSibsByCohort)), x$maternalGrandmothersByCohort, col = col1, type = "l", ylim = c(0, ymax), xlab = "Pedigree depth", ylab = "Count", lwd = thickness)
    lines(as.numeric(rownames(x$fullSibsByCohort)), x$maternalGrandfathersByCohort, lty = "dashed", col = col1, lwd = thickness)
    lines(as.numeric(rownames(x$fullSibsByCohort)), x$paternalGrandmothersByCohort, col = col2, lwd = thickness)
    lines(as.numeric(rownames(x$fullSibsByCohort)), x$paternalGrandfathersByCohort, lty = "dashed", col = col2, lwd = thickness)
    mtext("Known grandparents by cohort (i.e., cohort of grand-offspring).", side = 1, line = 5)
    mtext("Red: maternal grandparents, blue: paternal grand", side = 1, line = 6)
    mtext("parents, solid: grandmothers, dashed: grandfathers.", side = 1, line = 7)
    inp <- readline(prompt = "Press <s> to save current plot or press <Enter> to continue...")
    if (inp == "s") {
      s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
      savePlot(paste(s, ".jpeg", sep = ""), type = "jpeg")
      readline(prompt = "File saved.  Press <Enter> to continue...")
    }

    if (inp == "e") {
      s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
      postscript(paste(s, ".eps", sep = ""), width = 8, height = 8, horizontal = FALSE)
      par(oma = c(5, 1, 1, 1))
      plot(as.numeric(names(x$fullSibsByCohort)), x$maternalGrandmothersByCohort, col = col1, type = "l", ylim = c(0, ymax), xlab = "Pedigree depth", ylab = "Count", lwd = thickness)
      lines(as.numeric(names(x$fullSibsByCohort)), x$maternalGrandfathersByCohort, lty = "dashed", col = col1, lwd = thickness)
      lines(as.numeric(names(x$fullSibsByCohort)), x$paternalGrandmothersByCohort, col = col2, lwd = thickness)
      lines(as.numeric(names(x$fullSibsByCohort)), x$paternalGrandfathersByCohort, lty = "dashed", col = col2, lwd = thickness)
      mtext("Known grandparents by cohort (i.e., cohort of grand-offspring).", side = 1, line = 5)
      mtext("Red: maternal grandparents, blue: paternal grand", side = 1, line = 6)
      mtext("parents, solid: grandmothers, dashed: grandfathers.", side = 1, line = 7)
      dev.off()
      readline(prompt = "File saved.  Press <Enter> to continue...")
    }
  }

  if (lowMem == FALSE) {
    relatednessInterval <- as.numeric(names(x$relatednessCategories)[3]) - as.numeric(names(x$relatednessCategories)[2])
    midBins <- as.numeric(names(x$relatednessCategories))
    binLabels <- paste(midBins - relatednessInterval / 2, "-", midBins + relatednessInterval / 2, sep = "")
    binLabels[1] <- paste("0-", relatednessInterval / 2, sep = "")
    for (i in seq(2, length(binLabels), by = 2)) binLabels[i] <- ""
    plotbins <- NULL
    for (i in 1:length(x$relatednessCategories)) {
      if (is.na(x$relatednessCategories[i]) == FALSE) plotbins <- i
    }
    plotbins <- plotbins + 1
    par(mar = c(10, 4, 4, 2))
    barplot(x$relatednessCategories, ylab = "Count", xaxt = "n")
    axis(1, at = (0:(plotbins - 2)) * 1.2 + 0.6, labels = binLabels[1:(plotbins - 1)], las = 3, xlim = c(0, length(binLabels) - 1))
    mtext("Non-zero pairwise relatednesses       ", 1, 7)
    par(mar = c(5, 4, 4, 2) + 0.1)
    mtext("Distribution of relatedness across the pedigree.", side = 1, line = 5)
    inp <- readline(prompt = "Press <s> to save current plot or press <Enter> to continue...")
    if (inp == "s") {
      s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
      savePlot(paste(s, ".jpeg", sep = ""), type = "jpeg")
      readline(prompt = "File saved.  Press <Enter> to continue...")
    }
  }

  hist(x$maternalSibships[, 2], xlab = "(non-zero) maternal sibship sizes", ylab = "count", main = "")
  inp <- readline(prompt = "Press <s> to save current plot or press <Enter> to continue...")
  if (inp == "s") {
    s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
    savePlot(paste(s, ".jpeg", sep = ""), type = "jpeg")
    readline(prompt = "File saved.  Press <Enter> to continue...")
  }

  hist(x$paternalSibships[, 2], xlab = "(non-zero) paternal sibship sizes", ylab = "count", main = "")
  inp <- readline(prompt = "Press <s> to save current plot or press <Enter> to continue...")
  if (inp == "s") {
    s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
    savePlot(paste(s, ".jpeg", sep = ""), type = "jpeg")
    readline(prompt = "File saved.  Press <Enter> to continue...")
  }

  plot(
    x = 1:(length(x$missingness[1, ]) - 1), y = colMeans(as.matrix(x$missingness[, 2:length(x$missingness[1, ])])), type = "l",
    xlab = "Ancestral generation", ylab = "Mean indiviual pedigree completeness"
  )
  mtext("Pedigree-wide average ancestral pedigree completness", side = 1, line = 5)
  mtext("following MacCluer et al. (1983) J Hered 74:394-399.", side = 1, line = 6)
  inp <- readline(prompt = "Press <s> to save current plot or press <Enter> to continue...")
  if (inp == "s") {
    s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
    savePlot(paste(s, ".jpeg", sep = ""), type = "jpeg")
    readline(prompt = "File saved.  Press <Enter> to continue...")
  }

  if (is.null(cohorts) == FALSE) {
    gc <- gray.colors(length(x$cohortPedigreeCompleteness[, 1]) - 1)
    plot(
      x = 1:(length(x$missingness[1, ]) - 1), y = x$cohortPedigreeCompleteness[length(x$cohortPedigreeCompleteness[, 1]), ],
      type = "l", xlab = "Ancestral generation", ylab = "Mean individual pedigree completeness by cohort",
      ylim = c(0, max(x$cohortPedigreeCompleteness))
    )
    for (i in 1:(length(x$cohortPedigreeCompleteness[, 1]) - 1)) {
      lines(
        x = 1:(length(x$missingness[1, ]) - 1), y = x$cohortPedigreeCompleteness[i, ],
        col = gc[length(x$cohortPedigreeCompleteness[, 1]) - i]
      )
    }
    mtext("Average ancestral pedigree completeness by cohort", side = 1, line = 5)
    mtext("following MacCluer et al. (1983) J Hered 74:394-399.", side = 1, line = 6)
    mtext("Most recent cohorts are represented by the darkest lines.", side = 1, line = 7)
    inp <- readline(prompt = "Press <s> to save current plot or press <Enter> to continue...")
    if (inp == "s") {
      s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
      savePlot(paste(s, ".jpeg", sep = ""), type = "jpeg")
      readline(prompt = "File saved.  Press <Enter> to continue...")
    }

    if (inp == "e") {
      s <- readline(prompt = "Enter path (including file name but not extension) to which to save image: ")
      postscript(paste(s, ".eps", sep = ""), width = 8, height = 8, horizontal = FALSE)
      par(oma = c(5, 1, 1, 1))
      plot(
        x = 1:(length(x$missingness[1, ]) - 1), y = x$cohortPedigreeCompleteness[length(x$cohortPedigreeCompleteness[, 1]), ],
        type = "l", xlab = "Ancestral generation", ylab = "Mean individual pedigree completeness by cohort",
        ylim = c(0, max(x$cohortPedigreeCompleteness))
      )
      for (i in 1:(length(x$cohortPedigreeCompleteness[, 1]) - 1)) {
        lines(
          x = 1:(length(x$missingness[1, ]) - 1), y = x$cohortPedigreeCompleteness[i, ],
          col = gc[length(x$cohortPedigreeCompleteness[, 1]) - i]
        )
      }
      mtext("Average ancestral pedigree completeness by cohort", side = 1, line = 5)
      mtext("following MacCluer et al. (1983) J Hered 74:394-399.", side = 1, line = 6)
      mtext("Most recent cohorts are represented by the darkest lines.", side = 1, line = 7)
      dev.off()
      readline(prompt = "File saved.  Press <Enter> to continue...")
    }
  }
}