From 84bb060ccf6b1b53b820d0adcc77ba3180883c48 Mon Sep 17 00:00:00 2001
From: Voorrips <roeland.voorrips@wur.nl>
Date: Tue, 17 Nov 2020 10:07:49 +0100
Subject: [PATCH] Last changes (hopefully) for publication: demo data for
vignette changed, removed the "options" way to specify the ahcdir and to load
the ahc(complete)lists, other minor internal changes
---
DESCRIPTION | 2 +-
R/PolyHaplotyper.R | 851 ++++++++++--------
data/PolyHaplotyper_demo.RData | Bin 0 -> 15313 bytes
data/demo.RData | Bin 14821 -> 0 bytes
...otyper-vignette.Rmd => PolyHaplotyper.Rmd} | 154 ++--
vignettes/ahclist_6x.RData | Bin 4540 -> 0 bytes
6 files changed, 548 insertions(+), 459 deletions(-)
create mode 100644 data/PolyHaplotyper_demo.RData
delete mode 100644 data/demo.RData
rename vignettes/{PolyHaplotyper-vignette.Rmd => PolyHaplotyper.Rmd} (77%)
delete mode 100644 vignettes/ahclist_6x.RData
diff --git a/DESCRIPTION b/DESCRIPTION
index 5a57615..25a2a96 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -10,7 +10,7 @@ Description: Infer the genetic composition of individuals
License: GPL-2
Encoding: UTF-8
LazyData: true
-RoxygenNote: 7.0.2
+RoxygenNote: 7.1.1
Imports: XML
Suggests: knitr, rmarkdown
VignetteBuilder: knitr
diff --git a/R/PolyHaplotyper.R b/R/PolyHaplotyper.R
index 8834bfe..2b96e61 100644
--- a/R/PolyHaplotyper.R
+++ b/R/PolyHaplotyper.R
@@ -61,21 +61,18 @@
#'@importFrom utils combn flush.console read.table write.table
#'@importFrom stats chisq.test pchisq pbinom qbinom
-#the following line serves to stop devtools::check complaining about
-#some identifiers that do exist but it doesn't find because they are
-#loaded from RData files and/or created in the global evironment:
-utils::globalVariables(names=c("ahclist", "ahccompletelist", "ahcinfo"),
- add=TRUE)
-
-# documentation of the objects in Data/demo.RData:
+# documentation of the objects in Data/PolyHaplotyper_demo.RData:
+#
#'@title dosages of SNP alleles
-#'@description A data.frame with the allele dosages for 28 SNPs in 661
-#'individuals
+#'@description A data.frame with the allele dosages for 30 SNPs in 661
+#'isamples
#'@docType data
#'@keywords datasets
#'@name snpdos
-#'@format a matrix with 28 rows (SNPs) and 661 columns (samples) containing
-#'dosages (integers in 0..6 or NA)
+#'@format a data.frame with 30 rows (SNPs) and 663 columns. The first
+#'two columns contain the SNP name and contig; the remaining columns
+#'contain the allele dosages (integers in 0..6 or NA) of the SNPs in
+#'661 samples
NULL
#'@title pedigree
@@ -83,7 +80,7 @@ NULL
#'@docType data
#'@keywords datasets
#'@name ped
-#'@format a data.frame with 672 rows (individuals, some replicated) and
+#'@format a data.frame with 661 rows (samples) and
#'4 columns: genotype, mother, father, sample_nr
NULL
@@ -109,167 +106,159 @@ allHaplotypes <- function(mrknames) {
eg
} #allHaplotypes
-getahcdir <- function() {
- # for internal use only; get ahcdir from options("PolyHaplotyper_ahcdir");
- # if that doesn't exist, return working dir
+getahcdir <- function(ahcdir) {
+ # for internal use only; get ahcdir from (1) param ahcdir or (2) working dir
# result ends always with "/" so a filename can be pasted onto it
- ad <- unlist(options("PolyHaplotyper_ahcdir"))
- if (is.null(ad)) {
- ad <- getwd()
- } else {
- if (!is.character(ad) || length(ad)!= 1 || !dir.exists(ad) || ad == "")
- ad <- "./"
+ if (is.null(ahcdir)) ahcdir <- getwd() else {
+ if (!is.character(ahcdir) || length(ahcdir) != 1)
+ stop(paste("invalid ahcdir:", ahcdir))
+ if (ahcdir=="") ahcdir <- getwd()
}
- if (substr(ad, nchar(ad), nchar(ad)) != "/")
- ad <- paste0(ad, "/")
- ad
-} #getahcdir within loadAllHaploCombList
-
-
-loadAllHaploCombList <- function(ploidy) {
- # Function is called only by getAllHapcomb.
- # Load ahclist and ahccompletelist into GlobalEnv if needed.
- # Checks variable ahcinfo in GlobalEnv. If that doesn't exist or ahcinfo$ploidy
- # is not equal to ploidy, attempts to load ahclist and ahccompletelist into
- # the GlobalEnv from a file 'ahclist_nx.RData' and 'ahccompletelist_nx.RData'
- # (where n is the ploidy). These files are loaded from the working directory
- # or, if options("PolyHaplotyper_ahcdir" is specified, from that directory.
- # If that is successful, ploidy is stored in ahcinfo$ploidy
- # in the GlobalEnv, else ahcinfo, ahclist, and ahccompletelist are not created
- # or changed.
- # Return value is NULL.
- if (exists("ahcinfo", envir=.GlobalEnv, inherits==FALSE) &&
- exists("ahclist", envir=.GlobalEnv, inherits==FALSE)) {
- if (ahcinfo$ploidy == ploidy) return(NULL)
- #save any unsaved ahclist:
- if (ahcinfo$unsaved) {
- #the currently loaded ahclist is for a different ploidy and unsaved;
- #we save it (to the old ahcdir!) before loading the ahclist for the
- #current ploidy:
- save(ahclist,
- file=paste0(ahcinfo$ahcdir, "ahclist_", ahcinfo$ploidy,"x.RData"))
+ if (!dir.exists(ahcdir))
+ stop(paste("non-existing ahcdir:", ahcdir))
+ if (substr(ahcdir, nchar(ahcdir), nchar(ahcdir)) != "/")
+ ahcdir <- paste0(ahcdir, "/")
+ ahcdir
+} #getahcdir
+
+loadAllHapcombLists <- function(ploidy, nmrk, ahcdir) {
+ # This function is called only at start of inferHaplotypes
+ # ploidy: 1 integer
+ # nmrk: 1 integer, the number of markers in the haploblock
+ # result: a list with elements ahccompletelist, ahclist, complete_nmrk,
+ # ploidy, ahcdir
+
+ checklist <- function(lst, ploidy, namespresent) {
+ if (!is.list(lst) || length(lst) < 1) return(FALSE)
+ if (!all(sapply(lst, class) == "list")) return(FALSE)
+ if (length(lst[[length(lst)]]) == 0) return(FALSE) # in highest nmrk at least one mrkdid should be present
+ cl <- lapply(lst, FUN=function(x) sapply(x, class))
+ if (!all(unlist(cl) == "matrix")) return(FALSE)
+ cl <- lapply(lst, FUN=function(x) sapply(x, nrow))
+ if (!all(unlist(cl) == ploidy)) return(FALSE)
+ lnames <- names(lst[[length(lst)]])
+ if (xor(!namespresent, (is.null(lnames) || any(lnames == "")))) return(FALSE)
+ return(TRUE)
+ } #checklist
+
+ loadlist <- function(listname, ploidy, ahcdir, namespresent, message) {
+ # try to load:
+ if (message) cat(paste("loading ", listname, "..."))
+ x <- tryCatch({
+ fname <- paste0(ahcdir, listname, "_", ploidy, "x.RData")
+ suppressWarnings(n <- load(fname))
+ if (!(listname %in% n))
+ return(paste("file", fname, "does not contain", listname))
+ }, error=function(e) "error")
+ if (identical(x, "error")) {
+ if (message) cat (" not found!\n")
+ return("") # not an error but no list file
+ } else {
+ #checks:
+ if (!checklist(get(listname), ploidy, namespresent))
+ return(paste(listname, " has incorrect format!\n"))
+ if (message) cat (" done!\n")
+ return(get(listname))
}
- suppressWarnings(rm(ahcinfo, ahclist, ahccompletelist,
- envir=.GlobalEnv, inherits==FALSE))
- }
- newinfo <- list()
- newinfo$ploidy <- ploidy
- newinfo$unsaved <- FALSE
- newinfo$complete_nmrk <- 0
- newinfo$ahcdir <- getahcdir()
- #load the ahccompletelist:
- cat("loading ahccompletelist ...")
- x <- tryCatch(
- suppressWarnings(load(paste0(newinfo$ahcdir, "ahccompletelist_",
- ploidy,"x.RData"),
- envir=.GlobalEnv)),
- error=function(e) "error")
- if (identical(x, "error")) {
- cat (" not found!\n")
+ } #loadlist
+
+
+ ahcinfo <- list(ahccompletelist=NULL, ahclist=NULL, complete_nmrk=0,
+ ploidy=ploidy, ahcdir=getahcdir(ahcdir))
+ if (exists("ahccompletelist") && # in .GlobalEnv, loaded by user
+ checklist(lst=ahccompletelist, ploidy=ploidy, namespresent=FALSE)) {
+ ahcinfo$ahccompletelist <- ahccompletelist
} else {
- #checks:
- if (length(ahccompletelist) == 0) stop(" ahccompletelist empty!")
- if (is.null(ahccompletelist[[1]][[1]]) ||
- !is.matrix(ahccompletelist[[1]][[1]]) ||
- nrow(ahccompletelist[[1]][[1]]) != ploidy)
- stop(" ahccompletelist has incorrect format!\n")
-
- newinfo$complete_nmrk <- length(ahccompletelist)
- cat (" done!\n")
+ ahccompletelist <-
+ loadlist(listname="ahccompletelist", ploidy=ploidy, ahcdir=ahcinfo$ahcdir,
+ namespresent=FALSE, message=TRUE)
+ if (is.list(ahccompletelist)) {
+ ahcinfo$ahccompletelist <- ahccompletelist
+ } else if (ahccompletelist != "") stop(ahccompletelist)
+ }
+ if (!is.null(ahcinfo$ahccompletelist))
+ ahcinfo$complete_nmrk <- length(ahcinfo$ahccompletelist)
+ if (ahcinfo$complete_nmrk < nmrk) {
+ # we need an additional ahclist for the mrkdids with larger nmrk
+ if (exists("ahclist") && # in .GlobalEnv, loaded by user
+ checklist(ahclist, ploidy, namespresent=TRUE)) {
+ ahcinfo$ahclist <- ahclist
+ } else {
+ ahclist <- loadlist("ahclist", ploidy, ahcinfo$ahcdir,
+ namespresent=TRUE, message=FALSE)
+ if (is.list(ahclist)) {
+ ahcinfo$ahclist <- ahclist
+ } else {
+ if (ahclist != "") stop(ahclist)
+ ahcinfo$ahclist <- list() # new empty ahclist
+ }
+ }
}
- #load the ahclist:
- x <- tryCatch(
- suppressWarnings(load(paste0(newinfo$ahcdir, "ahclist_",
- ploidy,"x.RData"),
- envir=.GlobalEnv)),
- error=function(e) "error")
- if (identical(x, "error"))
- #create new ahclist:
- assign("ahclist", value=list(), envir=.GlobalEnv, inherits=FALSE)
- #update ahcinfo:
- assign("ahcinfo", value=newinfo, envir=.GlobalEnv, inherits=FALSE)
- NULL
-} #loadAllHaploCombList
+ ahcinfo
+} # loadAllHapcombLists
#'@title get all haplotype combinations that result in the given marker
#'dosages
#'@description get all haplotype combinations that result in the given
#'marker dosages
-#'@usage getAllHapcomb(mrkDosage=NULL, did=NULL, allhap, ploidy,
-#'writeFile=TRUE, progress=FALSE)
+#'@usage getAllHapcomb(mrkDosage=NULL, mrkdid=NULL, nmrk, ahcinfo)
#'@param mrkDosage an integer vector with the dosages of each marker in one
-#'individual, with the markers in the same order as in the columns of allhap.
-#'Either mrkDosage or did must be specified, if both are specified they
-#'should match but this is not checked
-#'@param did dosage ID: in combination with ploidy specifying the mrkDosage.
-#'Either mrkDosage or did must be specified, if both are specified they
-#'should match but this is not checked
-#'@param allhap a matrix as returned by allHaplotypes, for the same markers
-#'in the same order as mrkDosage
-#'@param ploidy the ploidy of the individual; mrkDosage should all be in
-#'0:ploidy
-#'@param writeFile The calculated haplotype combis are stored in list
-#'ahclist in GlobalEnv. Default writeFile TRUE ensures that this list is written
-#'to file each time a new mrkdid is calculated. With writeFile FALSE the list is
-#'not written to file.
-#'@param progress whether to print a message when starting calculations
-#'for a did; default FALSE
+#'individual; all dosages must be in 0:ahcinfo$ploidy or NA
+#'Either mrkDosage or mrkdid must be specified, not both
+#'@param mrkdid marker dosage ID: in combination with nmrk and ploidy specifying
+#'the mrkDosage. Can be NA or an integer, numeric or character which must
+#'specify an integer > 0
+#'Either mrkDosage or mrkdid must be specified, not both
+#'@param nmrk the number of markers (should match mrkDosage if that is specified)
+#'@param ahcinfo a list as returned by loadAllHapcombLists
#'@details Each column of the return value represents one combination of
#'haplotypes that yields the observed marker dosages. If any of the marker
-#'dosages are NA or not in 0:ploidy, a matrix with 0 columns is returned.\cr
-#'This function makes use of precalculated lists that are read from files
-#'if present: ahccompletelist_nx.RData and ahclist_nx.RData (where n is the
-#'ploidy). These lists ahccompletelist and ahclist are stored in the GlobalEnv,
-#'as well as another list ahcinfo. List ahcinfo may be extended with newly
-#'calculated results; if parameter writeFile is TRUE the extended ahclist is
-#'written to file.\cr
-#'These files ahccompletelist_nx.RData and ahclist_nx.RData are expected
-#'in the working directory, or, if options("PolyHaplotyper_ahcdir") is
-#'specified, from that directory.
-#'@return an integer matrix with one row per haplotype (corresponding to allhap)
-#'and one column per combination of haplotypes, with the number of copies
-#'of each haplotype in each combination (the dosages of the haplotypes)
+#'dosages, or mrkdid, are NA a matrix with 0 columns is returned.\cr
+#'This function makes use of precalculated lists ahcinfo$ahccompletelist
+#'and ahcinfo$ahclist; the mrkdid must be present in one of them.
+#'@return an integer matrix with <ploidy> rows and one column for each haplotype
+#'combination that matches the marker dosages; the haplotype numbers
+#'within a column are in increasing order
#'@export
-getAllHapcomb <- function(mrkDosage=NULL, did=NULL, allhap, ploidy,
- writeFile=TRUE, progress=FALSE) {
- # TODO: implement parallel computing of combinations for new did combinations,
- # https://artax.karlin.mff.cuni.cz/r-help/library/Rdsm/html/Rdsm-package.html
- # seems useful (using shared memory) but may need to delay saving ahclist
- # until all threads are finished?
- # Alternative: package parallel, but then no writing to common variables
- # like ahclist until all threads (foreach) are finished
- #initial checks:
- nmrk <- ncol(allhap)
- if (!missing(mrkDosage) && !is.null(mrkDosage) &&
- (length(mrkDosage) != nmrk || !all(mrkDosage %in% 0:ploidy)))
- stop("invalid mrkDosage")
- if (missing(did) || is.null(did)) {
- if (missing(mrkDosage) || is.null(mrkDosage)) {
- stop("mrkDosage or did must be specified")
- } else did <- mrkdidfun(mrkDosage, ploidy=ploidy)
- }
- loadAllHaploCombList(ploidy) #does nothing if list already in GlobalEnv
- # first check if this nmrk is in ahccompletelist:
- if (ahcinfo$complete_nmrk >= nmrk) {
- return(ahccompletelist[[nmrk]][[as.numeric(did)]]) #now as hapcomb
+getAllHapcomb <- function(mrkDosage=NULL, mrkdid=NULL, nmrk, ahcinfo) {
+ # mrkDosage has already been checked at start of inferHaplotypes,
+ # mrkdid are computed
+ # so we don't check either, only find out of did or mrkDosage specified:
+ if (is.null(mrkdid)) {
+ if (anyNA(mrkDosage)) return(matrix(integer(0), nrow=ahcinfo$ploidy))
+ mrkdid <- mrkdidfun(mrkDosage, ploidy=ahcinfo$ploidy)
}
- # else look it up in ahclist, and if needed calculate it:
- if (length(ahclist) >= nmrk)
- didix <- which(names(ahclist[[nmrk]]) == did) else didix <- integer(0)
- if (length(didix) == 1) {
- #the haplotype combinations for this mrkdid were already calculated earlier:
- return(ahclist[[nmrk]][[didix]]) # now as hapcomb
+ if (is.na(mrkdid)) return(matrix(integer(0), nrow=ahcinfo$ploidy))
+ if (nmrk <= ahcinfo$complete_nmrk)
+ return(ahcinfo$ahccompletelist[[nmrk]][[as.integer(mrkdid)]])
+ #debug check:
+ # if (nmrk>length(ahcinfo$ahclist) ||
+ # !(mrkdid %in% names(ahcinfo$ahclist[[nmrk]])))
+ # stop("mrkdid not in ahc(complete)list in getAllHaplocomb")
+ return(ahcinfo$ahclist[[nmrk]][[as.character(mrkdid)]])
+} # getAllHapcomb
+
+calcAllhapcomb4mrkdid <- function(mrkDosage=NULL, mrkdid=NULL, ploidy, allhap) {
+ # mrkDosage: an integer vector with the dosages (0:ploidy) of nmrk markers
+ # mrkdid: 1 integer (marker dosage id); if NULL, mrkDosage is used
+ # ploidy: 1 integer
+ # allhap: a matrix as returned by allHaplotypes, for nmrk markers
+ # in the same order as mrkDosage
+ # result: a matrix with <ploidy> rows and one column for each haplotype
+ # combination that matches mrkdid
+ #
+ # check if mrkdid or mrkDosage is valid (we assume ploidy and nmrk are valid)
+ nmrk <- ncol(allhap)
+ if (is.null(mrkdid)) {
+ # check mrkDosage:
+ if (length(mrkDosage) != nmrk || !all(mrkDosage %in% 0:ploidy))
+ stop("invalid mrkDosage")
+ } else {
+ if (length(mrkdid)!=1 || is.na(mrkdid) || as.integer(mrkdid) != mrkdid ||
+ mrkdid<1 || mrkdid > (ploidy+1)^nmrk)
+ stop("invalid mrkdid")
+ mrkDosage <- mrkdid2mrkdos(mrkdid, nmrk, ploidy)
}
- # did was not present in ahclist so we add it:
- #we copy the list to a local tmplist, modify that,
- #store it in .GlobalEnv and save it, and return the matrix with haplotype
- #combinations.
- if (progress) cat(paste0("getAllHapcomb: calculation for mrkdid=", did, "\n"))
- tmplist <- ahclist
- rm(list="ahclist", envir=.GlobalEnv, inherits=FALSE)
- if (length(tmplist) < nmrk) tmplist[[nmrk]] <- list()
- #calculate new matrix with haplotype combinations
- if (missing(mrkDosage)) mrkDosage <- mrkdid2mrkdos(did, nmrk, ploidy)
#TODO: the next algorithm is less efficient if the SNP dosages are high than if they are
#low (because with low dosages a haplotype combination will sooner exceed the dosage of
#at least one of the SNPs). Therefore in these cases (with average SNP dosage > ploidy/2)
@@ -305,21 +294,9 @@ getAllHapcomb <- function(mrkDosage=NULL, did=NULL, allhap, ploidy,
}
if (ncol(allmat1) == 0) stop("getAllHapcomb: 0 solutions found")
colnames(allmat1) <- NULL
- #reverse the row order and save to tmplist and ahclist:
- didix <- length(tmplist[[nmrk]]) + 1
- tmplist[[nmrk]][[didix]] <- allmat1[nrow(allmat1):1,, drop=FALSE]
- names(tmplist[[nmrk]])[didix] <- did
- assign("ahclist", value=tmplist, envir=.GlobalEnv, inherits=FALSE)
- rm(tmplist)
- if (writeFile) {
- save(ahclist, file=paste0(ahcinfo$ahcdir, "ahclist_", ploidy,"x.RData"))
- } else {
- tmpinfo <- ahcinfo
- tmpinfo$ahcunsaved <- TRUE
- assign("ahcinfo", value=tmpinfo, envir=.GlobalEnv, inherits=FALSE)
- }
- allmat1
-} #getAllHapcomb
+ #reverse the row order:
+ allmat1[nrow(allmat1):1,, drop=FALSE]
+} # calcAllhapcomb4mrkdid
getHapcombCount_1mrk <- function(mrkdosage, ploidy, nhap) {
# mrkdosage: the dosage of the marker (integer of length 1)
@@ -435,7 +412,8 @@ completeAllHaploComb <- function(ploidy, nmrk, savesec=1800, printsec=60,
lastprint <- proc.time()
if ((lastsave-lastprint)[3] > savesec) {
save(ahc, file=paste0(fname, "_", paste(allset, collapse="-"),
- ".RData"))
+ ".RData"),
+ version=2)
if (file.exists(lastfile)) file.remove(lastfile)
lastfile <- fname
lastsave <- proc.time()
@@ -451,7 +429,7 @@ completeAllHaploComb <- function(ploidy, nmrk, savesec=1800, printsec=60,
for (d in seq_along(ahc)) {
ahc[[d]] <- ahc[[d]][ploidy:1, 1:lasthap[d], drop=FALSE]
}
- save(ahc, file=paste0(fname, ".RData"))
+ save(ahc, file=paste0(fname, ".RData"), version=2)
if (file.exists(lastfile)) file.remove(lastfile)
invisible(ahc)
} #completeAllHaploComb
@@ -513,7 +491,9 @@ build_ahccompletelist <- function(ploidy, maxmrk, savesec=1800, printsec=300,
stop("overwrite must be FALSE or TRUE")
if (file.exists(fname)) { # in working directory
otherloc <- FALSE
- load(fname) # contains ahccompletelist
+ x <- load(fname) # contains ahccompletelist
+ if (length(x) != 1 || x !="ahccompletelist")
+ stop(paste("the existing file", fname, "does not contain (only) ahccompletelist"))
if (length(ahccompletelist) == maxmrk ||
(length(ahccompletelist) > maxmrk && !shorten)) {
# file in working dir already ok
@@ -534,7 +514,7 @@ build_ahccompletelist <- function(ploidy, maxmrk, savesec=1800, printsec=300,
if (length(ahccompletelist) >= maxmrk) {
ahccompletelist <- ahccompletelist[1:newmaxmrk]
file.rename(from=fname, to=paste0(fname,".backup"))
- save(ahccompletelist, file=fname)
+ save(ahccompletelist, file=fname, version=2)
file.remove(paste0(fname,".backup"))
} else {
if (totHapcombCount(ploidy=ploidy, nmrk=maxmrk) > 2e9) {
@@ -549,7 +529,7 @@ build_ahccompletelist <- function(ploidy, maxmrk, savesec=1800, printsec=300,
completeAllHaploComb(ploidy=ploidy, nmrk=nmrk, savesec=savesec,
printsec=printsec,
fname=paste0(tmpfname, nmrk))
- save(ahccompletelist, file=fname)
+ save(ahccompletelist, file=fname, version=2)
cat(paste0("saved ", fname, "for block sizes up to ", nmrk, "markers\n"))
file.remove(paste0(tmpfname, nmrk, ".RData")) #left by completeAllHaploComb
file.remove(paste0(fname,".backup"))
@@ -1033,27 +1013,27 @@ getHaplotypeNames <- function(haploblock, hapcount, sep="_") {
paste0(haploblock, sep, padded(1:hapcount))
}
-##!title Assign haplotype combinations based on output of inferHaps_noFSs
-##!description Assign haplotype combinations based on output of inferHaps_noFSs
+##!title Assign haplotype combinations based on output of inferHaps_noFS
+##!description Assign haplotype combinations based on output of inferHaps_noFS
##!usage hapcomb_from_IHlist(ihl, mrkdids, ploidy, nocombs=FALSE)
##!param haploblockname Prefix to which the (zero-padded) haplotype numbers
##!are appended with separator '_'.
-##!param ihl A list as returned by inferHaps_noFSs
+##!param ihl A list as returned by inferHaps_noFS
##!param mrkdids the mrkdids of all individuals
-##!param ploidy the ploidy (the number of rows in the retrned matrix)
+##!param ahcinfo a list as returned by loadAllHapcombLists
##!param nocombs if FALSE (default) the returned list does not have the matrix
##!hapcombs, only the vector ahccols
-##!details The results of inferHaps_noFSs, as passed in parameter ihl, determine
+##!details The results of inferHaps_noFS, as passed in parameter ihl, determine
##!the haplotype combinations that match each mrkdid. This function translates
##!these results in a set of haplotype combinations for each individual.
##!return a list with
##!2 element: $hapcombs is the matrix as before, $ahccols is an integer vector
##!with for each indiv the column number in its mrkdid element of
##!ahc(complete)list (or NA if more than 1 ahccol possible).
-hapcomb_from_IHlist <- function(ihl, mrkdids, ploidy, nocombs=FALSE) {
+hapcomb_from_IHlist <- function(ihl, mrkdids, ahcinfo, nocombs=FALSE) {
#haplotypeNames <- getHaplotypeNames(haploblockname, nrow(ihl$allhap))
nind <- sum(ihl$mrkdidsTable)
- hapcombs <- matrix(NA_integer_, nrow=ploidy, #length(haplotypeNames),
+ hapcombs <- matrix(NA_integer_, nrow=ahcinfo$ploidy, #length(haplotypeNames),
ncol=length(mrkdids),
dimnames=list(NULL, names(mrkdids)))
# we also return for each indiv the column nr in ahc(complete)list
@@ -1067,10 +1047,10 @@ hapcomb_from_IHlist <- function(ihl, mrkdids, ploidy, nocombs=FALSE) {
if (length(ihl$hclist[[i]]) == 1) {
#we don't assign a haplotype combination if there is more than
#one combination possible
- #hapcombs[, mrkdidind] <- ihl$hclist[[i]][, 1]
ahccols[mrkdidind] <- ihl$hclist[[i]]
if (!nocombs) hapcombs[, mrkdidind] <-
- getAllHapcomb(did=did, allhap=ihl$allhap, ploidy=ploidy)[,ihl$hclist[[i]]]
+ getAllHapcomb(mrkdid=did, nmrk=ncol(ihl$allhap),
+ ahcinfo=ahcinfo)[,ihl$hclist[[i]]]
}
}
res <- list(ahccols=ahccols)
@@ -1184,9 +1164,9 @@ checkmrkDosage <- function(mrkDosage, ploidy, indiv=NULL, markers=NULL,
} #checkmrkDosage
-single_cycle_infer <- function(allhap, mrkdidsTable, ploidy, knownHap,
+single_cycle_infer <- function(allhap, mrkdidsTable, ahcinfo, knownHap,
useKnownDosage, progress) {
- #function used only by inferHaps_noFSs
+ #function used only by inferHaps_noFS
#For each mrkdid, select the hapcombs that require the least new haplotypes
#(not present in knownHap); if useKnownDosages is TRUE, among these select the
#hapcombs with the maximum dosage of known haplotypes.
@@ -1209,8 +1189,8 @@ single_cycle_infer <- function(allhap, mrkdidsTable, ploidy, knownHap,
for (i in seq_along(mrkdidsTable)) {
# for the current mrkdid, get the matrix with all possible
# haplotype combinations (different combinations in columns):
- ahc <- getAllHapcomb(did=names(mrkdidsTable)[i],
- allhap=allhap, ploidy=ploidy, progress=progress)
+ ahc <- getAllHapcomb(mrkdid=names(mrkdidsTable)[i],
+ nmrk=ncol(allhap), ahcinfo=ahcinfo)
ahccols <- 1:ncol(ahc)
if (length(knownHap) > 0) {
# Based on comparisons of 20171023, testscheme 1 is clearly worse, so:
@@ -1268,11 +1248,11 @@ single_cycle_infer <- function(allhap, mrkdidsTable, ploidy, knownHap,
##!title infer haplotypes from marker dosages
##!@description infer haplotypes from marker dosages for a single haploblock,
##!treating all individuals as unrelated
-##!usage inferHaps_noFSs(mrkdids, ploidy, markers,
+##!usage inferHaps_noFS(mrkdids, ahcinfo, markers,
##minfrac=c(0.1, 0.01), knownHap=integer(0), useKnownDosage=TRUE, progress=TRUE)
##!param mrkdids a named integer vector with for each indiv its mrkdid; the
##names are the individual names.
-##!param ploidy all marker dosages should be in 0:ploidy or NA
+##!param ahcinfo a list produced by loadAllHapcombFiles
##!param markers a character vector with the names of the markers in the
##haploblock, in the same order as used for calculating the mrkdids
##!param minfrac a vector of 1 or 2 fractions, the second smaller than the
@@ -1340,7 +1320,7 @@ single_cycle_infer <- function(allhap, mrkdidsTable, ploidy, knownHap,
##!for each of the remaining haplotype combinations: (a subset of the columns
##!of) the matrix returned by getAllHapcomb for that mrkdid}
##!}
-inferHaps_noFSs <- function(mrkdids, ploidy, markers,
+inferHaps_noFS <- function(mrkdids, ahcinfo, markers,
minfrac=c(0.1, 0.01),
knownHap=integer(0), useKnownDosage=TRUE,
progress=TRUE) {
@@ -1359,7 +1339,7 @@ inferHaps_noFSs <- function(mrkdids, ploidy, markers,
while (TRUE) {
cycle <- cycle + 1
res <- single_cycle_infer(allhap=allhap, mrkdidsTable=mrkdidsTable,
- ploidy=ploidy, knownHap=knownHap,
+ ahcinfo=ahcinfo, knownHap=knownHap,
useKnownDosage=useKnownDosage[ukd],
progress=progress)
knownNew <- sort(union(origHap, which(res$nPresent >= minfrac[1] * nind)))
@@ -1369,7 +1349,7 @@ inferHaps_noFSs <- function(mrkdids, ploidy, markers,
#only be gained, not lost (except the origHap which need not be supported
#by the current data, but these have been added back to knownNew).
#Therefore the following warning was generated if haplotypes were lost:
- #warning(paste0("inferHaps_noFSs: lost in cycle ", cycle, ": ",
+ #warning(paste0("inferHaps_noFS: lost in cycle ", cycle, ": ",
# paste(haploLost, collapse=" ")))
#However it seems that it is possible to lose previous must-have
#haplotypes, although this does not happen often.
@@ -1384,7 +1364,7 @@ inferHaps_noFSs <- function(mrkdids, ploidy, markers,
#recalculate the result of the first cycle - very fast and avoids
#having to store it:
res <- single_cycle_infer(allhap=allhap, mrkdidsTable=mrkdidsTable,
- ploidy=ploidy, knownHap=origHap,
+ ahcinfo=ahcinfo, knownHap=origHap,
useKnownDosage=useKnownDosage[ukd],
progress=FALSE)
@@ -1404,7 +1384,7 @@ inferHaps_noFSs <- function(mrkdids, ploidy, markers,
if (length(useKnownDosage) == 2 && length(knownHap_ukd) == 1 &&
!is.null(knownHap_ukd[[1]])) {
res <- single_cycle_infer(allhap=allhap, mrkdidsTable=mrkdidsTable,
- ploidy=ploidy, knownHap=knownHap_ukd[[1]],
+ ahcinfo=ahcinfo, knownHap=knownHap_ukd[[1]],
useKnownDosage=FALSE,
progress=FALSE)
knownHap <- sort(union(origHap, which(res$nPresent >= minfrac[1] * nind)))
@@ -1432,11 +1412,11 @@ inferHaps_noFSs <- function(mrkdids, ploidy, markers,
names(hac1comb) <- names(mrkdidsTable)[onecomb]
# a vector with the mrkdids as names and the ahc column nr as values
# for the mrkdids with just 1 hapcomb
- hapcomb <- matrix(NA_integer_, nrow=ploidy, ncol=length(hac1comb))
+ hapcomb <- matrix(NA_integer_, nrow=ahcinfo$ploidy, ncol=length(hac1comb))
for (d in seq_along(hac1comb)) {
hapcomb[, d] <-
- getAllHapcomb(did=names(hac1comb)[d], allhap=allhap,
- ploidy=ploidy)[, hac1comb[d]]
+ getAllHapcomb(mrkdid=names(hac1comb)[d], nmrk=ncol(allhap),
+ ahcinfo=ahcinfo)[, hac1comb[d]]
}
#calculate the nr of individuals in which each haplotype occurs:
#select the mrkdids for which only 1 haplotype combination is possible,
@@ -1446,10 +1426,11 @@ inferHaps_noFSs <- function(mrkdids, ploidy, markers,
final_hap <-
sort(union(origHap,
which (hapfrq >= 2 & hapfrq > minfrac[2]*sum(mrkdidsTable))))
- resfinal <- single_cycle_infer(allhap=allhap, mrkdidsTable=mrkdidsTable,
- ploidy=ploidy, knownHap=final_hap,
- useKnownDosage=useKnownDosage[length(useKnownDosage)],
- progress=progress)
+ resfinal <-
+ single_cycle_infer(allhap=allhap, mrkdidsTable=mrkdidsTable,
+ ahcinfo=ahcinfo, knownHap=final_hap,
+ useKnownDosage=useKnownDosage[length(useKnownDosage)],
+ progress=progress)
mrkdidHaccountfinal <-
vapply(resfinal$hclist[names(mrkdidsTable)], FUN=length, FUN.VALUE=0)
#keep resfinal instead of res only if all mrkdids that had only one
@@ -1480,7 +1461,7 @@ inferHaps_noFSs <- function(mrkdids, ploidy, markers,
res$allhap <- allhap
res$mrkdidsTable <- mrkdidsTable
res
-} #inferHaps_noFSs
+} #inferHaps_noFS
chisq.probtest <- function(x, p=rep(1/length(x), length(x)),
@@ -1740,7 +1721,7 @@ testFSseg <- function(parhac, DRrate, errfrac,
##!usage solveOneFS(mrkDosage, ihl, FS, P, priorPhac, hacXmat, maxparcombs,
##!ploidy, minfrac, errfrac, DRrate, minPseg)
##!param mrkDosage dosage matrix with only the markers in the current haploblock
-##!param ihl a list returned by inferHaps_noFSs
+##!param ihl a list returned by inferHaps_noFS
##!param FS character vector with sample names of the current FS
##!param P vector with 2 sample names, one for parent1 and one for parent2.
##!There should be corresponding columns in mrkDosage , and orig_ihl should be
@@ -1762,7 +1743,7 @@ testFSseg <- function(parhac, DRrate, errfrac,
##!considered to be certainly present if it occurs in at least a fraction
##!minfrac[1] of all individuals; in the final stage for the "other"
##!individuals (those that do not belong to the FS or its parents) this fraction
-##!is lowered to minfrac[2]; see also inferHaps_noFSs
+##!is lowered to minfrac[2]; see also inferHaps_noFS
##!param errfrac the assumed fraction marker genotypes with an error (over all
##!markers in the haploblock). The errors are assumed to be uniformly distributed
##!over all except the original marker dosage combinations (mrkdids)
@@ -1798,10 +1779,10 @@ testFSseg <- function(parhac, DRrate, errfrac,
##! NEW: haxXmat not used or output anymore
# solveOneFS is not exported: internal function
solveOneFS <- function(mrkdids, ihl, FS, P, knownPhac, triedPhacs, maxparcombs,
- ploidy, minfrac, errfrac, DRrate, minPseg) {
+ ahcinfo, minfrac, errfrac, DRrate, minPseg) {
calcParcombok <- function(potPhacs, oldPhacs, allhap, Pdids, FSfrac,
- P.reject=c(0.05, 0.001)) {
+ P.reject=c(0.05, 0.001), ahcinfo) {
# Function within solveOneFS
# Calculates the combinations of parental haplotype combinations
# that can explain (almost) all FS individuals
@@ -1818,6 +1799,7 @@ solveOneFS <- function(mrkdids, ihl, FS, P, knownPhac, triedPhacs, maxparcombs,
# Pdids: vector of 2 integers: the parental mrkdids
# P.reject: two P values, decending: used in a binomial test for the
# nr of unexplained FS indiv
+ # ahcinfo: a list returned by loadAllHapcombLists
#TODO: parallellize; see also getAllHapcomb
# Phacs are column numbers in ahc(complete)list
if (all(potPhacs[[1]] %in% oldPhacs[[1]]) &&
@@ -1837,10 +1819,10 @@ solveOneFS <- function(mrkdids, ihl, FS, P, knownPhac, triedPhacs, maxparcombs,
if (isnewpc) {
#this is a new combination, check:
FScombs <- getFScombs(cbind(
- getAllHapcomb(did=Pdids[1], allhap=allhap,
- ploidy=ploidy)[, potPhacs[[1]][pc1]],
- getAllHapcomb(did=Pdids[2], allhap=allhap,
- ploidy=ploidy)[, potPhacs[[2]][pc2]]))
+ getAllHapcomb(mrkdid=Pdids[1], nmrk=ncol(allhap),
+ ahcinfo=ahcinfo)[, potPhacs[[1]][pc1]],
+ getAllHapcomb(mrkdid=Pdids[2], nmrk=ncol(allhap),
+ ahcinfo=ahcinfo)[, potPhacs[[2]][pc2]]))
FSmrkdos <- hapcomb2mrkdos(FScombs, allhap=allhap)
FScombmrkdids <- colnames(FSmrkdos)
missedmrkdids <- setdiff(names(FSmrkdidsTable), FScombmrkdids)
@@ -1886,7 +1868,7 @@ solveOneFS <- function(mrkdids, ihl, FS, P, knownPhac, triedPhacs, maxparcombs,
} #calcParcombok within solveOneFS
calcParcombstats <- function(pco, parcombok, Pdids, parcombstats,
- didNhapcombLst, allhap, ploidy) {
+ didNhapcombLst, allhap, ahcinfo) {
# pco: number, parental combination (index to parcombok rows)
# Pdids: vector of 2 integers: the parental mrkdids
# parcombstats: list with statistics for each parental combination, see
@@ -1908,10 +1890,10 @@ solveOneFS <- function(mrkdids, ihl, FS, P, knownPhac, triedPhacs, maxparcombs,
# $expHCfreq: vector of the expected hac frequencies in order of hapcomb
# $dupdid: logical vector, TRUE where mrkdid[i] == mrkdid[i-1], else FALSE
parhac <- cbind(
- getAllHapcomb(did=Pdids[1], allhap=allhap,
- ploidy=ploidy)[, parcombok[pco, 1]],
- getAllHapcomb(did=Pdids[2], allhap=allhap,
- ploidy=ploidy)[, parcombok[pco, 2]])
+ getAllHapcomb(mrkdid=Pdids[1], nmrk=ncol(allhap),
+ ahcinfo=ahcinfo)[, parcombok[pco, 1]],
+ getAllHapcomb(mrkdid=Pdids[2], nmrk=ncol(allhap),
+ ahcinfo=ahcinfo)[, parcombok[pco, 2]])
# parhac does not contain NAs!
segdat <- testFSseg(parhac=parhac, DRrate=DRrate, errfrac=errfrac,
FSmrkdidsTable=FSmrkdidsTable, allhap=allhap)
@@ -1947,7 +1929,7 @@ solveOneFS <- function(mrkdids, ihl, FS, P, knownPhac, triedPhacs, maxparcombs,
tmp <- calcParcombstats(pco=pco, parcombok=parcombok, Pdids=Pdids,
parcombstats=parcombstats,
didNhapcombLst=didNhapcombLst,
- allhap=allhap, ploidy=ploidy)
+ allhap=allhap, ahcinfo=ahcinfo)
parcombstats <- tmp$parcombstats
didNhapcombLst <- tmp$didNhapcombLst
} # for pco
@@ -1978,7 +1960,7 @@ solveOneFS <- function(mrkdids, ihl, FS, P, knownPhac, triedPhacs, maxparcombs,
Phacs <- vector("list", 2) #Phacs: Parental haplotype combinations (column nrs in ahc(complete)list)
for (p in 1:2) {
if (is.na(Pdids[p])) {
- Phacs[[p]] <- integer(0) # matrix(nrow=ploidy, ncol=0)
+ Phacs[[p]] <- integer(0)
} else {
#get the parental haplotype combinations based on ihl:
if (is.na(knownPhac[p])) {
@@ -2032,7 +2014,7 @@ solveOneFS <- function(mrkdids, ihl, FS, P, knownPhac, triedPhacs, maxparcombs,
parcombok <- calcParcombok(potPhacs=Phacs,
oldPhacs=PahccolsTried,
allhap=ihl$allhap, Pdids=Pdids,
- FSfrac=FSfrac)
+ FSfrac=FSfrac, ahcinfo=ahcinfo)
PahccolsTried <- Phacs
step1pc <- parcombok$setaside_pcs
parcombok <- parcombok$parcombok
@@ -2045,16 +2027,17 @@ solveOneFS <- function(mrkdids, ihl, FS, P, knownPhac, triedPhacs, maxparcombs,
maxP <- max(maxP, parcombstats$P, na.rm=TRUE)
#print(paste("firstcombs, firstsel: P=", parcombstats$P, "length(sel)=", length(sel)))
if (length(sel) == 0) {
- #No parental combinations from inferHaps_noFSs, preselected with
+ #No parental combinations from inferHaps_noFS, preselected with
#stringent P.reject, fit
#1. See how many parental combinations are possible without first
- # inferHaps_noFSs (i.e. with getAllHapcomb).
+ # inferHaps_noFS (i.e. with getAllHapcomb).
# If p1*p2 < maxparcombs we try them all.
nwPhacs <- vector("list", 2)
for (p in 1:2) {
if (is.na(knownPhac[p])) {
- nwPhacs[[p]] <- 1:ncol(getAllHapcomb(did=Pdids[p], allhap=ihl$allhap,
- ploidy=ploidy))
+ nwPhacs[[p]] <-
+ 1:ncol(getAllHapcomb(mrkdid=Pdids[p], nmrk=ncol(ihl$allhap),
+ ahcinfo=ahcinfo))
} else {
nwPhacs[[p]] <- knownPhac[p]
}
@@ -2073,7 +2056,7 @@ solveOneFS <- function(mrkdids, ihl, FS, P, knownPhac, triedPhacs, maxparcombs,
parcombok <- calcParcombok(potPhacs=nwPhacs,
oldPhacs=PahccolsTried,
allhap=ihl$allhap, Pdids=Pdids,
- FSfrac=FSfrac)
+ FSfrac=FSfrac, ahcinfo=ahcinfo)
for (p in 1:2)
PahccolsTried[[p]] <- sort(union(PahccolsTried[[p]], nwPhacs[[p]]))
step2pc <- parcombok$setaside_pcs
@@ -2249,7 +2232,7 @@ checkHaploblock <- function(haploblock, mrkDosage) {
# FS solutions (i.e. not all the ones present at the end of step C,
# which might be conflicting between FSs in the same group)
-hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
+hapOneBlock <- function(mrkDosage, hbname, ahcinfo, parents, FS,
minfrac, errfrac, DRrate,
dropUnused, maxparcombs, minPseg,
knownHap, progress) {
@@ -2257,7 +2240,8 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
# mrkDosage is now checked, a matrix with only the rows for the markers
# in the current haploblock
# hbname: the name of the haploblock
- # all other parameters as in inferHaplotypes
+ # all other parameters as in inferHaplotypes, plus ahcinfo: a list
+ # returned by loadAllHapcombLists at start of inferHaplotypes
# parents and FS are checked, and FS has length 0 if no FS families defined
# Return value is a list with elements:
# $hapdos: a matrix with one column per individual and one row per
@@ -2268,16 +2252,16 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
# $mrkdids: the observed mrkdids for all individuals, in same order as the
# columns of $hapdos; names are individual names
# $allhap: a matrix with the composition of all haplotypes
- # $FSsFit: logical vector with one element per FS family; TRUE if a (or
+ # $FSfit: logical vector with one element per FS family; TRUE if a (or
# more than one) acceptable solution for the FS is found (although
# if more than one solution they might not be used if unclear
# which is the best solution)
- # $FSsMessages: character vector with one message per FS. If not FSsFit then
- # an error message; if FSsFit either "" or a report of mrkdids that
+ # $FSmessages: character vector with one message per FS. If not FSfit then
+ # an error message; if FSfit either "" or a report of mrkdids that
# were not expected or that were not assigned a haplotype combination
# because multiple possible (and likely) combinations would give the
# same mrkdid
- # $FSsPval: numeric vector with a P value for each FS. If one solution was
+ # $FSpval: numeric vector with a P value for each FS. If one solution was
# chosen, the chi-squared P-value of that solution; else the maximum
# P-value over all considered solutions
@@ -2476,31 +2460,31 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
} # function newSolstats
# hapOneBlock START ####
- allmrkdids <- mrkdos2mrkdid(mrkDosage, ploidy=ploidy, check=FALSE)
+ allmrkdids <- mrkdos2mrkdid(mrkDosage, ploidy=ahcinfo$ploidy, check=FALSE)
if (length(FS) == 0) {
- ihl <- inferHaps_noFSs(mrkdids=allmrkdids, ploidy=ploidy,
+ ihl <- inferHaps_noFS(mrkdids=allmrkdids, ahcinfo=ahcinfo,
markers=rownames(mrkDosage),
minfrac=minfrac,
knownHap=knownHap, useKnownDosage=TRUE,
progress=progress)
hapdos <- hapcomb2hapdos(
hapcomb=hapcomb_from_IHlist(ihl=ihl, mrkdids=allmrkdids,
- ploidy=ploidy, nocombs=FALSE)$hapcomb,
+ ahcinfo=ahcinfo, nocombs=FALSE)$hapcomb,
nhap=nrow(ihl$allhap))
result <- list(
hapdos=hapdos, mrkdids=allmrkdids, markers=colnames(ihl$allhap))
} else {
- # one or more FSs
+ # one or more FS
#first calculate ihl without using FS info:
# - for elements of the returned list
# - for a pre-selection of likely parental hapcomb
- ihl <- inferHaps_noFSs(mrkdids=allmrkdids, ploidy=ploidy,
+ ihl <- inferHaps_noFS(mrkdids=allmrkdids, ahcinfo=ahcinfo,
markers=rownames(mrkDosage),
minfrac=minfrac[1], #without final inference
knownHap=knownHap, useKnownDosage=FALSE,
progress=progress)
ahccols <- hapcomb_from_IHlist(ihl=ihl, mrkdids=allmrkdids,
- ploidy=ploidy, nocombs=TRUE)$ahccols
+ ahcinfo=ahcinfo, nocombs=TRUE)$ahccols
ihl <- ihl[names(ihl) != "hclist"] #sometimes big, not used here, delete
nhap <- nrow(ihl$allhap)
# we calculate the FS sizes per haploblock because of differences
@@ -2596,7 +2580,7 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
# No break, go on with new cycle
solstats[[cycle]] <- newSolstats(knownHap, FS) #only the specified haplotypes
# 1: for each FS find a set of acceptable solutions
- ih <- inferHaps_noFSs(mrkdids=allmrkdids, ploidy=ploidy,
+ ih <- inferHaps_noFS(mrkdids=allmrkdids, ahcinfo=ahcinfo,
markers=rownames(mrkDosage),
minfrac=minfrac[1],
knownHap=solstats[[cycle-1]]$knownHap,
@@ -2609,7 +2593,7 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
knownPhac=solstats[[cycle-1]]$Pahccols[fs,],
triedPhacs=FSdata[[fs]]$PahccolsTried,
maxparcombs=maxparcombs,
- ploidy=ploidy, minfrac=minfrac, errfrac=errfrac,
+ ahcinfo=ahcinfo, minfrac=minfrac, errfrac=errfrac,
DRrate=DRrate, minPseg=minPseg)
# we update FSdata[[fs]]:
FSdata[[fs]]$message <- sof$message
@@ -2789,7 +2773,7 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
# knownHap for this cycle
# - mark the hacs of the parent(s) with the same hac over all
# solutions as "known"
- # - mark all FSs in group as FSsFit
+ # - mark all FSs in group as FSfit
tmpapc <- allparcombs[bestcomb,, drop=FALSE]
# limit the solstats[[cycle]]$pcos to those allowed in any
@@ -2816,8 +2800,8 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
# for each parent we determine the haplotypes common to all selected
# solutions and add these to the knownHap for this cycle:
for (p in seq_len(ncol(tmpapc))) {
- phacs <- getAllHapcomb(did=parmrkdids[p], allhap=ihl$allhap,
- ploidy=ploidy)[, tmpapc[, p], drop=FALSE]
+ phacs <- getAllHapcomb(mrkdid=parmrkdids[p], nmrk=ncol(ihl$allhap),
+ ahcinfo=ahcinfo)[, tmpapc[, p], drop=FALSE]
for (i in seq_len(ncol(phacs))) {
if (i == 1) {
phaps <- unique(phacs[, 1])
@@ -2944,10 +2928,10 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
# we assign the FS individuals from the FS calculations based on these Phacs
# we assign the rest (including the non-fitted FSs and their parents) based
# on the known haplotypes in solstats[[bestcycle]]$knownhap
- hac <- matrix(integer(0), nrow=ploidy, ncol=ncol(mrkDosage),
+ hac <- matrix(integer(0), nrow=ahcinfo$ploidy, ncol=ncol(mrkDosage),
dimnames=list(NULL, colnames(mrkDosage)))
- FSsFit <- logical(length(FS)) # all FALSE
- FSsPval <- rep(NA_real_, length(FS)) # all NA
+ FSfit <- logical(length(FS)) # all FALSE
+ FSpval <- rep(NA_real_, length(FS)) # all NA
rest <- setdiff(indiv, c(parents, unlist(FS)))
imputedGeno <- matrix(integer(0), nrow=nrow(mrkDosage),
dimnames=list(rownames(mrkDosage), NULL))
@@ -2970,7 +2954,7 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
# else: length(pcorows)==0, nrow(parcombok)==0 and message not '';
# message already contains "no solution found" or something more
# informative(?), we don't overwrite it.
- FSsPval[fs] <- FSdata[[fs]]$maxPseg
+ FSpval[fs] <- FSdata[[fs]]$maxPseg
rest <- union(rest, c(parents[fs,], FS[[fs]])) # note that some parents
# may still be filled in from other FSs
} else {
@@ -2989,19 +2973,19 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
# one single selected hac for this parent, assign:
Pdid <- allmrkdids[names(allmrkdids)==parents[fs, p]]
hac[, colnames(hac)==parents[fs, p]] <-
- getAllHapcomb(did=Pdid, allhap=ihl$allhap,
- ploidy=ploidy)[, solstats[[bestcycle]]$selPahccols[fs, p]]
+ getAllHapcomb(mrkdid=Pdid, nmrk=ncol(ihl$allhap),
+ ahcinfo=ahcinfo)[, solstats[[bestcycle]]$selPahccols[fs, p]]
}
}
if (!is.na(solstats[[bestcycle]]$selpco[fs])) {
# one solution for this FS selected (but there might be more
# than one solstats[[bestcycle]]$pcos[[fs]])
- FSsFit[[fs]] <- TRUE
- FSsPval[fs] <- FSdata[[fs]]$Pseg[solstats[[bestcycle]]$selpco[fs]]
+ FSfit[[fs]] <- TRUE
+ FSpval[fs] <- FSdata[[fs]]$Pseg[solstats[[bestcycle]]$selpco[fs]]
} else {
# no selected pco for this FS but >=1 FS solutions (pcorows) found
FSdata[[fs]]$message <- "multiple FS solutions"
- FSsPval[fs] <- max(FSdata[[fs]]$Pseg[pcorows])
+ FSpval[fs] <- max(FSdata[[fs]]$Pseg[pcorows])
}
# Next we assign all FS indivs with a mrkdid that has a unique solution
# over all (1 or more) possible parental combinations:
@@ -3016,7 +3000,7 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
knownPhac=FSdata[[fs]]$parcombok[pcr,], #known Phac both parents
triedPhacs=list(integer(0), integer(0)), #to have sof recomputed
maxparcombs=maxparcombs,
- ploidy=ploidy, minfrac=minfrac, errfrac=errfrac,
+ ahcinfo=ahcinfo, minfrac=minfrac, errfrac=errfrac,
DRrate=DRrate, minPseg=minPseg)
# This sof has only one row in its parcombok because only one, valid,
# combination of parents was tested! So sof$parcombok is different
@@ -3059,7 +3043,7 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
colnames(mrkDosage) %in% FS[[fs]])
if (length(NAind) > 0 && length(NAind) < 0.5 * length(FS[[fs]])) {
expmrkdos <- mrkdid2mrkdos(dosageIDs=expdid, nmrk=nrow(mrkDosage),
- ploidy=ploidy)
+ ploidy=ahcinfo$ploidy)
impdids <-
matrix(c(NAind, rep(NA, length(NAind))), nrow=2, byrow=TRUE,
dimnames=list(c("NAind", "mrkdid"), NULL))
@@ -3096,7 +3080,7 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
# the imputed mrkdids should be among the expected mrkdids, check:
if (anyNA(impdids[2,]) || !all(impdids[2,] %in% expdid))
stop("error in impdids")
- if (impstats$stats$P > 0.1 * FSsPval[fs]) {
+ if (impstats$stats$P > 0.1 * FSpval[fs]) {
# with the extra imputed FS mrkdids the fit of the FS family
# is not too much worse, so add the imputed mrkdids
allmrkdids[impdids[1,]] <- impdids[2,]
@@ -3104,7 +3088,7 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
# of the return value:
FSimpGeno <-
mrkdid2mrkdos(dosageIDs=impdids[2,], nmrk=nrow(mrkDosage),
- ploidy=ploidy)
+ ploidy=ahcinfo$ploidy)
rownames(FSimpGeno) <- rownames(mrkDosage)
colnames(FSimpGeno) <- colnames(mrkDosage)[impdids[1,]]
imputedGeno <- cbind(imputedGeno, FSimpGeno)
@@ -3152,7 +3136,7 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
# mrkdos2mrkdid(
# mrkDosage=hapcomb2mrkdos(hapcomb=hac[, colnames(hac) %in% ind],
# allhap=ihl$allhap),
- # ploidy=ploidy,
+ # ploidy=ahcinfo$ploidy,
# check=FALSE)
# }
} # for mdhcol
@@ -3163,8 +3147,8 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
# solution or parental conflicts (but not the FSs that were left unsolved
# because of too many equally good solutions)
rest_ind <- intersect(rest, names(allmrkdids)[!is.na(allmrkdids)])
- rest_ihl <- inferHaps_noFSs(mrkdids=allmrkdids[names(allmrkdids) %in% rest_ind],
- ploidy=ploidy,
+ rest_ihl <- inferHaps_noFS(mrkdids=allmrkdids[names(allmrkdids) %in% rest_ind],
+ ahcinfo=ahcinfo,
markers=rownames(mrkDosage),
minfrac=minfrac,
knownHap=solstats[[bestcycle]]$knownHap,
@@ -3176,14 +3160,14 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
allmrkdids==did]
if (length(ind) > 0) {
hac[, colnames(hac) %in% ind] <-
- getAllHapcomb(did=did, allhap=ihl$allhap,
- ploidy=ploidy)[, rest_ihl$hclist[[md]]]
+ getAllHapcomb(mrkdid=did, nmrk=ncol(ihl$allhap),
+ ahcinfo=ahcinfo)[, rest_ihl$hclist[[md]]]
#TODO remove DEBUG CHECK:
# nwdid <-
# mrkdos2mrkdid(
# mrkDosage=hapcomb2mrkdos(hapcomb=hac[, colnames(hac) %in% ind],
# allhap=ihl$allhap),
- # ploidy=ploidy,
+ # ploidy=ahcinfo$ploidy,
# check=FALSE)
}
}
@@ -3285,9 +3269,9 @@ hapOneBlock <- function(mrkDosage, hbname, ploidy, parents, FS,
getHaplotypeNames(haploblock=hbname, hapcount=nrow(ihl$allhap))
result <-
list(hapdos=hapdos, mrkdids=allmrkdids, markers=colnames(ihl$allhap),
- FSsFit=FSsFit,
- FSsMessages=sapply(FSdata, `[[`, "message"),
- FSsPval=FSsPval, imputedGeno=imputedGeno)
+ FSfit=FSfit,
+ FSmessages=sapply(FSdata, `[[`, "message"),
+ FSpval=FSpval, imputedGeno=imputedGeno)
# TODO?:(if rest_ihl with two minfrac values results in additional haplotypes,
# see if these can help solve the remaining FSs? And then re-iterate?
} # one or more FSs
@@ -3350,7 +3334,7 @@ getLinkedFSgroups <- function(parents) {
FSgroupnrs=FSgroupnrs, Parentgroupnrs=Parentgroupnrs)
} #getLinkedFSgroups
-getNewMrkdids <- function(mrkDosages, ploidy, haploblock, maxmrk) {
+getNewMrkdids <- function(mrkDosages, haploblock, maxmrk, ahcinfo) {
# non-exported function only for use by function inferHaplotypes.
# Parameters as the corresponding ones of inferHaplotypes but already
# checked, and mrkDosages as matrix not data.frame with only the selected
@@ -3364,21 +3348,24 @@ getNewMrkdids <- function(mrkDosages, ploidy, haploblock, maxmrk) {
hbsizes <- sapply(haploblock, length)
maxmrk <- min(maxmrk, max(hbsizes))
result <- vector(mode="list", maxmrk)
- loadAllHaploCombList(ploidy)
if (maxmrk <= ahcinfo$complete_nmrk) return(result)
# we need to check all mrkdids for new ones
knownmrkdid <- result # list of same length
for (hbsize in (ahcinfo$complete_nmrk+1):maxmrk) {
- if (length(ahclist) < hbsize) {
+ if (length(ahcinfo$ahclist) < hbsize) {
knownmrkdid[[hbsize]] <- integer(0)
- } else knownmrkdid[[hbsize]] <- as.integer(names(ahclist[[hbsize]]))
+ } else knownmrkdid[[hbsize]] <- as.integer(names(ahcinfo$ahclist[[hbsize]]))
}
for (hb in seq_along(haploblock)) {
if (hbsizes[[hb]] %in% (ahcinfo$complete_nmrk+1):maxmrk) {
- mrkdids <-
- unique(mrkdos2mrkdid(mrkDosages[rownames(mrkDosages) %in%
- haploblock[[hb]],,drop=FALSE],
- ploidy=ploidy, check=FALSE))
+ # mrkdids <-
+ # unique(mrkdos2mrkdid(mrkDosages[rownames(mrkDosages) %in%
+ # haploblock[[hb]],,drop=FALSE],
+ # ploidy=ploidy, check=FALSE))
+ uniqdos <-
+ unique(mrkDosages[match(haploblock[[hb]], rownames(mrkDosages)),,
+ drop=FALSE], MARGIN=2)
+ mrkdids <- mrkdos2mrkdid(uniqdos, ploidy=ahcinfo$ploidy, check=FALSE)
result[[hbsizes[hb]]] <- c(result[[hbsizes[hb]]],
setdiff(mrkdids, knownmrkdid[[hbsizes[hb]]]))
}
@@ -3390,6 +3377,67 @@ getNewMrkdids <- function(mrkDosages, ploidy, haploblock, maxmrk) {
result
} #getNewMrkdids
+addNewMrkdids <- function(newMrkdids, ahcinfo, progress) {
+ # non-exported function only for use by function inferHaplotypes.
+ # newMrkdids: list as produced by getNewMrkdids
+ # Parameters as the corresponding ones of inferHaplotypes but already
+ # checked, and mrkDosages as matrix not data.frame with only the selected
+ # indiv.
+ # output: a list with one element for each haploblock size (nr of markers),
+ # from 1 to the maximum haploblock size in haploblock (or to maxmrk if that
+ # is smaller);
+ # each element has all mrkdid (as integers) that occur over all haploblocks
+ # of that size and for which no allhapcomb is yet available in
+ # ahclist or ahccompletelist
+ newdidcounts <- sapply(newMrkdids, length)
+ totnew <- sum(newdidcounts)
+ if (any(newdidcounts > 0)) {
+ maxnmrk <- max(which(newdidcounts > 0))
+ # add the haplotype combinations for the newdids to ahcinfo$ahclist
+ # and save the new ahclist in ahcdir (overwriting any earlier version)
+ if (progress && (maxnmrk + ahcinfo$ploidy >= 12 ||
+ sum(newdidcounts) > 100)) {
+ cat(paste(sum(newdidcounts),
+ "sets of haplotype combinations need to be calculated;\n"))
+ cat("this may take quite a long time\n")
+ }
+ cum <- 0
+ for (m in which(newdidcounts > 0)) {
+ if (length(ahcinfo$ahclist) < m || !is.list(ahcinfo$ahclist[[m]])) {
+ ahcinfo$ahclist[[m]] <- list()
+ # initialize also all lower nmrk:
+ for (i in seq_len(m)) {
+ if (!is.list(ahcinfo$ahclist[[i]]))
+ ahcinfo$ahclist[[i]] <- list()
+ }
+ }
+ allhap <- allHaplotypes(1:m)
+ for (i in seq_along(newMrkdids[[m]])) {
+ did <- newMrkdids[[m]][i]
+ if (progress) {
+ cat(paste0("\rcalculation of new mrkdids: ", cum + i, " / ", totnew))
+ flush.console()
+ }
+ ahcinfo$ahclist[[m]][[as.character(did)]] <-
+ calcAllhapcomb4mrkdid(mrkdid=did, ploidy=ahcinfo$ploidy, allhap=allhap)
+ # save after every 500 mrkdids so we can resume after crash:
+ if (i %% 500 == 0) {
+ save(ahcinfo$ahclist,
+ file=paste0(ahcinfo$ahcdir, "ahclist_", ahcinfo$ploidy, "x.RData"),
+ version=2)
+ }
+ }
+ cum <- cum + newdidcounts[m]
+ }
+ cat("\n"); flush.console()
+ ahclist <- ahcinfo$ahclist # because save will only save an entire R object, not a list element
+ save(ahclist,
+ file=paste0(ahcinfo$ahcdir, "ahclist_", ahcinfo$ploidy, "x.RData"),
+ version=2)
+ }
+ ahcinfo
+}
+
#'@title infer haplotypes for one or more haploblocks
#'@description infer haplotypes for one or more haploblocks, for all individuals,
#'using FS family(s) (with parents) if present, and infer haplotypes for
@@ -3397,7 +3445,7 @@ getNewMrkdids <- function(mrkDosages, ploidy, haploblock, maxmrk) {
#'@usage inferHaplotypes(mrkDosage, indiv=NULL, ploidy, haploblock,
#'parents=NULL, FS=NULL, minfrac=c(0.1, 0.01), errfrac=0.025, DRrate=0.025,
#'maxmrk=0, dropUnused=TRUE, maxparcombs=150000, minPseg=1e-8,
-#'knownHap=integer(0), progress=TRUE, printtimes=FALSE)
+#'knownHap=integer(0), progress=TRUE, printtimes=FALSE, ahcdir)
#'@param mrkDosage matrix or data.frame. Markers are in rows, individuals in
#'columns, each cell has a marker dosage. Names of individuals are the column
#'names, marker names are the row names or (if a data.frame) in a column named
@@ -3422,15 +3470,15 @@ getNewMrkdids <- function(mrkDosages, ploidy, haploblock, maxmrk) {
#'considered to be certainly present if it must occur in at least a fraction
#'minfrac[1] of all individuals; in the final stage for the "other"
#'individuals (those that do not belong to the FS or its parents) this fraction
-#'is lowered to minfrac[2]; see also inferHaps_noFSs
+#'is lowered to minfrac[2]; see also inferHaps_noFS
#'@param errfrac the assumed fraction marker genotypes with an error (over all
#'markers in the haploblock). The errors are assumed to be uniformly distributed
#'over all except the original marker dosage combinations (mrkdids)
-#'@param DRrate the rate of double reduction per meiosis (NOT per allele!); e.g.
-#'with a DRrate of 0.04, a tetraploid parent with genotype ABCD will produce
-#'a fraction of 0.04 of DR gametes AA, BB, CC and DD (each with a frequency of
-#'0.01), and a fraction of 0.96 of the non-DR gametes AB, AC, AD, BC, BD, CD
-#'(each with a frequency of 0.16)
+#'@param DRrate default 0.025. The rate of double reduction per meiosis (NOT
+#'per allele!); e.g. with a DRrate of 0.04, a tetraploid parent with
+#'genotype ABCD will produce a fraction of 0.04 of DR gametes AA, BB, CC and DD
+#'(each with a frequency of 0.01), and a fraction of 0.96 of the non-DR gametes
+#'AB, AC, AD, BC, BD, CD (each with a frequency of 0.16)
#'@param maxmrk Haploblocks with more than maxmrk markers will be skipped.
#'Default 0: no haploblocks are skipped
#'@param dropUnused TRUE (default) if the returned matrix should only contain
@@ -3456,12 +3504,15 @@ getNewMrkdids <- function(mrkDosages, ploidy, haploblock, maxmrk) {
#'by printed messages
#'@param printtimes if TRUE, the time needed to process each haploblock is
#'printed
-#'@details inferHaplotypes uses lists that for each combination of marker
+#'@param ahcdir a single directory, or not specified.
+#'inferHaplotypes uses lists that for each combination of marker
#'dosages give all possible combinations of haplotype dosages. These lists
#'(ahclist and ahccompletelist) are loaded and saved at the directory
-#'specified in options("PolyHaplotyper_ahcdir"). If this option is not
-#'set (or invalid), they are loaded and saved at the current working
-#'directory.\cr
+#'specified by ahcdir. If no ahcdir is specified it is set to the current
+#'working directory.\cr
+#'If an ahclist or ahccompletelist for the correct ploidy is already in
+#'GlobalEnv this is used and no new list is loaded, even if ahcdir is specified.
+#'@details
#'First we consider the case where one or more FS families and their
#'parents are present in the set of samples. In that case, initially the
#'possible haplotype configurations of the parents are determined.
@@ -3474,19 +3525,27 @@ getNewMrkdids <- function(mrkDosages, ploidy, haploblock, maxmrk) {
#'needs less haplotypes, that one is chosen. If there is no clear best solution
#'still the parents and FS individuals that have the same haplotype
#'configuration over all likely solutions are assigned that configuration.\cr
-#'For FSs where no good solution is found (because of an error in the marker
+#'For FS where no good solution is found (because of an error in the marker
#'dosages of a parent, or because the correct solution was not considered) the
#'parents and individuals will be considered as unrelated material.\cr
#'If several FS families share common parents they are treated as a group,
#'and only solutions are considered that are acceptable for all families
#'in the group.\cr
#'Finally (or if no FS families are present, immediately) the other samples
-#'are haplotyped, assuming they consist of unrelated material. If FS families
+#'are haplotyped, which are considered as unrelated material. If FS families
#'have been solved the haplotypes in their parents are considered "known",
#'and known haplotypes can also be supplied (parameter knownHap). For these
#'samples we consecutively add haplotypes that must be present in a minimum
#'number of individuals, always trying to minimize the number of needed
-#'haplotypes.
+#'haplotypes.\cr
+#'InferHaplotypes uses tables that, for each combination of dosages of the
+#'markers in the haploblock, list all haplotype combinations (ahc) that
+#'result in these marker dosages. In principle inferHaplotypes uses a list
+#'(ahccompletelist) that, for a given ploidy, has all the haplotype combinations
+#'for haploblocks from 1 up to some maximum number of markers. This list can be
+#'computed with function build_ahccompletelist. If this list is not available
+#'(or is some haploblocks contain more markers than the list), the ahc for
+#'the (extra) marker
#'@return a list with for each haploblock one item that itself is a list
#'with items:\cr
#'message; if this is "" the haploblock is processed and further
@@ -3513,18 +3572,20 @@ getNewMrkdids <- function(mrkDosages, ploidy, haploblock, maxmrk) {
#' marker dosages match different possible haplotype combinations.
#'The next elements are only present if one or more FS families were
#'specified:\cr
-#'FSsFit: a logical vector with one element per FS family; TRUE if a (or
+#'FSfit: a logical vector with one element per FS family; TRUE if a (or
#' more than one) acceptable solution for the FS is found (although
#' if multiple solution are found they might not be used if unclear
#' which one is the best solution). (Even if no
#' solution was found for an FS, still its individuals may have a
#' haplotype combination assigned ignoring their pedigree)\cr
-#'FSsMessages: a character vector with one item per FS family: any
+#'FSmessages: a character vector with one item per FS family: any
#' message relating to the fitting of a model for that FS,
#' not necessarily an error\cr
-#'FSsPval: a vector of the chi-squared P-value associated with the selected
+#'FSpval: a vector of the chi-squared P-value associated with the selected
#' FS model for each FS family, or the maximum P value over all
#' models in case none was selected\cr
+#'If for new combinations of marker dosages the possible haplotype combinations
+#'have to be calculated, an ahclist file is written to ahcdir
#'@export
inferHaplotypes <- function(mrkDosage, indiv=NULL, ploidy,
haploblock, parents=NULL, FS=NULL,
@@ -3533,41 +3594,22 @@ inferHaplotypes <- function(mrkDosage, indiv=NULL, ploidy,
dropUnused=TRUE, maxparcombs=150000,
minPseg=1e-8,
knownHap=integer(0), progress=TRUE,
- printtimes=FALSE) {
+ printtimes=FALSE,
+ ahcdir) {
mrkDosage <- checkmrkDosage(mrkDosage, ploidy=ploidy, indiv=indiv)
tmp <- checkpops(parents, FS, mrkDosage)
if (tmp$message != "") stop(tmp$message)
parents <- tmp$parents; FS <- tmp$FS; rm(tmp) #a.o. all as character
haploblock <- checkHaploblock(haploblock, mrkDosage)
- if (maxmrk == 0) maxmrk <- max(sapply(haploblock, length))
- newdids <- getNewMrkdids(mrkDosage, ploidy, haploblock, maxmrk)
- newdidcounts <- sapply(newdids, length)
- totnew <- sum(newdidcounts)
- if (any(newdidcounts > 0)) {
- if (progress && (max(which(newdidcounts > 0)) + ploidy >= 12 ||
- sum(newdidcounts) > 100)) {
- cat(paste(sum(newdidcounts),
- "haplotype combinations need to be calculated;\n"))
- cat("this may take quite a long time\n")
- }
- cum <- 0
- for (m in which(newdidcounts > 0)) {
- allhap <- allHaplotypes(1:m)
- for (i in seq_along(newdids[[m]])) {
- did <- newdids[[m]][i]
- if (progress) {
- cat(paste0("\rcalculation of new mrkdids: ", cum + i, " / ", totnew))
- flush.console()
- }
- ah <- getAllHapcomb(did=did, allhap=allhap, ploidy=ploidy,
- writeFile=TRUE, progress=FALSE)
- }
- cum <- cum + newdidcounts[m]
- }
- cat("\n")
- }
-
+ nmrk <- max(sapply(haploblock, length))
+ if (maxmrk == 0) maxmrk <- nmrk else maxmrk <- min(maxmrk, nmrk)
+ if (missing(ahcdir)) ahcdir <- NULL
+ ahcinfo <- loadAllHapcombLists(ploidy=ploidy, nmrk=maxmrk, ahcdir)
+ # see if there are new mrkdids for which all haplotype combinations
+ # must be calculated:
+ newdids <- getNewMrkdids(mrkDosage, haploblock, maxmrk, ahcinfo)
+ ahcinfo <- addNewMrkdids(newdids, ahcinfo, progress)
nmrk <- sapply(haploblock, length)
result <- vector(mode="list", length(haploblock))
for (hb in seq_along(haploblock)) {
@@ -3583,7 +3625,7 @@ inferHaplotypes <- function(mrkDosage, indiv=NULL, ploidy,
hapOneBlock(mrkDosage[match(haploblock[[hb]], rownames(mrkDosage)),,
drop=FALSE],
hbname=names(haploblock)[hb],
- ploidy=ploidy,
+ ahcinfo=ahcinfo,
parents=parents, FS=FS, minfrac=minfrac,
errfrac=errfrac, DRrate=DRrate,
dropUnused=dropUnused,
@@ -3993,8 +4035,6 @@ getGameteFreqs <- function(parhac, DRrate) {
HMx[1:2, cols] <- parhac[d]
HMx[3:(ploidy/2), cols] <- combn(parhac[-d], ploidy/2-2)
}
- # dimHMx <- dim(apply(HMx, 2, sort))
- # if (is.null(dimHMx) || length(dimHMx)!=2) browser()
HMx <- apply(HMx, 2, sort) #sort each column separately
HMx <- HMx[, order(HMx[1,], HMx[2,], HMx[3,])] #reorder all columns
# HMx is now sorted, and comparable to that from getAllGametes
@@ -4552,9 +4592,9 @@ overviewByFS <- function(haploblock, parents, FS, hapresults) {
if (hapresults[[hb]]$message != "") {
messages[hb, 1] <- hapresults[[hb]]$message
} else {
- FSsok <- FSsok && "FSsFit" %in% names(hapresults[[hb]])
+ FSsok <- FSsok && "FSfit" %in% names(hapresults[[hb]])
if (FSsok)
- messages[hb, 2:(nFS+1)] <- hapresults[[hb]]$FSsMessages
+ messages[hb, 2:(nFS+1)] <- hapresults[[hb]]$FSmessages
ovw[hb, 2] <-
sum(rowSums(hapresults[[hb]]$hapdos, na.rm=TRUE) > 0) #nhap
totimputed <- 0
@@ -4588,9 +4628,9 @@ overviewByFS <- function(haploblock, parents, FS, hapresults) {
}
if (FSsok) {
# fit: was an FS solution found?
- ovw[hb, basecol+1] <- hapresults[[hb]]$FSsFit[fs]
+ ovw[hb, basecol+1] <- hapresults[[hb]]$FSfit[fs]
# P: P-value for best segregation, even if rejected:
- ovw[hb, basecol+2] <- hapresults[[hb]]$FSsPval[fs]
+ ovw[hb, basecol+2] <- hapresults[[hb]]$FSpval[fs]
}
} # for fs
#data for "rest": the non-FS, non-FS-parents indiv
@@ -4628,6 +4668,83 @@ overviewByFS <- function(haploblock, parents, FS, hapresults) {
list(ovw=ovw, messages=messages)
} # overviewByFS
+calcPedstats <- function(pedchk, indiv, haploblocks) {
+ # function used only in calcStatistics; for input and output see there
+ if (missing(indiv) || is.null(indiv))
+ indiv <- dimnames(pedchk$ped_arr)[[1]]
+ indix <- match(indiv, dimnames(pedchk$ped_arr)[[1]])
+ if (length(indix) == 0 || anyNA(indix))
+ stop("names in indiv not found in pedchk")
+ if (missing(haploblocks) || is.null(haploblocks)) {
+ haploblocks <- dimnames(pedchk$ped_arr)[[3]]
+ } else if (is.numeric(haploblocks)) {
+ haploblocks <- dimnames(pedchk$ped_arr)[[3]][haploblocks]
+ } else if (is.list(haploblocks)) {
+ haploblocks <- names(haploblocks)
+ }
+ if (!all(haploblocks %in% dimnames(pedchk$ped_arr)[[3]]))
+ stop("not all haploblocks occur in pedchk")
+ hbix <- match(haploblocks, dimnames(pedchk$ped_arr)[[3]])
+ if (length(hbix) == 0 || anyNA(hbix))
+ stop("no selected haploblocks present in pedchk")
+ pedstats <- data.frame(
+ haploblock = haploblocks,
+ stringsAsFactors = FALSE
+ )
+ pedstats$mrk <- apply(pedchk$ped_arr[indix, "mrk", hbix, drop=FALSE],
+ MARGIN=3, FUN=sum, na.rm=TRUE)
+ pedstats$imp <- apply(pedchk$ped_arr[indix, "imp", hbix, drop=FALSE],
+ MARGIN=3, FUN=sum, na.rm=TRUE)
+ pedstats$hap <- apply(pedchk$ped_arr[indix, "hap", hbix, drop=FALSE],
+ MARGIN=3, FUN=sum, na.rm=TRUE)
+ pedstats$match.NA <- apply(is.na(pedchk$ped_arr[indix, "noDR", hbix, drop=FALSE]),
+ MARGIN=3, FUN=sum, na.rm=TRUE)
+ pedstats$noDR.TRUE <- apply(pedchk$ped_arr[indix, "noDR", hbix, drop=FALSE],
+ MARGIN=3, FUN=sum, na.rm=TRUE)
+ pedstats$noDR.FALSE <- apply(!(pedchk$ped_arr[indix, "noDR", hbix, drop=FALSE]),
+ MARGIN=3, FUN=sum, na.rm=TRUE)
+ pedstats$withDR.TRUE <- apply(pedchk$ped_arr[indix, "withDR", hbix, drop=FALSE],
+ MARGIN=3, FUN=sum, na.rm=TRUE)
+ pedstats$withDR.FALSE <- apply(!(pedchk$ped_arr[indix, "withDR", hbix, drop=FALSE]),
+ MARGIN=3, FUN=sum, na.rm=TRUE)
+ pedstats
+}
+
+calcFSstats <- function(ovwFS, haploblocks) {
+ # function used only in calcStatistics; for input and output see there
+ if (missing(haploblocks) || is.null(haploblocks)) {
+ haploblocks <- rownames(ovwFS$ovw)
+ } else if (is.numeric(haploblocks)) {
+ haploblocks <- rownames(ovwFS$ovw)[haploblocks]
+ } else if (is.list(haploblocks)) {
+ haploblocks <- names(haploblocks)
+ }
+ if (!all(haploblocks %in% rownames(ovwFS$ovw)))
+ stop("not all haploblocks occur in ovwFS")
+ hbix <- match(haploblocks, rownames(ovwFS$ovw))
+ if (length(hbix) == 0 || anyNA(hbix))
+ stop("no selected haploblocks present in ovwFS")
+ startcols <- which(substr(colnames(ovwFS$ovw), 1, 6) == "parmrk")
+ FSstats <- data.frame(FSfam = 1:length(startcols))
+ FSstats$bothparmrk <- colSums(ovwFS$ovw[hbix, startcols, drop=FALSE] == 2) #na.rm should not be needed!
+ FSstats$FSFit <- colSums(ovwFS$ovw[hbix, startcols+1, drop=FALSE] == 1)
+ FSstats$mrk <- colSums(ovwFS$ovw[hbix, startcols+3, drop=FALSE]) / length(hbix)
+ FSstats$imp <- colSums(ovwFS$ovw[hbix, startcols+4, drop=FALSE]) / length(hbix)
+ FSstats$hap <- colSums(ovwFS$ovw[hbix, startcols+5, drop=FALSE]) / length(hbix)
+ startcols <- which(colnames(ovwFS$ovw) %in% c("mrkrest", "mrkall"))
+ reststats <-
+ data.frame(FSfam=c("rest", "all"), bothparmrk=c(NA, NA),
+ FSFit=c(NA, NA),
+ mrk=colSums(ovwFS$ovw[hbix, startcols, drop=FALSE]) / length(hbix),
+ imp=c(0, sum(FSstats$imp)), #only imputation in FS families
+ hap=c(sum(ovwFS$ovw[hbix, startcols[1]+1]),
+ sum(ovwFS$ovw[hbix, startcols[2]+2])) / length(hbix)
+ )
+ FSstats <- rbind(FSstats, reststats)
+ rownames(FSstats) <- NULL
+ list(FSstats=FSstats, haploblocks=haploblocks)
+}
+
#'@title calculate some statistics of the solutions of all haploblocks
#'@description calculate some statistics of the solutions of all haploblocks
#'@usage calcStatistics(pedchk, ovwFS, indiv, haploblocks)
@@ -4680,71 +4797,29 @@ overviewByFS <- function(haploblock, parents, FS, hapresults) {
#'the rows above.
#'@export
calcStatistics <- function(pedchk, ovwFS, indiv, haploblocks) {
- if (missing(ovwFS) || is.null(ovwFS) ||
- (is.vector(ovwFS) && (length(ovwFS)==0 || is.na(ovwFS[1]))))
- ovwFS <- NULL
+ if (missing(indiv)) indiv <- NULL
+ if (missing(haploblocks)) haploblocks <- NULL
+ if (missing(pedchk)) pedchk <- NULL else
+ if (!is.list(pedchk) || !identical(names(pedchk), c("ped_arr", "parents_arr")))
+ stop ("invalid pedchk")
+ if (missing(ovwFS)) ovwFS <- NULL else
+ if (!is.list(ovwFS) || !identical(names(ovwFS), c("ovw", "messages")))
+ stop("invalid ovwFS")
+ if (is.null(pedchk) && is.null(ovwFS))
+ stop("no data")
+
+ result <- list()
+ if (!is.null(pedchk))
+ result$pedstats <- calcPedstats(pedchk, indiv, haploblocks)
if (!is.null(ovwFS)) {
- ovw <- ovwFS$ovw[rownames(ovwFS$ovw) %in% dimnames(pedchk$ped_arr)[[3]],,
- drop=FALSE]
- if (!all(rownames(ovw) == dimnames(pedchk$ped_arr)[[3]]))
- stop("haploblocks in pedchk and ovwFS don't match")
+ FSstats <- calcFSstats(ovwFS, haploblocks)
+ result$FSstats <- FSstats$FSstats
}
- if (missing(indiv) || is.null(indiv))
- indiv <- dimnames(pedchk$ped_arr)[[1]]
- indix <- match(indiv, dimnames(pedchk$ped_arr)[[1]])
- if (length(indix) == 0 || anyNA(indix))
- stop("names in indiv not found in pedchk")
- if (missing(haploblocks) || is.null(haploblocks)) {
- haploblocks <- dimnames(pedchk$ped_arr)[[3]]
- } else if (is.numeric(haploblocks)) {
- haploblocks <- dimnames(pedchk$ped_arr)[[3]][haploblocks]
- } else if (is.list(haploblocks)) {
- haploblocks <- names(haploblocks)
+ if (!is.null(pedchk) && !is.null(ovwFS) &&
+ !setequal(result$pedstats$haploblock, FSstats$haploblocks)) {
+ cat("The haploblocks in pedchk and ovwFS do not match\n"); flush.console()
}
- if (!all(haploblocks %in% dimnames(pedchk$ped_arr)[[3]]))
- stop("not all haploblocks occur in pedchk and ovwFS")
- hbix <- match(haploblocks, dimnames(pedchk$ped_arr)[[3]])
- if (length(hbix) == 0 || anyNA(hbix))
- stop("no haploblocks selected")
- pedstats <- data.frame(
- haploblock = haploblocks
- )
- pedstats$mrk <- apply(pedchk$ped_arr[indix, "mrk", hbix, drop=FALSE],
- MARGIN=3, FUN=sum, na.rm=TRUE)
- pedstats$imp <- apply(pedchk$ped_arr[indix, "imp", hbix, drop=FALSE],
- MARGIN=3, FUN=sum, na.rm=TRUE)
- pedstats$hap <- apply(pedchk$ped_arr[indix, "hap", hbix, drop=FALSE],
- MARGIN=3, FUN=sum, na.rm=TRUE)
- pedstats$match.NA <- apply(is.na(pedchk$ped_arr[indix, "noDR", hbix, drop=FALSE]),
- MARGIN=3, FUN=sum, na.rm=TRUE)
- pedstats$noDR.TRUE <- apply(pedchk$ped_arr[indix, "noDR", hbix, drop=FALSE],
- MARGIN=3, FUN=sum, na.rm=TRUE)
- pedstats$noDR.FALSE <- apply(!(pedchk$ped_arr[indix, "noDR", hbix, drop=FALSE]),
- MARGIN=3, FUN=sum, na.rm=TRUE)
- pedstats$withDR.TRUE <- apply(pedchk$ped_arr[indix, "withDR", hbix, drop=FALSE],
- MARGIN=3, FUN=sum, na.rm=TRUE)
- pedstats$withDR.FALSE <- apply(!(pedchk$ped_arr[indix, "withDR", hbix, drop=FALSE]),
- MARGIN=3, FUN=sum, na.rm=TRUE)
- if (is.null(ovwFS)) return(list(pedstats=pedstats))
-
- startcols <- which(substr(colnames(ovwFS$ovw), 1, 6) == "parmrk")
- FSstats <- data.frame(FSfam = 1:length(startcols))
- FSstats$bothparmrk <- colSums(ovw[hbix, startcols, drop=FALSE] == 2) #na.rm should not be needed!
- FSstats$FSFit <- colSums(ovw[hbix, startcols+1, drop=FALSE] == 1)
- FSstats$mrk <- colSums(ovw[hbix, startcols+3, drop=FALSE]) / length(hbix)
- FSstats$imp <- colSums(ovw[hbix, startcols+4, drop=FALSE]) / length(hbix)
- FSstats$hap <- colSums(ovw[hbix, startcols+5, drop=FALSE]) / length(hbix)
- startcols <- which(colnames(ovw) %in% c("mrkrest", "mrkall"))
- reststats <- data.frame(FSfam=c("rest", "all"), bothparmrk=c(NA, NA),
- FSFit=c(NA, NA),
- mrk=colSums(ovw[hbix, startcols, drop=FALSE]) / length(hbix),
- imp=c(0, sum(FSstats$imp)), #only imputation in FS families
- hap=c(sum(ovw[hbix, startcols[1]+1]),
- sum(ovw[hbix, startcols[2]+2])) / length(hbix)
- )
- FSstats <- rbind(FSstats, reststats)
- rownames(FSstats) <- NULL
- list(pedstats=pedstats, FSstats=FSstats)
+ result
} #calcStatistics
#'@title produce a table of nr of markers vs nr of haplotypes
diff --git a/data/PolyHaplotyper_demo.RData b/data/PolyHaplotyper_demo.RData
new file mode 100644
index 0000000000000000000000000000000000000000..324cec872e681f2e61fa6edfcc8ec637c4d9bbe0
GIT binary patch
literal 15313
zcmb2|=3oE=wzswHGjvZEhP}V?PKnb%Old|+PI9E%rlWT~cYCv)PFUn9q2p7rZ`!=3
zwc=myG2d)>qBtW*xYr=z0Lx@!we)v;t(j*e6rDTyWRC8uH-7$qV}oU%y{WAF?z#JZ
zc=<D*7yDQ5{`LOVtXVtr>n=yX3g7!Q^xV!<7c-w~J-Yk2>%V2u@pn`A*w+4gnSTG<
zna}4gX8%~k|6yW(+Rqz5-?-O)`&(XfD>%RE`|bG_yEA^5)oyvGe?Hv(_vX*L;$_$V
z{#8@ECG*_*YwqXues8XREwt-j+;*#~^_H*o!`yHGTNE2#eY13J;QeK<uc^kwd=Wpp
zD=juVqvrjZ*ZSApcUPysUi>w+_I$?N>aR-|7T%7TBPsQL`s-hD^Wv|(UR-;*BwKga
zv$)vxru%Z<t$lqv{@&{Arn}1v*T-F7v-gGc<Bamkv@K88O%8WIuW!AY&vyInO=mA}
z-}j-f`un+We?s5?Ty1$`{WpdEyKL*fR{yiQ|Kaq`{HnX**2{KI48Ae%%<Z`{=Z~+h
z{4cWZV|~T?=ITFNuh!I>p7`3&z0df=Mn3h=AMeT^tc+iG%)Y7kpVjk^j_t=;|8txw
z(9^kpT5yZ>x8C)?Iq!J8MgEic_Fk>V;`^op@`pCrJ^E75^Py1Z=?C3^mTWrDuO(l8
zG+%fJ-|N+1s{=o@f7oUHhu?7D?e+gxbJ*RgJoRtxuJVKO*LJL6eiZ$+`(c7~1<(5L
zd8?1_pI`3t_xbPq#rx&{b5FZ|KZE^>eCgrO5$tP}&z-)uH!A$~jg8u~|NB4qQRMe>
z^NqU?&e(0wF}nU=^B?c8`>k)~UuUnCf0h5SzUU8c)7DG(FI-(E{Al^H`<xFfOa1z}
z<JGhOoqiboYxynvWA<WUJ(Gie{kDp#%KGhouRC$+oAZ{t|9;<SpZ(|gFZ=EPZY-)l
zu)pB0=|2HE?}yQk?@K=4@^9(m^-KOH|5<(8b^HE%Nju#y@a_C4FP*-8;^U1!*>ewl
z`oBt1=kTTY$I-uzADf>l9Q{`RLH&+pF8?QY)LH(E-q5!7pV<E&xt#xXzL|SSXZ`N!
zv%hfs{_IcoUbo_3FSzaSTYlq5@z?)X{1-mB{YU!C*^&Q~{>fjvKIO-L%m41*yVp;?
zdi~&!*{>&m-e2(Vpw+QI`>*_&F7-e9c6|SSng3DW!cYCb_QU_@f6M>3v;N=sGyhTi
z>c4Y0m)e$oo3)4K{w3pM+(&<D*F4Yp?f-SZ>A%0<Zmj*d_S1K{|L5L`i(kAS^`GUR
z`kebl|F-{_e{F`wM|-b7{zv+M+FzTq;r)-xulr@{IR8(5VEtI;(Qm8N=kl%dKkfJV
zd-s;!q(A02o<6@ef9b#F&Dlwh?^{iqRbTY)y4c_1$M(nSrQ}a{KRli!yYuJqUCT@U
zWd8F$8o%wwf5X3~|6f%H|C@jHUznuW|ET9%{-s4X|GOb)`QLoozwWQ~UH`xQnci*B
z`e*lQ<rf<^{a3uB_m=<rl1J|i|9bzmXaBwQ*Yr>RXZ~;ObN_u`!+!g}%fFJPUw&=Z
zvwm71FE5?9E&5q}&A;i7_KQvn{<~k>{o(XA`(vKWY0ckSpY?b8SNrx2y&r6@op#-@
z!8h!2e#G>Jx%y`R&i<3XmYR9*>+!tJ|2lrkZ~S-P>i?5J@j3q|eRb!OKizk0|E}`x
zztM9-jQ$0B{pp|bf9|ILLI3&BN50*ENIvge^pD?X_OnOVpZTg)`7d<Szl;CkPkz>y
zsW<%l{p0=(JV(5a&JQ%t?s@(}e9z85LjSqHp0}$nWRLr=@z?vO`{VU*@(UWzZ=S+>
z{M%Q1uRrsh*X`GxcR5#h^Zz%$u0Oqf|5E(p>96{qZ|IrxZ|R@lL;Ej&lYjJM<)3DW
zx^w^c+L<N4yZvutp`i5t_6PT;{0}txpYiYVxA>;@U;EG6YyQ6+_224o+};06Ki6Nl
z7WJQ_Zsq^ww~tBH|Gsr!?7#c({g%JC|60=WzvTDhy@r38*F>HV{lDAvzv|C&j>1K6
z?c>vX_uu|i-;w<B5$A*KzlA>QuGvfene;RL;{BWd!gGA`D*r7#KEL(7*N26^ULQ)W
zqkla6#a__#t^dJdzI|^*{zv?M|K7LP{`CLuv-VQQ=Uey5@IRAZ`rq#7@{1qsr|%8_
zfB5hAUi)Kj>gRurKlg9jnf<Fbn%(%7d$Qi^<@}@N2LGL7zdikb<&W@<|2<pJgHpWZ
z|Jz;rm;L#F<Im)-{l~7why7;%ll(&d_5YB6kv#iL{$2mK|Je0k)}`m9>gDP$|B%o5
zbKu`q)0JQQ^S9sl(S3IR&Y$OZGX0)i$T#zUfYE>VxA$A$ia)l$qT9Cq?fob76?M3e
z-|tSW5Bq=o*!*vtUVr9`{Ezx0e);OE=#~FdAHM!5{z3mm$Je{*|D1K+{$KtkKl@*5
z=D*s$`t`3TUaa^2z5G{w_V4$?|21RlSs&{+_aD!DTrYn1e?jf3|GU@z<6F1!$NDAz
zb^pxx`}$-5_i2w7|9!vM{+9f!{H^up{&{`AU-WwajQ_Dl|08}Ke>Okn|MVB}m;PD)
zx&0yk<e&MU_6PsXf0=aFbGvN)w7=02JHP42{*Tx<(>`zO#@6@!Q*Xrc*r)wJC;Dsu
zx}N%?f6G(<|9)nF^w|E^<g?tG|BJ8wE0*n^|9^FM(&O{<^*i@Z*|vY_{%7;ECR81{
zWo7z*>a+hxfAGJxKVBd4Z{N4&$6NQSoBm&R?7rPU`M&*g{s$iVAN_3ql>ajy|6luW
zzv+LMkLhn#{0rC~^FO<0ZGFz4=x_bCo8JF&zg@rOzvJ<#|Nh;mIlJ)s{pdgM=lo6X
zDZToS>+pN`ukvrQ^MwER^!~rQ?a5DX&i`+u|65)8EWgF{^nd<0^&O9=@2Q^~T_5!?
z-sgYqz5g5E#6M2|yMFGxwsrESzOomcZOXs$((Ld2C-U33nRd_q_R#)F|IPia^Q-?S
z$p1es^S?QEe(3-6*JWS*z3-6keZA+O;eYMx|Cc`I|7z9pf8wV9cQ?NI9bR1hX7>LT
zZ}VsVe_rJM=)UUv`LD{({hWW!^;SK{?Gst`-v1@`&HHb<Z|gtq56^$yQTvy5v!17Z
z{>SOB{?GZ>exhFc<Nt4x^$~S>|IYrDpZ0G>HdlT2KjUi`zu5N}ZP@pF-~TeE>&Ktf
z*ZX(=e|uV^^52rT`&0f3ZmwV5x4rhC;Nk0kqLb=(ovo2RdaYip?#lm88S5**uXp61
z-8wCIO3(k8f4*n_Yt5|BpIpE6r}>5ZH$Shxs<!>){$oX7_g;(Ap88+-^Zx1I<a7R(
z?wxb`e0|;Q*#En>{uiHJ@BPO9`hU;Au~q*r|J{H7^Ze<D|8zgiy=wpcU;nY|Z~Lp<
zZr4w^zx^oxkMxMg(O*A(^;m1a_MiXW|ACkOAN<1{;dgibCHt5EHU55;>Hod)WB%0o
z*?-z+{11NeU#sT%wypc0J#DBDx&7i}f9jeq@`tL6{?^Cdt2u9UT;hl7`}yjB*Z-_q
zy~p9dW1afP*nNMQZB6aw9Dl6&Km2_C()06!e?0pUp1b(Bz1e@Yk9Vv7-9BQU_y31C
zkNvfOhW~cv{<$dn;(g4W|9-arwJQI`tlFPi-+oE{^|q=1o2+=vZ$J9+Ss<>t^gz|y
ze`o%g_vP09zkR6xWqrq=>bL*J&+peW{LlMvW9@x=&!6?nHEYsu&kuDwbMD*L|Fs{2
zzs=A3zxFYIY4^|nepAJ_)`uV1&-%mi-_1?`PCY#SXZ@8Q{a^QY+8;38{r6|N;eV|^
zSALz|wEyZ(rkcCw#QsnJbpF+hL;WA?^Zp7S){m@z_9(ycfBu8Te_wyQFIBHw`JeH&
z$L;yw?D?nGOMN(>F#mAb%};+{uKmdV{Mxtlt@ST|*q?o{{+a#K{YUFJe)T`2w(ZRS
zUH^l>ZQo=XyZq<(Yd^1NNgpX+cjy0wKiOL&d+LMU&i}mMr{=NT-Ty1<uly6e7BsD~
z|66-%ba9MN{n>iC<kz3~Px{0B_y7JQ&;AGfN%*(@md?WZbN_Ei{a^3(;i1&se+T!v
zMMf;w{Udxsa*Fin`Xa&DwcqZ~`u|X=)cH@&zqwh9AN{z#;7fb%g8wzXC%>;Re!Kp4
z!pnH8|Ia_=Z+P?KSNO~MkM2+XpZ{R`kNukeHvYPw&%d<(!G8zIi2Xu;?)(-v{OA8u
zUhMz0AM;lnUiRnyLA&-}?tkAa?nygZ&$jQF=i2>S{?FXBU-wUN*ZTLqKkb+PTK^|s
zgZZ}lh4`0Q50fR2o|iZN_598Gtl6&r9Zz{Yv@SjOWd6=IwV(9EbC-SD@`wNQzv9RD
zqy8WJbN<S&(_8i*TzcKUVBM3S!RP+>{Irk!e_z(@pTq;{h@C&#Py4Oi->>vE^Y`{Q
z`YRTnzJ34Hxwa27+q}c-m%pjs_ss0Qz01$z7v#UsZH)i&=jes~H~uO85<MQze$m+U
zU#abeHB0{IU&%N9$NN|JkHr6|r_blrvw!3N@%xATu3UlrDStOU{k>oHU%KZ%?{ojV
zZqC<R{p7#>#t+5c^0)0@nEmc=`;A}9pn|c!^MB}~`piGR(_<g~U;HnkEabo99`jrO
zsvn4dt#|q#^E-JiM_=x#{{g>_U$fu%Q-ABfd%ylAp803<@7nkJjX9@}pRujK_#^uB
z{h%Mwx9;EHeg1d2<-e=H+<%I{>ip_^biUUY{mB0o;xFrWn4A57X*K=D{;dC!%a{Bw
ze!YLnzl}eayVSGq&lj%0f5?94*ZymY*Z&#*`+nQ^<Nw&$lHU@g8lC@Lf8W2qdfV(z
z^&)kr{@?laCx8CG{)f?jrhm5I@MgvT*$2*_t8du<{#Um5@92N^UG>TTzOUS6zx`kD
zBmP&jetuW`XZk8^*T2&r<4^rt`|tQ?srxbY4gWTNy+7@n|IhgwzxrqWZBKphQ{=zG
z|JPs7Uwa;}`Y-<EpWsLJ5&!z->eUN+pULO^e)>b<s|6bm`}gqO*?S~@82;7$Z@+Z^
z$@*g(*Z!NHSDzzZZ}aVsrsV%7GtW=@F~X1MPx{(^^*?*<|2r2<kMB49cXI82pMS<5
zihnJCcK^no^l$sS*VhOA+`jQ=e)j+0Yc2m>*?H-!__z3rzYbrpKW^{(alWF?`A17{
z{MddizWu{~rr)2xHdihCQ2Iapr+(x6D|@G{J7`)gYreOhvuf$B`$yWVlfR{JmF)Yv
z^#AkRoWG&J{wRL?V{Z1}s`Q-Hf6f2N^S`Dq*dMrk#{Z?C|9k(j?|pyT^XLE7-{#-`
zIN|;Lr2mn{lY_q>`(eIxzuw_j_1*uy{^q}qSoweN1OA75pZ&NU%Y1$6tN2I%Ev8?;
zZ+ZKl*N3%M%YPjIQvWQQPqhB_t$+V&PX79T<>^oJXZ1(k*r(Q8{=R?wfA39j)mr}4
z>E_aZ>8Ji%y8OG%b^k(rt9{GYEgv@a-2Tvi_20Rt*X^U9&A;@Azczo`|7YKB?D|)(
zxA=cdp68#PF8>$IU)ZzcFaPUJo9^E?KbKt?G0#fw@Beb%zp{VqzI>i9^nZ8Bzl~qc
zud`SFwf@zTFZy5iZ`@P+^>}jhr~R={VixSb74>7<fBDb*zx>$zOYK6vUis_)^)o;H
z*Rf6h`+j@-hW#=B-Jjn-cujTc|Hnc<!hXLu{{L*=|H#t5v->^&8vpnw@hJUsJezd<
z?RcI4>r(#RH~A;xHtj#}U+;VNr+fd;_#3@1E$`p`r2iFiU;Cfke{FU4dhXf!;=H?N
z9rY)Esz0q?^5?&>bNy<6uSfA!e;!*N`F`B{ntl5P`)mJKem#EekGuE(Nng*~{mVc3
z=lPvox!Zd7-}#gMQoj1XNOFRd{j&Z4>_0O7XYQN*XZ5QY(f{9k{$qXp|2d0rtN&|r
ziq<{Py=UKd&c5xRq2I2vzw0%BpWe3r()~l%i)MeE|G7Tsul(csEAhYQNB#@{vGn}@
zm?sr^|2F@4zJ;NF^WU9+{v`jjk1)F@{QCcuU&`kHA8piqKR<ute&;pzLEpq5-xvC=
zef0cU`{?8M7e4>L_Obr=zpiWduXE<QzjH&*>1*3=-H-cqpQ&E^|M5TbCF)N7+t@SX
zpWofxtoC8cbw8?i*Yj=tck|Z&?GN|sK3>1&zc;AWnLqb`(O>zi^&I`Vr~dbE;eECK
zb>hqXhtpr=Z{2@P|6{$*pSAzahyB}pZt9KtGw1tn&M&U~|M=I9g3Wiu<A2w8|0~|P
z^7;LxKccq%kiYr$xX#h{_7cBlzqNn8_5A043A=6ed4B`n{FnZ5ZU_5D0rj$)_&-9w
zwPpWr+-YX|FX&(YQTyoU`|tdTzw~eE=lC<%bbsIa^>1h1?Al-VFV)ZgDgM@8;(znq
zY0+=bZ~A=WPrRe_Jd4@?w>`6ee(e9XZ~V8)*S){??f9pXy<X4dSO1+``0v86`&R!i
zALu_@ulhsv_5T@vk00Bg$G-fX|IPY!|7GV`mo5Ca^QU>{|K_!a?fz{qF8qC4{lC)l
z`dOdfui^i6;gs6*;@S7=14^C$aepwKc02#i^Edn-r~jLO`tGkkhcCsSUU>ii{pUZG
zUkT^duXtYcHzB@Q`rpP^p_l$wtWo(V8~uODkM!5S^9_G@znnR%-sjJWjkQ1ZP5%|$
zk1xL<7yZaS`r3T$b%w7(X4Qu`t^c`S<zK_w_c8bWeVDc4U-8f1HT5>j>d*Y&`eE+v
z{J)7W=S%(H`SrNsJp0Jk_OE3giGSYTxV}{5p#Hr7ZKXd|AKM?bH~iQA^}pTE@Tx!4
zf9;?4zx?*s@&*4DcCr2V{BiC5`Xfd`KYqWOAN_yovHkps|Igk0HGNHe$=#%H`7h_s
z`hW7m{YUaP%PW7fY*T&gzkR>yzePXKud=cIKYfk;wte+49>{EKyk`IWzvYg!%YV32
zkGx6h{-^qQ|CayLtDjp<TeIJ5|No?~|HJ+k|4E-Z@sQ`u|0@^nPn~V?FZo~oyPw@J
z>!p5#-F^b9&DQS!9`)b4bN<pV_Cep5eb~P4zxVb15pgAd(!cCK`^(t;|Igf$e>Ru?
zbN(|Ub_uAzD}K2?`=9TszrjEGm%cu~=D*?Vsk8oj{rx}X|HRk-Ywy)xNsV}{zy04O
zE2U}I^u_=Cto@(#*uHIk^?z;C|I=RVe{kL+djI~gPk+uw{m=YeE*Kkbc(Stk_V@iI
z`(6IpKfaO4_|NWd{ICD}r~Z&%`ThNte=q*On_R#Cb<WoMygy5C-Oswe_V4`c`ucrO
ze(Ha#=lmZHst-0FpI@4z{ePj|wEs`Hw|Ci3`87Xco8>i=e}8^&`*~jb-^|m#52ydS
zf9v0tf7-w0uZW#C`rtip|I%yr8$V9pTwlND$#3hbKZ%uw|4;YNp7Zzass8~_Usikn
zSwHEc_~-q>kM+&2FSsB1>EGJcdg*`OAM<(Y<$s?(Yajh=|D6Ai4Oi{7dA{H;`wIOx
zzjI&y4qy4RocsUtYw?CPoBy*Py!0dN|MNfV*FKiNeeB@Jgg^Wbzy3ZJzhVC7*Zq>4
zpZr}v>Ery_|2M`h{l`9S|M%a^H@N<p^v8aGZQ#EB#_`j)?q~X2{5Jn%e{tczo6mo`
zpZkCP@B8<$XYCi5Z#z8g>-=Ne-zXj3ukH8c$LEa_c8_+P_;}}+{_nrtOn;nv_e<2}
z{k!;U`Jee$KCAD)yX)WPXYt<umHub{TYLMz#QR3G?ydhft_t6(U9R_=<=fxC$1Cm~
zIBwYTH{;>&ZJZw>HJbLn;5}aSSpMT*`J(m5kAGVpS<ts{-Tq5=KkW~w=f3}?enHH8
z_R!CN_vv0R`)j%S@5cXYzr~+xl-n1l{r9fbx^=fxR{u@BJK?{0O$PTB(*y0FcmDhS
zIHG^E#=hw9`{lc(<iCEezxB^_-+r}>`YYj6_g}idJ#gW_mEX2ktUUFip6AE3e~o_+
zmTvs`e7*ZUSLw_5L%X)_&#Y(LXZruakM+;0YxXDpj-Sq7sxA5d#{Ux&Z~SNc-%!`M
zZUJ{uLD3)Wzsxn_-?VT2H@p*maNm6;_1vfR-1jYCrXPua^*H4JVYAG7J;Q%78TJ3>
zDXYuKNY_{B9IqGq|KQfC?c$&Q9s0jJ``PPhuRlhx-M=LNlmBo1>i=JV#6QzmQ9Y@C
z!Mu$hmmjx(wm;&Zdtbfxt@@+;z3QLHuU~!K-TK5&`;9;4Kc4;g{IzX|A7i|O^pj`u
zr~iAs;!lnKaby4Qp8re#EB}f=Z(Fb~XvcBAEe}BrIR5AL+5arte*bwNyyN<H&foX{
zPd&3gutsBh(EjGH=37oq`{r(IyX~3$mwJ<F_oIGI|549(ZegyX-LL=AA2#On3O=rn
zcenjl6DwA)($xR<{qz5G{>y%N|Degfp7+=Ps(;;I_y0fr)86v$-n{Pr%&ANMyl;MJ
zy}@35W<Bp8iT}U0b)T)j_Gf?3f5(!jWA-oR%hpT&E&R9j*nLi!<8|AAtWWu;8&`jJ
z#_F%@^QA53&90yPE&sv$f78D1=RO{v^;h_`{iZ+pZ)Q*ZZz``|+WYPL*FX6$*Z*BF
z`k#N%?x+7*eiY<q{oUUBfA8PNdt^TBkM6gh{GToQ|BSy+|C~SYZ0f%o?)qJScYpn}
zf0Mn%wC?(6*Y4l=*ZtSJqqkCaf$o~G-|c>M|2Lfe{Pm_E^H=<T&}0!h@vpP(#eXJs
zJ^yw8b@tj%{;&K)?0@v1<e&08KZYCqfBi51;qD3b%>SdWuDtPIp2hw4{=EOK+m(M!
zKNdgv^L#Dy1N|?m?urTiHD7Q4yg2N)^7FNi>?Qwt{|G;CfACkjo!saBg+4i({yRQ&
zmb_8c_<!YxufNJC|6OTy?3TUw|M(yEcKcfYdwu1<X0P|<|AznGfBd!oEjiv6SFiGu
zeS_|A?T`Bd>uvR4{9GUOBVxJhe?=YXukP{p9e18MSh4<j*mI?AEI*>^xA-54cxI{o
zw|dp7*q1Ms*ng?;J^Fvk-=D8J3Owdsf7>6k`Z)i(`Y%5&|H>|0^-j$0+GCfgTi4k0
z{;j=f@vP*3_W$_LO|9#{`k$>2`@eW@WlfCmtAF1s{HC@(5Bej%x&Da%&H3-@>+eL>
z=UlzNv-J1+WAX`#{0Vc#>TWgJbMS53pZPza?!y1wkK=f1>d)-lP-FjZ?YHZ{9+z?d
zx7fa|-uB<#dAt4##`kYjF1<VL-*ubaNB^xebU5GeU+LHKAJ_X+ZT|)RU4Oe?`@iD8
z{QCX-J{&Lo!+w1Kw*R3%*Pb4?f0WQ;aP$A?{MCQ^%cGtD-}%%0IR3l+z12VcKj-t@
zpT0KNY{UMa`$Ych|M<LhR(;ZE|D1p4d;h2W-d*}%-u;8FeWmGt(HHVJ?59Xx;q06K
ziT{@6yz4)9T>Z9x%|xD>)qj)!EdCn*(ck6&(S1rkqS(3Jf7=`W-6}8BbCUO0vb9C}
zdHcZs{x9Dz@caEkUhMy@U*~t+yEd_YgZ#D+{x5$x^VIjh`#<CV+)ws`@we;G{NJ8r
zY4|_A`)~fo|IB~W@67l=Rr>$5|IvH?G5GX!PWzOf^QW4pcl}?7HNv-2&;8r;<EOp(
z|3q7ze;*&8-~LeS|LOSTpHDlgd9Ei$>y>uBsPFjd|KR^ZJ>xXJe*)S6XZ)YM$v*eA
z=D+fDdG$y1CI8R(zWRCn^V^C|>%VTlwm<rr{95m>|8q7U-=CMB{a1c=z3!jamMV+t
z&;9>&eaW@>UH`xSV!!QpJHGG#+|9r1pVpuI75VyowzT1w`$qq#9+N-w^?IK4!R=qJ
z|9KuW>wnPa^N%LGe+Yi9e}CGvzd!fi{xHA$LGthQkMf)EJU@N7YR$9zKWEgNZD0C7
z_|NLE{oDS#{Z0O%`g;G2e;YrxTdnyb|74F`{mvi!j~;JcIq8Rc&OhyK^(&5l-Ff%@
z&l#rwqkqkR^*i>~|Gsne0n@p@ZC_Xa<&MPc`cm71|J(k>m;E)&U;J-<^!EC`|7J($
z$K3gUuxd(Az5cxVt`D_8og>1|{r5au&;5UuQS!I;pEok>pPc#`|9XGwpM~oe|5R^_
z7y7U9X#UUbt@6?TEt%!tt~dE7^FjZ?gIOlWKtr@r^|$^XUi-}c{O;&C|6T6ZKe@Ii
z<W>EH8yXkoKL-DK|G&h~N#|dP>lgc$zsc4tANwUgTP&*YtLOZpv+=*izT0B;j(_jE
z{ayQSec4xuC%bh2zJI&t$UmEX&9C>|`)yNytZ?(u|F;j>YySOPX!~!{XZZ&|tlqw=
z<2jz6{a;G{wPRjxee~o19r?EP_mk?+os@|X`mY&Tf9z5Gbt}Fio9X{UAJzZ({jJNs
z>Q3~7`0p1t*yaDK|6J-9TmC5LfARO}-St0yF2C~sLaOc@iHr88{|f%DecfL5$V{r<
z<ilD^nUCVXXTP^?{3ICvx*=8fygm2V>mT*I>wW(-8wk{A|6%;<eorr|{_J(O$M>(F
zJaBwh(%i`(*|*EpJN`K?`G-GuvUI)n|N03<g7Me?R&dnMUzL1S?#HgVUw1wEwc}Xp
z{p-mmUgte{eCyihmSg3|EdTzAUvDiIZFcL)k)yxmkMCV-e`~?M>oxzLUHaF4?S5AM
z5yRW}r}I92{ZoJ0e!2RM-^wd$3J%V{`rh&1Z-KbkQuUgDE*~kp+h+CupTw($Z{?rQ
z7vBEu>fiWPaoopOR;@YK*0}!sYqNrL+jx(x50Sj}Z>#^DSM}V-Kk*9Q-@bjFdyjO`
zTlv?|W&X4Nd|t5DRQ=d3`=$Hs#2@gVlc|s1zrXy^qrd$>#QsM{e_>7j^}2dZX<uA@
z%JKMx|5X3A|NmUI{+GNVoA?)J&i@(zjklib{%HSf?(;ADng3P(XVvtszLgLud~oBG
zzq-FvK5YH_+GyIQv@=&e_h0?d{@Ucf%waRo+~l?o>(^DEI^XljJ}L3<#~BOjbN+38
zqkr|^>au_;`L6$Gf8M_k{ks20)t}c@%kJ5qd0D@9d*T1wU-x&G-dA}pukkq~Hh57z
z=l_Vu`46j>)<1pp_I>H)Qa8Rs*^3_Y7c$y^|7-j_=U?x`=^xJL-oF1yep>OW7;kgi
zYqry){zOaus+PX-=zrqB?xXc9Kc@eTmy-W#{#Ran=Ih?%?tjiJX3sNTQ|5QBJ}$h}
zo;OjF{djr8|F}P;AF5m5Kl&q|{a?rKU)a+{d^?}Vi`51Fm#jW&zhVEO<JD9CU9-^s
zRv%*#x9G<@oqvs0ZrA@8%Nkr~{Aq6bfBKsGFF&sTslIFH{Nr}~`s*nx>U-)-ev3bp
zj<Eb!c+RfB^h2`C=kxDo^uLtf`9X4e-~Bm<O#6lHj)l$TzI|uYW{Lm5o=cyOd3^2t
z{zZKK-1|)bZCUas|7>oR>^7O=g}nP${8+VSy6wLYj~?Ye{6FQtbivv~_c#6f8SuzD
zeNE&-zIo|3*Aw=xdEW16yWrpTJ3L>#*6ru1d48zuV@&_|teP+KNw5D){I7_eUM6?w
z^Q`G6|AOKV9Q}6I_W!}H9sdto9sP0si)vuZuD^f(Y`Ap4wKnM?dl^H7_iTwDQT6?e
z+uF`+6xhWVvQO!cs9ODT^K1VDJs*{>X=c>#sBc~0_uHuH>-NI@tUuZ7y0@NO_I&<l
z`IW^7f6P9;?MM3ea+@O`xBff$OZG!)jL`e9pKH(j(^zCbt+({pw=2JXS6S?EfBxdd
zkJKsu*(~>;@n7+O&;B3A&D?KZa;&o~dr<EnsrsY+!T*pyD^30@Z2$N5?Rf8N($
z-eA`zGvkwe(#Q9~s(I`7Z<l@Joc&L`_Fvw^<myGei~imCxAez~+56u~tc(64|76e6
z|J?r`9NH;Vmpn82ll|6zhJQk3j<s%^amZcz(?|Q&b4~Zz_q|@9@W(m&$fJDKU)H5Q
z*H7mco-K~`&k^R9f4%%QUuxvH=Q20;Yu$e|zxFZz4fkt%N>-iXxK(xU(VYpK>fhXF
zO&6$HGh42HUVZWX`F3nI_WKq5*FAb4FfH%@f_Kfocz+yU<)``kdbR-j^#!$h+g6uc
z&5hKo(f+sa{l+iXi=BI{*XKsWZojeO(*Be`ihPfvzeoM+e=2<J!}{Oy@2>yber~x#
z%$)jw>!;S{8a=O%er#|0Pv*tmP5&$HTi0KYmpk4c-SfO+p7S5&BjtH-|BJp_|IzgP
z@mohNAJiZC^tbija?$@6{w?qOob`9@SN><)v;JMq4chh3_hbEo`sQ=X{;%vidd>dc
zcg`>Ezc=K6xOd?Hmx$&owHNE{#4g|8R3D>!^vn0#5q~6WkHq|cY|-~|{ZjjLvR4_N
zKdYS+|Mi!7M7@SiJ+EE$);Ii*&+qzY`9Jvf{Y(FySEO(HfA^Q|pE-Yz|MAcMlUU32
zueoUbdESNVBo=SKb^p-zKX0bS3mj&T-FVmg_Wh{z=#MPfZ}V)kqvoAG__6kxe2(;)
z?&H<lkJ;DUyp^{0`Td<Am(Q<f{k-*h?92S4@muCU{n&rCQeHgzS6;*01${ene_ofm
zzw^SsH6PbskFWX@{bL&c)<^Xg>3#M$)+I$cSM<u&Z~U?O>iJ!Nb^k@%{=fP5{C?fP
z8*`3C6tFGb&;8bDh4rgf-)-%#9G+Eg^8d)Ocpbe*r(=^s&auATuli5m$o>aXUrsXp
z=l@#v@!$PT|M=!@Zs1$_G5?%t!5X6vb2tAy{wV(J{~P5m%2oE)zd7`;|7ZOV9l6&X
z#lo-Td;h;S+xfG$=5XxscT;y4``@pBb@F-T^Tj*o@SkD#pTFz>9dVnVbN+w5mOiiU
zVeI=E^|JQ&<6?Av*H&{LKmW8s`S|&H843Ga|MGpS?*H~DvR3fULvP#4f1loTZhpGy
z*`Ir#UjGW|wV(CCs(w{f+_zb&>tny|N?jlS?cDUe-Kpob-OtaTdooXVcki2Rn{J=X
z)8GAkQ?=r5roj7?V;xtst`*x7bs-DE$Yxn9xW(;)lp`u5>_U`dHtSl^EhvoFENg|g
z_+7AaL}kQXsB(PGx>kG(3d34_OR*z1Vz<B+Zbz(y$QJL$9U>hl1k0pd<^12yC+u&l
zU9hy;A?16+iFaQ&)bsro-@v~ky`h>7jfhBZ2xSiAjWA8P3TG@x=8fP@Fm1TX5XOwe
zc+oN0mwSiqgKZ7Bnevbr20!PPHr!^+W8WeBpbUk<zr*&yw}#uyc_@r|%<n88U?buf
z-myMFCiE=y(jVNJ6U*>U^+6mHCc&hC^?<i&K_Ks+VCMf<4%nN1@W;uV&!|@r$Qu*F
zEPnNXD@adF2ux2bvlxaHlB%J>G(F(T`@@6xLlCoF5VOox26xi~fxHDFBr*R5G1pvW
zh&MeD&if%0$|TtkGVR`Tg+br+fH%kk5T+|{K_GL?6$WvT&q!l($6R3&-+I7Rw;+%z
z%%=V;O#inWu-E<I&;92zgjwLnJ?Aph^DPG~bwBuU+gxTkzvY0j?gyW4N6!Z{#VVGu
z&dY2}k3wbcSjw80*%%%5V44_$xy@bYg5aO|7+Mbg&ulD?dN4=KVma$MbSAO^NM?bV
zFB_t*9!%rgv6MNFI?NBnRv(J39?aphSk8RzWdnkF;Q7mjdaDQX`F<>C{`V5Z>xT1?
zBruqN=J4%U!L0YHq1)=gA}AAK4u*2D>DWaNJO^7l-Ri+ozKWI1ajzOMne1_|8mCu1
zSSnX5ClPP^OWJ^~x<}QnknQ)rGqc%?(huryQPkXiT!Q^qitf4fEBv{C2k;f7^2gj{
z4*!<ELV61;jB{YB?T2LkKR1}|za8kuWU7C=A%7$NK&Nd%691eV%<A9Z%<o?hl-qs@
z<gdBj{Qm0p8`mS$UNF^MXTA>-{t(L#VHVmJ#PILA#=QRPnaEcszeV4#{r~gl$!|p;
zLi2b3KlImp-`}0T{cTt3)jushzvp}Jzr-K!IBQKm*vtQU!2X{zf1hKl&w0C_X8*pM
zTECcmdbNJ}JJIiRs=uaNeK)MXxwF5P^|bYz-#g|XFxvRX__^J`@7CWrw-(R;(H-$`
zM)CT;Un2iq=&iYAJ-1r+{`Xz?kAM1luGZ-O?_c@nKW#m?w=w&3+x>@kmjCg3ZnrgY
z=lM^eI{VM<eCB&@YyHKY?;g*s{;Y2OJ<@;UpUmfW&y(j=t55&<Q@-S#<~_T=@^h=7
z<y(D!Q(pYe_WtFa`hVQc?R|gV`uo52#qYfB=C_{yv1$9Xli_>Hb>n}|&VRCc-_PC}
z)7M{2_xHAbt^aU1vX13+`o=#V)8vnSx~XG7>*?_ge=MfSpZ~;i`uU~}b^LbytsCVE
zqU|62tNKtTmbaKa{=g@b)5bdIH|=SEYAN<^7rXw34eR)K*zxb1@U%OkPUy6MTBZEy
z-;w)fwBG)F&_DHO<~iFxUen^we#(xh)8*d&+Bxrf`~G=P%{TmUo3=jfr?1YM`-krz
z`V<`*vAuo&+}7VuAHGld8LP8@_S5eh{=|vBJKlchIsg8)r~DiKIEv*>Zm*y5w12~&
zK(TGj?e%k>`dfVO@cj{ZzUW-h!J?C&Vk39#m#ZzG_Wt0Zzu#U|Z>-6m_Wtmv+DN<J
z*0(P%eoz0&Gfn=$C+5gHfz$kHKNWTCr?lpNzQ~_kY0UlosPemR_dDYKJM`rW<87+A
zxxXJ@d_M7^UFDbLxw{`Fhu6I6>W}@narIrpcdNDWpl$!P->qL=zHvisY~lXcy|Y{2
ze(w7H4u;-#{l<b`E#FXZKeiAR{Vka`|KKO%$T~&t?^i{?Uk0PYqTBzAZvS|7eq!Z;
zPaRL~H~dL|eoR?R?waz(4;FI%5q^;sF+DbX?sBglua4Mf$ld?(<Kl>YlHA6pe!Pmj
zaqjr_GoN%L_NkuUp7isYj?KSg`G-H{M(&g4?tl8xdP9xQ>GG7HJ0dIY^vyr;X>Nqw
z#HU+z<{$ZFB_?<Nad}#$?P=)<yVj?{k@rq~ijA02e9U-SmDsbkjj4&zQa7i$ZCg4q
zH!}r4^Jhv*wbXuj?%P*>NN(N`+nW=veEaH&GSf{6<|hgHo7>#Bt(};gopM_W!hD;N
zVl8#E&u!b{i3sM$;`E*N+_x{E_-3-nwl~LL89Vdy!4qbtn|yn7f|WlXgjx02KP3@j
z`09yovr=|T-CXCk4U@UQr8PJEPo$3h?5EE+{E3-1-~9CX^q)b~^pk(G>f|&lfByZW
zJF<@Rw0`Qtor!-^b>tI%KGlgo_(?lr#_8knM-TlmHu~Rl#$NgJ`4dP?Gh3toJvrBv
zZ+}G=MN)9b*ywvt&Tr-7^C!T}Gq!N1jh@@Hc@uZ0|Cy{)-_?5j<d5{oiZwm_k#c8?
zjv8$;nBHIC*ZREoqyB~ucT7I*G5n`0Ht+kf{f7>HKKUa&@}Jb{`1GIBIyUcn>bsxj
zZ~h}YZU5;{-y{CnowiT-*{<`y<Ei}S4^yoF1&hryK3;$BlYYcM)zkZveqPtHG4K1|
z_VoUyKi1Rg&mT&zdh#>(kH60Djep}K{%N1CPyYE{=YR9l`b~fIr~N<ssXpSL{^|e8
z4`0?kv78iIEL?o($Q%Phb)PewNX*OH>4|YOW{5A#H{$+$_{bjvLwlby{KCbDkK8dZ
z)b}~VE&Tb=kv({r$c7=Ai<?<vU?}f%hF$pc!6T$IS%r%a9*Hq96!$qpCX-pX_`sLS
z%?9i~XBdTx4;+y(FeHs>v)|_ov+(BwM{EoX`H5nZVoU%3@{bV}{NJA!+myvU`1mDG
z?|*H4m9qK1r{_fPo;=sT=jOg&#t&~i`Tq3AIsU2Hc2_oiPCR^n{r_9M-|O~@S#JLv
zdjEgv!>{-5*T@|G{Y`&sMRnDjCrk4t2i@HM-urq`d67B<wNKrezqOnDT~fa_XYJ|Z
zo@HOB{JmaRwWIc~_JQ!4*3$YtcVf>Ni&$>&p7`#_aYOm?$VbQZH%4rHFg+^nm6_($
z+TZ_{=H0j6IT?h0o|rn@Qs>Lrvh#Dwi~idEoc@09qlTBtj#cc-dJ>`>ChTf_X)s|^
z!^<t=-*&%x`1_4`nHAHrB?(gO%S;nau`XMfkj1?$J0Xf~nMJ}YPtPv4Wx)x%*p`JQ
z{PGUu;$AlCfoX)=qoh}C-&mK)cE&k+2`syw&}F!9)5=K`lj6l2#NW=dX#cCXLppfg
zh9w4nJ?kE(vmff`K4R}AP_t9v%ucSNjRx-A?3~YDektoaAA4@w{rj^v_g3C&_-*kd
z=l)Zt^S}Ee^}fXVz1j2s)9KE=r<1?`mXEA_!(RQL^?s=3kKG2}zIVp;tQUIwSySmx
zQ`6r_r5gXAg+E*b_r?E7{JvA^55s|D{^EOM_J^17yw|x`yE{X2>W^LlnfE_B7wz_!
zy8Gt1UFMg>FOtfCR4#P?^IYs)d+f!#%dY0XH~Nd${-|HaVXs}b+h6={?Yb9f{|kS|
zT)eB^mG_PRbEnduId`_57d@%&Z_gsH_iMwih>LgGyViYwZx>+le*YZ%0G9b1_Wll=
z_~w4?YOvNYLAwW2KQ))1`Dj~GV=i<5<8cd3{)w-@?^?~jYxeac{m+ZnR~(=0Z~ys_
zxb5#Je|~>inI2cMyhh;f+wTvrYtP^R|H|g~b)Pzy&)=>2ra~{L@>lH>!TvqJ+~dA0
zum3Ue`mYnKYo5HWc{0`NlWvUeuW;}6Vczq@y!*qv>%+YB!@T3eyjO>K+lP7Uhk46~
zdGm*P{}1(kAL@NS)cbtsX^|k^y<xhS-M9M4U43*X>~u<?Zm{cCACap@tRdd-L%IKk
zbN>wI{u<7GFudD0gnMT=_tx<4uOZy2;oLjJxI@FZt;4#N!*z8cb!8%Tc_MZHMCiVW
z(7h9(`z~s0LCRKx6xF_4hx)D_>b-iX=jx&EsV)=5{nm)Ntr4?ZBWAWnOmB^t+8Qyr
zHDY3G#Q4^Tv8@sNw_5DmYO!~##hxj1`b=b5;PJ>`<)Qkd7s8$&I#s^82vsR{+Dv4b
z=kdr-B~g9S3SrL+ohr{<go>0pEk3SW;HeYi<gYSOebND8&pVweZ(M{*l{l>?t||9C
z<iAHdL+JA~3lUF=E)_0Up%$f1$B8WqJd^@d7&RvS5cd4hsq)K3s79&Nc4Ess4<&z<
zgKCo^ggsYus$6mr%24Vwp4c+eW08*vx9p?}Vb2|%Dz{vOa+EqvC$`M;Sd^^ds6J_i
zu;+<Rl}9c@1xlUf6I*6`Eb_%=h@;vh31QEUPL)G0LJ3Np`V(8Gdo1#1@h??5=OUD*
z)M+%4Wu`}xkBXw&Bo1NEhEA1zE<$lioq7{lrh6oLt1MKVbVA6pNW!zHvt0USZtC{Q
zcjL;=r+#-TWA)wP;^fV8u$0qi;*l91PChDvYLgg*J!?8u_P7YeDsk#fJTlG0$$O7>
z#vxUqWF=0+i7qodn0!<|s!sYK<XO_Gvcp9vT8UG4qRTW7CU2EO)kzzKJVQEFR=5ZS
zD{*R0beZD8<fZaaWm19A!Oi-dIunmf^*H3E5~w<9fsm(9r^*r+p+F^0jfqDldmQpq
zxu`NJLC905Q)P;akf#!-%ETiRJq~%O1gcC@5c1^dRB3S$a#Z3}n0Tb$<B+?`Mde8i
zf}VRiRAO9&bd@*-C%SZaFuACFRGjodz_X-7#ll5MTB#yTCs1jUf}m$lhl-1fkfIW&
zz{Df%9*3M&E-FrXAmEwPp<?19B&o#7Gx12P#~~+`K*c8YGeYu8oct4A+C7?_RW>S4
zx**`0(xGDDA|$TF$vx4f)uYKt#ZYn51OZQ%4iyC#A$}!J_K7ad9!-ua8x<x62zaV=
zl!JrTNkvd`QiFh}O^1q%ix968C)>m$O&(5;DiakZIS6>lbg1yS2(c=0GEF?v;Nj$;
zBB(HlLBR7*yUHJDp}&fp|0W!%^Ki0PnJ7Q$0Key*c9l2ILSGd*e@!@2<KbkdA}Bwp
zf!}eDdy>r+j_Dp-+Es2i3%yq4{5rv<%7e)k5@#P3Ie$zzQtffbR^_7Hqy&D?HSH=_
zoP}O0a(<a`q|)P%jY^>0Bn5uYo_3WZ&O#3rIX_G|QtolcTIHhbqz8PSIqfQ2oP};G
za=w{xq}1b(l}e!Oqy>DQKJ6+?oP{nba=w^wq}bz-rOHK_NeO(OI_)Y`oP|y*az2@G
zq|oD#g-W2zBn3WCo_3WM=LtI%56VpXlRjxm`lOt+NhxWQV$vpsq)qZko8*!<$tLZd
zsL4Yi&kVJb>1v+oYMz$ro;hkKv(-Ga)lO!qd1k4d%vAHtR6CiW23K=3P0ce+?PRK&
zXR4a#Vs+0HwUhby)SOIF1F6}o?73I@<Q`?uJ<2C{D|_x%KDkBNbBprH&B~sel}~O`
z_S~d=a-*{6L*+`|mM6LVllEj!+LAqKP4=WE*^}mEPnwcFsV93<OZKFi>`5irlX9{r
zrDRWv$(|IFJ;^6~l1ug^o9sy@*^_j#C#hsllF6PVl0AthdlF0bq(51czGO{$lQro{
z)}%XGldfb<I+Hc&NY<o1S(CP8O*#`Z=}63^Ju#EE#7tTfGigc8q&YE@ro>F@iJ8<A
zGpQzKQc29DoR~=|F_U6qCWXXI@`;(`5;MsrW|B$FB%PQ^Dly>jJ9$ji^O)+%qpF@q
zRZkvK^*o|_^02DsVbzm|R6P%=o;;}Pc~JG_0aecfswelWdhS;}xlh$|pX$lIsy|%}
zJfGQkKC<z=XXAOx#`Btu=Or7@b2gr*Y&?(IcpkFx+-Kvt%f@q?jprsC&viDQt86@%
z*?2Cp@tkMlIm^a#nvLfq8_zx)&n_F!HXF|-8_zl$&ng?wG8@k#8_zr&&nz3yG#k$(
z8_zf!&nO$uFq@OnrF#$q*-dIOO6y-w+|%j)Nl518?j)7<ize>*<ML@w+mpR^ldgaF
zxYwmznb`9r)@)MvO^<si%9V+|Pv)9|JzbeN@ky=gr0b%d_rg>vHT#~-HJx<*w@2I+
zrAosI;Cy(GMWr&#^%IxG$?OJ|e8x$8SUf)6lRvq-LWTdhM_iFgC4c`DPpOm9)hha?
zllDY=e6mt}&o!%{bI)wgPurwV&hAw?|GXmnLW_6RL!rv&y&yA}i0i#<^{#q2DQ*Wc
zb2_`3_v|q5=8*1PR~LnEeG#Plb^V`Yw@>|&{m+aKUSClzKeO!qhyRSa4Ba6!PB1b6
E0G!=-5dZ)H
literal 0
HcmV?d00001
diff --git a/data/demo.RData b/data/demo.RData
deleted file mode 100644
index 0deeba01a92e59aa76a5bfe51fed67a7695c1f1a..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001
literal 14821
zcmb2|=3oE=wzswHJ9KUqUfUV_J$4II&dp7&4AJ*uw`eDx{503ESYJ}7(UYyk(W&LD
zSJc{+?B+7@6y|Nx8%~;C<dEi>kx+C@ZRVeQv-6B*%qgD0IMG<_>i4qpzp=6EU(VRP
zom}$#cUJuUna0oGTP=TOUz+vm%f6qz`Df#QU7K$8+{itB=IIW7f8&4WJFQFkKkQFG
zZ~ym=_58|vmQ|MKzZmv^Hk#%yZTGijyUo8R&uzXQ`@Zk>#<Q-G)4R*{yI(y@`j}T5
zc=dGo`PEsQ-z09fHLu!n;&txEoY~i>-%h=K`q7`FkMF*1+86QdKzX#zeci6N%i_L8
z=WIW1z1sYC&CPB5Y|fvb{x|Pc)n=LYs^p$0KaW=K%h|p1v+vd`Z}(_T%l&M1U+;ZR
zbnd68Gp0@#e_C`#Z&gwDoygL>J-c_VG57BMyvi>A-g`^_BSk8wyJu~kar^d7>$r+f
zfA-AzyLIN@P}YB1O!a4f-%<M>#Qsk8O8Sd+&brR!bKn0ApZn+fwDtQwU#<DVb8yLp
z`b)vbKD@vAp-%tr>uVq8HJ|_bJNftKhkv6BZ0po_^B2s`GT&w^c}wwH=kuzXht<k_
zNxOeDe0NEVsy%w_C(Es5`Cak5<`>THH?EC$xg-32;=#R7U!HEau46y$cEA4rPt7g6
z5B~je;B3h3k5YxlOJ8xVp7^~qT7KQp^E2J!uPBDb|G9H$?zMvt?}hgkZ?3WTyq~si
zJ^yvd+RWPE%H%iU>wo`}^R1PZHcttc`p!~d{DS{2SG@7=pz9ldUftzWx8%!($2A=9
zRy=><ZgDQ-{Jui7!iROfyQS>*KYjDF_vggm_iq=kH$5l$BKf0Rj%?4zt1~ATy2ZTF
z=G;EB-}%hE+K-##zQ1R+|6SnU`z>&{*O9M3^=~a*`u@uFtmW&Se<jV|SSs-Q)649e
z`&wr{a$hm~nqu4Tu=vw_L4PuG-beSReY1J-=<(_cQpqoOAM01Eo&B>S^4p#?*)KOY
z9}VAFs@RhKbB<|xcg5kO-!)fEU#fp@ee3Jl@o(L)S+u{=GCOBy?^Jv3)#|r5FKz#S
zJ*vL;SKoS(+K;I(vu{l`RJ#=x^xXW&`CIFf3ad}qi5x$cetz!*I}_D&H_Pwx|C|{7
z_S%OOuRVgZ;<e8m49}TU`_1O{=AHlUbge%ZGF#-%l+uIG_kX`MPtaU-{UoXQrMEVx
zhKKEGv}$>@S|IO4eeFk;Z&P3G|FKi<@s(2f_?3lgj$N&N^e1V*hUTwj4}SeIZ&gqJ
zm~x_j%i3R$zP+$8KmC33x8H}q*Sx-%Ge?^F`UUa$m5+PwfB%@0G3R?%@f5r6pSjZ_
zvVU@YOy9#;tN!uW1-bpt7py&BZ|a{_d;L>rPOL@ozdvnT`CVndeT`q+UtD%)-?jR+
z+phhI`jvJ2V(re4bs@W7)|#5F{Ji>W>}I_sA5{91JJhcq3s<Z)U3czy)z4Ee?jPzs
z8a{2`J{xDQ{3CC?K3}}`=Skor&(_yh(z74!61n%`=w!{`i!WHTKl*m#*X}IatY?oW
z2<1n7emy5Q{D+rpy8OF}w~y1k@0-4EI>vie|5fv(Rrx381sAR@Rcx6pbp7YE>t}?u
z?l^zrR<z&r=lZMDm&~{QoAq((=i}$Y;_dD%vACxD{QDK%_xe-Koj$Mr`t9j$U;8W;
z<u>Wd<;SGY?)htHGjo?e=i2>>_dh(eTCD5kweQ>`&n?GIzsbzHf8$^I!hGZ79uJ?-
zu>D=Ic>m)S)gK?e`g-Ef*H0<N_4cWMzn)#6TEFu5JCW+RBa_{tPpb!CSB{Urp;sEe
z?LgME@0-)+)E4sn@_JPzB9)#Sy>L&4$vMVa)8&t<Z%teK^X~Qjq>r2X%(v{jHfPD3
z{w@19{;*W-lfABZ<@EA=u63s84Cnb=6~FHsGv9W0Wt&=Gch!%AkM9F&V>8ax>`<(&
z{nvee{>$tysj<t~xQoUA{PRMha*3q6$@SctN|o$4eUsEnUVr1hS2MS9cmABm($8<+
zV0tHYkK^Zh{dLwQvOX)DZ?&5zeSZCUaqZ{HhX1er^4rl{+yDJYJ^LEd&)xT<t@mrb
zz8JYgzgK!w?Uk=w;V)MoPm{OzI#l}lt(M9DxZPP#^6y*i*uZGMeX4WacCNMDzdyO@
z^lkm+6WgEkYuSUSUu-{DMJBENc)u{N_E+0A*5`j7+idf{YQDujYV*HkA4Fg6O^e^A
zwl69>+_2Vp$99|Cjq64J$L{v`lzeXbZ+4UQ&Ctoq@=xu0fA_yuZSrU9m&SLZb7BvK
zYsJ6Z_vO#plkrC%%@Nt{6<OaedHv{{pBHx6F3Uf&E+BukzWlp|zgfNS@4a4o^}XFc
z>BGD)r8lh(u&wtvZfq|9JfOx@`uVe6wQKghllXe?n#J)qmTBRNF>hyQ*Db$2-Q46~
z&w8odf3If$cyr;ig_u9v`Q<!24D0&d_9y*{Dtse#$MkpW-XD>7vtRAgc{-<~d7sm-
zV-?Dx<!`H|2pymG`Sr8U3w~;e`LCOKf2OVZ>5q|5-Fs*66{|C;kz4h8^D)bo*^BZ!
z)cJeWgTHehKQ#ACM$Vm^yX?2z{9P02aNqe`s#N%O-l)&r$&U|T+423C%mtHshyThK
z$4{GLUi|v!&35b0rTu)9b0vSho9*6pR61(<g<bpG%#VFjtgZb#yJ(;FRzJ7b(mnGX
z@7G@ct5R#auKRfD%Ce<jx0ZBXJ^o00+q>Td>5pw5tT~x~<csgu_}<-C`<h-)sYu@A
zcXIa9XRfs}n;-Bpn?G?s_tt*D;_gH5v*V>pZ)|^8btGo<yXAewUNPI#>OM~G>Yo4l
z=hcrXYtH?-Q1RHOR(+~^@Fx$xTWx<%q=(y8e@k4t=J(zuPv-Z(+#fc7Wm)M7D~oU2
z^Z(4*#g`S6d;Z+f$Gk_rzMQe@eTA&m^j&Fp`)_{!vA1we;QcATZa9Dbv#;Xt^Xy&v
z@_DQFX@2`NF@n8z&bF}lKQ^q_qu2}I@Xd{%yk08veYbk>byn{3N4<>4H>FQ8{k?ea
zC!Vu=_*bvLV)Nlvh4Q!EpB3w#Y}VfMe(R3pYI(h+>f<{n?lZ}|RJCN^{sXdh{;Q-j
zFBMFFB-6V3e5T&7rt_wGTK;|iUlq?<`0GaUyYGc{{A>4n#rAKpJ^n3qaoobslU9Db
zI?HPNzfzIiMW^@7e45SotLv!r#=lLi(qELf=h=3>oGr9{{Z0NyH-0=bE7|w^%;oT`
znD+Mz<@W8Ujko72+<A=m=v%RM#m5iwo?UNx?(*}qbvruj54Zm7d+V<F|LC@TKP-<;
zuIQh(@7d0xtLIMFp53GSsO^vW5$Ukq$xkHHYE18ZpY?0kgxIA=Hj2gnoEL1eFR3tx
zBmDN~S@-7#U;8Uz|NGT^=WlaL@-}WPbpBYJ^WXZW=!5Kg2d10uPrI}3SM=e@x9W9c
zrij!XvoO25F6r}Zndd4Ge&xiU|HxDR;`@YB1Cw)`zh)fhdTK87V)uKijNXTa=NMmZ
z6MNSpdr|y&$`RWo`|^ZZU(ea=@3~R?v{=yd>)(Grkh*iB>WOsv^PZoww~oB=`mJA=
zcQbSg->Qjq*Peu4d$8=yX+FvH`@7~(V(<HS?BaL#voY%>-bqz|`*_Uo+TR7X=AZw)
zYn7k*dG><$iXS)U&zXC%?ChUoKUV#}eyV#`ndz3!r|W|rN5t2NM%O#<d3ZGU*>|z-
z+Z*o-nwzdqdii>Dda-QCj`mM&t<r0Z+c$DOe017+ztr}Px0nA4JW=l+`6tgbA?mBg
z_p1NWv*YLYRoH93UXU7o@7HRzHThmOt-oGXT)JeQI>$86AkV#4di#UAxBI0|m%qJi
zK6SpcS=Pg4Yq+ny-%wR|{ezoNeDD5eKhxTOyw_D-vEI(z;PUpJ|DFd_yTwjAchsdo
zx9{V_SB2kf>pw;OOrQV%*SYKke826~71PdNyRqbmrRJ|+kIMf3ejjDuwRZny3++43
z*Z!LRlU{v#>yOI&_vfaUi5~d+^7Umy)pH*oF0z{@fAmZ5T6e|T@4f25UuFOMo;z-~
zAm8}C`2NS2%eV1q-8s|z+GL&k{Ry_&F=x-;`<XlK@A`}3kIxrO-FN&5?^<`wpAo|1
z@621eYuo3k7w6P{`xwghy>|UKO^fpEbGpAuwpgfz-k0$IwV^;}(x)^o@%`^Fm*4-p
ze%tEIIlrqjmh9cPJu?1+-|s*nf34c){rWpfKkvR)|2>vH`1nop{TuguS1i0cv-`N^
zjcbhW6u<TSYWd`G?M7>b^0DxJ6+-*h=PI<#wmPTz`jF54hjQO-zne%uZ`fyNZ{>3T
zLe-Y{|IE|E745sO=K6haefywl+4@)7BF|YqTP}ES(QkL|NVN>>-|+s;{N}3n=N+1?
zaq@hyY(erzw@I<QwLd?3RM}dY%s2kHg!@k%+xLPuKfFFK$+z|T7%v<DcDDG}pZTj|
zA68GxU3>DlT=M<5tDCJ<>inLjpOtPuwsGMPuh+9T`4@8jS~NFmhvB?C<+4|HWPh5x
z;=EniRGV{Wg|0t38FFWWZFObP*WDer+v{s9+iqAy+ZFZPUz_xJ&0DwIQ%%($Yp2dh
zuG#u4Ci?O5OV{mcYahRfNtko>xs><qeg5}1-i<4r-S#dlI{Tc@y@r1^f96=7JNVa$
zp*%7CY>YnlcZolH%51$PpNs6i<$ixgtYX=TdCt`g-uFfJMV(oHDScAfy7&9`KR)wC
z@|V}z*W3@=j-?r%;|hPB{qgF?PX`v>Kizk7@~!@BbJF7%=YO(0s{2~wo$H=Po7o}a
zOG-YxlzyBZ9A7@e_I^;VbD3z=cU@V&v#WbHeY+E!^+%>V+4)aSq2sF9)`GQ#`>ua}
zEuEfScRajH{%66~bDF=Ft#RL2cx}hR;vHh|_I+M`{8l;N3G;I+e?`^){POkhzWIyt
zd*&bi7gYCuYeD};{;kI>&Gs)Wye2cbOL|RydHctdj`byVpW;4FZLMAryS7g~HmWw~
zmBr)dH?B+7Z7bPO;9Mv2__^~Bt+wv}a=$G@HvZjby5F$cuXBCUufD>Z7ymcce4U%N
zF8@`X;)%&^`xNKdwC&4&W`1>E@}9=Wd!AWJWPd8nex&_;+KgLn??Yv(J!F%0?m6tO
zwb@x_>$Nfc&Gqg4Uwc0*zGk{Lb>97mSVxf6QT*NOlU`p`{8)Tq-I}i&uJ;RTYIfYp
zOZ#Q^x%&E!X?;@WV#j|S`zf%y>!<#=mGh0oHm?qeS-bui-%;snzg|l}7un6#yV+D;
z^JDQVi{p>V?pw4M)P9$^Y+m+`^Y?x6>oN7xnfuG?lIyml)#=(>E%|qI`@diFq&KA5
zOC0+8G1RNB=>Aua<Cb={Tl@QKm#*W^zINaF``v{<dDn3tt3CTGi$%Y!d;ak^T2FTh
z{QecFG4IOji&FbGSsZ)5dzSf}=yQjqc6)7nzPfJqYqOPKTiTzPo9xqH|7iAo=SSLa
ztX1>omt-ul=Pr-0zP88k_3D}5m#z<fbMv_L#{1&eZcO!$&uQP!|La&vjqKB!uFvcD
zSJZP?hrN~A^8I1fy@=>~(MP>gcGU9buK`s-a)(VT|JDA~dgA_9drN$T@wx0DuO!YF
z9-ZF${PgNK+2I>4!p_gnzsvtC(Bb}>9mf~Qn|}WN?dW;K-G{<g{ZBjl=ic&tvv%y$
zEaZFC{%Wu8?px*i<v%IUf4ltI#3Y-So3)=m`dt0OVzJ)J--~C+-~F^Dwx_uGZkliH
z-Gx78;_vBfi~sO5@Kjx*kn{Zr8_C1HM>j5RxmMdNeX0D7wd?&Ew$YzoZ`S^6-LkL#
zNYC8zYdaKcgDaBP{+G+%-;lhw%{=JQ4Ut;im~1<>-C0X|`OS|=$NkTJRvYnHX7+uF
z_=kVjK72kyUiGaw+uxwi2iiZhy}P_!e#;;8CHXnFvwp5`lz!aneD31!Ys!1xZ~eU6
z?B0i~tA7;iwUNAd{nOcy|Gu5KgFk6WO}|qUxAv!RZ2y*hn@jia{AUwpzf<j8p<A?!
z{r#QGZ&Xd`TK#qI%`KPAx4kdRtGB<DJmbxNUm5cw(#i8ROLt~J+GTvzY#RHwhkC~%
z?@nI&epbcdoj)v=Kbe^HtIw!!^ReEmt7Rtbck1hV<~sZ54Oy%Gnw6WkK9+ghJ3&6O
zxTy7Z!@gJlf(mP2J9)*G+<)$rA^-8|*>4~3PuX#9Rr9BhOO`h8^V|5`bJlCKjem|A
zep?tdU1ayBUtKS|=SS}mJmtTxv)HM4(c8+mJ7-jHwg2?<;>JhcR>m*YudA2e-F38d
z-;TR0%M|rK#2zm=`}KNp{qe5XHQ&~}&01@F{85B$+rC}Dww6sjH2K<}V--fJu>}{d
z9s2sQ_(c8O&#$fP)Q?DC{$qB1{`!^9c28vPo3wF>e?3@j^UUIR!TZ#EIhQmSy|?=;
zeaQ8HY;=0viOSxW`Ofzz{9pE<Z~cv3>;Impj=OZU()s@4rR(+QZMpSGd+K_Vbx)P$
zPZkSryf2<+bJHXK`=dYk*UtRjaO`K?A3yQa#x2+4Pk-C{dG(jCYxn>7|HHCHJ@DF`
z|FOdUY>%y`Z{t6DePf+|^zL`ZzxZ~(7W>To>(`@i&YwLDWA6Wv`}I3t=lp}q=1<%;
zKR$W&_u_?b6KZy!l<uvuPdVRytkm${hJyZy^=5y=A4X-r3OsjIqjv7UI|qEPnrEEn
zx_&ymH`?a=<5SC}4t~7JbnkzL|KsnI3m)DttULDm+rxQH5BBuE4Js^rbL@3=-plLR
zU(eR;cxUzR^0%Cqt8*PsuPQ7Hy5I6w|C#tdx#JydlKZ#w>%For`kHb0nnkqD&Rgr+
z&Tpxmv+L=ee*G=01!vz8x-Oo-V4p?YiIXe8<vg6XBlS(*EAH>*tBzf1m)`Vz!oJ5p
zBFj==t$r7JKYvr~!}Q)4!kbpN$~(%RQ~ja;qV-!1_r3Vk@Ut@ITlcZ$dBn>+54jV_
z|7CG>-b;78I|BJ@J}N(dR=qnU=5l|{j@}o0Hl%-*t;zf!YrFq!4EvSxw-&-%ZnJ-9
zd1oCN|DyiwpL1ViZ=F9De|Pb;+UMol{_QgGPm<sKHRFZ;s`Zx=|HMs_cf5VFU*z4b
zuh)ClU$Oo5=lYfKhu`&$&l!ICdUkzWzfa8C_aF9MD}QtP`1!ayX??q|^<P`Zy|2*k
zVC;vfRr2w+a*p#?AHNxUe*Lxbvli{w7V9)rKYCMP^u&E(q2Dv@uQT<}aUOenq5Yiy
ztirj`>zXC()juy`H=n64{j&V5#r>x2=}oU~YEmDt{mpo@KXAvY{t5EEp9Q<r#eX%$
ztt>PBCAe?9*}0GJobT4WH@m!j;=W`5bZZOeoILxzd7p)P>-@rbKGke7r|ZvJuz%V3
z#_M(V+tX_HP1W3=+=Q=-7uNpkULR0ndii?$>uG;%wsD`$UuC=ejd<E)*2fbTt^XJ2
zb$8h}<@x{H-rQ_|9aD2(_UrvZp0@h*kD^~TYUMB6mw#0H2iIBsSJjW0_B`y_aaZ^8
z^MIPCvd=sA&5}>5+?@QB(_Qw<uA+=NkIz20sFvOGMl!y;^kBR6Mf21-AE(-!)5{A#
zxubV>VM)jOYwO&z|Gp5|edyo$m>c~Q)+arlBl29M_UqrC^%wS8d&V!nc7NgBms~OH
zZO!w{_Mb5Pap|$mO^H40udmx^zlZ(xcSqUoxeoES{-6I)d%*R7+=o==Uu*?)H$RpM
zzMp1imK-m!TmO&8)}wV_Z?t}qExq=*r+)d5{Aml;Z)5G@UbE5fVC<RojncE@&xmI2
zE;{ynTK<w3)^&N0_<m*9uKgKneXgo*tKIQKvDW9V{fVkG+4tpF_RVd<?|XkOiar1A
zdb9PqpHc4f_np&y%-;K&C+5$euJyB)^6jfGIA?z@cz^Gq(oJcZ^Lid;?Yp$&*PB`A
z>-|G#Pl$c_T{q@z<gNOX8IMnYZhuz(C+7C&tL1I~#D2$LsD9)9|DW`(Zx0vQ?JNBD
zbG_*{#s9X!e<Q#3d}O;fzg7BEy!^W9^L1;bw_pErtROJ{`Lr6dPi`4?YYY26w&lye
zdfYpu^xz}jx9Wy>M5^O9^Cj)RP&r4W_Umu;53+CF#nu(iE0r$@ejdKwaAB$ZXRRuq
zq}LPIuif_e!k?&bHv4L7o<*~N-Ieh3eM<PMc<nykm#^<!3QPL9*;EfSa^m*4{DO9r
z{ml+JJK5wK(XQ78br;|5-EUo8a{YV9uSd1*=WKglJ89&}-<MBz=w4s4->~{x*;cbY
z)8mx~IG={P$zD5h-miAu&jV8Fx$mp~`p&!`VDs{^=jyvRroXdqvEG;dY00(iGq=nC
z-t?$!kJ<J6$JV#=hu^-x%EnpwM#tjHOFK*d%sC`27Js&4!^T4KZ+AG>rG{&Mn{#S%
z)cu=(R{fqG%Rm3+R}HcJAAdBCNxyi#(84Sw{MNkWpK}6Y($B3w<huJ(&DLZ6GBtIr
zsV^S$)H>HqIVAnKw(9MvqsJp(zy2KGYkm3pfvdj*&fd4`^Zgiq#k^B}>bw>2Z&|n>
z-)LOhc3!OPs?P7Lzq0NIUvpgle^s^5*Z9BxRupfNx3tdlkL>?)gXQO{kFUGGzxuoP
z%;NR=e~;SLmbF$L|CIJiE+GC(yi9uW{)9WXXIUKIxm3RO(drAocE7K4eBAcK>wf&T
z4mr)bC%2dXcyr3K`hWJ@-TgDqOa9dm@o&3!-}KL%Lz6Wt+h)WreNwD*yZdc@a*gTd
z`;6(l{>y%>-x|+xZ1S&tF;@8*wZ==E`(BIHvAo@Q_J>zzviP?<Cv1=Z@Omt_OP~AM
z?YVp9mj1{;WmU24R%B_s75|<x`xnyBo8Lr!HJ*L{*tPI}<3~!*HdY$VNH_R$)h091
zwKAyJ|J6R9?^6$j8_P`J_Ws(poWJkZ6mH$O`D^&=)djURTW7>y<-gZo-1cAUmqK1z
z%{%Gm7rtIyJ@NO(N4(j8xBk#-k-pWwXkW*Bk#$G)H4A^2+~wGo8ZNf(>+UBO=fdXt
zc*Os?6D7W;H2${jKBwPPf^#gsw3fzaS3DN`V)!ljWk*$y{#)@I>;L^<mbvcIXX|^r
zuZiEb+CJz0gnd71|M0!C{<Paa{;zaL-)oLrk^f#tvj<;mU9m53-TFSknjMq(8Squr
z`QH)JHaQoz=k=YUy5t=eY#V<bEoJ?@T*Cj?8?VPvl@~8Fp8p=harW_w{1f%Zt4m((
zZI#%u@8{LppKne*celtFduP$NS?X<C=jyLlzwLSTJ@1A1&fGo!CD!Giz4H3Zqn!f1
zTdsfA+8o2V<@3X~KN)k>@4qjxP|G;KuhgK;`sJc8uJ*HThHjRzuC@JUdF%5m3$gdo
zw=@5LUF-hBJ?w9p$NeRLFMcTRQP+(*J#TMmyx+~uS?!Mv_O1Cj_nPuj{XX$q#a#YL
z(a)b{@0@piO2y+Vy2o!m(w<ZPeA;81%dU@`?Dkb{p1H5UUO4aZ@uMHjzRDhFE;;<q
z`LWC;^CRbX-3{IEwLV~ua7?|U{Mqx<{_8#3m3A+`YhSWWcHidU*Hag+&ySV3{p(u%
z*`Hu@8mhm>_U}6pzVGMJkDva=rF~R6mz*)j_yxD&MZLwpYcGiB-ua|lqx*yF*ZjBQ
z5B>_o|Gd^Zb01f^n7@fvyt#b5ZMyQh^Zf;XukP5qm9g;m>xcPYCa*63-TS3{ncer7
ztFO08+FMJcA8D;Kv41S_rnmc-Q?%{=mA?;4<SWRZy%x`($M&o3k9kVCv7AfoY|C}|
zclqB27w|@Im6M#`zq2`@j`{udRfdi+IpV9$H(OkPq<!W3t?3W<YyQrAyMN{9*Rsza
z{=QheGOtojEcgEOe%sPwzkN&heNK9ucCNVYLZ|QJ?{3x?z8mk;nOE}ReRps5qFX+>
z(g*h~-S_^n=O3B36(TY0dv`4O%==ETK>B6!wDs9b)<^v5DZHz6ZpZlo*)4}Z?)+uD
zbY}lu&>YjtuN?L9kv}HBs?K@t(7)P7Zi{s3@qg=X9lEh_pV+!<Cy&=Y{`>fr$=B;1
z(YfxI_n(y!|IY2Wuj@tj!*j=?A1<>IPm(eI^KV5p&v&sGveEr{EBB?nu2cS=(z|(E
zx=j9+@a?YEU!R^`6Op|~yxUyt_u^w`pM|eiRA12jE#7wjj9`~LQ)*<r-jsXpwTt?-
z?rZ$%b!&bet(J%nmQUV)VZVO<vck0%$EPdDKG?NWR{x9XvFESkq_1B8@Ll)E8L8*N
zciFyN+A(YY$4HLF_blS=YHOd*^O;w2s@=g~=K0C;FBT2!@(w6%-;!**|C(Luv9HIz
zA3v{q?|_9^?RBa2x4Tyt-o5l~t*Jy|ZE3C8_oJUIrM`Q;|5zs-eZN&+Gv-S>tG`op
zxBA+@PNzOkT(fcO53kqOFLPeJzHR-8`*`r{*U`HR?oCK{-e10WeMYUXz0|$Sg|bCq
z9p)eYmCP%s;r|)={Px<~rT+2y_R>{}bDrF8e^uyKdtu$0&3CW;d1bqx`TdVay*Gco
z&e>;G{JHw_j?23k=0}$;Jud$H%5x5<-+wc0*kkTnMAs}mUbyYP{O{oJnP>Ag&Tag1
z^>rZ6b@Tlh=QH{xm|xZTEu5==UVhW_+D8$0m-o&0?=?%lU*R7vv)G^}H~G(;MelXy
znH*y*)UDOM7xDSeF6}<4Tj{TAm*09+^)_1E>HPj#Uu4(YzJFXDb^oZ{(ux1XKeoL3
zAF|c?+!1lj*a@3U_WAvvm-h3bMgOezU*^2{;niRL!sgA5@M}67-S*bL-Zyo=O`&(~
z^51#StaQFL|B_qq_eR+574x#Mm~ZiyjIWA~=&4lB)!n`6jpX(B0`b1ZR(sZe=M#I^
z(fj-Tl|RRpT)(ky;`s^oUv0DYNgdZXGI`%uj>=2o`Dd=Nh_ASRVMjxIjp@<w3+`O&
z&iB`t?y*wa_J-?kR6k?i=SSt*+wWNW6|+@6wrGE3b@%$2b*6d=5C0y0#QXfw;~q=B
zxBGt?%YToIFuW7#v+`?b^5ban3Gagn(<@n8XP*jRp11E`-}mI33q+r{)TPINxO&0t
zYq?eJ|7-KIb>m+a`klS*`%|Vh`yIP&?ZfY1fAf`Gut~i+)jvKyR{q(tExcLxIG#;h
z!2dRTvvK}o(>wk;`wrD!{P6ek>%VVTRxkMf=jHMFRgWt7RY-3CdEmpdo5B6|-(CG}
z%>RAuOP^Qw(AJ-?ey&}7T#V24(n`+v`A`2k-MeoScA%8=x@gRG<y)U;|4I7RyYYJT
z`S?7!+dgNz-v0b@{i?g=zrZuMdgW8kpAI#Qoq9RcH1_J{t7W@(%~<F5^u^D!ZhFN#
zb)IaRa@tETd*_`~Pg`TST<Y|uH%;S?(#zPn;3<L;#u25M(Yipi36-&Q!PKTOt|-k6
zR0e01YR1$Bp-o{ZjI9f<HidCV>1LoXwux?8yC9b}S34ux@iy13unXIoz(TPcw^|oq
zCAe-ab%^0Qf<iEP&YgewE8Du_9c)on>`K<`O4jQiR&ie!FW}uFec&#`JGKJzhTV*J
zcpv03lru+2A2`LR&AvuvLlK;D!H|6o(*~{sS`4chkQgOLjxJ-lBUIphK$melf)Nw(
z$4vA=SVJ^p9_tQN26G<o4$}wM8lq7c(Fg7_y%Q_g-mshXj_!jz#&Yf*U?J8AaSZQ}
z38x2}8Xu$;Gwzt)P|b-+aPOJItpD_YtK|oOzCV+h|35ikkCVBNG3LP{nH?JJdZ!w@
ziy+K@CmZXF9?Z94*Mabnq=p<*tD&2#V<OjsrA$dnm@HN?^l3FL=6bM#B<7=~Og~mJ
z?9*yk&-Gw6lu5E7WZGS^k|9p3VLHeI5N0>mgGEd`Rxs#+d`22mXvYc`J*~!WuLp~$
z!fcUP$?{LLvEJ*!e4!u9q0Bi#7Ry=AX*L#nJ(wv}v5aM&W@9>tIeF3nW7Q8noHoHs
z=ck}D3%oh+1Tk%&azIxV!A$LqSfG5w9z%=MpJ1lvQw~_F7Wi?_L1!WxfMgbg$+$h_
zfUanPH|w1ssxSqfhX_0mIbbPT;Kw>Am=VErJQK|LKjeVD=m&q+KfxfLcqojAB!R&^
zVkue>$Ql#EC?0yi70N`IgP|O3I(AXVGhl1=Ll1b1eh6aS6T*ne<lPg(tbg@@x9KsP
z2l1xIq#rb!J`<9^ec<bj=U@4YX5Wm@aK4n^YsuGh(kmuDaK6yGh4K$hN$uFgwr-92
z3h5r!f>hp^NapZs2d0{SNap<$!EAr+KtCo^{hEcmMf!nG(}E=4IT6h2*Wk?WR}Yk%
zo(bfw33q-II^E*>4fPjHHQ~(nLBbzmc_GX~(}Eb@Jz>o2ubzqg9KTJz;@!;e`xF0^
z?*6~hHv5OY{J#(F|68*vUdz6({qyso`iGZir|<jxQCzP6@a}n!-+zw3ziyw?dApB(
z|K6Xj&Y1pnZ~FaLcf8Nd{k=r~=^gF=tBTiGah{%h&h|%eq?MRk?F)X(?*;Gftge?l
zz1`ybiP>yFm9Kw#v3qVcZ~c4I`Uy{)&+VPI|M$%K$39J;Q!RV{&+9sK?)vv%_aEDN
zUMoT_{U_VB`!DXK|CgL*|2pQQ+w)WV-^6?sIL)8(@ZW0r$9HDOzrOR_{^+5%_r7^O
zzjwB_{$cyv>aPE_zwE8PpSfT5?&kZqcYg1iUmX8OcFylh`ybqy{<rk`zE9kdJNB=y
zbDvf7_V@HpU-!R_+V^ns`lr&7*YkhP)v?jr{_|em|A=Yl(|;D~*w1^at7G5zRDIK*
zlxg!1eUdmWzG1_;?f1_eir8Msm-qS0ojP&u{=0AWlPYIDm7JCydGGiq8E*e~Z^a`b
zuHUY)-@bpwQ}2j6tJCKbeg;j`PySgXCincdefQIJ>+<KXBkN>OuTTHUt7CJmykgD!
zKWfwNpZOFWv2Wtj?VC2#mDgFHexHzd_uU^ov3cvQ_qRX&zUhzMwEO2ju}4(Yz2BgF
z=cn4V{Is8kb!^ni|7o0-PyYE>$L4(bKbF(;OWx_t{&sg(>gl3`MJGS)ji}i7u1<d1
z`_qR$mwo@e`45BGyzjg1JDwItRQ!AYPhi^l)SnA=?0cTJZ~ntEP5$JkN}Zg|<@F|~
zXGd1(z55h%=jqNJr;7`xKY#QtI<ms<-M-$|=TF}G7v&v)|MtgT+5Ouee=qx6Z4+)i
zZ?)|Ct9{S2ju&rzyz|z_n70+{${@nOkM=#!{w^oBZN2RI+wb=3<jl98w_Eo7-M;5}
z$3Jg-TzUKBp0^eE%53tj=WUlgf4lE_?swzJig#r;<p_1f+kac?#2@@r99gHy4YDWm
z_~$K;D{p-K@xG#`{EhMM+dAcuc1J$VpPTI;apP`ul4a`7qeaIK)m_|Fzs-m7t#{-;
zW$xc~zd3d0pE$&FsqT_O&hxjY*?)>m^G~g0Jzbud_&<BU_UZEEpYL?uJbt@;^PaA!
z#u0KdEwwvpXFRRmuxHj&!D;T1afc86yHn>s?R-k5?CIB$cHK{Tb>^Qv6mumn;*;YH
z{ace-b62KDXWg6@wQcD%{LE+5roCPH_n=Nr@Z0Bnhsv&QirtzMu6sN5G!oP9-lp8G
zInlbe!%lCz3S|bKF1xZxc59Bi?rk{ptvJu89|}4Cy0?Q*f4j2Db}M$~8Ft;<A*ao*
zZt~p<VnVI5?`_Rpoq9X#=DMhDtEatPmAV@n^RJiKwzc1$Z~Rj;E&jkK?#Q}}3j1d~
zT^+G+*3;=5Ht=pc^ZTuS`cF}v{ga+%N7RX)W{;>dJ$-#sjV$-)>fhOsGwP$B&7byj
zdFsE+nF!|1Wzo;(PkXy7b$1?;C^EA;>e;+$x%*S!WzGaMcP@uB|1JVEZ`mLHv^?^L
z|1JCDhxT2%m$=jVMBmcXM7{j~Hr(r_-^QmrG+q6)bn#zpojr2f>U&$SPo8%D<hT3{
ze|)FyKYHl3-|zB>e|o3yC;$Ad^S|M#{H8z3)9TNB@{jnZdD=et=X0Gs)3*Je^fZ3s
zpX_P%hd<><R-Db*-~06b#y{25{vSAGclXn;YyZ@z-P``3`*eNk&wQQ#bD!34_~SqA
z|IttTBmeQA{-64@UT4pn?R(rdWvxs0doe@U>RC&l<w7K8b%MaN*1pPxzx|pG>`Tt@
z%YJq~{^x?BdC3`G+0RbL_gpa4FFC`Fn>pu#p?b*~PT69o<3uvqWj{L}ueo3-Uvh?I
zX3Pac@scyFvc-<a$z<AGFk~;umv6CLkj(dD2D8<(hCb4mkH0QWmU}US-RfCmU*!U#
zm?YS8<G=lhjknq(rt@#qlTV7ZJF;s}``^d4HpkCzwBLKc?pxilKbO|;nl^pI+&xqN
zKAU_q-hS%ahoU+0^Uu|Nkj#0`zscA3gxLp$-+!9lIP-5@Yj>phyTJF~0=NGk{`39t
zo%e@xwD~{(us@Uh(ctx0gZuXnpONn0aMtF~=l>l2mzCcn_kSvvfBL&>%KMjr&w8$J
zsNFy9`OmFqSl6e|KKJwfocHUKqJP*to3Q=Eoqes(zpcIVF89ye?VIl#_U=}nxqRM}
zviHe7-;LMTu-^W0y*41})jzx6zxQ^{zklWU^;(ho+Wgv>zmvkh*7q05)_*TrdieH-
z&B>>~R(wmo+#hRy?(YwqKc?5`f3^CSeC_p({(AY#@x9k)uRFe?|7CS;@RIx78845o
zoNnj#e(8K>x$WALkKg<bJht~-x9;Necl<vk^?R;*WBq+)^$F>n6Co(sTVDUW<w?8s
zH=n!zj@a6t`CDRZzvu58JI_!2_VZ?V{pr0u^7faM48LBDdoX+bcdptq|GWHm_C@}d
z*w_DC{^z#JIk(>{SG}~aKdXN4a`m2D+rLy#D2%^!{^(-=*!pYFPyPPu{6ByH{1?Sf
z|L$9Ozx>PfzjjLgQ|ykfeVTRtQrz<+)smmDec$b~UwO}oq^JMlb1F9;R|)+#_qTX;
zP5u+NJP>*u{!8Y)Ms3}S_0AI4x7t0w{r~p!%HO-1x9|OQw5PFOZhcz#dEXz`_Ig)t
z@W~20*0MA(;-8Ur;Gs?W8fOn~`pErzCeP>RW-&KX%_FK>pL#{?YhLOak=Lpk6EUsj
zse3@3h4$-%MXRT5TEg_6_ZXvh_l6SA(>)8`Y&>?!<Bk03qm#V$ZTP}<J7IoBe#EPm
z;;z&xBel&!Ryyv7C)L<ITGREt#{Y3c^UD6*hr4FgN_(?xTwb|z!JpRtmow~p4*$Q<
zcXP7cGt;ikw_Mj`x}+c8we8lr&fS)e4dT|@K5lp}#qhn(;&DN7S;CHe=g;xoxiW$O
zx_qta6Z7YG$6oiZFP3+=)^AJS-B(w3rf%N#y%k;8?fTyxQt!Lb|M5`q-|vT&_kOVX
zcl^%dy&oDSwdKElz4)tUeXzyj>d$M|pVa<1)9S_hA0F#xyEnwf|CVF?X87fXeof0Y
z+j^E88|8kNeqOQub@u(m`|Q5Szi;h-X#FNI-n{?&u6_IN%e^oC?6Lm!>&G^ixaIfn
zc=s{xQ^$k;1I_#FmdU*<{p_%Q_3PQ}RlL{b&+>n1a<{v0bK84e{=B!PwYB?gZ*vO;
zUHZ$GY@V|H=J~`Q(wl3_zge)?JYlwLW0YUss2TQ6@44I$gZ+8`_TK;ZRz7#)*ZTcy
z=C4jX7x(Fju)qE1L(^@3ANlb8X7KqvAM$?IJMDR0_~UisbNTvLH^19`ym9$_Y`L3r
zdT`>AlNp@qi(TjUJgYI1dU!J9g!*F9`90f8jHGU#d|`4kW%|u0D{ekncJs-Sn@<+q
ze6rx?lX*8&*4}(F>*kXgH=j(q`DDt?CzEbInQ-%o`^_gkH=lIf6p>Gv?41^zGr@PF
zqgiEwnJWiKD7a*T??y*6SGJoX;wh8grv>{=@D+44t4uL-WxgrGof51v!S{8o+0RO|
zwO?;;xs)!foN}3Ko^N)p+1fcbw?w1`$ISQbPBYV;d2>rb+U25|zTK&2x>Ik8>`lKs
zX{IlCnp)xSEfHy#gJ$}2r<&<bp802if$qecTO86Zf12UjoouGtdvl9H+T~3%e7lp(
zbUSZuQAoR7G{d($(M-4X<`#iTPc%J0O;ma5F)2sI(`b^4xaXuUm6I-$RFphFO;E}7
zn53iPDKtri-E&ft%1MVwDhi(0BqxQac&bfO5%iqYuHxxDDMit9(gYPtk4azHJ<TTR
zv-;LCP28g588%5p({s{f6;IDeDPS{xdQ2)&@id#HBI!A)TgB6LQi_tN&_tD;?vt)4
zdp?_}^3h{bxr(Rtq!d}tNqs86<0gGk_WU+c<)_D_Y86l0Nhxxkla_+45mE7Eo20_%
zIjLU7(|%Hlyyv8TT!vV>Pnx3ad2FJ}Lyt-MDxT((Qlvd6^|IK1^1L-s<u%AXCoLx}
zk?}O@Q{i-<RHN*<Yof|+P*9vSowP*C)2LUa({0i$=SfqP<Hd?1Po=6(-x`zm^plew
ztM3jMqh6N%TRqNARC($#sZhn!V$u_7&r7{3H{B+sD0{A(sB+a~Ql<)Lm1j|(il+Od
z9%au%6IJeeOv+O^X*Nkk+Ow!vWu@CBA7#%)6IIT8OiEKZX#@(Iq8=4Zw@Ez8o=p=~
z_Ipf<Q#q+8^hZRi(ruE7vgfpkDknWAC98NEPD+vR6zWm=={o6)l4sdOm7N}wqE$R~
zC#8sa3U#Ygx=z}n<QX<mWu?cYU=>fzNhu<pLR~69T_%+%d74dBndvdfUd2;>QVKsP
z4keu@wJ3V-o1hZwF-cd&Q*6=`F3+G=mZi*-OjJ%vPFlj_Y1F2|={%`M(R0@Xm1t0!
zI4L@5374l)t4gQSBpXG~Srb&eJtnEBXjn`rQS{t4K_${-lD3Mc=%f@bPoY+oN~cLC
zik{OZsCarzQdaR4oRq@hDb%9!({a)l1<x{Y+L2cA<eikl<|)*yVqGvvMA5Ttf{LTZ
zBzYB2{z)n9o|BqYb~;YFqTrb}LB-HxlDLW|_oNh7P@=SSoHRwj(`|x^qQ@kD6;Jj_
zDa@Xe8dY{WObSu(RGXk8=rM_1#glnb3Zv(w1{F((Ng@iKY!g%%JtqBE_W194s-Nqj
zkmtXNPyTrX)vL_3pL9gt^IpF`D1L<+RDRk|`XcZ7tzYG*`=sB>o_{B%{PPg1SE;n0
zv_;<YTEEIm_er0XJ%3J2`Qsr}r}EQoQi;6hwtkhH?vq|Cdw!jm^2<Z0R;AKzl8L<M
zw0@P7?voxXdw!gl^20-@M&+mNq%U%wW#F`XTiNsN#FTFyLe(mjwv)EVd4~0?taP7r
zS=saD#FQ@_T$Kvj)t`P)5&fWY^sUO)w<=5Ds!V;W()w1V^sP$jTa}j%lbFgU{VJRE
zsch1#vPr@mDy8pKQs1eBzEg31r(*g}MfIJE=sOkGcPd}ssyuzGa`i2WIm@*?UfuFM
zb=Ontjwjb0&sVoSPu=!hb=$M+wrAFDPp{jaTDLv9ZhOAE<>_VUsb%QNW$5|J!1I)W
z=PCnFE`85edY-5BJXh&?cIkO$>3MqTd1_gC`iW1JQBu-8tiJM)y5=GEl?T-|52~*`
zpssm9edT_2&Hd^t_o-{{Q(w7PU30Jc%023ud(>C%R@dCEzH*nk<}US>JJmIJs;}Ik
zuDL^f<#u(=?dmJHscUXiU%6FXbF2EwE$W(E)K_j+*W9eWa+A8|6}6R@)if`wt-PeB
zc}Z>MMK#TfYAY|OX<krUd0tKPyxPihYMSTNR-RSUJgc_yjGE>dwUwvUG*7FoJf)_2
zN^RvyHOs}WL34eB=J*E9_6?fl8#L25Xohdlbl;$9zClxcgQoZfP4*3%<Qp{6H)w)y
zP`__bpKnmFZ%~hKP`7VTm+z;R(n-xCD$g}M`T{4+S)^d;;#kzy^n`=;Bq!G-W)YR=
z>K=Xm6XwiQu(WY3Dr<W3g5~5%&PmO}D$mtC`g|wMnWbQ9;#ic{^yCD~$xhBm%)%<q
zRXzH=C(N0qV5#F+6xQ@)1<OfK&PmNeD$i9s`aCDhnWSK;;#d?hSzmkSnx{qD)5VK4
zGY+X1aZP#BpnmeE*QCYGD(b?X=MGK$!{%NTGv!H$#z{u+NuQfl+@(GJ%qGd?s95gz
zC^|9`oY2o5oFvn%YRT$d#4+tjkIG3ww@JZ^Ri1zII5&He%x@LTn_fi+rabxXK1sMy
zMcvDjZ_1=O>Z+DnUPVTepRoH++PqoixRGa_tg7V^&!VttPb{=fa)wM|UagXz@7ecf
zl1!b7rHc2b+1e*JhE5XRs=^-a+2`<S=>pE_R}QIeE%4Nh3F5}X{F7zu>uT1UY}OlO
z#_Mjz>u9!euBYvBHm-feJM<pxn(}(hq<d<6w?A4aUR0@R-=BZ(p_JDy^MJCQjfa-=
zhDhf;_;Qb}uIPTw!(UA67|q%5$P^SGIE%*MKA_96o#75s0V)GNJ0Wt%=hSR32VKUk
pELVgAv>nz~U0Hp9r~h&Jv;J-8#r|LC`2P04yl0qZPvve#1^~s2>LmaG
diff --git a/vignettes/PolyHaplotyper-vignette.Rmd b/vignettes/PolyHaplotyper.Rmd
similarity index 77%
rename from vignettes/PolyHaplotyper-vignette.Rmd
rename to vignettes/PolyHaplotyper.Rmd
index d65b86f..045afc5 100644
--- a/vignettes/PolyHaplotyper-vignette.Rmd
+++ b/vignettes/PolyHaplotyper.Rmd
@@ -9,8 +9,11 @@ vignette: >
\usepackage[utf8]{inputenc}
---
-This vignette shows how to use the main functions provided in PolyHaplotyper and
-explains their output.
+# PolyHaplotyper vignette
+### Roeland E. Voorrips, `r Sys.Date()`
+
+This vignette shows how to use the main functions provided in R package
+PolyHaplotyper and explains their output.
The PolyHaplotyper package contains tools to derive phased haplotypes from
unphased bi-allelic marker (SNP) data in collections of individuals of
@@ -34,7 +37,8 @@ are available, or if the parents have missing data.
Load the package and demo data:
```{r}
library(PolyHaplotyper)
-data(demo) # contains data.frames snpdos (SNP dosages) and ped (pedigree)
+data(PolyHaplotyper_demo)
+# PolyHaplotyper_demo contains data.frames snpdos (SNP dosages) and ped (pedigree)
```
Two inputs are necessary for the assignment of haplotypes: a matrix or
data.frame with the dosages of the markers, and a list indicating which markers
@@ -76,7 +80,7 @@ hblist <- haploblock_df2list(snpdos, mrkcol=1, hbcol=2)
# number of markers in each haploblock:
sapply(hblist, length)
```
-This small example has data for 7 haploblocks, each with 4 markers.
+This small example has data for 7 haploblocks, each with 4 or 5 markers.
### Marker data
The marker data must be supplied as a matrix or data.frame, with markers in
@@ -141,7 +145,8 @@ dim(snpdos)
```
Function mergeReplicates has converted the snpdos data.frame to a matrix,
with the rownames taken from the column "marker". For each individual present
-in multiple replicates only the first column name is retained. We see that
+in multiple replicates only the first column name is retained, and the
+merged data are the consensus scores over all replicates. We see that
snpdos now has 31 less columns: 1 was the removed "marker" column and
30 were columns for the duplicated samples, corresponding to the 30
duplicates removed from the pedigree.
@@ -183,8 +188,16 @@ results <- inferHaplotypes(mrkDosage=snpdos, ploidy=6,
haploblock=hblist,
parents=parents, FS=FS)
```
-(The first output line mentions ahccompletelist. This can be ignored for now.
-More information on ahccompletelist is mentioned in section Remarks.)
+We assume here that you run this command in a directory
+that does not contain file ahclist_6x.RData or ahccompletelist_6x.RData (more
+about these files in section Remarks).
+The first output line says that ahccompletelist cannot be loaded. Then you
+will get messages about sets of haplotype combinations that need to be
+calculated. These would normally be available from the ahccompletelist file,
+but failing that they are calculated as needed. For this small example this
+will take about 1 min. When this is done the actual haplotyping is
+done, and again progress messages are shown.
+
This function call returns a list with one item for each haploblock.
Each of these items is itself a list with several items. We'll take a look at
the results for the first haploblock which are in results[[1]].
@@ -195,7 +208,7 @@ allhap function.
```{r}
names(results[[1]])
# show part of hapdos: the dosages of the haplotypes, for the first haploblock:
-results[[1]]$hapdos[, 11:18]
+results[[1]]$hapdos[, 1:8]
# show the composition of the used haplotypes:
usedhap(results[[1]])
# show the composition of all possible haplotypes:
@@ -203,16 +216,16 @@ allhap(results[[1]])
```
Haploblock 1 has 4 SNPs, so there are 16 (2^4) possible haplotypes. These are
all listed in the matrix returned by allhap, with a 0 indicating the non-counted
-(reference) SNP allele and a 1 the dosage-counted (alternative) allele. With
-more markers there will be many more columns, so mostly function usedhap is
-more useful; this shows only the haplotypes that are present in at least one
-individual.
+(reference) SNP allele and a 1 the dosage-counted (alternative) allele. In
+haploblocks with more markers there will be many more columns, so
+function usedhap is often more useful; this shows only the haplotypes that
+are present in at least one individual.
Matrix hapdos has the dosages of each haplotype in each individual, in the same
layout as the marker dosages in snpdos. Only the haplotypes that occur in the
population are shown: in this case haplotypes 1, 3, 7, 8 and 11 (appended to
the haploblock name). For each individual the haplotype dosages sum to the
-ploidy (6).
+ploidy (6). Non-haplotyped individuals have NA dosages for all haplotypes.
Another interesting item in the results is imputedGenotypes. For some FS
individuals that have missing data in the marker genotypes it is still possible
@@ -221,8 +234,8 @@ dosages. Below we see these imputed marker dosages compared with the original
dosages:
```{r}
-results[[1]]$imputedGeno
-snpdos[hblist[[1]], colnames(results[[1]]$imputedGeno)]
+results[[1]]$imputedGeno[, 1:8]
+snpdos[hblist[[1]], colnames(results[[1]]$imputedGeno)[1:8]]
```
The other items in the results for each haploblock are mostly intended for
@@ -245,8 +258,8 @@ knitr::kable(ovwFS$ovw[, c(1:2, 27:31)])
#
```
overviewByFS returns a list with two items: ovw and messages. Both are matrices
-with one row per haploblock. Matrix ovw first gives the number of markers and
-of assigned haplotypes for the haploblock over all individuals, and then
+with one row per haploblock. Matrix ovw first gives the number of markers (nmrk)
+and of assigned haplotypes (nhap) for the haploblock over all individuals, and then
for each FS family specific information: parmrk (0, 1 or 2: the number of
parents with complete marker dosages), fit (1 if a polysomic segregation could
be fitted, else 0), P (the chi-squared P-value of the best fitting segregation
@@ -283,14 +296,18 @@ and parents_arr. In both cases the first dimension are individuals and the third
dimension are haploblocks; the second dimension are the columns as shown above.
Array ped_arr shows for each individual in the pedigree whether it
-has complete marker data and an assigned haplotype combination, and if its
-haplotype combination is compatible with that of its parents assuming Double
-Reduction to be impossible or possible (NA indicates that the individual itself
-or both its parents have no haplotypes assigned).
+has complete marker data (mrk), whether some of its marker dosages were
+imputed (imp) and whether it has an assigned haplotype combination (hap),
+and if its haplotype combination is compatible with that of its parents
+assuming Double Reduction to be impossible (noDR) or possible (withDR);
+NA indicates that the individual itself or both its parents have no
+haplotypes assigned.
Array parents_arr gives information for each individual that is a parent,
-whether it has complete marker data and an assigned haplotype combination,
-and how many of its progeny are compatible with it. For
-details see the help (?pedigreeHapdosCheck).
+whether it has complete marker data (par_mrk) and an assigned haplotype
+combination (par_hap), how many of its progeny have complete marker data (mrk)
+and haplotype data (hap), and how many of its progeny are compatible
+with it, assuming Double Reduction to be impossible(nonDRmatch) or
+possible (DRmatch). For details see the help (?pedigreeHapdosCheck).
### Summary statistics
The results of overviewByFS and pedigreeHapdosCheck can be summarized using
@@ -298,13 +315,15 @@ function calcStatistics:
```{r}
cst <- calcStatistics(pedchk=pedchk, ovwFS=ovwFS)
knitr::kable(cst$pedstats)
-knitr::kable(cst$FSstats, digits=c(0,0,0,2,2))
+knitr::kable(cst$FSstats, digits=c(0,0,0,2,2,2))
#
```
calcStatistics returns a list with two matrices. Matrix pedstats gives per
-haploblock the totals over all individuals from pedchk\$ped_arr. Matrix FSstats
-gives the totals or means per FS family over all haploblocks from ovwFS\$ovw.
-For details see the help (?calcStatistics).
+haploblock the totals over all individuals from pedchk\$ped_arr.
+The total of columns match.NA + noDR.TRUE + noDR.FALSE and
+of match.NA + withDR.TRUE + withDR.FALSE is the total number of individuals.
+Matrix FSstats gives the totals or means per FS family over all haploblocks
+from ovwFS\$ovw. For details see the help (?calcStatistics).
### Number of markers vs number of haplotypes
Finally we can get a table of the number of haplotypes vs the number of markers
@@ -312,13 +331,13 @@ per haploblock:
```{r}
calcMrkHaptable(ovwFS=ovwFS)
```
-In this small example all haploblocks have 4 markers, so the table contains
-only one row.
+In this example there are 5 haploblocks with 4 markers and 2 haploblocks with 5
+markers, so the table contains two rows.
### Segregation in one FS, one haploblock
With function showOneFS it is possible to study the segregation of one
haploblock in one FS family. Both the segregation of markers and of
-haplotypes is computed.
+haplotypes is shown.
```{r}
# show the segregation for FS number 1 (FSnr=1) and
@@ -334,70 +353,65 @@ first two columns represent the parents. Their column names are the names of
the parents in (brackets), and their first row does not contain a frequency
but the marker dosage IDs of the parents. All following columns have as names
one of the marker dosage IDs occurring in the FS, and the first row has the
-number of FS individuals with that marker dosage ID. In all columns, the next
-rows contain the dosages of the markers or the haplotypes.
+number of FS individuals with that marker dosage ID. (A marker dosage ID or
+mrkdid is a number encoding the dosages of all the SNP markers in an
+individual). In all columns, the next rows contain the dosages of the markers
+or the haplotypes.
The third element of the result is a matrix listing the composition of the
haplotypes present, in terms of their marker alleles (with a 0 being the
reference allele and a 1 the alternative allele).
# Remarks
-PolyHaplotyper creates and uses
-two or three variables in the global environment (.GlobalEnv): ahcinfo, ahclist
-and possibly ahccompletelist. These variables are used to manage calculated
-lists with all possible haplotype combinations that result in the same marker
+PolyHaplotyper uses one or two pre-calculated lists, ahclist and
+ahccompletelist, that contain
+all possible haplotype combinations that result in the same marker
dosage combination. These may take a lot of time to calculate and are completely
-reusable between haploblocks and runs of inferHaplotypes, so it would be a
-waste to discard them.
+reusable between haploblocks and runs of inferHaplotypes.
-The ahclist and ahccompletelist lists are also stored
+The ahclist and ahccompletelist lists are stored
as files with names like ahclist_6x.RData and ahccompletelist_6x.RData (where
-6x indicates the ploidy).
-
-The ahclist lists store the haplotype combinations for any marker dosage
-combination that is encountered in the data. They are modified "on the fly"
-and store only results for haploblock sizes and ploidies for which no
-ahccompletelist is available.
-
-The ahccompletelist lists store the possible haplotype
-combinations for ALL marker dosage combinations, for haplotypes from 1 marker
-up to some maximum. These lists must be pre-calculated. They
-are only used but not changed by inferHaplotypes, and speed up the
+6x indicates the ploidy). These files must be either in the working
+directory, or their location must be specified in the call to inferHaplotypes
+by setting parameter ahcdir. If these files are not present the haplotyping
+will proceed normally, but the computation of the ahc data will add very
+much to the processing time; a message to that effect will be shown.
+
+An ahccompletelist stores the possible haplotype
+combinations for ALL marker dosage combinations, for haploblocks from 1 marker
+up to some maximum. These lists are pre-calculated: they
+are only used but not changed by inferHaplotypes. They speed up the
calculations enormously, but it takes considerable time to calculate them
and above 7 or 8 markers (in tetraploids) or 6 markers (in hexaploids) they
become too large. They are useful if haplotyping is a recurrent activity.
-The package does not include ahccompletelist files because of their size.
+An ahclist stores the haplotype combinations for any marker dosage
+combination that is encountered in the data, and that cannot be found in the
+ahccompletelist (because that list is not available, or is limited to haploblocks
+with fewer markers). These lists are created or extended "on the fly" and
+then (re-)saved to file.
+
+The package does not include ahclist or ahccompletelist files because of their size.
The demo in this vignette doesn't need them as the number of haploblocks
-is small and the calculation of the ahclist data requires very little time.
-After running inferHaplotypes you will find file ahclist_6x.RData in the
+is small and the calculation of the ahclist data requires little time.
+After running the vignette example you will find file ahclist_6x.RData in the
working directory.
-The ahccompletelist files can be computed using function build_ahccompletelist
-completeAllHaploComb:
+The ahccompletelist files can be computed using function build_ahccompletelist:
```{r eval=FALSE}
build_ahccompletelist(ploidy=6, maxmrk=5, overwrite=FALSE)
```
-Apart from the file ahccompletelist_4x.RData this will also generate
-separate files for each number of markers, with names ahc_4x_nmrk1.RData ...
-ahc_4x_nmrk5.RData; these can be deleted.
+Apart from the file ahccompletelist_6x.RData this will also generate
+separate files for each number of markers, with names ahc_6x_nmrk1.RData ...
+ahc_6x_nmrk5.RData; these can be deleted.
Alternatively the following ahccompletelists are available upon request
-from roeland.voorrips@wur.nl
+from roeland.voorrips@wur.nl.
ploidy 4, up to 7 markers: ahccompletelist_4x.RData, 41 MB
ploidy 4, up to 8 markers: ahccompletelist_4x.RData, 769 MB
ploidy 6, up to 5 markers: ahccompletelist_6x.RData, 9 MB
ploidy 6, up to 6 markers: ahccompletelist_6x.RData, 515 MB
-The ahccompletelist files and ahclist files must be copied to the working
-directory to be accessible by inferHaplotypes and other functions in
-the PolyHaplotyper package. Alternatively they can be stored in any other
-directory, say "MyPath/MyAhcdir", and an option must be set before
-calling inferHaplotypes:
-```{r eval=FALSE}
-options(PolyHaplotyper_ahcdir = "MyPath/MyAhcdir")
-```
-
For those interested in some background on these lists:
An ahclist for a given ploidy has one item for each number of markers per
@@ -407,7 +421,7 @@ for which the set of haplotype combinations has been calculated. The name of
such an item is the mrkdid (marker dosage ID) and the item contains a matrix
with (ploidy) rows and one column per possible haplotype combination,
with the numbers of the haplotypes that are present (e.g. a column with numbers
-0, 17, 17, 24 means that haplotypes 0 and 24 are present in simplex and
+1, 17, 17, 24 means that haplotypes 1 and 24 are present in simplex and
haplotype 17 is present in duplex, in this tetraploid combination).
The order of the mrkdid items in the sublist for a certain number of markers is
the order in which they are calculated, and they are accessed by
diff --git a/vignettes/ahclist_6x.RData b/vignettes/ahclist_6x.RData
deleted file mode 100644
index 33c9e61a1cd7e9a06b693a16045d316326d877b9..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001
literal 4540
zcmb2|=3oE=wxhA9yHt*g`j++|(^b6PX=AoIqb()I<g5WtWdAyLjyr7ijp2v1IinLA
zm-Ozr;HCXI-|&~zWBrff5$oihpY5=;R6pN&YRyE)eU9f>N~^!qj++@;|L5@X)%^B<
zH>_IK{p!=j)Y8Jr9J^f>|L#{`UHqzf@t>Q|>;HZZkFWc^J^6Du|NpW-@Atpl-`&3d
zZ^4xP7N?Hy|Ly<vSH0r<^~Y38N`Fq*-*4A7e}7tW$oF6K@3q9pD&DX8_f`CUS&3f5
zIcbxS@A;Q@{P_5`HvXT*wduEy^VHgIK7RF2yjRWLldt{4mi@~X7XP%eliSYUN-2N)
z)!%Jv?KbT+ZN0P3eP`^ML%MHP?_gXmxuZ7c_Ub9php)5P+fDO6ey+@Nmtmpe{QjcN
zNeAA&K6d2n@w4kQ%@@~3nx1{P+O)7#J$~koj_uLH_xD(qS^qnp`}E3|wP)ArK6&=Y
z*?w~_@298-k$aQ(IqY5cUTVUb@)h$R-n6iduYYhQF>-g_y(ydK*@!2n|M*xH6w;pB
ztNqk0?|EgGY{Be+*J5SyYctzBs;9=Az9|1}_#v-O)1Hn0ihbVg<88V2o1gcGJ!pEY
zE!=PF@laaA`Y=2HmZ^u7>z}@UusDhH=Z+s2&L{C#KCYNqXU4uEBGpv(<B7v(TiH}j
zvtBv>!Svwtr+&68*_4*QG%eVmCV!^*?!EW_T06NdAJ3Zejq@|#(Wz$_el9A0YbL+F
z_3YcL7VD0??^(HV!Cdd->;6>m{&?QxJ~MZ2<=g7p+umBNJFah~SyLF__%G=GUJ;8g
za}Ru<-~DTKjqD%q<jllJ*$@0b*!Hj7Se5;#d4Jd2<qPINsQwpRobFhlQ_ivPe++ZQ
zwAlyun@qoG#}hC7ui<^=oG|UWw`U%2a;vEIopqeO)^L9~OZ~B(MHj@oD~;;s7q64~
zvh<Ve8|}~KuWrng&u$9;7@roMzV6H4b6KV*PnN$bH_cyk|LPv&(l=+`FR=^SKR+${
z`1`}z?w5{RzLI;({pQ;3hbDP`J-p|6IF~2oaUbxNKAl-^zDB?7-_PJ%k#o;>ZCtVW
z#{47a0?b3~<}XgWwC}{dmRsV-7rlwGD*Iv`S<Sa3Li${UaVgJ!i^4bh4Cc8yemwba
zHy*g7wR7Um8D|eNM~61`XGa@dUv0VLy?Vu&x4T&%h|DjHxBi~AZ~a2G6B}1ZF>O7(
zPI2kJG*hpu6J8u&tGIOEKE^qbm(_m!y1#9G^nsNd8RFI+P0dM=J#b?+OTBs6$(wV}
z?B3zbDOdcQ>00Pi?mStw^iHF1hw}6ER#@?;_m%LoRu@*sX{^=Fy0`MA1>^nCGgg;P
z+#Ff9%{^PTJAF&=hxx0gEs72Ov!QaQVuhLh^lN41_tvxi4}RVhE!<yMcAax+G4DEO
zsYZTDn`LKH_v>0No-pUwZ!N3E6KWsstCQN^9l7J6Wb4^n!<n~@Z|4=5X&)E*_@<Bb
zJY$UgHMP1s{ntKU-L`hawmJ64YV05Fa=I!0`rDo4{<<@@S9t2a<hEw}iSD*YKbyPy
z#c|u{nR%;^=gAviDcHaM<1Tam_<+gxPwt5GdG}g*_xkB-e6J2HDR;QDPsDtV{_^7L
zJJ;t=&_DiEMLsk|X3hEKEdF-p2l;M0t`#nmc*7juasKFyjci6Xx7{_PH&&$AocK5C
zlcJ%<&m9%Z<n7EKuE{j~F!#gPTGr&kk9j!;w|DKcGwtVCulfB`^G5Ax!@Sh5$4`EI
zQvQ^GTJyU5p?Ze<wQYXgJv!Uc?%?~(+?6>;a-RNLWU?;izG3v#qYJ<9%Q9UTU1Xix
z_wAaN{~kNGJ3jX3Hia4PJiBRLNw!h+;wev~ZdP47%UXVJ-eLE3d!K!o80~y#>$y30
zMW0X1Wq$VM;?K`7X6CkUyTI2G{_Anvn`1}XkNjS6Ip^tF?spMMZ?>z|_01?O&Z#<M
zv#r+HowxA#)oF!Y6|+-{x34eTB-@t#INa&ct%!$@a|E*2U+Pcq;i>1XiLu$)@A^k4
zJ2_Q3P5Zj|54ok;+T!`Yl7e_^9(^@@kap~<PImJ3{u90r)US4BujjuwJ$(Mj-HXy+
zXs=y-=Z2f#^X~SAvggu2Ka2`LDzZ8|_3PP<ht(%Va|abV+eANqJcIfCwVyoae<av9
zs-Ne$Sh(5L;InKr^X`3Ne+~Yz_$NI+Jo(YJv-isb&)@yoUuYd~Tc1{(&ApSgPO%{8
zjh0ws&)hWa_v^0Bm43G3a9+IbjYO&Kn{01wx%2TF!#f+3_Oo~QMID&A(W?9Geb+ZT
zYg=XS#LPX+mnS<@`P-d+w<LcZdM?yD=Slqb^~)yz@Ohkh=hQd5u3*jW&wB1mn0%V|
z!{#@q&&k`@#Xo*LYq|dUw#&b=?&+yNUGVFj=ACqHnK^TPcJAM{b8AELj}6j&OXpR-
z`C%J%YW~J!K{htGtIq50NUltM^z8F)d;KXzv;Y0MVLUtfwr9Cq=f)kMJzn03pS9zG
zah`Vm{6|Hzlh?lZc2(?F;;*Q{;_iJS+u5GT+_|&s{W0k`;^`~@cg@};UBtb|+P`Iu
zuid(}1?%6O6Mnp{^ws`zXMcXvxB3@rZ`LqhE?NHe*F{sL51X^t+lf8j&D3^3|NNt%
zwdprz_Mg($w_E%E-GVzeBF<EmT-(($=ay4mWBUEqr;f$Vu9>p-#?85r{pX^lpNKdY
zG(E=NjK$tPc9nIkVzeT&_?$e}eTo&Ic5?3gbFy}M;s+b$$8GO3bXe@Z9FK@=a;fqC
zVPNy=;PVL*c0V2z`}OcY?w{}c_sz}~Vg)_@hq|AyQ@b-g_r116uy)FW=Gk3KrM^r*
zHt&@3v+SenKX-3jlC++2WfK2KrD+}O=Y9{cle*&jWYLa9TQ<u~`H;guV)GW|ojdUM
z_rcw4UD4N-Yxgmz?#o{A;m5kh)9OFC?fG8qJ9NCt^J&nDU6U*Fk9|m#pIUuF{`~|=
zyH8F}vra6VeBRgQnvtmZLizcp6MjflXxMMAobqsD|DnfKVMkqGq@VXPwy6$_de}PG
zV{O4c^`HK=YaRG}!t$RdJfC&w`dpp0u6eS@VuCBqn6~Cyy}Y(-*W{&J`M$?>^+!ub
z>S_GVd&^amnQK@)|H&c!wVQu=6*xxDPUd^LPO$LJ(v9iA@-7vh__`zSr0)auwP)|t
zN)-kFUh(pdSo3%GR#~?9{JSc)&x>Kcw{_}nu`Gu2sDD-EQ{K(<EmZyy`HoR^`2<6;
z@_DZIyOPTvYkf$I@m;0+?%upY;UAvoSDQTEb#~vL=T-&lD{jrZ_%hqNB0a2^|9*04
zq13OtrL1-FGrG&Ki*|$;);fRa5&z|L)%(c&ZO^vX`4>O1e$q7QX4k{35e)V(@3HS(
zxzQkx#eQYjj7n?gf3M1uXUa<4vp)X$TJozuJEqJ|H_x|yHm|d#<h+Hz`*p|f7`INy
zt9AUdJEDK*^WZM&gXcf}>v$w`Y3_5UU7yAJ?)GYz-+T3B$=4ay{swtW_TIgwJAI@2
zXKwwpruMA;x`)>PPEC8a^R0#X`Tcpab1w5`{mEe8tac`eee&%kH;Rwj?_IzDez_Fe
zq^G8jUiod3I{IOan<C#=;~Ec}$A=3Q_A&nDk+=D<&ii4Y{jHU3oxcS8e>rd3DP5>|
zbJhkX$v-!D$^JC>ZCvvry;|nU**E7utdKu0TXUnBMMLoXk9L=PJ@&8i9hB;q-9M`R
z_woJpzolcN_5PXu=<S@KJF$b?^VM&M`kvmaSyxuSJ0>9i>!gk8gT?Qp1NQIocHXnY
zTlzz$)yr#Hh2PD7^zwiG6aHea){*wQuXmU~<y)QnCbc5oR=x4p<Ki6Iri#+?6F)ZZ
z<$vt|fw^9>-fVlE=%YaU$orRmgk(2Qc(iWHvE`TgKfYWSdZk$X#Qv8mw^uK!`#9%R
z!k&fkYv-=x`nWLOmc6;K`D5n~DV5)Tp~uvJY&KTC(*N__ukN{`3i4lN_#dA?w?m@7
zN&ZjU#NB3fK^^XYE5gqIFgXzavGC(U`NzT^XFSi!jc1lG*Kc@ybhpHn{u+sSFOTn*
zUON8)`+mjNUpMbPxpidmUs1asn=2A3#Q1yanJ3*nyRfd6|F5Y>*}ULJp1*J22!7rC
zv-_d5y|uCXy*2X|*M-(g$o*81tK3lI@|V}6Y<kx(QMu>4C%n-SEL1%CJ>}xhxo`h!
zKY04c**-eHHRbc}YO#aYykD(pxO;Iy-L|--KhE~n(vE*~{w_YJRIqp9vD5En1lULC
z??_i}SpIpd)$N17)_jPVFFEh*bsqW0A+IGt-cIMwI<WaG|E<l8C(mta;-Bq5wN}LL
z$Lg=5()!WyivQVEs>N8pYTH~sKQA`%m!#ZlwMT~Uo%cumSCISq_?P2R^<&2Wo-I2i
z{zupH$>Di#9slahdwJZdrSkNT`$96Ok55=J|6}^E$Unw4FOFA)RrL5jxi4Vz$?2Y+
z+}Gq9Et^;Ug$n#_@d7fBjqA$(&MnmccK)~DVTWxlO6%S&sovJJaqWk*vy-=;H;c9T
zq;D%-cz#n%&5QE$eK+p0d|iF((1#Pp@19=0xruGb?(L-?HE!Np!W8qO{Q9@0)r?Qh
z|Jaco_~hZbzfu$4MV&Zm{MRDSZ2Hb(`NzRgtMfVg9v;trzcsLW|E~8YUnX5WB358p
zQ222T+tvo>zpJ*bp4Ie&?Tc~EhX7O7k1OWiT*-E_XvW8M%VR%wS3EiYA@7{RyqE4Z
z&L6V$G9KKY85`Pr_`IPVi+|DQ<hiro|4lAPzs?oeC~`iF_uP(%zJ>o*JJq}pKR<`}
zoWq@{CGS27_E(D@UH)=O^y?$BcD9qlf9#t4Do^ewOZ``wf*18!kDTjk`tF`pmixQ<
z!@>D?&l>;Zs%M>_Y%2b>UF7^EzDMq7YeUtK##;Bvm*<6+KR==VWA*br&i2oyKa&4-
zHAeo|)mV#1&ikwWFJ!Wrx6<HmMtv*)-#70pmsFeIVd|TE)Zp(4_m|=2dpr35zR8>V
z@aEktkMHj`v&wubyq@;voxumu7svlN)x1dkdHTb``FB~_;?|y(ko#&}<NRUGyQF=q
zx8~jTEH!WYabf*EeJ7h|`j77Jo^~L4-wsf?yo+_^-Z9~oZRnBBo2!q>AMY;hdzt=Y
z){FGYm-FKHFv&b_uGq5wR`!Cruk&;t&9puKU24DiIYztP?@hjJ5;ZC?EeQK@X>09e
z&+^=uUDMqDEzM>w*M0KHd4JUH6}QfPSTX;hbW?b7>Vwd=WyQfC3qLBEd%nBkdf)N+
zK7Gm02`|6fO|9#n{QF)y>z}I@1^?H7C{@3j&)9b4Z|9-MRpAc{|E~Jrt3H4Bj|=g?
zWXrSuH`L_$)bQF&nSb-&oQKW)w<bTlS>^ZAy;{rW*5#wz*1kXZcCr4|wSOWaRJ-QG
z8TL2DXO{7t&sp`sXMcM1bJ;%&b0>03SHJu&J>}QAwUhJqDBb&VcgcIX^g8X=ulIV^
z3R<7JCH^w`+CCTg%fG6AY*&7E@$mP@A@A!y|GXYw6TLU?M2LM!{hkQfMO*K=WpD1g
zw(tLK#eZcx@2m;h{Qq^;cCWL}_uq)d<Xu}i^WU!Mdtc99)iSN--Mu8|{{OERH~HSp
z6`KDx@p8Rh>2`0O^4Y#Y8}l}QGs=Fy>6_haqt|uw?RMW=yZM`1_Wh-0efQP>&aOMd
z|8`0J=IFGyJezN)ziqj%yFIdc`c~n(?Ou1vZtq?*x7O_I)~xq@EBEGZE<2xFYT|Ou
z@5fW4-tyTu)A#$%;g@{weC+3hLc0_je#!HP%clSMqWI;)JHyq=ombBNSa0<2b<cmR
z$NT;M^uNu$=J$8q`%UN5wneF*J6?9Ia(?gM2_<QU+p^4WM!cJ(A2z!sQg+iu86&&Z
zYwm}B`~FocJKe_orUV$J2ivT__WMlUZ=dIOn-AMvu0QeP%j3N<-!sc|^n#72_LpjD
z)oRJj{(Z&g>r(l{Ph<UG-gsl06D+RD*MD}OU(W9{HfF+Rn${|B6i>e>+j{O+WT<Cc
z#JX*olQy5Q+_rgty5J_iz`b|Ue&lv_JvDPZYPdUif4j)T<1f78vTvteJGaTY``emF
zvdRT>w;J1)N3!lXU6K~LVVC!hP%+WfI$B|oPuAbQa=mJ*>n*G9Ev0fNBbTJ!xfXeC
z&E_>HwLWXDKCs$E+jfib-*%q)Hwz}`?1;IS`p@@I|Jx_syHnbqn%(?!<6O?|XE&!8
f-g_cbly~z6e|t}Kbm`~6{~2Dhs_CCR&A<Qvv@fYT
--
GitLab