compute_GRM: Compute genomic relationship matrices.

In DominikMueller64/dmisc: Dominik Mueller's R Code

Description

Compute genomic relationship matrices according to different methods.

Usage

compute_GRM(..., scaling = c("VR1", "VR2"), offblock = c("canonical",
  "chen", "melchinger"), p = NULL, weighted = NULL, check = TRUE)

Arguments

...	Genotype matrices of the differen populations or families, coded with (0, 1, 2).
scaling	The scaling procedure that should be applied to the columns, one of 'VR1' or 'VR2'. With 'VR1', all columns are scaled by the same constand computed from the allele frequencies and with 'VR2', columns are scaled individually.
offblock	The scope of the applied allele frequencies, one of 'canonical', 'melchinger' or 'chen'. If 'canonical', ovarall allele frequencies are either computed or taken from argument p and used for centering and scaling in all populations or families. If 'melchinger' or 'chen', population-specific allele frequencies are computed and used for centering and scaling.
p	Allele frequencies, either a vector if offblock is canonical or a list of vectors if offblock is 'melchinger' or 'chen'.
weighted	A logical value. Should the allele frequencies be weighted? Only regarded if offblock is canonical and no allele frequencies (p) were given.
check	Should checks of the input be performed? Not doing so will be faster.

Author(s)

Dominik Mueller

Examples

 Xa <- structure(c(2, 0, 0, 2, 0, 2, 0, 2, 2, 2, 0, 0, 2, 0, 2, 0, 2,
                  2, 0, 2, 2, 0, 2, 2, 0, 2, 2, 0, 2, 2, 2, 2, 0, 0, 0, 0),
                 .Dim = c(4L, 9L),
                  .Dimnames = list(c("PopA_Ind_1", "PopA_Ind_2", "PopA_Ind_3",
                                    "PopA_Ind_4"),
                                  c("Locus_1", "Locus_2", "Locus_3", "Locus_4",
                                    "Locus_5", "Locus_6", "Locus_7", "Locus_8", "Locus_9")))

 Xb <- structure(c(2, 0, 2, 2, 0, 0, 2, 2, 2, 2, 0, 2, 0, 2, 2, 0, 0,
                  2, 2, 2, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 0, 2, 0, 2, 0, 0),
                .Dim = c(4L, 9L),
                .Dimnames = list(c("PopB_Ind_1", "PopB_Ind_2", "PopB_Ind_3", "PopB_Ind_4"),
                                 c("Locus_1", "Locus_2", "Locus_3", "Locus_4",
                                   "Locus_5", "Locus_6", "Locus_7", "Locus_8", "Locus_9")))

# We assume that the there are bi-parental families from completely inbreed parents AND
# that all loci monomorphic in the progeny are also monomophic in the parents (which is by no
 
hf <- function(x) {
 # Little helper for getting the expected frequencies according to our assumptions.
 x[x > 0.0 & x < 1.0] <- 0.5
 x
} 
p_spec <- list(hf(colMeans(Xa) / 2.0), hf(colMeans(Xb) / 2.0))
p <- Reduce(`+`, p_spec) / 2.0 # Weighting or not makes no difference here, as the population size are equal.
 

dots <- list(Xa, Xb)
scaling_levels <- c('VR1', 'VR2')
offblock_levels <- c('canonical', 'melchinger', 'chen')

library('tibble')
library('purrr')
library('dplyr')
library('forcats')
library('ggplot2')

d <- cross_d(list(scaling = scaling_levels, offblock = offblock_levels))

dat <- map2_df(d[[1L]], d[[2L]], function(scaling, offblock) {
  
  if (offblock %in% c('melchinger', 'chen')) {
    tmp <- p_spec
  } else
    tmp <- p
  # tmp <- NULL
  G <- purrr::invoke(compute_GRM, dots, scaling = scaling, offblock = offblock,
                     p = tmp, weighted = TRUE)
  if (!is.null(G))
    mat2df(G) %>% as_tibble() %>% mutate(scaling = scaling, offblock = offblock)
})
dat$row_names <- forcats::fct_rev(dat$row_names)

ggplot(data = dat,
       mapping = aes(x = col_names, y = row_names)) +
  geom_tile(aes(fill = value), colour = "white") +
  facet_grid(facets = scaling ~ offblock) +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       midpoint = 0, limit = c(min(dat$value), max(dat$value)), space = "Lab", 
                       name = "Coefficient") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
 geom_text(aes(col_names, row_names, label = signif(value, 3L)), color = "darkgreen", size = 4)

DominikMueller64/dmisc documentation built on May 6, 2019, 2:52 p.m.