implement pairwise_grm

NAMESPACE

@@ -126,6 +126,7 @@ export(loci_names)
 export(loci_transitions)
 export(loci_transversions)
 export(pairwise_allele_sharing)
+export(pairwise_grm)
 export(pairwise_ibs)
 export(pairwise_king)
 export(pairwise_pop_fst)

R/pairwise_allele_sharing.R

@@ -1,7 +1,7 @@
 #' Compute the Pairwise Allele Sharing Matrix for a `gen_tibble` object
 #'
 #' This function computes the Allele Sharing matrix.
-#' Estimates Allele Sharing (matching in `hierfstat`)) between pairs of individuals
+#' Estimates Allele Sharing ([hierfstat::matching()])) between pairs of individuals
 #' (for each locus, gives 1 if the two individuals are homozygous
 #' for the same allele, 0 if they are homozygous for a different allele, and 1/2 if at least one individual
 #' is heterozygous. Matching is the average of these 0, 1/2 and 1s)

R/pairwise_grm.R

@@ -0,0 +1,31 @@
+#' Compute the Genomic Relationship Matrix for a `gen_tibble` object
+#'
+#' This function computes the Genomic Relationship Matrix (GRM). This is estimated
+#' by computing the pairwise kinship coefficients (coancestries) between all pairs of individuals
+#' from a matrix of Allele Sharing follwng the approach of Weir and Goudet 2017 based
+#' on beta estimators).
+#'
+#' The GRM is twice the coancestry matrix (e.g. as estimated by
+#'  [hierfstat::beta.dosage()] with `inb=FALSE`).
+#'
+#' @param x a `gen_tibble` object.
+#' @param allele_sharing_mat optional, the matrix of Allele Sharing returned by
+#' [pairwise_allele_sharing()] with `as_matrix=TRUE`. As a number of statistics can be
+#' derived from the Allele Sharing matrix,
+#' it it sometimes more efficient to pre-compute this matrix.
+#' @param block_size the size of the blocks to use for the computation of the allele sharing matrix.
+#' @returns a matrix of GR between all pairs of individuals
+#' @export
+pairwise_grm <- function(x, allele_sharing_mat = NULL,
+                                    block_size = bigstatsr::block_size(count_loci(x))) {
+  if (is.null(allele_sharing_mat)){
+    allele_sharing_mat <- pairwise_allele_sharing(x, as_matrix = TRUE, block_size = block_size)
+  }
+  Mii <- diag(allele_sharing_mat)
+  diag(allele_sharing_mat) <- NA
+  MB <- mean(allele_sharing_mat, na.rm = T)
+  diag(allele_sharing_mat) <- Mii
+  # estimate the pairwise kinship (coancestries)
+  coa <- (allele_sharing_mat - MB)/(1 - MB)
+  return(coa*2)
+}

tests/testthat/test_pairwise_grm.R

@@ -0,0 +1,26 @@
+test_genotypes <- rbind(c(1,1,0,1,1,0),
+                        c(2,1,0,NA,0,0),
+                        c(2,NA,0,0,1,1),
+                        c(1,0,0,1,0,0),
+                        c(1,2,0,1,2,1),
+                        c(0,0,0,0,NA,1),
+                        c(0,1,1,0,1,NA))
+test_indiv_meta <- data.frame (id=c("a","b","c","d","e","f","g"),
+                               population = c("pop1","pop1","pop2","pop2","pop1","pop3","pop3"))
+test_loci <- data.frame(name=paste0("rs",1:6),
+                        chromosome=paste0("chr",c(1,1,1,1,2,2)),
+                        position=as.integer(c(3,5,65,343,23,456)),
+                        genetic_dist = as.double(rep(0,6)),
+                        allele_ref = c("A","T","C","G","C","T"),
+                        allele_alt = c("T","C", NA,"C","G","A"))
+
+test_gt <- gen_tibble(x = test_genotypes, loci = test_loci, indiv_meta = test_indiv_meta, quiet = TRUE)
+
+
+test_that("pairwise_grm compute correctly",{
+  # expect error if the tibble is not grouped
+  coa_hier <- hierfstat::beta.dosage(test_genotypes, inb = FALSE)
+  grm_hier <- 2*coa_hier
+  grm_gt <- pairwise_grm(test_gt)
+  expect_true(all.equal(grm_hier, grm_gt, check.attributes = FALSE))
+})