## ----echo=FALSE--------------------------------------------------------------- if (capabilities("cairo")) knitr::opts_chunk$set(dev.args = list(png = list(type = "cairo"))) ## ----message=FALSE, fig.width=5, fig.asp=1, fig.align='center', fig.cap="An epistatic network generated between 50 QTL."---- library(epinetr) # Build a population of size 1000, with 50 QTL, broad-sense heritability of 0.4, # narrow-sense heritability of 0.2 and overall trait variance of 40. pop <- Population(popSize = 1000, map = map100snp, alleleFrequencies = runif(100), QTL = 50, broadH2 = 0.4, narrowh2 = 0.2, traitVar = 40) # Attach additive effects pop <- addEffects(pop) # Attach an epistatic network pop <- attachEpiNet(pop) # Plot the network plot(getEpiNet(pop)) ## ----------------------------------------------------------------------------- # Inspect initial phenotypic components head(getComponents(pop)) ## ----message=FALSE, fig.width=7, fig.height=4, fig.align='center', fig.cap="A graphical representation of a simulation run across 250 generations."---- # Run a simulation across 250 generations pop <- runSim(pop, generations = 250, truncSire = 0.1, truncDam = 0.5) # Plot the simulation run plot(pop) ## ----------------------------------------------------------------------------- # Get the allele frequencies af <- getAlleleFreqRun(pop) # Get the phased genotypes of the resulting population geno <- getPhased(pop) # Get a subset of the resulting population ID <- getComponents(pop)$ID ID <- sample(ID, 50) pop2 <- getSubPop(pop, ID) ## ----------------------------------------------------------------------------- head(map100snp) ## ----------------------------------------------------------------------------- nrow(map100snp) ## ----------------------------------------------------------------------------- length(unique(map100snp[, 2])) ## ----message=FALSE------------------------------------------------------------ pop <- Population(popSize = 200, map = map100snp, QTL = 20, alleleFrequencies = runif(100), broadH2 = 0.9, narrowh2 = 0.6, traitVar = 40) pop ## ----message=FALSE------------------------------------------------------------ pop <- Population(popSize = 200, map = map100snp, QTL = 20, alleleFrequencies = runif(100)) pop ## ----message=FALSE------------------------------------------------------------ pop <- Population(popSize = 200, map = map100snp, QTL = 20, alleleFrequencies = runif(100), h2est = 0.6) ## ----message=FALSE------------------------------------------------------------ pop <- Population(popSize = 200, map = map100snp, QTL = c(62, 55, 92, 74, 11, 38), alleleFrequencies = runif(100), broadH2 = 0.9, narrowh2 = 0.6, traitVar = 40) pop ## ----------------------------------------------------------------------------- getQTL(pop) ## ----------------------------------------------------------------------------- dim(geno100snp) ## ----------------------------------------------------------------------------- geno100snp[1, 1:10] ## ----message=FALSE------------------------------------------------------------ pop <- Population(popSize = nrow(geno100snp), map = map100snp, QTL = 20, genotypes = geno100snp, broadH2 = 0.9, narrowh2 = 0.6, traitVar = 40) ## ----message=FALSE------------------------------------------------------------ pop <- Population(popSize = nrow(geno100snp), map = map100snp, QTL = 20, genotypes = geno100snp, literal = FALSE, broadH2 = 0.9, narrowh2 = 0.6, traitVar = 40) ## ----message=FALSE------------------------------------------------------------ pop <- Population(pop, broadH2 = 0.7, traitVar = 30) pop ## ----message=FALSE------------------------------------------------------------ pop <- Population(pop, popSize = 800) pop ## ----message=FALSE------------------------------------------------------------ pop <- Population(pop, alleleFrequencies = runif(100)) pop ## ----message=FALSE------------------------------------------------------------ pop <- addEffects(pop) pop ## ----message=FALSE------------------------------------------------------------ pop <- addEffects(pop, distrib = runif) ## ----message=FALSE------------------------------------------------------------ effects <- c( 1.2, 1.5, -0.3, -1.4, 0.8, 2.4, 0.2, -0.8, -0.4, 0.8, -0.2, -1.4, 1.4, 0.2, -0.9, 0.4, -0.8, 0.0, -1.1, -1.3) pop <- addEffects(pop, effects = effects) getAddCoefs(pop) ## ----message=FALSE------------------------------------------------------------ pop <- Population(pop, narrowh2 = 0.4) getAddCoefs(pop) ## ----message=FALSE------------------------------------------------------------ pop <- attachEpiNet(pop) pop ## ----fig.width=5, fig.asp=1, fig.align='center', fig.cap="A random epistatic network generated between 20 QTL."---- epinet <- getEpiNet(pop) plot(epinet) ## ----message=FALSE, fig.width=5, fig.asp=1, fig.align='center', fig.cap="An epistatic network generated using the Barabasi-Albert model between 20 QTL."---- pop <- attachEpiNet(pop, scaleFree = TRUE) plot(getEpiNet(pop)) ## ----message=FALSE, fig.width=5, fig.asp=1, fig.align='center', fig.cap="An epistatic network where 7 QTL have no epistatic effects."---- pop <- attachEpiNet(pop, scaleFree = TRUE, additive = 7) plot(getEpiNet(pop)) ## ----message=FALSE, fig.width=5, fig.asp=1, fig.align='center', fig.cap="An epistatic network with a minimum of 2 interactions per epistatic QTL."---- pop <- attachEpiNet(pop, scaleFree = TRUE, additive = 7, m = 2) plot(getEpiNet(pop)) ## ----message=FALSE, fig.width=5, fig.asp=1, fig.align='center', fig.cap="An epistatic network featuring 2-way to 7-way interactions."---- pop <- attachEpiNet(pop, scaleFree = TRUE, additive = 7, m = 2, k=2:7) plot(getEpiNet(pop)) ## ----------------------------------------------------------------------------- inc <- getIncMatrix(pop) dim(inc) ## ----echo=1------------------------------------------------------------------- inc[, 1:5] inci <- which(inc[, 1:5] == 1) incs <- rowSums(inc[, 1:5]) incm <- which(incs == max(incs)) ## ----------------------------------------------------------------------------- rincmat100snp ## ----message=FALSE------------------------------------------------------------ pop <- attachEpiNet(pop, incmat = rincmat100snp) ## ---- fig.width=5, fig.asp=1, fig.align='center', fig.cap="An epistatic network derived from a user-defined incidence matrix."---- plot(getEpiNet(pop)) ## ----message=FALSE------------------------------------------------------------ # Include the 20th QTL in the first interaction mm <- rincmat100snp mm[20, 1] <- 1 pop <- attachEpiNet(pop, incmat = mm) ## ---- fig.width=5, fig.asp=1, fig.align='center', fig.cap="A user-defined epistatic network featuring a single 3-way interaction."---- plot(getEpiNet(pop)) ## ----message=FALSE------------------------------------------------------------ # 20 coefficients, 3 of which are 0 coefs <- sample(c(rep(0, 3), rnorm(17)), 20) pop <- addEffects(pop, effects = coefs) getAddCoefs(pop) ## ---- echo=FALSE-------------------------------------------------------------- foo <- which(getAddCoefs(pop) == 0) ## ----message=FALSE------------------------------------------------------------ pop <- Population(pop, QTL = 1:15) ## ----message=FALSE------------------------------------------------------------ pop <- attachEpiNet(pop, additive = 1:5) ## ----------------------------------------------------------------------------- getIncMatrix(pop) ## ---- fig.width=5, fig.asp=1, fig.align='center', fig.cap="The epistatic network generated with the first 5 of the 15 QTL being purely additive."---- plot(getEpiNet(pop)) ## ----message=FALSE------------------------------------------------------------ coefs <- rnorm(15) coefs[6:10] <- 0 pop <- addEffects(pop, effects = coefs) getAddCoefs(pop) ## ----message=FALSE------------------------------------------------------------ pop <- Population(popSize = 200, map = map100snp, QTL = 20, alleleFrequencies = runif(100), broadH2 = 0.9, narrowh2 = 0.6, traitVar = 40) pop <- addEffects(pop) pop <- attachEpiNet(pop, k = 3) ## ---- fig.width=5, fig.asp=1, fig.align='center', fig.cap="An epistatic network consisting of 3-way interactions."---- plot(getEpiNet(pop)) ## ----------------------------------------------------------------------------- qtls <- which(getIncMatrix(pop)[, 1] > 0) qtls ## ----------------------------------------------------------------------------- interaction1 <- getInteraction(pop, 1) # Return first interaction array interaction1 ## ----------------------------------------------------------------------------- interaction1[2, , ] ## ----------------------------------------------------------------------------- interaction1[2, 1, ] ## ----------------------------------------------------------------------------- interaction1[2, 1, 3] ## ----------------------------------------------------------------------------- components <- getComponents(pop) head(components) ## ----------------------------------------------------------------------------- mean(components$Additive) ## ----------------------------------------------------------------------------- var(components$Additive) ## ----------------------------------------------------------------------------- mean(components$Epistatic) ## ----------------------------------------------------------------------------- var(components$Epistatic) ## ----------------------------------------------------------------------------- mean(components$Environmental) ## ----------------------------------------------------------------------------- var(components$Environmental) ## ----------------------------------------------------------------------------- getAddOffset(pop) ## ----------------------------------------------------------------------------- getEpiOffset(pop) ## ----------------------------------------------------------------------------- mean(components$Phenotype) ## ----------------------------------------------------------------------------- var(components$Phenotype) ## ----------------------------------------------------------------------------- cov(components$Additive, components$Epistatic) cov(components$Additive, components$Environmental) cov(components$Environmental, components$Epistatic) ## ----------------------------------------------------------------------------- cor(components$Additive, components$Epistatic) ## ----------------------------------------------------------------------------- geno <- getGeno(pop) ## ----------------------------------------------------------------------------- geno <- geno[, getQTL(pop)$Index] ## ----------------------------------------------------------------------------- additive <- geno %*% getAddCoefs(pop) + getAddOffset(pop) additive[1:5] ## ----------------------------------------------------------------------------- getComponents(pop)$Additive[1:5] ## ----------------------------------------------------------------------------- head(getEpistasis(pop)) ## ----------------------------------------------------------------------------- epistatic <- rowSums(getEpistasis(pop)) + getEpiOffset(pop) epistatic[1:5] ## ----------------------------------------------------------------------------- getComponents(pop)$Epistatic[1:5] ## ----include=FALSE------------------------------------------------------------ geno2 <- matrix(sample(0:2, 100*5, replace = TRUE, prob = c(0.25, 0.5, 0.25)), 5, 100) ## ----------------------------------------------------------------------------- geno2 <- geno2[, getQTL(pop)$Index] geno2 ## ----------------------------------------------------------------------------- additive2 <- geno2 %*% getAddCoefs(pop) + getAddOffset(pop) additive2[1:5] ## ----------------------------------------------------------------------------- epistatic2 <- rowSums(getEpistasis(pop, geno = geno2)) + getEpiOffset(pop) epistatic2 ## ----message=FALSE------------------------------------------------------------ popRun <- runSim(pop, generations = 150) ## ----------------------------------------------------------------------------- n <- 10 pmf <- 2 * (n - 1:n + 1) / (n * (n + 1)) pmf ## ----echo=FALSE, fig.width=4, fig.asp=1, fig.align='center', fig.cap="Probability mass function for linear ranking selection on a population of 10 individuals."---- df = data.frame(Individual=as.character(1:n), Probability=pmf) ggplot2::ggplot(data = df, ggplot2::aes(x=reorder(Individual, -Probability), y=Probability)) + ggplot2::geom_bar(stat="identity") + ggplot2::xlab("Ranked individuals") + ggplot2::ylab("Selection probability") ## ----message=FALSE------------------------------------------------------------ popRunRank <- runSim(pop, generations = 150, selection = "ranking") popRunBurnIn <- runSim(pop, generations = 150, burnIn = 50, truncSire = 0.1, truncDam = 0.5, roundsSire = 5, roundsDam = 5, litterDist = c(0.1, 0.3, 0.4, 0.2), breedSire = 7) popRunTGV <- runSim(pop, generations = 150, truncSire = 0.1, truncDam = 0.5, fitness = "TGV") popRunEBV <- runSim(pop, generations = 150, truncSire = 0.1, truncDam = 0.5, fitness = "EBV") ## ----fig.width=7, fig.height=4, fig.align='center', fig.cap="A graphical representation of a simulation run using the default parameters."---- plot(popRun) ## ----fig.width=7, fig.height=4, fig.align='center', fig.cap="A graphical representation of a simulation run using linear ranking selection."---- plot(popRunRank) ## ----fig.width=7, fig.height=4, fig.align='center', fig.cap="A graphical representation of a simulation run using truncation selection and a burn-in period of the first 50 generations."---- plot(popRunBurnIn) ## ----fig.width=7, fig.height=4, fig.align='center', fig.cap="A graphical representation of a simulation run using truncation selection based on true genetic values."---- plot(popRunTGV) ## ----fig.width=7, fig.height=4, fig.align='center', fig.cap="A graphical representation of a simulation run using truncation selection based on estimated breeding values."---- plot(popRunEBV) ## ----echo=FALSE--------------------------------------------------------------- allGenoFileName <- system.file("extdata", "geno.epi", package = "epinetr") geno <- loadGeno(allGenoFileName) ## ----eval=FALSE--------------------------------------------------------------- # popRun <- runSim(pop, generations = 150, allGenoFileName = "geno.epi") # geno <- loadGeno("geno.epi") ## ----------------------------------------------------------------------------- geno[1:5, 1:8] ## ----------------------------------------------------------------------------- ped <- getPedigree(popRun) ped[512:517, ] ## ----------------------------------------------------------------------------- qtl <- getQTL(popRun)$Index af <- getAlleleFreqRun(popRun) af[, qtl[1]] ## ----------------------------------------------------------------------------- geno <- getPhased(popRun) geno[1:6, 1:10] ## ----------------------------------------------------------------------------- geno <- getGeno(popRun) geno[1:6, 1:5] ## ----------------------------------------------------------------------------- ID <- getComponents(popRun)$ID ID <- sample(ID, 50) popRun2 <- getSubPop(popRun, ID) ## ----include=FALSE------------------------------------------------------------ pedData <- getPedigree(popRun) pedData <- pedData[,1:3] ## ----------------------------------------------------------------------------- pedData[201:210, ] ## ----message=FALSE------------------------------------------------------------ popRunPed <- runSim(pop, pedigree = pedData) ## ----fig.width=7, fig.height=4, fig.align='center', fig.cap="A graphical representation of a simulation run using a pedigree data frame."---- plot(popRunPed)