Merge pull request #12 from mrc-ide/feature/tutorials

Feature/tutorials

NAMESPACE

@@ -30,12 +30,13 @@ export(mipanalyzer_file)
 export(pca_wsaf)
 export(pcoa_genomic_distance)
 export(plot_coverage)
+export(plot_distance)
 export(plot_map)
 export(plot_pca)
 export(plot_pca_contribution)
 export(plot_pca_variance)
 export(plot_pcoa)
-export(print_full)
+export(plot_wsaf)
 export(rbetabinom)
 export(rdirichlet)
 export(sim_biallelic)

R/analysis_distance.R

@@ -231,12 +231,16 @@ inbreeding_mle <- function(x, f = seq(0,1,l=11), ignore_het = FALSE, report_prog
 #------------------------------------------------
 #' @title Get identity by mixture
 #'
-#' @description Get identity by mixture distance.
+#' @description Get identity by "mixture distance". The mixture distance between
+#'   two samples is the proportion of loci that have identical within-sample
+#'   allele frequencies (WSAFs), or alternatively have WSAFs within a given
+#'   tolerance. This extends the idea of identity by state (IBS) to continuous
+#'   WSAFs rather than categorical genotypes.
 #'
 #' @param x object of class \code{mipanalyzer_biallelic}.
-#' @param tol tolerance on mixture comparisons. Defeault = 0
+#' @param tol tolerance on mixture comparisons. Default = 0.
 #' @param diagonal Should the diagonal of the distance matrix be changed to a
-#'   given value. Deafult = NULL, which cause no changes.
+#'   given value. Default = NULL, which cause no changes.
 #' @param report_progress if \code{TRUE} then a progress bar is printed to the
 #'   console.
 #'
@@ -247,6 +251,9 @@ get_IB_mixture <- function(x, tol = 0, diagonal = NULL, report_progress = TRUE)
   # check inputs
   assert_custom_class(x, "mipanalyzer_biallelic")
   assert_single_numeric(tol)
+  if (!is.null(diagonal)) {
+    assert_single_numeric(diagonal)
+  }
   assert_single_logical(report_progress)
 
   # get basic quantities
@@ -277,6 +284,7 @@ get_IB_mixture <- function(x, tol = 0, diagonal = NULL, report_progress = TRUE)
   if (!is.null(diagonal)) {
     ret[diag(ret)] <- diagonal
   }
+  ret[lower.tri(ret)] <- NA
 
   # return distance matrix
   return(ret)
@@ -292,10 +300,10 @@ get_IB_mixture <- function(x, tol = 0, diagonal = NULL, report_progress = TRUE)
 #'   at every locus.
 #'
 #' @param x object of class \code{mipanalyzer_biallelic}.
-#' @param ignore_het whether to ignore heterzygous comparisons, or alternatively
+#' @param ignore_het whether to ignore heterozygous comparisons, or alternatively
 #'   call the major allele at every locus (see details).
 #' @param diagonal Should the diagonal of the distance matrix be changed to a
-#'   given value. Deafult = NULL, which cause no changes.
+#'   given value. Default = \code{NULL}, which cause no changes.
 #' @param report_progress if \code{TRUE} then a progress bar is printed to the
 #'   console.
 #'
@@ -340,7 +348,44 @@ get_IBS_distance <- function(x, ignore_het = TRUE, diagonal = NULL, report_progr
   if (!is.null(diagonal)) {
     ret[diag(ret)] <- diagonal
   }
+  ret[lower.tri(ret)] <- NA
+
 
   # return distance matrix
   return(ret)
 }
+
+#------------------------------------------------
+#' @title Plot a distance matrix
+#'
+#' @description Simple image plot of a matrix of pairwise distances.
+#'
+#' @param m square matrix of pairwise distances.
+#' @param col_pal which viridis colour pallet to use. Options are "viridis",
+#'   "plasma", "magma" or "inferno".
+#'
+#' @export
+
+plot_distance <- function(m, col_pal = "plasma") {
+
+  # avoid "no visible binding" notes
+  x <- y <- value <- NULL
+
+  # check inputs
+  assert_square_matrix(m)
+  assert_single_string(col_pal)
+  assert_in(col_pal, c("viridis", "plasma", "magma", "inferno"))
+
+  # get matrix into long-form data.frame
+  n <- nrow(m)
+  df_plot <- data.frame(x = 1:n,
+                        y = rep(1:n, each = n),
+                        value = as.vector(m))
+
+  # produce plot
+  df_plot |>
+    ggplot() + theme_void() +
+    geom_tile(aes(x = x, y = y, fill = value)) +
+    scale_fill_viridis_c(option = col_pal, na.value = "white")
+
+}

R/main.R

@@ -97,6 +97,48 @@ get_wsaf <- function(x, impute = TRUE, FUN = median, ...) {
   return(wsaf)
 }
 
+#------------------------------------------------
+#' @title Plot within-sample allele frequencies
+#'
+#' @description Simple image plot of within-sample allele frequencies. The top
+#'   row of the plot corresponds to the first row of the input matrix.
+#'
+#' @param x matrix of within-sample allele frequencies, with samples in rows and
+#'   loci in columns.
+#' @param col_pal which viridis colour pallet to use. Options are "viridis",
+#'   "plasma", "magma" or "inferno".
+#' 
+#' @importFrom methods is
+#' @export
+
+plot_wsaf <- function(x, col_pal = "plasma") {
+
+  # avoid "no visible binding" notes
+  locus <- wsaf <- NULL
+
+  # check inputs
+  assert_matrix(x)
+  assert_bounded(x)
+  assert_single_string(col_pal)
+  assert_in(col_pal, c("viridis", "plasma", "magma", "inferno"))
+
+  # get matrix into long-form data.frame
+  df_plot <- data.frame(sample = nrow(x):1,
+                        locus = rep(1:ncol(x), each = nrow(x)),
+                        wsaf = as.vector(x))
+
+  # produce plot
+  df_plot |>
+    ggplot() + theme_bw() +
+    geom_tile(aes(x = locus, y = sample, fill = wsaf)) +
+    scale_fill_viridis_c(option = col_pal) +
+    scale_x_continuous(expand = c(0, 0)) +
+    scale_y_continuous(expand = c(0, 0)) +
+    xlab("Locus") + ylab("Sample") +
+    theme(axis.ticks = element_blank(),
+          axis.text = element_blank())
+}
+
 #------------------------------------------------
 #' @title Produce ggplot map
 #'
@@ -171,25 +213,3 @@ Pf_chrom_lengths <- function() {
                                2895605, 3291871))
   return(ret)
 }
-
-#------------------------------------------------
-#' @title Ordinary print function for unclassed object
-#'
-#' @description Calling \code{print()} on an object of custom class, e.g.
-#'   \code{mipanalyzer_biallelic}, results in custom output. This function
-#'   therefore stands in for the base \code{print()} function, and is equivalent
-#'   to running \code{print(unclass(x))}.
-#'
-#' @param x object of custom class (e.g. \code{mipanalyzer_biallelic}).
-#' @param ... other arguments passed to \code{print()}.
-#'
-#' @export
-
-print_full <- function(x, ...) {
-
-  # print un-classed object
-  print(unclass(x), ...)
-
-  # return invisibly
-  invisible(x)
-}

_pkgdown.yml

@@ -12,6 +12,10 @@ navbar:
       menu:
         - text: Filtering data
           href: articles/filtering_data.html
+        - text: Structure analysis
+          href: articles/structure_analysis.html
+        - text: Comparing genetic distances
+          href: articles/genetic_distances.html
     - text: Functions
       href: reference/index.html