mgdrive_run.R
In MGDrivE: Mosquito Gene Drive Explorer

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  echo = TRUE,
  hold = TRUE,
  fig.width = 7,
  fig.height = 10.5,
  eval = TRUE
)

# set seed for reproduibility
set.seed(seed = 42)

## ----eval=FALSE---------------------------------------------------------------
# ####################
# # Load libraries
# ####################
# library(MGDrivE)
# 
# ####################
# # Output Folder
# ####################
# outFolder <- "mgdrive"
# dir.create(path = outFolder)
# 
# ####################
# # Simulation Parameters
# ####################
# # days to run the simulation
# tMax <- 365
# 
# # number of Monte Carlo iterations
# nRep <- 5
# 
# # each Monte Carlo iteration gets its own folder
# folderNames <- file.path(outFolder,
#                         formatC(x = 1:nRep, width = 3, format = "d", flag = "0"))
# 
# # entomological parameters
# bioParameters <- list(betaK=20,tEgg=5,tLarva=6,tPupa=4,popGrowth=1.175,muAd=0.09)
# 
# # a 5-node network with 5% per day migration rate
# sitesNumber <- 5
# adultPopEquilibrium <- 500
# 
# # auxiliary function
# triDiag <- function(upper, lower){
# 
#   # return matrix
#   retMat <- matrix(data = 0, nrow = length(upper) + 1, ncol = length(upper) + 1)
# 
#   # set index values for upper/lower triangles
#   indx <- 1:length(upper)
# 
#   # set forward/backward migration using matrix access
#   retMat[cbind(indx+1,indx)] <- lower
#   retMat[cbind(indx,indx+1)] <- upper
# 
#   # set stay probs
#   diag(x = retMat) <- 1-rowSums(x = retMat)
# 
#   return(retMat)
# }
# 
# # fill movement matrix
# #  Remember, rows need to sum to 1.
# moveMat <- triDiag(upper = rep.int(x = 0.05, times = sitesNumber-1),
#                    lower = rep.int(x = 0.05, times = sitesNumber-1))
# 
# # batch migration is disabled by setting the probability to 0
# batchMigration <- basicBatchMigration(batchProbs=0,
#                                      sexProbs=c(.5,.5),
#                                      numPatches=sitesNumber)
# 
# ####################
# # Basic Inheritance pattern
# ####################
# # Mendelian cube with standard (default) paraameters
# cube <- cubeMendelian()
# 
# 
# ####################
# # Setup releases and batch migration
# ####################
# # set up the empty release vector
# #  MGDrivE pulls things out by name
# patchReleases <- replicate(n=sitesNumber,
#                            expr={list(maleReleases=NULL,femaleReleases=NULL,
#                                       eggReleases=NULL,matedFemaleReleases=NULL)},
#                            simplify=FALSE)
# 
# # choose release parameters
# #  Releases start at time 25, occur once a week, for 10 weeks.
# #  There are 100 mosquitoes released every time.
# releasesParameters <- list(releasesStart=25,
#                           releasesNumber=10,
#                           releasesInterval=7,
#                           releaseProportion=100)
# 
# # generate male release vector
# maleReleasesVector <- generateReleaseVector(driveCube=cube,
#                                            releasesParameters=releasesParameters)
# 
# # generate female release vector
# femaleReleasesVector <- generateReleaseVector(driveCube=cube,
#                                              releasesParameters=releasesParameters)
# 
# # put releases into the proper place in the release list
# #  This performs the releases in the first patch only
# patchReleases[[1]]$maleReleases <- maleReleasesVector
# patchReleases[[1]]$femaleReleases <- femaleReleasesVector
# 
# 
# ####################
# # Combine parameters and run!
# ####################
# # setup parameters for the network. This builds a list of parameters required for
# #  every population in the network.
# netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, nPatch=sitesNumber,
#                              beta=bioParameters$betaK, muAd=bioParameters$muAd,
#                              popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg,
#                              tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa,
#                              AdPopEQ=adultPopEquilibrium, inheritanceCube = cube)
# 
# # set MGDrivE to run stochastic
# setupMGDrivE(stochasticityON = TRUE, verbose = FALSE)
# 
# # build network prior to run
# MGDrivESim <- Network$new(params=netPar,
#                          driveCube=cube,
#                          patchReleases=patchReleases,
#                          migrationMale=moveMat,
#                          migrationFemale=moveMat,
#                          migrationBatch=batchMigration,
#                          directory=folderNames,
#                          verbose = FALSE)
# # run simulation
# MGDrivESim$multRun(verbose = FALSE)
# 
# 
# ####################
# # Post Analysis
# ####################
# # First, split output by patch
# # Second, aggregate females by their mate choice
# for(i in 1:nRep){
#  splitOutput(readDir = folderNames[i], remFile = TRUE, verbose = FALSE)
#  aggregateFemales(readDir = folderNames[i], genotypes = cube$genotypesID,
#                   remFile = TRUE, verbose = FALSE)
# }
# 
# # plot output of first run to see effect
# plotMGDrivESingle(readDir=folderNames[1],totalPop = TRUE,lwd=3.5,alpha=1)
# 
# # plot all 5 repetitions together
# plotMGDrivEMult(readDir=outFolder,lwd=0.35,alpha=0.75)

## -----------------------------------------------------------------------------
####################
# Load libraries
####################
library(MGDrivE)

## -----------------------------------------------------------------------------
####################
# Output Folder
####################
outFolder <- "mgdrive"
dir.create(path = outFolder)

## -----------------------------------------------------------------------------
####################
# Simulation Parameters
####################
# days to run the simulation
tMax <- 365

# number of Monte Carlo iterations
nRep <- 5

# each Monte Carlo iteration gets its own folder
folderNames <- file.path(outFolder,
                        formatC(x = 1:nRep, width = 3, format = "d", flag = "0"))

folderNames

## -----------------------------------------------------------------------------
# entomological parameters
bioParameters <- list(betaK=20,tEgg=5,tLarva=6,tPupa=4,popGrowth=1.175,muAd=0.09)

## -----------------------------------------------------------------------------
# a 5-node network with 5% per day migration rate
sitesNumber <- 5
adultPopEquilibrium <- 500
patchPops <- rep(adultPopEquilibrium,sitesNumber) 
# patchPops is optional. If all populations are the same size, parameterizeMGDrivE
#  can take a single number. However, if you desire different population sizes, it 
#  must be a vector of length equal to the number of sites

# landscape
# auxiliary function
triDiag <- function(upper, lower){
  
  # return matrix
  retMat <- matrix(data = 0, nrow = length(upper) + 1, ncol = length(upper) + 1)
  
  # set index values for upper/lower triangles
  indx <- 1:length(upper)
  
  # set forward/backward migration using matrix access
  retMat[cbind(indx+1,indx)] <- lower
  retMat[cbind(indx,indx+1)] <- upper
  
  # set stay probs
  diag(x = retMat) <- 1-rowSums(x = retMat)
  
  return(retMat)
}

# fill movement matrix
#  Remember, rows need to sum to 1.
moveMat <- triDiag(upper = rep.int(x = 0.05, times = sitesNumber-1),
                   lower = rep.int(x = 0.05, times = sitesNumber-1))

moveMat

## -----------------------------------------------------------------------------
# batch migration is disabled by setting the probability to 0
batchMigration <- basicBatchMigration(batchProbs=0, numPatches=sitesNumber)

## -----------------------------------------------------------------------------
####################
# Basic Inheritance pattern
####################
# Mendelian cube with standard (default) paraameters
cube <- cubeMendelian()

## -----------------------------------------------------------------------------
####################
# Setup releases and batch migration
####################
# set up the empty release vector
#  MGDrivE pulls things out by name
patchReleases <- replicate(n=sitesNumber,
                           expr={list(maleReleases=NULL,femaleReleases=NULL,
                                      eggReleases=NULL,matedFemaleReleases=NULL)},
                           simplify=FALSE)

# choose release parameters
#  Releases start at time 25, occur once a week, for 10 weeks.
#  There are 100 mosquitoes released every time.
releasesParameters <- list(releasesStart=25,
                          releasesNumber=10,
                          releasesInterval=7,
                          releaseProportion=100)

## -----------------------------------------------------------------------------
# generate male release vector
maleReleasesVector <- generateReleaseVector(driveCube=cube, 
                                           releasesParameters=releasesParameters)
maleReleasesVector[[1]]

## -----------------------------------------------------------------------------
# generate female release vector
femaleReleasesVector <- generateReleaseVector(driveCube=cube,
                                             releasesParameters=releasesParameters)

# put releases into the proper place in the release list
patchReleases[[1]]$maleReleases <- maleReleasesVector
patchReleases[[1]]$femaleReleases <- femaleReleasesVector


## -----------------------------------------------------------------------------
####################
# Combine parameters and run!
####################
# setup parameters for the network. This builds a list of parameters required for
#  every population in the network.
netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, nPatch=sitesNumber,
                             beta=bioParameters$betaK, muAd=bioParameters$muAd,
                             popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg,
                             tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa,
                             AdPopEQ=patchPops, inheritanceCube = cube)

## -----------------------------------------------------------------------------
# set MGDrivE to run stochastic
setupMGDrivE(stochasticityON = TRUE)

## -----------------------------------------------------------------------------
# build network prior to run
MGDrivESim <- Network$new(params=netPar,
                         driveCube=cube,
                         patchReleases=patchReleases,
                         migrationMale=moveMat,
                         migrationFemale=moveMat,
                         migrationBatch=batchMigration,
                         directory=folderNames,
                         verbose = TRUE)

# list folders to show that they have been created
list.files(path = outFolder)

## -----------------------------------------------------------------------------
# run simulation
MGDrivESim$multRun(verbose = FALSE)

## -----------------------------------------------------------------------------
# first and last repetitions
list.files(path = outFolder)[1]
list.files(path = list.files(path = outFolder, full.names = TRUE)[1], recursive = TRUE)
list.files(path = outFolder)[5]
list.files(path = list.files(path = outFolder, full.names = TRUE)[5], recursive = TRUE)

## -----------------------------------------------------------------------------
# read in male and female files
mMat <- as.matrix(read.csv(file = list.files(path = outFolder, full.names = TRUE,
                                             recursive = TRUE)[2],
                           header = TRUE, sep = ","))
fMat <- as.matrix(read.csv(file = list.files(path = outFolder, full.names = TRUE,
                                             recursive = TRUE)[1],
                           header = TRUE, sep = ","))

# look at male file header
colnames(mMat)

## -----------------------------------------------------------------------------
head(x = mMat, n = 2*sitesNumber)

## -----------------------------------------------------------------------------
# look at female file header
colnames(fMat)

## -----------------------------------------------------------------------------
head(x = fMat, n = 2*sitesNumber)

## -----------------------------------------------------------------------------
# First, split output by patch
for(i in 1:nRep){
 splitOutput(readDir = folderNames[i], remFile = TRUE, verbose = FALSE)
}

## -----------------------------------------------------------------------------
# first and last repetitions
list.files(path = outFolder)[1]
list.files(path = list.files(path = outFolder, full.names = TRUE)[1], recursive = TRUE)
list.files(path = outFolder)[5]
list.files(path = list.files(path = outFolder, full.names = TRUE)[5], recursive = TRUE)

# read in examples of new male/female files
twoFiles <- list.files(path = outFolder, full.names = TRUE,
                       recursive = TRUE)[c(1,6)]

# read in male and female files
mMat <- as.matrix(read.csv(file = twoFiles[2], header = TRUE, sep = ","))
fMat <- as.matrix(read.csv(file = twoFiles[1], header = TRUE, sep = ","))


## -----------------------------------------------------------------------------
head(x = mMat, n = 5)

## -----------------------------------------------------------------------------
head(x = fMat, n = 5)

## -----------------------------------------------------------------------------
# Second, aggregate females by their mate choice
for(i in 1:nRep){
 aggregateFemales(readDir = folderNames[i], genotypes = cube$genotypesID,
                  remFile = TRUE, verbose = FALSE)
}

## -----------------------------------------------------------------------------
# changed female files in first repetition
list.files(path = outFolder)[1]
list.files(path = list.files(path = outFolder, full.names = TRUE)[1], recursive = TRUE)

# read in examples of new female file
fMat2 <- as.matrix(read.csv(file = list.files(path = outFolder,
                                              recursive = TRUE,
                                              full.names = TRUE)[1],
                            header = TRUE, sep = ","))

## -----------------------------------------------------------------------------
head(x = fMat2, n = 5)

## -----------------------------------------------------------------------------
# show non-aggregated female file split by patch
#  This is for patch number 1
fMat[(releasesParameters$releasesStart-2):(releasesParameters$releasesStart+2), ]
cat("\n")

# show aggregated female file
#  This is for patch number 1
fMat2[(releasesParameters$releasesStart-2):(releasesParameters$releasesStart+2), ]

## ----fig.alt="Example of how to plot a single repetition from a set of simulations."----
# plot the first repetition
plotMGDrivESingle(readDir=folderNames[1],totalPop = TRUE,lwd=3.5,alpha=1)

## ----fig.alt="Example of how to plot all repetitions from a set of simulations."----
# plot all 5 repetitions together
plotMGDrivEMult(readDir=outFolder,lwd=0.35,alpha=0.75)

## ----echo=FALSE---------------------------------------------------------------
####################
# Cleanup before next run
####################
unlink(x = outFolder, recursive = TRUE)
rm(list=ls())