Abstract

We provide two examples of parallelizing flexdog to genotype many genes.

Analysis

Before we begin, we’ll need to load the parallel, foreach, and doParallel packages. There are other ways to perform parallel computation in R, but we find these to be pretty handy and easy to use.

library(parallel)
library(foreach)
library(doParallel)

Of course, we’ll need to load the updog package to actually do the genotyping.

library(updog)

We’ll work with the uitedewilligen dataset, which is a subset of the data from Uitdewilligen (2013).

data("uitdewilligen")
ngenes <- ncol(uitdewilligen$refmat)

You can read more about these data by typing ?uitdewilligen in R.

Suppose we wanted to use parallel computing to genotype this individual. Also suppose we only wanted the MAP (maximum a posteriori) estimated genotypes for each individual and each gene. Then we can use the following code.

nc <- 2 ## number of cores. 
        ## You should change this for your specific computing environment.
cl <- parallel::makeCluster(nc)
doParallel::registerDoParallel(cl = cl)
stopifnot(foreach::getDoParWorkers() > 1) ## make sure cluster is set up.
gene_est <- foreach(i = 1:ngenes, 
                    .combine = cbind, 
                    .export = "flexdog") %dopar% {
                      ## fit flexdog
                      fout <- flexdog(refvec  = uitdewilligen$refmat[, i], 
                                      sizevec = uitdewilligen$sizemat[, i], 
                                      ploidy  = uitdewilligen$ploidy, 
                                      model   = "norm", 
                                      verbose = FALSE)
                      fout$geno
                      }
stopCluster(cl)

Note that I used model = "norm" because we have a very small number of samples.

However, this doesn’t return any sort of diagnostic statistics for each SNP. So suppose we now want to save all of the output of each SNP and not just the MAP genotype, then we can use the following code:

nc <- 2 ## number of cores. 
        ## You should change this for your specific computing environment.
cl <- parallel::makeCluster(nc)
doParallel::registerDoParallel(cl = cl)
stopifnot(foreach::getDoParWorkers() > 1) ## make sure cluster is set up.
glist <- foreach(i = 1:ngenes, 
                 .export = "flexdog") %dopar% {
                   ## fit flexdog
                   fout <- flexdog(refvec  = uitdewilligen$refmat[, i], 
                                   sizevec = uitdewilligen$sizemat[, i], 
                                   ploidy  = uitdewilligen$ploidy, 
                                   model   = "norm", 
                                   verbose = FALSE)
                   fout
                   }
stopCluster(cl)

Returning the diagnostics for each SNP is pretty easy. For example, below we extract the variable prop_mis from the list glist and make a histogram of these values.

library(tidyverse)
prop_mis <- sapply(X = glist, FUN = function(x) x$prop_mis)
qplot(prop_mis, 
      xlab = "Posterior Proportion Mis-genotyped", 
      bins = 30) +
  theme_bw()