From ca61ea168209c6a7c0e89e9e343d8d0a1815ce0a Mon Sep 17 00:00:00 2001 From: jendelman Date: Wed, 24 Jul 2024 14:43:27 -0500 Subject: [PATCH] v1.11 --- DESCRIPTION | 2 +- NEWS | 3 ++ R/Stage2.R | 24 +++++++---- R/blup_prep.R | 13 +++--- R/class_var.R | 4 -- R/dominance.R | 42 ++++++------------ docs/Vignette1.html | 18 ++++---- docs/Vignette3.html | 92 ++++++++++++++++++++++------------------ docs/manual.pdf | Bin 123222 -> 123003 bytes man/class_var-class.Rd | 4 -- man/dominance.Rd | 9 +--- vignettes/Vignette3.Rmd | 10 ++++- 12 files changed, 111 insertions(+), 110 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 3d7a26a..441e1f1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: StageWise Title: Two-stage analysis of multi-environment trials for genomic selection and GWAS -Version: 1.10 +Version: 1.11 Author: Jeffrey B. Endelman Maintainer: Jeffrey Endelman Description: Fully efficient, two stage analysis of multi-environment trials, including directional dominance and multi-trait genomic selection diff --git a/NEWS b/NEWS index afb1a17..4601c36 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,6 @@ +Changes in 1.11 +* dominance function is now scaled to population, not just var component + Changes in 1.10 * Added support for multi-population/multi-ploidy analysis * Fixed error in Stage1 related to multiple experiments per env diff --git a/R/Stage2.R b/R/Stage2.R index e48f6ef..d08453c 100644 --- a/R/Stage2.R +++ b/R/Stage2.R @@ -118,8 +118,8 @@ Stage2 <- function(data,vcov=NULL,geno=NULL,fix.eff.marker=NULL, vars=new(Class="class_var", geno1=geno1, geno2=geno2, resid=resid.vc, B=as.matrix(B), model=result[[1]]$vars@model, - diagG=mean(sapply(result,function(x){x$vars@diagG})), - diagD=mean(sapply(result,function(x){x$vars@diagD})), + #diagG=mean(sapply(result,function(x){x$vars@diagG})), + #diagD=mean(sapply(result,function(x){x$vars@diagD})), vars=array(NA,dim=c(0,0,0)), fix.eff.marker=result[[1]]$vars@fix.eff.marker))) } @@ -134,7 +134,7 @@ Stage2 <- function(data,vcov=NULL,geno=NULL,fix.eff.marker=NULL, data <- data[!(data$env.id %in% data$env.id[missing]),] } - diagG <- diagD <- numeric(0) + #diagG <- diagD <- numeric(0) dom <- NULL if (!is.null(geno)) { stopifnot(is(geno,"class_geno")) @@ -147,7 +147,7 @@ Stage2 <- function(data,vcov=NULL,geno=NULL,fix.eff.marker=NULL, id.weights <- table(factor(data$id,levels=id)) .GlobalEnv$asremlG <- geno@G[id,id] meanG <- as.numeric(pvar(V=.GlobalEnv$asremlG,weights=id.weights)) - diagG <- mean(diag(.GlobalEnv$asremlG)) + #diagG <- mean(diag(.GlobalEnv$asremlG)) if (non.add=="dom") { .GlobalEnv$asremlD <- geno@D[id,id] @@ -158,7 +158,7 @@ Stage2 <- function(data,vcov=NULL,geno=NULL,fix.eff.marker=NULL, names(dom.covariate) <- id data$dom <- dom.covariate[data$id] - diagD <- mean(diag(.GlobalEnv$asremlD)) + #diagD <- mean(diag(.GlobalEnv$asremlD)) dom <- "dom" } } else { @@ -696,8 +696,8 @@ Stage2 <- function(data,vcov=NULL,geno=NULL,fix.eff.marker=NULL, } out$vars <- new(Class="class_var",geno1=geno1.vc,geno2=geno2.vc,model=model, - resid=resid.vc,diagG=diagG,diagD=diagD, - vars=vars,B=B,fix.eff.marker=fix.eff.marker) + resid=resid.vc,vars=vars, + B=B,fix.eff.marker=fix.eff.marker) #random effects if (n.trait==1) { @@ -783,12 +783,18 @@ Stage2 <- function(data,vcov=NULL,geno=NULL,fix.eff.marker=NULL, } } - if (!is.null(geno)) + if (!is.null(geno)) { + out$params$scale <- list(add=as.numeric(pvar(V=.GlobalEnv$asremlG))) rm("asremlG",envir = .GlobalEnv) + } if (!is.null(vcov)) rm("asremlOmega",envir = .GlobalEnv) - if (non.add=="dom") + if (non.add=="dom") { + out$params$scale <- c(out$params$scale, + list(dom=as.numeric(pvar(V=.GlobalEnv$asremlD)), + heterosis=as.numeric(pvar(mu=dom.covariate)))) rm("asremlD",envir = .GlobalEnv) + } return(out) } diff --git a/R/blup_prep.R b/R/blup_prep.R index a802e53..7b174d3 100644 --- a/R/blup_prep.R +++ b/R/blup_prep.R @@ -32,6 +32,8 @@ blup_prep <- function(data,vcov=NULL,geno=NULL,vars,mask=NULL,method=NULL) { stopifnot(method %in% c("MME","VINV")) } + if (vars@model > 0) + stopifnot(inherits(geno,"class_geno")) if (vars@model==3L) stopifnot(is(geno,"class_genoD")) @@ -43,9 +45,9 @@ blup_prep <- function(data,vcov=NULL,geno=NULL,vars,mask=NULL,method=NULL) { stopifnot(n.trait > 1) } - if (length(vars@diagG)>0) { - stopifnot(!is.null(geno)) - } + #if (length(vars@diagG)>0) { + # stopifnot(!is.null(geno)) + #} data <- data[!is.na(data$BLUE),] @@ -65,8 +67,9 @@ blup_prep <- function(data,vcov=NULL,geno=NULL,vars,mask=NULL,method=NULL) { } if (!is.null(geno)) { - stopifnot(is(geno,"class_geno")) - stopifnot(length(vars@diagG)>0) + stopifnot(inherits(geno,"class_geno")) + #stopifnot(length(vars@diagG)>0) + stopifnot(vars@model > 0) id <- intersect(data$id,rownames(geno@G)) data <- data[data$id %in% id,] id <- rownames(geno@G) diff --git a/R/class_var.R b/R/class_var.R index 02124d0..33b0e91 100644 --- a/R/class_var.R +++ b/R/class_var.R @@ -4,8 +4,6 @@ #' @slot geno2 second genetic effect #' @slot model 0=no markers, 1=add, 2=add+g.resid, 3=add+dom #' @slot resid residual -#' @slot diagG average diagonal element of the G matrix -#' @slot diagD average diagonal element of the D matrix #' @slot vars variances for reporting #' @slot B var-cov matrix of fixed effects for gain #' @slot fix.eff.marker names of fixed effect markers @@ -17,8 +15,6 @@ class_var <- setClass("class_var",slots=c(geno1="Matrix", geno2="Matrix", model="integer", resid="Matrix", - diagG="numeric", - diagD="numeric", vars="array", B="matrix", fix.eff.marker="character")) diff --git a/R/dominance.R b/R/dominance.R index aafa44c..a5cfe7a 100644 --- a/R/dominance.R +++ b/R/dominance.R @@ -2,22 +2,22 @@ #' #' Report dominance parameters #' -#' The dominance variance (Vd) and baseline heterosis (b) are quantified relative to additive variance (Va) and std. dev. (SDa), respectively. For single traits, the estimate and SE of the ratios are calculated based on the delta method (Rice 2007, p. 162-166). For a multi-trait/loc model, the result is just the ratio of the estimates, with \code{index.coeff} specifying the coefficients of the standardized true values (see also \code{\link{blup}}) and \code{gamma} specifying the relative weight for non-additive to additive genetic merit (see also \code{\link{gain}}). +#' The dominance variance (Vd) and baseline heterosis (b) are quantified relative to additive variance (Va) and std. dev. (SDa), respectively. As of v1.11, the variances are scaled to the population (previously, it was just the variance components). For a multi-trait/loc model, \code{index.coeff} specifies the coefficients of the standardized true values (see also \code{\link{blup}}), with \code{gamma} indicating the relative weight of non-additive to additive genetic merit for the standardization (see also \code{\link{gain}}). #' #' @param params list returned by \code{\link{Stage2}} -#' @param digits number of digits for rounding #' @param index.coeff merit index coefficients #' @param gamma contribution of non-additive values for genetic merit #' -#' @return data frame with estimates and SE +#' @return data frame with estimates #' -#' @references Rice JA (2007) Mathematical Statistics and Data Analysis, 3rd ed. Duxbury, Pacific Grove. #' @export -dominance <- function(params, digits=2, index.coeff=NULL, gamma=0) { +dominance <- function(params, index.coeff=NULL, gamma=0) { if (!is.element("heterosis",names(params))) stop("Dominance model was not used") + stopifnot("scale" %in% names(params)) + id <- grep("dominance",params$vc$name,fixed=T) ia <- grep("additive",params$vc$name,fixed=T) vc <- params$vc[,-1] @@ -28,7 +28,6 @@ dominance <- function(params, digits=2, index.coeff=NULL, gamma=0) { stopifnot(!is.null(index.coeff)) if (colnames(params$heterosis)[1]=="trait") { traits <- params$heterosis$trait - #Vd <- Va <- matrix(0,nrow=nhet,ncol=nhet,dimnames=list(traits,traits)) Va <- f.cov.trait(vc[ia,],traits,us=(nhet>2)) Vd <- f.cov.trait(vc[id,],traits,us=(nhet>2)) } else { @@ -49,31 +48,16 @@ dominance <- function(params, digits=2, index.coeff=NULL, gamma=0) { Va.mean <- as.numeric(crossprod(index.coeff,Va%*%index.coeff)) Vd.mean <- as.numeric(crossprod(index.coeff,Vd%*%index.coeff)) b.mean <- sum(params$heterosis$estimate*index.coeff) - - ans <- data.frame(ratio=c("Vd/Va","b/SDa"), - estimate=round(c(Vd.mean/Va.mean,b.mean/sqrt(Va.mean)),digits)) - return(ans) - } else { - Vd.mean <- params$vc$estimate[id] - Vd.var <- params$vc$SE[id]^2 Va.mean <- params$vc$estimate[ia] - Va.var <- params$vc$SE[ia]^2 - SDa.mean <- sqrt(Va.mean)*(1-2*Va.var/Va.mean^2) - SDa.var <- Va.var/(4*Va.mean) b.mean <- params$heterosis$estimate - b.var <- params$heterosis$SE^2 - - ratio1.mean <- Vd.mean/Va.mean*(1+Va.var/Va.mean^2) - ratio1.SE <- sqrt(Va.var*Vd.mean^2/Va.mean^2+Vd.var)/Va.mean - ratio2.mean <- b.mean/SDa.mean*(1+SDa.var/SDa.mean^2) - ratio2.SE <- sqrt(SDa.var*b.mean^2/SDa.mean^2+b.var)/SDa.mean - - ans <- data.frame(ratio=c("Vd/Va","b/SDa"), - estimate=round(c(ratio1.mean,ratio2.mean),digits), - SE=round(c(ratio1.SE,ratio2.SE),digits)) - - return(ans) } -} + + Vd <- Vd.mean*params$scale$dom + b.mean^2*params$scale$heterosis + Va <- Va.mean*params$scale$add + + ans <- data.frame(ratio=c("Vd/Va","b/SDa"), + estimate=c(Vd/Va, b.mean/sqrt(Va))) + return(ans) +} diff --git a/docs/Vignette1.html b/docs/Vignette1.html index e78c39e..e5be5e2 100644 --- a/docs/Vignette1.html +++ b/docs/Vignette1.html @@ -609,7 +609,7 @@

Stage 1

library(StageWise)
 ans1a <- Stage1(filename=pheno1a.file,traits="total.yield",
                 effects=effects,solver="asreml")
-
## Online License checked out Sun Sep  3 11:48:17 2023
+
## Online License checked out Wed Jul 24 12:47:50 2024

The Stage1 function returns a list with several results. List element “blue” is a data frame of the individual BLUEs per environment. Element “fit” contains the broad-sense heritability on a @@ -673,7 +673,7 @@

Stage 1

different models. A figure showing the 2D spline and spatial distribution of residuals is also returned when SpATS is used:

model2$resid$spatial$advanced
-

+

Before proceeding to Stage 2, the BLUEs and variance-covariance matrices from the 2015-19 analysis and 2020 analysis need to be combined. (If other software is used for Stage 1, it can be incorporated @@ -707,7 +707,7 @@

Marker data

## 2 solcap_snp_c2_36608 chr01 508800 1 0 2.00 ## 3 solcap_snp_c2_36615 chr01 510745 1 0 2.03 ## 4 solcap_snp_c2_36658 chr01 527068 4 3 3.00 -
geno <- read_geno(filename=geno.file, ploidy=4, map=TRUE, min.minor.allele=5, 
+
geno <- read_geno(filename=geno.file, ploidy=4L, map=TRUE, min.minor.allele=5, 
                   dominance=T)

 
## Minor allele threshold = 5 genotypes
 ## Number of markers = 12242
@@ -741,7 +741,7 @@ 
Marker data

 
##        F.G        F.D 
 ## 0.03414348 0.03352837

 
ggplot(x,aes(x=F.G,y=F.D)) + geom_hex()

-

+

 
The above code shows there is relatively little variation for \(F\) (sd \(\approx\) 0.03), which is relevant when
 interpreting the results of Stage 2.
@@ -954,9 +954,9 @@

Stage 2

##   estimate       SE
 ## 1 59.52934 17.83285
dominance(ans2d$params)
-
##   ratio estimate   SE
-## 1 Vd/Va     0.38 0.14
-## 2 b/SDa     8.00 2.49
+
##   ratio  estimate
+## 1 Vd/Va 0.4374054
+## 2 b/SDa 8.4843420

The standard error of the baseline heterosis estimate is large because of the limited variation for \(F\).

@@ -1053,7 +1053,7 @@

Pedigree data

ggplot(result) + geom_line(mapping=aes(x=w,y=h2)) + geom_line(mapping=aes(x=w,y=y2),colour="red") + scale_y_continuous(name="Genomic h2",sec.axis=sec_axis(trans~./axis.scaling,name="AIC",breaks=y2axis,labels=y2lab)) + theme_bw() + theme(axis.text.y.right=element_text(colour="red"),axis.title.y.right=element_text(colour="red")) + ggtitle("Blending G and A for Yield") -

+

Based on the above figure, the blending parameter w=0.4 is chosen, at which there is no longer much dominance variance. The optimal value for w will not be the same for all traits.

@@ -1300,7 +1300,7 @@

Marker effects and GWAS

gwas_threshold(geno[[1]], alpha=0.05, n.core=2)
## [1] 5.103417
manhattan_plot(gwas.ans, thresh=5.1, rotate.label=TRUE)
-

+

The GWAS peak on chr05 is near the gene CDF1, which is known to have a large effect on potato maturity (Kloosterman et al. 2013). The following code extracts the most significant marker diff --git a/docs/Vignette3.html b/docs/Vignette3.html index 08ad541..3aba9aa 100644 --- a/docs/Vignette3.html +++ b/docs/Vignette3.html @@ -379,7 +379,7 @@

Potato dataset

effects <- data.frame(name=c("block","stand.count"), fixed=c(FALSE,TRUE), factor=c(TRUE,FALSE))) -
## Online License checked out Mon Sep  4 08:35:40 2023
+
## Online License checked out Wed Jul 24 14:18:21 2024
names(ans1)
## [1] "blues" "vcov"  "fit"   "resid"

As with the single trait analysis, Stage1 returns a data @@ -443,12 +443,12 @@

Potato dataset

vars=ans2$vars)

For multi-trait analyses, the blup command was designed to predict total genetic merit, not individual traits, so the user must -supply index coefficients. The gain command allows breeders -to compare the expected response for each trait under different indices. -The argument “merit” specifies the contribution of each standardized -trait for total genetic merit. For example, if genetic gains for fry -color and yield are equally valuable, and maturity is ignored, the code -is

+supply index coefficients.

+

The gain command allows breeders to compare the expected +response for each trait under different indices. The argument “merit” +specifies the contribution of each standardized trait for total genetic +merit. For example, if genetic gains for fry color and yield are equally +valuable, and maturity is ignored, the code is

merit.coeff=c("total.yield"=1, "vine.maturity"=0, "fry.color"=1)
 
 gain1 <- gain(prep1, merit=merit.coeff)
@@ -597,6 +597,16 @@ 
Potato dataset

 ## 4 AF5445-2 37.71679 0.5367704
 ## 5 AF5450-7 37.05010 0.5219305
 ## 6 AF5484-3 36.03935 0.6062087
+

As explained in Vignette 1, the dominance function +expresses the inbreeding depression and dominance variance estimates +relative to the additive SD and variance estimates, respectively. This +also works for multi-trait models once the index coefficients are +specified. To be consistent with the above gain +calculation, gamma=1/3 is needed for tetraploid breeding values:

+
dominance(ans2$params, index.coeff=index.coeff, gamma=1/3)
+
##   ratio  estimate
+## 1 Vd/Va 0.4506758
+## 2 b/SDa 7.6863051

One last comment: groups of unrelated traits can be analyzed independently through blup_prep, and the outputs from this command can be combined as a list in blup to make the @@ -617,8 +627,8 @@

Genomic prediction with secondary traits

drought environments are distributed with the package. As with the potato dataset in Vignette 1, including the Stage 1 errors in Stage 2 lowers the AIC substantially.

-
data(wheat) #load the wheat data
-head(wheat.blues)
+
data(wheat) #load the wheat data
+head(wheat.blues)
##       env      id trait      BLUE
 ## 1 drought 6569128    GY  3.395834
 ## 2 drought 6569128    CT 36.950000
@@ -626,11 +636,11 @@ 
Genomic prediction with secondary traits

 ## 4 drought 6688880    CT 37.902333
 ## 5 drought 6688916    GY  3.459078
 ## 6 drought 6688916    CT 38.706000
-
ans2a <- Stage2(data=wheat.blues, vcov=wheat.vcov, geno=wheat.geno,
-                non.add="none")
-ans2b <- Stage2(data=wheat.blues, geno=wheat.geno, non.add="none")
-
-data.frame(vcov=c(TRUE,FALSE), AIC=c(ans2a$aic,ans2b$aic))
+
ans2a <- Stage2(data=wheat.blues, vcov=wheat.vcov, geno=wheat.geno,
+                non.add="none")
+ans2b <- Stage2(data=wheat.blues, geno=wheat.geno, non.add="none")
+
+data.frame(vcov=c(TRUE,FALSE), AIC=c(ans2a$aic,ans2b$aic))
##    vcov       AIC
 ## 1  TRUE -16.97154
 ## 2 FALSE 539.30801
@@ -639,7 +649,7 @@

Genomic prediction with secondary traits

was lower without it (non.add=“none”). Genomic heritability was 0.45-0.50 for yield and canopy temperature, with an additive genetic correlation of -0.81.

-
summary(ans2a$vars)
+
summary(ans2a$vars)
## $var
 ##                 GY   CT
 ## env          0.621 1.58
@@ -662,37 +672,37 @@ 
Genomic prediction with secondary traits

 marker-based selection (MBS, see Vignette 1). Since the goal is yield
 prediction, the index coefficients are 1 and 0 for GY and CT,
 respectively.

-
id <- unique(wheat.blues$id)
-N <- length(id)
-folds <- split(sample(id),cut(1:N,10))
-MBS <- NULL
-for (i in 1:10) {
-  prep <- blup_prep(wheat.blues, wheat.vcov, wheat.geno, ans2a$vars, 
-                    mask=data.frame(id=folds[[i]]))
-  pred <- blup(prep, geno=wheat.geno, what="BV", 
-               index.coeff=c(GY=1, CT=0))
-  MBS <- rbind(MBS, pred[pred$id %in% folds[[i]],])
-}

+
id <- unique(wheat.blues$id)
+N <- length(id)
+folds <- split(sample(id),cut(1:N,10))
+MBS <- NULL
+for (i in 1:10) {
+  prep <- blup_prep(wheat.blues, wheat.vcov, wheat.geno, ans2a$vars, 
+                    mask=data.frame(id=folds[[i]]))
+  pred <- blup(prep, geno=wheat.geno, what="BV", 
+               index.coeff=c(GY=1, CT=0))
+  MBS <- rbind(MBS, pred[pred$id %in% folds[[i]],])
+}

 
In the above code, the “mask” argument for blup_prep
 only has the variable “id”, which means that all Stage 1 BLUEs for those
 individuals are masked. To only mask grain yield and use CT as a
 secondary trait for marker-assisted selection (MAS), a second variable
 named “trait” is used.

-
MAS <- NULL
-for (i in 1:10) {
-  prep <- blup_prep(wheat.blues, wheat.vcov, wheat.geno, ans2a$vars, 
-                    mask=data.frame(id=folds[[i]], trait="GY"))
-  pred <- blup(prep, geno=wheat.geno, what="BV", 
-               index.coeff=c(GY=1, CT=0))
-  MAS <- rbind(MAS, pred[pred$id %in% folds[[i]],])
-}
-
-ans <- merge(MBS,MAS,by="id")
-
-library(ggplot2)
-ggplot(ans,aes(x=r2.x, y=r2.y)) + geom_hex() + coord_fixed(ratio=1) + geom_line(data=data.frame(x=c(0.2,0.8),y=c(0.2,0.8)),mapping=aes(x=x,y=y),linetype=2) +  ggtitle("Reliability") +
-  xlab("MBS") + ylab("MAS")

-

+
MAS <- NULL
+for (i in 1:10) {
+  prep <- blup_prep(wheat.blues, wheat.vcov, wheat.geno, ans2a$vars, 
+                    mask=data.frame(id=folds[[i]], trait="GY"))
+  pred <- blup(prep, geno=wheat.geno, what="BV", 
+               index.coeff=c(GY=1, CT=0))
+  MAS <- rbind(MAS, pred[pred$id %in% folds[[i]],])
+}
+
+ans <- merge(MBS,MAS,by="id")
+
+library(ggplot2)
+ggplot(ans,aes(x=r2.x, y=r2.y)) + geom_hex() + coord_fixed(ratio=1) + geom_line(data=data.frame(x=c(0.2,0.8),y=c(0.2,0.8)),mapping=aes(x=x,y=y),linetype=2) +  ggtitle("Reliability") +
+  xlab("MBS") + ylab("MAS")

+

 
The above figure shows that using CT increased the reliability of
 genomic prediction, by 0.2 on average.
diff --git a/docs/manual.pdf b/docs/manual.pdf index 1b55e10f44adcbcb66c991d5169facc4078e18c5..a60bb1e0b6ca0fc82429201e45beda0ffa38d738 100644 GIT binary patch delta 15986 zcmajFQIcJ@fm6cCKMr74Z z#gmaytAmgq{gCnGz#K`;c+?4sSpcs04R|7j>4%B38c}kp7SPdYxMOB1nPN!rac5h2 z`nKuHlDV|PoDwcz?%8L|w&H;OVCdpJX?AfDE-OOvB=^<2>n2>C5%Xk8Xt>*qGcU=G z#=K;sVz)at#;g2-?prN0%EZ=YmFiPcQ9H;}9z(3K$52j#)0S9Vvhw;~0ARs>C?3Wt z+M6WSC18@-2n^|gMf-_hfp5_t;X-5VxJX(!QSouuTfrEcd4CZ z`!2GL!(5W=VEdxr#WPt9aAd+-l0M>2_Q?6{1zs>+?!j9q&yfS7oJPRvkZvE76U;;P zy2l4-N{SJ#TfUMNh+}rlG{8cX_a!uD4@X3gnVx%%Z1Q4t;*#qg$GfH!Wq2LUYb!l< zmegi!%yR{S{V*nHRZdGdX|-qm*J+66Q^p{D6G{lRgCX8L0Q&J%`_I{dl%WJFm`ZXq zc28sU8dInPpyd!ui2vCmgtc+(D?vS3);s^&h`HkKb7@T5xBMQ?9}ss^`SViV^lVYr zt(H-1Vi2vy-QTZxeCMN9)_u`7@?4_Yb~OIr<;6k!>BV}+A1l}?&zRiy%L&>MJB}<( zjP`D8s_{H>x;4~g(Kb_li}uPKf5k_gzRW=cw#1Ef0cD4;2FE@mi;m+FiWrRn^ComV zeir_;lzkR(8*xeK2f!H>Z`~CfRfg=#W|Sbk?fdj&Z#WMwD74o=6%QccWJIhA?lgpx z7cKbxivqf+l?DnIT+}|uy9#1~&oIU#S8Xq{2meqYD zWb?_Ng_39ocrYzJ%Y7YBqMAbo_y=?<*GOHqST7@H!r=z`6|j%BM?1{u_*#l5CT4|0 z5b$~cwcL#Q5o}K)KiR&k(Zp{Sp;w_Rd~CP7>d$w`8Klyb^yc(es)t+#xtFEpa_oX_ zwOVD*(T^i2Lv7QSrc(#e*QZ?yK+%g3Lc6+1iUmgge1WC4mh~~(^Ik;TlR@Q|OjDpp z_anmX+9^s{#3+aPZmA8I-IF+yLqup#t=x1y%DP|L3$J8w_K|wy=ojepdj4f{@9On{ zZ%&K+(1Zqtg9XaP&HeBDp#{i09C9G_d}!Wb6|tuo78Z*#2nz<}hqpQ@ZjFc}ScaC< zcO)NB{PdD0sK{*?ABi{K1+vM{c5`N93!tiz1o`W4S7YY+C!myYrN|XNLnw`DU@tto0~=s0-@*Ywr$nd%G$i{Wh8Cw+4`0w9-YmJOUoQ?{9+!lZpNw z1Ttr<1U~#Zc*G2XdAr9LP_zpAzHGaR&qo^z8uE?0!2m?Rl$t08OOCNo4X7xgB|%A& z37Iql`>MawG$egoANIbwa6qser=$8y z!(F8$uxe6}YRQIxe*%)Mor_*;YH8z6kpCQB7D4MQpvg*D++DZ#j=BWdtzer|nTA?y|b)dww-{=&tFJk z(LQH#ZTnuXLtn~&Bie%W5K)%}c5kkz!atJ7Y8TxD1}N3AxD{m?UF{L!KcYv%5T#16 z%WJl_0@=W}m;mB034R-Tc6M^gpm)*?HMHdNMc{rF!*owg3r>b8$pju+PIo40xkIF) ziI_9Ua#^efyP^&*%)F5;bZ!|OyM>=J^;^aGXNpnB0~yD!mW3~%uWf5M_G0{bEnNxF zWKVOCtZh+-JLVSb>O@K(<<7|{>W@|`d+DPSO2U{54IbO7BpPGCYqnIuw&d(CU%xm6&~$2Yz>xCOinN2BQwV59 zI0cOP@k8(Jh=#|o`#a7lJa+CJ@QtD=bw}-h>_cRX3(7e3-0^s;WZf1t{5@e2H0c=y zu~4siOai!jR6LyN4*0^2+3ocZ>1nkebyk)k{(3CZFIdx7Zs)!{5lqI)_2RF-_rHEp z_rf)GFcTQ}|1fFIxwSwG5+)uJVo7?-hRQZ){q!l9|47wydJ#~;g-c|lg7>{TMw!02 zP#FN?fvuV`4_k~BMu}1BoOBm`GB@Iv?rFjY>V!D>@q3-mNrEGy0%hXlNVKE|XvoHG zaUk`6Xr57lvQLFeNrw~Y6w@PtfoF1ckO>RRHfx5Gb&|Y3-}0hS>n>XN>~PS#5{-;a zzs^Jj7m0=x8EOjJY3}}zNF^8PDIwdvA`f6Sh>2D7N3i19O`W;Wts1sD?B`?Qo05p8 zmlBlo9R+{xP}IQ=CxG;`ZT~wxkxM)l1XyK#GC$pqDzur^~4xa&J})qo1&73 z}leGy4W5NJ%CCie)k<&@Je)2Z1d=N|Hr8#>OtUUCSQ(lX<75 z;oeh#=L1a_qfT#EqPrxkOGfzsR-)A%mg50p;1~o%NtW%Y<&aF>d2#jxP#lkzlE=#s z$4xAXJ_V(PsHTV>wKa$0xZw!}r2f4mjNvO6?@@G;<7e3$pPhXR;_gk#dAIM`#~NE5 zVDbMBCiMS!yU|})|LCE)O45%Pt_A)k^+Xv9-;!uq)g!wDafI|F5H|@f0||;i5H9ur zyp-qEJPT|GBbq6DM1-9KJ}7mxBf(}MsKx1H-a*+6RkDjClsWV^JUTlJ@X7fH9AAUM zDVM}axrbPiPrxXsNwlFR_d2U)oRif(SY!klWYUec!H`i$h-qpywkIXz(W|~PA3Wn( z5J{yD4Nkd|cc^HqHoIv)8j1a~U-+Wb%^cKDPnlej*dTrq>4wVzkgUu?xunDjDKED5 z8441@go5BkZA8#Ch@mvGrySYRAdA}Wkuc+;5Y)@omkT$|J;ImrgQ)rHE#ct&)Z!IJ zbz*w7+L=b!*OKS#1ZEbrnw8nD{zV?Yn25*01xvUKy+AZ5sUR~98d@9OcTNtyPchvm zcWBZ5uVQZSk4j$wwFYN*pyO*~t+gu32+JkPCA$yKDE;ORpOLm^UzRigjp6u?wirmb ztQvpdaFalbyYM%z@Jj3@#cDIPgZ{G z_s#!i#$rtb_-84=RZ<)r7`p%3VE^Q6`>CLuvM_Lt(T+2=oT+S8LoAlU<`f$Xk-IqYni^VS@Uo%Le9 z7jqGYu=B@b;*N>!sdvYBARS``HQ3;)qTBX%ys7|=yED|fNmmtb7h!IrgQB(b&v;{q z=2l>*i;!saZU?yLz6knr9FQDBrFtvPXB=gxse;r?AF!KSCP7jo85Jle6KfI?G&Mj= z`#({{f7g)sO&Y>B3{4zgNR`Nr130o!xGU?47#c?VU=z-CVxb}HXw4%moNhzrzz!B| zU=)R0bmn~quA=R^{Pp3YY^#crv2a{g8V5YV!{?M)1{#r+F_sFM@6I0fwx&-XMo@L- z;p>`FXjbBgI7SK&+%cUZI^igC6R;GR8f&>0>xn6(Sd%^R}|Dsy$3L-YOG|!&8ki*sifObaBr;8 zUBRI_+T@~uq1=2Fe9>O1@JNPOBv=|}xkU4S@t+8Bh`Y4@Xem@W?Lws%2KdZ-&R|K0 zGHqIAq(uYwr3jO0IC@%k1<|c#!}i6W97(--lP_m;&atmePOK4_3Zxc%@-adEXxkvv z+)m%}lTJ6dZ;rQfC9W;vMnx<;JXR-bu+*@B^Nfo6oh|UWne8Zdi;}Ln$m~{iI>gxC z`OShDy7i}&eh(!E!tM=a15|!Qcim4We8FGJ6vk@0cL*G$oSl@|%p^Q}&iBqPx zFWleMbVx}fRD$^c~qtpCaU69WXlN@XfZXeGS0r&}1?#Zk@xzDYa zT$PwfT(+rXqBb_x*qsarCVIkCGsG6(3gfF?~nbPt>fVRL(%>McT z)|>qN6KxfoEM|42)w z+9lwxx_*rXispKA7-0APPffLHne$^~T1?~Vj`6(JOHK)&x1AJ%w`<(@VA1vvs6-GK zSLPhi%keky>N@KDD3X~~wJ>kt?*jiOIZ%AHxF9?G^r-5oYWYlC2ldeW^$f}lck4Vi z2$OPAzKdPG%Q*dA1HbNHn{c+}8YIX>c#jWXD_gJ2v_`1;p#Vq6r(lXdKEzyn2+Ve} zwx-2VM{-uhk&8P|h%gezrw?c`-#j7rA^5+BVz;TxZ2II}gZ`=M&SJ zdD;@96;XB7!vVdz@-Q*r+F2PO@!u8i@7K?DeYGdnD#W@gMD+$Da4D@eYCHNebyW&J z*fqez4oN!Co1NWhGa>FyUyoH;BUb|N+OQ3FPqU>C{6@#Zyl)7qJXk)GROLR9kb z>+QNx*>urPp-gGt2!xs-BoldcgnbsTbjuBN*YE~3t|#l*QK zihoL$OsUz=TVwM8tk%}w0+XI~Tx^wLg`l#>i~!D!=O*eb>L)hhr0_9&hlBTw3_!K; zFLQ7x%&yj=kL*}auJKY($LCOkbBIAQSP%R8$P_tvBl}d|(Q(y-%9Kx>Tu}k%e>E_; zz8s@tXRGZl>?P~jTwD?YpOBDz*jR7k!o_>6C~()g`|SF!AM=ORRUrm!(c~Bvth=V2 zr2yjoa@+)^JTq@x6i&)nNbC`zKlQ35&9*vZ>yUvw^jxD;Vvc-H784d=gheci2!7;k z>An<`0TWF%n_vY_v1?Fv;Xf8;Ru1sCYi6YYkF@#*O=wAQn3?pybr51pYY_7nWUVV1 zu)R=S&&G|JL0?{&$f72<<{R$LjK)NSm;k(^Lw=~iGHf&K@QJJ85-V5+Z_r7=(Uf(`=0fV7tJ{N8mq_r0FQy6i%qg~BlXf@pd5tetVX-t%1|lK4dfRp2 zM8O(O`n{&E>mn3)^WCZYg>zwxf6=X-V^EvhqmZ7IKnIO->U_u~BaKBq%Tqvx0D9W! zD9~B)`u^ZSBxfFCoSYC}N@Py7vEv=x+5J)IQF8`xucV8&5T~Cfh$IZN-FpdUC&0l@ zu}9A~2{>}YiBs#SbgP=21v*2nK+2vmhhSCkV58Ly5BlYy56g((#hc+c4#gxLL~^t4to z-4O@4;PKClr;5SLV@@1tSU^Q51{&ZKG;MH`o)%h|+1u+~J`pOUcG~Oh;t>>NmTcqA zxSg)3&Qyt<%{FWkL2KN}tiolC&^LG~B#l7ssN)Axsp+u$tC8#x7KaJW1u@|eFLw2q zmmb{P#ADCjULObS&*@G+T8k+qqR(Tr$V~xoW9MpMlfCO)^qAT6SXa=ZGfzHk-IeP)fC@0tq+cPL0&;kK}AfB`qkz;a579nuZwIfllcmIFQ_b=n9R zechSr=jP&~hLW-@y78HVRdMp8vnWa@n@DyWUvTCZK)Y*Q>@;QVxt!%ONG`GEA&)Y_ zdB^DIGq5zuRKZuLF(Ol-EP;y!=O6CuZ|a8So2ip^LSrIM4ka}jXo*RbpQ45JrG2)a z16D%8RgcTO@8rUg?Fb?Z4;V{KmMsVP{mB*evRV>r(;K9}c)xId_8So*zNoM89#a-u zRjPfoe!gdQgJ}=xCNfUtX#=%77K_ihgrcX;iCQs4Jcyaq9WMN2k&s4VrVN9@72DVP ziCg)$9UzFOGb@4xsH!i;VgEfLWGuLuf^DKTO{1;$<}B>N<=^z)-^?k%Ys~|Ugh(s^ z8g`pSyz|i|sL8F&ulZ2*cP7otAr+;T{A{%E`nMZ4XXSc1Uopv#Di61Yjt>iRr#SJS zO~Z3hbqO@g%N>|Um~7p*=?v2m$=O=d4kQg;LW$mM)mMVX-D&hz$uw)Uwi^r1zcz={ z=p{lj;5&;-zo*4M&@AU0~pla0T8HmaW=m+duDs*!7ag(o2iG zLl7R1fgj#RQn8Wee<>$LcM;Z5%hZJcXkl^9bGBEV%c;9(ud&fODr(Yv8ZrNVXuK7c zFzmi9WWpa+s;Y*sOumE{Kq0UY)FBBiQ!Qm^`1o`2OkSS_YXd8nYTCYfiH)2SNysX*mrH5 zFyCtvLtB%Sg{&gIzWV@iqW%@F(S-&v*rkN}(udPhZ$x|KcMcH&#a@e-v-~K5+;rpH zs7(z9%vpjvnT*!FJmcBjbu(>26^0xocP(el_(y3@fohBXuP{O(sx_13ePvg!hHv2n zDw~4f^I4^Rf_v?n@Z*)mAZ44SXb^LUrq^Sa^NW}=0^uD%-sTi{wzN3XATU4o2%Yde zFMS9E#quy__NFeb&Su7T|3MDFtYKJ^knGW#R+ve^fLWT}*u6kOIR3Nurvl7q%R1n4 zBK6+Z?KSi`PIx7NMFV;D+XZ4kUEwX_aC61tfU%Pe#gTlyewiPyN7A=7lY?`QPvv}Z zeR5@EBb@6|701IuF2u`=JhtID4^LsEOzp(QWI(?YVl zE*Vt_(+wow5Bmg5xY7r`ssKbP1Vj*8-aVq@Dg2O5NG$U?PgzWfo zgUEYmRXibTCV5fh)u0?g$W@T*4lo9LT;O(?@_0tIAx|(z7be>Q_{(7&GHyGh1T{;d zATW*a;d!RpkmH=vbXuH9r2?ZI@hF9{Gc$_F5vqcQcpk}!l@Dsx6snx z&HJyt=gXDgz=pFcgJIw+=Hpa-8|J7VXpq&09U<606Nug6 z<>Sc2{uK*e=>2om7gKMBUoRXdXwUnpZ5ZMhp~cV}O7QGS8fHYTRlG4@x;|a)tfD!l zYKUUtt~sFqKA=y_7eP`Wmm?-5Lt~<$)dTDL4UeXuLh&x8Bh8!lb+xo2ZBRL$=4QhR z{bNKSl&Z7Ome5v<0r6F!ekEf3a>zWGP>E10nG|K}r2%{003|x*{%H#`6$MxDs-tSQ z^I`#|Z&aK#1p&*OG#?m@@lzI{;kX>rcgAYvKY+25JCsWy-Ra*tg=u3E)9c;kb?|(j z%Vbguar?lh&nD6iFYhp2I$O!?#(muDFyhjy&|&7fsv*_jv}~GctNKP2Zkxa#{m?D$ zRenY%U}FXN0sKmA9+|OF-?BNSqUY$iewykS%AILz&VF(R!QeMX?1Es7kJyZP9k#4> z4S+MXRW0OYR4)CZ3u~;w%E9BuVB*j3rwd!3-(py~F`#_0iONn_XtJI`XoH4V!YG7s zf-`N|Gkf^aGn$5Qz+rDpYzXz>lx4*$i?Wlx3Za8~z-JH#BjVCq(4M6LWZ zJ&_A82c4oW4ux@sA7EyXTdMR)j_`29T-9goC5IJJ0g57jHYrT)qqbOr9pRXLtTT-9 z%|`o>*RMtXU-*QRGl`f8wP{LF8xolPzmQ8394d*32(3w3kq+{|AVH%S0+^#oNk0Py znCrhbUlwQ(z^qLu_6?xGoK3?{L!iK{X$mXAxJfZca7lG;aG+dF-2a=GChYMzk$UDd z@6z?Sszb$5kl;K6lYpR&)&u ziIF8Q0+SJR(qCFY&c+>2r30&mvx$=aJf zuriJMY31{CRwKwF+7r?F4Z!CtZQ@=XL}3D$MWvEo6k| znro+_9mXo%RIt)3_uq(2GC`1TjQQbv2-BeipVVSTc(H+rBA$PY5Q*h9H^Kv~M#5#3 zoF~-oDQj4a)P+uPKQR05&YM~j()&AHU_YcL*&jcBrg^!L7LU?-#7 zNp6Lw`^kZ;CIMG4FeJeZ)`2dp7?F@h;ysdgdo-j#O~3HmCv>_&9V|*iX%m~g6rneK zV^r?r^`YU=wRDgiIw;2Wo=+-mU{d)Cjc82Pw9pY9dx@UlSM~dkDeO7LhhaV$>Vqo0 z6-Pv$b1~?mUWw{iG>X$lsyioCso+|`X=#JVP_*(bX#r8e!6OgST&AQ)n~(LM&uzQN z2h@MXLa7Injs2W_KdtZg27Ha0@DlV)v<@ui{h+c_3N4B`0xfiDT6g9B8yxVl>2u&7O{(eqlj~g1& z<<-q{V3s*(SfEZR==W9c^a^@6zat^|*;5hZOwo25ZqZJ=yTTyb2gconUp_;iAX1ik zq}0yv#W1EaG8p73_5u9+zFiu8o4NhEcyi38>Hw*N3;-Pgf#UP<%O{Aa0Rl7|H`IQ{ z+4`WnU}Ub5Y^?h5q<4dmW5leqETBuk>t~#suZyp%J4T)%-k3+#^NA1Ve5y~Reh1`t zZeo?xawlm=YlCJ+IkWM<01Gc%SQ%^-l85?%q!lGm#0N3zj?#2j)6j{Zw-xA<*P!zl z1JJ5zw5Zm|W!JEQ*fK_Bjh&fOQbAH(59Vu?y2Ywi9}9NF^^snmDI~`RiyAJzCg|rP z93JF@!l?hqpIn4mtu-;q=KR$lOb_LstnzAlHP!z^MtYM8K4Whxyi_*n2U ze#0wRn*S>z5;Pk_3-xw^-HSy!95R!gXRLI3Km5hT(=FmWkZtnE~2P0gV-3c#5p zD$<3%aCLifKV|4)G?XD@Y(#A&NGK?62wYMw7?b$)>MeeMBBDVK49-Eqi|9d|aYZdK za_CZzS{uWpOj1@*(4>5%9)iL#m&>zdK`4rUDuNK55PR;^kgQKzaV?2OQSkmWbE{Qw z2jJ%zY;Attn5QpTuGxoj&1R4F12&WdGok}K9IH9XAfZa4g+ zQrBOK$#0&9a&onNvMfqS)tYDG&dF2HA<)NQvL)zOEHioqnqN#w8Js%LWEUK`)OE=` zHTO+d=a%}jlOsSXrIgUx{j% zYMxu4T3=c%s@=T4WW%-6EcA2TyjG93{|`uEsB`ry0Zh^A8~jufl}uf>t*Uf(FCPEX zyL+`vU2^uQbanghQnpN0=^I=>Sf^=rb*~<`Fx359??z51ZLQ9-*v+eE>T-aDKcIOG zHnMle*`cJ36(BBbo@zuzJ%Pv;CeQGZHXot-pzzHtELRAW-OxRp5rQm_ouU`K0woFk zqS%kW$0YykPVOmz55}X^wFvAbYS~vlRD%Bvtm;)@o~Aquj0Vh|Ht-7=wF%k}K@6Dr zKNzivZlD48AEW*kKWx)O1GE1Jqc^disIdWa{uj69;Q$j!fv|A5{hP6*{vS__Hlb)> zGA@bXSX7(1YmPD*2tu)-dTB zr&c-caMoBX$Cqgu9>|mIVQv<}g*j-9lo4$f>hdO4xk)6SO)DuC7!RUFtqYrs=0F}^ z07XynZ0>{+(JDAaQJke%X)zcQW4cNhP9pU+F=Tk-;8XKjf9ZWEr<|9_c_0UpTpVmc z;G9n>Igw>t7tJrFM0+%QNEOsKw|F*!Z|F#bKTE)jV&%2L#u8xNFzpi37gy;E7s(gz5!0Ta6Ut1emUvJtB1)pn0NR)fD@QYQqYKOQP)nR}Z5tv=9&>!P1!O$P zebjvp>kHHj)cj$nf$Zqz_-QLD%aC8QTz^fkESNy2tEyP0Pth=%iP%6JVPXofk(GW(LB6Z%Ssk$+qm!!#;zA| zX(l17U~fT%vkUJ#!I1pl54V?z_u)kgo~wstoP_4*e=V8*E+58?COtYoC=lt4_#nw*^f6S)~RjhSxO04?14l>Y)tIB00#Mg?2leZmW zoHcq^I!?^ix(-Hci!1^J1}GlJ5N!r8Igi0Q;;=(b!2HaHA{R$az~BtN+X5|!27!XJ zX2T0ibfo*`tfIfp&ONyu%*=KWXdVp}jwFRCfuXRGW2A(^W@O%sQV|(RQ(PUPy-}&f z)klgGU!)YQwGS(v(mTm;j%M6;7u%Pfx~nzn{>|1h>P3Q!YRHU2%>2SsrMA*$bl6;g!xcrjumfVglq zihFTUG7?{z76WzlQk8}g=n>ja<3bu{x0*nJrPSJnh2yXuVqArewk`OpWIjMCLaK8+ z!Z=FIC_`g@a1|>#7637Y88crTn*W?G?2d*J!_3B_QtXU@P+-|1t~wf`i1BvL$jiVzx<`ur)c)c}Yep15hI8@6>7p^qdKV zknE_5m6oq&*dsvWv(j1>lIE8dYS>yP3W z>C6^T4WATM8Obduj~V#54c_3CMkC7TKL5wm5vU=alO74`qm46jF5kY&Ix^V&rXiU* zB62!c0jxN4=zid9;qZ$ZQ1xp{4b06_`oW7q=}o!nE#_;=MeyQ#%%uAbE9g*c!2Mgm zecJIFHb54%rwPP$zOT7JwvsicH+zWR(judE(ruG#tPjB@LTb)ZHvLA4u;HNZNpEp% zwNk#TYK9AVO&}?6w3y#0v*hC2mEwJdH!D^GVZsujlMwxye(EB?-=N|pbJLZ~cO3jF z`-#-oYxg93bnx+n=eg*keeEvMv*?!lY-!J$36S|z(79V_a%Wx){z4m~M z9g;BXTFjh8mo?k`D$pht%r`HfJ&N^d;g3vA!D(ETD(|;P1r!G88{>jhz`J=$beo`)&NStFs5FS-*8*{L>`=8#)!;O%80Y?v$XbO=}T8EA?mBj z>srEiu~1qy@+r;@gIVlhKg7YlnUWji5?hYj2 z_gqTCXI`#M+~#7n`iR5yVUw&?p?Trr z-{)6mAKv+g&tbPhHu^;Sg!Z;+F|M+B?D)x+r9C%>QhmjByKgXQ%>gm5kN^GL(UpTh zer)Kw`e&~HDkm9}b4b7q?}f2+TFrsao9%m2iTXx*vjG#iTm_$xpH?}98`s;n4*>YG z;cjmqF2SHWp1i=`35Xuj25@otaCKs(!XYYJ9Ip@xZ6DTMdPF}XIc?-AVZVR^#tHBJ zzSDmO5#Lb#hvVK|{09}^WxsiEAHSmkdlsLzzHTe0r`<=#r&-_efFnXkLmB@8MhY{4 zhQb?p-{;5QM!b9ApCrJ6r~q9VK-GUl(sz2SDADlo^1B!hek!i1lAW0DvojlR^z@VU z$CouX<&-eYcd*QN)FLgZ@u;b9^fPQXT8T$-ZZs(oun)Me`*iqLwIU+Ix(QK?}Rni3PKv=?ijGpIXm6CoghM`dq1>_2TRQ7Rhr863FA)FGA8 zuw~OEIbC2rCZxdj9^CBFQh+Z<`ijXO@4hW%ZlB6N^`1DOj1W=)zM2$|=@!)6Xa{j_ znm$u*(uy9xpBfd+kJp-ve<~YgWNZk>IFopfTNi*Nx@3d82zKNi>C!gX*vmiU}G=pf4jo~svyw@Gw&*<$#%0X=}vI)qlK0(CxP1LfRE=8oy4C|JG)Gmg5m zW3Q1`7oYJ6mPXpt+~T*4u~|N5x3G~Ksd1>@l={x@Z&l3n@Mdwk7@V5d~+U2 z*MR7Ou{iQDF50rI91ww0GXj39g-j*;gXr`0a2-sIXdQ8Y#fW%~#Wk1uHu7vS{XiX> zbT?U;iYJ9QEsUy+uFbE>|=|X_+wa>6WKlP?-Gk2RgbzMgK7K|bL@6DAsodLHaM-QBjNMu zwPx_+#xdJhNHHsZ4X%$Fd^_AAT`&}ls4ymdbftj}e4fAd!%Y9`iqq8kHcKmBrW(r- zMr)v*K7VsGef5Y0s(RFn1L#(%wmEd3hR0Z0Yv3)wo_&Z#rXhN99j`PRTdoN{Ju7u@ z%Dfj`DJw(AFKTtVdsi)v@8nX8-ob4#7O&Dom8F26*zqm6Q&nP2UN@?SS2o;6*i%hC7( zWEnvKpKO06R{L*yit3&al+Sy zkZU*9q(9a);nj-Imz8SK6l;<6Xv^!#K|g1XKSy!8;d*DQll#n*>o(hTzNg*Bcj?A= znQgH&&}8qmdg-abJ=Ym*vGP#8dSC6F>jpVjd8}P?sC~>=m3sOE#0r^c@IObeNvFQ> zAe_u??2Et#AW$r9%xr8-%-m^N%fPt*SL$5`j{BGR{X160!)AhEl(2Gkas5vQ_y1&b z%p5FS|IOxN(&ZhtIgr9`ztQ#)a=IFTq!JJTzy=`fj`R zZZb*8H-tcJ!iV-*8L1JdS0W^&XtPlge~_UWli7ko9>}D03o7>e6tL_muR~qla6@dn z(S~DGWIKrxMQ~~CS2tfR0o8m z$}d6HXcE%ytuk2fB2}Txi5c|JCgj;MObAA+FsRVsNij`jBM*+b)Onz_&w_R}&+vgSgk*Ca%St!+90F zlrbxcgENj8uZ1;&6c$@aE-z*BTyUZ_FfLN*`^E{_4*+Jj5DXwJG;3YCF(jPlHbZ5!)D`xPPr|`Cqt&@Iv zso9dB13k`W5^uBRSlit?j2Qar<=wt*r9G06MPL^5d|7`oeARE;e6qQbF{f!jHrudb zsXKLdEuf1EZ|!pkFSw3_{nd47Cx#8^I6@LNWgV!pZ0QgctfHmOesuidZ~3a<__bP{ zCu1MYbHkdkRly#cJTua3qgjMMIoCWAOYGpg?sxg?alZkdj-mpDc#2em7x`o~i)N zIchYTll;JIyh&5s&6dbD%}Xd4-Z;c~zsNZ%So{;AQM7@ZT#FsCt{sanEL;O}ZB zUf3HMVB3>PW)`1d68fratoyJU62YOY-+RRA3X*(a=i!neqgoa)1hQ;6xrYaoJk`0coj|H7(exEze$?~tPeM0v zP%E1rq8f%Z#GR&1C@Z&)awr;Q*12njv7N_lFsrJW5F+}Nf1pAnq*?8sziO>dkc_#1 zNx=X3%X1fjSpnob{OjM4Y^scjYihzmn_IIX zK&|Xj5Q^PO_8@4uxnU796J9~Mifd$`-R~mjMRYV7&~R^2$lF4jdyiU`CKw&#Wa3s} zw=otq#79MR4F5g;U*A;vPq&Z9pb33P|IZ410Itb=dT9LO8k)o?LID~3vBithBX^wD z#4Y4%`&6IiqV^J%TSxP%%A%%vjSz&>uF16_op%Vl<0)H5ntEdbAGM3L)HVmC-sF+P z-0~`6;_6yo;%Wv%Bt+w?$wT5rk>mR!mD@C#p72KX<9nsnN#Ik#M9jN(<&e(+bbH`l?##<$VnF_I#uxANzm&6d-zIlI?*j9if;TfnKVOVGXQ z=I>(xFSkqESnUrIv=|ql>4z}Gag!-&-xn2xlO>IM8<+&Zst?1cYUOQ4#Ky@A!>CE5 z%S^;d#QaaG;^5%=Psv2YM5F@4C~a@<@K0j?-;x-SE}tkjGaEM-hd3)2E4PR^>p${i zVSQRr!#rce+akNHk#PgWTS2el zGw6!BE`Kf1GA6erENPb+vk9}LP#?ZYVI7x#1O8ENW{sHP9{9D~zK59Z9t^ttx8je5 zHxgbUhuk&662ZLeHA_pE+y_uy9%eZ`E8rSZvYfxzJ^WnRG5VaUea00C4>EX<4XZU- z=XpdiAjv2KU9(@^5tPo+7%V!7b&n1U^t9j05m*oG&@KUI*&V$#iqE~UQT5?I<55A^nLTUF2kRDBv`xd-g@=nMJ)%Ffq7a&Y| z21?nF*FKi>4pAT$U(oMQVS`ScBrBcj_X~OQtFfe1g1?ERq?)=~hWQ?$76O+fxqw*K zM15u4Wm)^evxZuis#$-XtU$g{f4~6h8hnvJh8G(V1k*0{#qA8jl#d69I+i~md(>YD z6+cYm4R8WOYUhH8&_-$^QV}lx1r?u71gFfp{!w-#5*7pp zjtz5wI!qa%D0eqv7Ip+yYO=W4?@~i9DjXK{26hikg7lxs_+X;00n9ku)UX|!YkDD$ zAX~5w7!FmWrhc=0=6S5e@;4$zORB(b&qyl24|Q2mLCYBJ@&FzsL>wT1Jd*>2yb-Kr z>(hf_%{s&xpkril7;>dG7~hC;9eg=_nQD`MCg$AZLiA~38+_TmGca>nIG@8A?@3EX zYni15`8lR2z&-);?V$PD7a}@9Sa#~iazxq@CL6+H0E-Lqtknx8T1RJNJ`#411Pc;s zMx6=Xyq_+g_H@f=al2-4ua|bkKO5EBduZgQAT5E7v@gJ&srT=M9Z6L7V z|84yLGyKn8{D<%VSp0vDe`f@0DVi9Y&N=tmQU9dT?*kLUaBy(JkduomO2GUtft#Ng delta 16240 zcmajGb8u$C*De|-6Wg})#MzbY-w=v2gGonHrO=m1|p;3I2gnG_7wRk7LP>|3j!yJ+=#%sqAB zK0g}&4YlTTGHOfm46?dV0%s#7oo#_9`@ zY10l{g(Oy4ao7n_i5n^gB*=g=dWS4aa&u02`@F-KjYjI8rV-1PV4qD74!4la^~n$V z@*!od-;v`E-YQ zL+1t8SfWcxmHtAn9yyiF3JWMsnlBcjUVGOBi9zmHNSp)jeRh$CnHFf0OBUV;hrT|N zLVJ@j`W%HEz4_BWB0K_aDUQnsIKok?V>SA}`&&4&oz|1w%S{X=3Zk9mM_Wua|BczB2_}_ z5}(85i>->%sS z6c7(WMM|jjT%s57eqB=S(&~O2c&&!)OGhV9yj+s;vyfs3N-PbVImqr)L7torz}4s* zI9S!7(_rXPCvwr0Ys`EBw3LDP8{`jk$MqeSYsC1ByZm z3Oy84W1OgbVFN_KX%1w_NlaIObHI57!in;C3*&YNr!1jCh2Mchdk0Nyw(s`N6`o)* zbAnTVICc&icsdJ95(^D?8lP(buwE?{x3}I+`_YPUfnR|rc0Ot}W8M+UaO#OU58RgS zXGW+Na(-~OowjCi1YyDx>jr6b4rngI(&{0~9v94J%~2Kg1(8p@x@A?x6EFl|9@cr& zspmyUN6Q=?+T6OE#`Dg6lnu=VtGiab^=|On|3s`yHGFuln}&#fR2jV$sDLO-$%tB< zO!|aB6^kdYVsVu9NsIS;p8Ta`*f zr!PN>h`4y&C!!8I{JWwG*yMIW+u*9i5Uj>2^_E&5T`y57hxOOhHn*H`IJQp{XNTfA zn?;q{JV10|yDM2s2V54Gg!M$9%b_j(dWr_$eH0pN#=<8x4RfwAR29!pFk*2SG6U+& z_j&L*`lme|cD?vFQ2l3rb=5|$2ww|FR3YEVE+Zb8s@K>QVHkb>oATIxR!=a6W$nno)CzV(fr-G?^N@`(U6))ffgA7AiF03uR;^6T}=%F#W)L z&+V~H11N6Jer5gzOucHwF_ulUy`r0~4Q6KI6O0L1Q4eOU{ zmYL7W3rsZiJay|H&2C-tII>QNP5n!pedRIneCcb(Gf{u?jAwAQ`$S^IyX!65s@@-t!L?30cYH58B*YZgXPs{EqK55>Q`KgE%SYjCMZ!wO>#6ZX z7Em4S@`OFERl|{uBfWK2aCD17Abs#=mK>KlArT-7Qsq@+J{AI5U7dnHirJ zC>y`Uf!z15dHN5U3a=O|49Ap)gAfj++|f%aD4^J^7DCle+!_0Iqf}xhec{$(jdv-` zM#Q>1&H2F<9V{9UysP=kUQ_TBGnrUaTkg=P8wqAKLfk1~7~4rA*bc=P_V6kFN_QPQ zo}-CNJcH7xY)hr0%T_E|w%aYV_mdO2wX#bGCW@cXY(0srAf2#Qt!U_|)PJxAlGeo5 zanrUcECwS(3Ux^P+VJMm!1?X|ibHfrg*}4u19e(iq!T>7x-EAvNULdrB&sRpvC(!d zC;TewRw>1zw-BEfj{b)^gwstqRevT9R!Xw23H|8xQ=0OIjeK~#2n#u{#&a35% zY9J;k3y=wgmqeoANrVwxL5UJuoJq32e#j=S{U?6c$|(WsQ(%PI+$?N}=PQQ053Sa} zb;mJA)YRCT^J~NJtHbMMwtn){;P`lgcAZ!=P;z#csh=DXsa%Ic!VqDPdr#mS^R_2y z^oODeHR&c+;-=q~(w1IzuOKuKNR6;5I*+s~csGfPF8bG#&Wobu3y#mw1)m~9d6j6R zy~^M^4mIU|U>h2od?{MGDbk#J0tVhB#oPRLnqhWQIpr!L`Pye;5aM7o8~goZ=v z%!d1k?&as?(9UE}MALXDKGxz-YccJ!n@%6jP4BryrL<#>HhXw7UMUfvQF=+A2tkJQ zhj!S3C-SCCpM?MpB;{Wv|0_t@0$y}6{5cn{W-->3GEw<{1~bP%Tyy-*T$XmL-bhaf z`atkIgJ!xgDe4L}cK`*|j)fA12ns)m_vW0TwOE5+0Fe^;q)gO}`b6v_3 zYH4-?8oZp1$;SpV-$TEFRdK|;J?%DyA4o!dl!c8v;@lyOE%d*Wn$bWxTrMoc$f_ju zyLz2qUY^x>`y19P9QGQNjgWRPW|W9dR8|gell>*nuWPf5?aMSmYk0Cwo7qs_yAmAM zM8}qHc`!=_G!B&Kh~>zd-B`@lS+v#tf2Di#c52!=h+&VB1+FCng&P$;`Nrb)7`iUZ zMuzSLqv3SYGmtcitSZFcJ!u%{XGq@k2HMJQ-BaIfxSf{~Y7$N}gPwj$mb!IT|IssI z>v~lP%rRzJ%80*+5ouDm7pFahE5j^fAj)%)*Ahb+xcs*3!H(@0ajvE7+41X5mt1)Z zfAHO*(CVnGNt18~*1x3+#007%a(yKWOi#-G_7M4QLh2^WkG>jnHv@;=N^`@08?B*6 zdKis)D{Z3r13rmSG3pY#CETByMSnoF)6(gf#tExl$$Fj8MiB=k2!nL9419a{w4|!R zon!S%h-EXs_U{m4x`P{hL?Xi-woEo8rv_(ZX8CWHarpg@Wp3(`os&SAL{)A?k?gOq zVHV3;bN{hS0>w2Wcx?*`jo{k_ukXZSok@Lc5nzuXX=tDC7ml&RmGO-C!q4pC*xvIo0R=#N(& zTLWL3_=j^CA5>w*RG)i8S-2a`*N=&s9w*v)eC<-!H>r>#+2YbRh- zg%$&xsD0+|mNZ2i{|VVCF=ho}H<`v?;G))ff5{7FlgKAoy{Y`z^tT;J-*xUw+dl2K z(<=bo|a}f;)7g z2>8-IsEXsn=$15}(x}cvja5xan7;#7@$FSZN8Ig9op;F~rLIy~#-a-fdOxqT%@uOMP7LjKNgqT%`|NKexrF868^@!kAiEL!zkw%Pb zKAV=^wR;6V3&9jn(sv~^YTy?pb6d{!w^uhLW=<9~^`l z>pJGiS~L+{QFAbSu`v8mUU;RZJ(gjgaea_zfZ(m>@wWA9EtarKSt1xPce7jP!)*`1 zmnT>@E&-4sJrL&>C*is;*G^YtR8FxRF5MpKfDP+zZWF_i2>HX*I5I8gX<-pnguO@O z?h;GQ+Xr6ZV35qJig)=@idCl=WaK3F=8rn!-B*Ulo+V+TNOmNZ+KlUmzD2nAV6irp z+-YDHNIkm#J@WB!cQynJ*IteKu~XUhP2QY~UGP2`mwvDPRY*~3eavML{4LoS)z$|Y zg8k6P{e5+~lC5nYk}zEqLzAvR@=xv*=PnCoCbBs;r7if+J z5?fI0N3-6m_$h(tm_>0Ktt$wlv2gDL*z!W>fBLA+?dXEzS~yE!tC~fDor5=Z_E$93 zKUb#1CdB-dMZis%5XS*pM3zPy9DDYx{d~tkaHAiWwt9E33x#B7zAUVyzu2TOZw=KL zReJ}jh@I~i3mgA}Cp{UwaFJ;QAcMSwzRj~8naZ=sk3PZ0fM14tofyjKc$g(?(outR zvIG9{GBvOjd&_={VM{ZCaP?2v31NQ-4~UZgLHho$@@YU303KdP$PI;V%~{ znM8ahc{q4sVuRUFg{2XDVJs$vJiQz~&%D>@oIlbD^o(tA`6{Rq6QuIgkwZpK%8WZ+ z8}~yaHy{~uAVbttBr^dqpd0#3e&@_j$~U}q&rLvDl#>u1i2ymO@e5ENf)YkUJAikG#18bis#Gzz7NSWXapL$*LaeR)8H2>qPw=B*#tq)<`P%pk=@Diww zvUy2SDnzNQ_#YH0rbb89nHJARTaBAPJP?c>dU4-K?jFGFF`LdWMk_GA%O&Pie?tHX zwPtc@8-FS+WFq^mC5Zk8!{#IjVB6%uMNWOjZM;Om11J#X)+qo(zYU?B2^@qC8=|e; zF*AT?3-%qudm=4BPby8;7w&z$jT*eTLgPUrTv{Q)0I3F<;pKofl}H-+fr|2-Eq;eM za~J$Xu0BbQa&oycLSc8taB!@L_DY~6O`)=8w~x77B5VNtSu=@+4CWA(c${A+!pX zV_pjyd%2B~%O5H?{IO@;9H@|pDn1gZfrwQ7 zpq0|ERow0&_dl(;;A;a2HPgUZl$VKdoZ?CK1x=L)lGRK?4N)!o_cF)Xu@E)bpjmzH zgD6d3Fw(5WKaZN-?E7H5RF>($s+EM9IQj1O^MXpdfLRF`e9)kf*UD_dF9Q=}a66!h zKNs30^N3VJL56xI}V-WB^!ec=bI>W5esA?*cpdkc`D>PqYo&+fZ9go{a! zwvE?7zAn^T1g`3F@j*+7P* zMv%48*GKWU$_BJkL^D*JiK1Eakvt&p4;#f(wDbmu5sX1Gc~v3dy|(W1eYfUeJd9bnP3gtPKhT z*4*CA#nsu|*zP};gNY3+3mZEzfcQTaA73IaZ1Wa?1Ok+``IF5H9E^?Qzf^HXTUSw! z6WL$z*I@irbo|LdBSWE7XGEI{op?fc8#PjhZE=Ma7K2jhci(PKX+vNsb3-=e;_|-XZ!kqmV`I`Y-tz7Y#6Op8CkqZ+O;8FOz`meZ#6KxQpdsACu;*<3 ze0av(`>lT&fkBv5A5r)H1rv~1CgJMF@D1IVh65#y^#O9tsE+2+eKNo(VRxy~qI}t~ zKpt67NJ{7i^pNa&fGY`cBDoaDE~E-%w-jjIYy1WVq%})i^0Cplh|V1bLY8mZT?-Gu zRA?3z4{8Gt;9V}OIwI5vrl8KEL9DL69`#Tr@dX=ZdRzpyf1u(zDZT}=U>)eLCDhSH z3>otC2?N5@R14 z&CwnL+zousNFmw<0U#CuW&trY2K&e4mGkln;jr3K%rOofaS-@=;qYtdk19Hr{n7fK z_`}rFQu|P|<@XQD;?o{qs2#&;rokNp>=zO+R0F#$PS}o@6)wn(GaLBNj#&L(N5y7q zFq^}Rqk}(=?kI7CuWu{<*t^pL`XLv=J0E@=L(zvw&BnfwqLvTQE}m^ebJQ<;&02fA#qepZTBMCvJbAkHJCMkUxWWlgJInjNDH_!7b(C&`sHWD&^Ym# z)@@&Nq2Zg~TDSFR%FupCFlBefoW?7VKZjM|$wt9yIpG4_IPdIQVNjs_01z7Vb!Qw4 z0&~8(W0)|fUYd_wO}u%(xW3r>_=>C#`11CgZ7DR68gxvB>x7dc!W@fnsxr8I1STE; zF@p3v^kB<(+XF-op3dIB{)9-3iid_{Q#blXTgK%r>jj(`aR34{>gW~L7?wXJhS>yW zUlkSTDKbDg)!Xz47BJ~)6Z@zA%MNDFb0qPL_3YCk=FBS)+Z7WL{|oSRcpOKn3|+ZW zGnP0hOYx=fKX=bv2?KNe}E;u8k|JKP$USMGF{*z9X#l$hy%IB8^@Q zqsa=&ucJZtjAZ+2NA*{q6{U%z^_?h*8Ag1D2gkr$Fz-Ev(P2GO`Chds)e)_I*-M5W za+d`Nj@lQDSvlV>-$J;asEP-QnClJ(=lIKci4OE6xlf;zb&+Kpg1< zf`Crt76oe549H{ zCUwA_ds2raJU*G6tQ)sokMeX&U3k(;@%XPYhC!bm@c58f^j!C8XC^eFIKkO={0Sq1 zn6myTU`>To79(nG=byM&|F)dm?KrFDH?mEhXC4_Eg_GaV9LjvHNIBKSE-eaxQYjN& z5W3;n+a9z1~W>jz-VYrh>dd))0LG&YCIboGn@(JLUE$^Jdy z@dZ8K_}3|Tc0Z23^3|}-;VF@z2Z#iNVm8Z@grkVS_U()bnN}0lm!H;6ex)j=+*gVK z00%a4>q=ene!JAgy6;erR?(I5NS_UPo>=Db=gDN4ROz$SiOQ8no!g~6j%^|*&N480 zpklVjRb<3yxf3a;x9E~Q=P56x@kh<+sR|P&3XA$Ss8;{TXqO)uw%zWiiZ5Em&(M^g zL>W5bKFiXX^KPVn{8Wjs<*tXpi>XF~drT(2|1kgCsFMny&l$O+5!nnppzQBQ=3}k= za!*T8ja&J+iP}TlfVdK_AvdgM*P5V#*bITcZpzF^)WCm&?1)0RPWe$sHvVA9mFtlPePJ$A0+&ILy;2D8|d_ z0y!)r$p>u-#-M3>lRk&Fu!PpFNEfn=G)+b&mAMqne0>uX4kjBLSqw&M27+rDSSsux z!-xo6i8^>V0A~9D1yDOB0@%>OPOgBT3c80xRb$Kp0IMfS3i{bZkVpeyFb&CL%EXbT zl|{`!weFy%iw>^DVN|LihgDnWgW2>#xJ@#v(viX9^%jcverlBy*Yk^M%19)nRgu$6 zJ%(7LpMoWG3>*xl$%#13*V9Nrw%KeJl~VGUV+h3ozdbMbeFsD;m`yke9EV%aY z)3`+CUG}O{5NNZ&e*WYZ8OHBup?^j-s?n`;QR4Cr%Tx&JG|#xj$F7UV@FE_au<*+@ z=#`m5gHRvNxYDIUi852l>DL#!`$Yxg>X`1N3Pe3n(+jryyL5E?nZ4r1BR`ez&nG?A zUTs*nv@PWD?b0|H+#Fv$at{sPrDJGV+F;$&r`2bmZ)eZATVK}}=(+W@wi2#rx$EFB zUCzt_F?UyC{Gbu!3ABV?R zl)~>o&(AqGf1tm+TdLu1&Pjar({Z-a{Oo&^0KQH9 zP0rNLZAZ(a(F3$ww~aK}C{eHF7YKH?o%^wL196vP6n$>N6$i8gYT7nsDRE7)CojNc z7w4;<@1WCQ{xcEl@r~U8WX=B+Qa}C-ve7-Ute|N;3SW_6GliPYs7F5UF^F>wm_aMp z@@4<9`YXCOzK=wKP-l{iEVk$vN5Wh1 zrJ~XTVPA8-+ejTDGIu6Wx#3S6RT)0$H2(|Bk1(48%x`{3;=-|?vd~Zk6%L{$)D8~1 zgtalx6`-)8esN8^dlrO&Rq4MAfjAe`!7+L8O6k?jNIo8CqKBnsU9~$df#wy~aW1eF zmp_N|2zXA{>UT*YUjL9U+-F`Bek`0?PSJTwK%i~5op%A0QYPv^5u^E zFLScxql7WOw)?!lHxOib@6J1oBj1Y~hW(AS4p96(L#XPCTZ_ZU)LQr#UFYfu*}dG+ zIW%#`rp2W?xS$)4+LFmuO*+)0Pe1I6%m~uw?t&= z(L!YpChz>xId$!S`ncVH_o~=c$t?*VUZ=+WMvm9fbi~F1%iuYvd;|eJUyfJX+rN7c z11G_A(D_KfzB_ec-`iN#pLhR|dW1fX*UddVs@_GYU9Owt|4e@V2NigI+q(#^N9uEO zr}VpxCBr=ncEFB=A6hV`qJEZ`Do$NEw~GJ3fXq9aP*PDLXA6_BqcokAvUnNIhg&Cy z?I%l1AD70~P*OijYT_?X1wLm3?CiIoup z>a1@jBNKBUL}MN>EL+x69%eEL%Q|lWo(%$GB4= z%H)7;KV!#ybt85}5ezC4MRM&Ah^V|U%0LO^(X~h!jM&RqNaOLbki$jEa6%;$F&o7r zl>wlg@*~?o&XCdqElLx?5$bBd36nym`=Gu-la;NyZTVA_?FbH;N(&?I1v zQxdGCj$6vYOlGfa#5RF6SL@z zL^K=iKtc#kb)ZUReWDqvl~)V=M5ueExI<(Cvl3^P2C*mxo3t=O1DH$Dvynstf5EC4 zJ+98HaFY#?g;p{bfqqEEvRZ0knK)I+@?_j-{&v2ia5j)Wm)lGYc$Waaya2*$@Y!WcPA{6H=3UDW1J^LY1gHq`gYic;25q|W)VYj=8h$dfacqPe zXHNoe%f?CAx=w?U+h9mQz?y`kNj6DI(O*mwwwA@k^loKIQ_#i5#IB(mN2DTYp^0x* zij84F0zFZZ6VC_jBb1L3QwsNyQff%y!Ekm-ow&?Bf`ql_10HQDL@?CdVYHR|iPh+H zc2N)_lq)D}ZbOS?3LB-*;Fvr_MbV0s(9nTdY6~IM#&9?8xP>|}@*(Ge4WA$@q*h>R zNI|@}GH9-a%IO>i){ZVK9Q(u@u~c(epaB|GFkBDI0LS^l~{??5j{Luv{4O2pt+<(Wy5rHp2j{5 zbPIZ#l&*Cm9K92##6}1-Np!_D3!hYI3x!8;;W22nQ^^?xlGgT8lLzP2w)e!gNk#f^ z-Vr6HPF?s`MsH?Yw^m(q`zs}ujUyn7;yXxg$TmI0IjYPMR&bo%kJ5o~jSG_Dt7ee7 zm;naufK=>8GMI>{#|b64P&G(QdaZ)F&Re-=8WFoEa5=jq+J1;yTX~fM@~g5c|8C+c z>{8o)ZhGMkr4D)XWdh#-YvvNvGOU6BSSZUJ8_~j(SDGyh5#FIINQ3r3;7WZAT?*1w zmEi#6?@>|-J3(V=8H0)f@_wn!v3|c_*0E$y1fCM$MR6&g#Fx+M*$fd~ZP(CF(0_GF zcorvrDeTHrl2Ui#LED*_BNWE1^&hu%(Yg_6dacjVIMy`)Wgauc+=capla*e+yi zG%Wm*;Hne$Ip}W9VW$bZP;0DbjXIPsTQ7F(X6-4p1k}}C>kKC(By~2vELj}vHiBbN zl>Ho~7Q9iQmeF++r5b(eApLY5Cn4o8_OtF1^&8Fto^_VQ*Q#$-oD;Z0=)KGNdsEqo zsisxEFJ=DWG<}PcK%B5b>>|1|pfYze@;9h_&f0V(_l0eSo8g!?`=HY zkGdfaw$0YGAeZn`Puy*YHC+lGlY-Tr*0#`$z0henoz0%tdRBr7E2K$Rot7&56&VXr z+WM8Ei&sTEp!aGPh(711P4mQ2N#q|97u`nTx-uSkg*KJ^hxeF)?J{#%CsY(s0L!5O;0+QJkU`S&q+=PZR8O0K2xgp1 z3x2FP`8~PWY2xnlR=IQ%0*@Q~@1 z88`tv5f#w8MfwQImt$fK`$$Re`S4}Oe~AA$Je8f8d)=tPXPm?lJUIax=B1+v902bX zeW%BON4~R6uOW(mT~+R7{kd~elb!xLDYsgLgc_zbFFd9IkMrxxqVFGvGUi@$m|hQj zHPU?V)=~`bs(;9w{GD*sxM%uc`bzzt`KGa4BBJXej0GrGohTXAX&hXv+`ThFJM8wn zS%071;nDD%?W!hTEig?<(^Lk`vI+fS%lJtQJkOuVGdcE3PP?Ri=M$&vZ-2Vf})Ovw9^FVr5lTzocdOFukKF5(sodf8% zv8t)mH#QojKQzmnJoq&7&^$0Cd*6PX&LnF-l^Dn|#g*MT-dKC_=&&+jt>@pz0z(>YwC)ax_u_#?SO?uk>L- zQ-BC<^Q&PU1$(G{ePm%^|LbcDv3s4D5Q!=pSgh7*3bzs64oD0K1WqskiVC2`upDkmfkcnuJebQ*Mq@?OrUX@&sgB`ndActwq4c>&rnw$p)979sJW^$N zZw#NCggd)3WSjKXoF)isWziNjjXYcy8dPDKXt^IaRq~qDW;=JH9pDSV^u+c_*yyjV zGfyBZxmmvF8RhTm0Y`%}q*BiFsZ^MxudI;i96{WfV>o9Gs1pHv#;tmZ^5_Szw};eO z$Ck70pFjIAh%eoch0Uo?i$axGe%2nGwHx#ite5dz^KbwObP?rgGmWFS`r9fAyW$w^ zhb1DbMljm~Yni6>hRUCSS9fQ2d)@jojs&GWaf|EtWidE%-3S?evbPsDeV|I&tsJKg zAeR-R`ABn<+({3QQAwX5YpgO?kL5bnkI4%(MOJb*p&vHyVm!LJeywMu7bh|$?>E0r zFsklx%3K-f{H(_|D$g&g?MQU+{r~-@(Rx@c?w9t|StTiK(ZsI-bP%@l4u_9bddqHy zs6ioh9ow_ijAnMe?O7$Lc1k=yb=X@Jr_io+{{!G8#- z)495^ji8DHq4r$XUo6{z#C1F|7N3Y?{EjhwDRAsvBW{CS^2k{%9DjQf1@BKs!o=Nm z8UC8e1&8v59HXQgaNOF^6<0R9yV zG62}A001f)SSEP~GcjXVb7E>S9u@!#8-VkFyQE(%fyO5bbK&E%F~c%RSUbD8{tNB6 z|0m32=3ry}uP`rGOWxs-1G)D^^R~mFAW{tRSHK8Z^zo{Dw^Vz%3DTyP8w8iL4c+_Q z4R3szT~ta+j@%p{;x3;Geo~Scg&9Q+dMgThVnu1fZ&wi>Iw@*92papw+!VA;lyC2S zY#Fu4uGRCX3d~mL81mstY#Etlg(6}6-TuhYL3>0c97Bg(?72#?XN_KeyAF|9(L|$#Mg#YGDx)k6k!zbrq7$`?|aS_M_tK=4bhBAx*P*3tVkEWJKJan7*4DkqqX)+$?cD z(I{MD9^8IVI90?9o#!9StKesjR#F~Rxa8Q7#K{hd*%t%8&`0PY3fUs)gL5eA%7a7V z$m9h%;;)219**5TzNZXDn~5|Ap`i|Dwts!HXLA7UdQlguD{7^UsHU?+;7digXF>5{cPjL4{}+0 z=5@_+GuW-o9vWI07t-v7ir7ErQDm1tI08w=jq&*%ggk@NHB(uwNdoS{56K_NcQFn1 zc+L*&2jt~1R!Kb{Nijt{*KU}$D#(oX<^1x3M0_;QI<=oO9iNPo2TfaF?H;RhczaKq z8M`tQc<<-AV9z|=AP`9YbFsyyKL;+2n+-;vIWzNe%yn?v;DChGeTCb}044$h6n6bw z^V_QH_8-!pkfY!*$@IvYzBrJ?l3Dh*eG3O#U}t@G5V!H}NR*_>1|XAW{xA@~xF-CB z$IuA)m!kmA$qTybPG@E6HJer)3M6;t_86Clo%?OoNb;s+--X5rxCvbMwGQ}s+Jx0@ zAnKt2l4pA(V0nFuZubA)HEw+KF$$MfRi-S#M=2z8Qw20IQ&|hXY&&75($<3Jl)M2y z3+=)#+*sLSs92AL(w9a(zUTQ6uD%#-cN9vKP-K13BV9kcKArTPUOO`KI<7$enrO5pwkse+=KJ=oV%$kqW^%u>Q=r zF5-zDFWK|jL@mlgOu?dyGEvu}s1azVRhy)3N3H2xRvDtT^a`GV=@t4-tX2msY-u>g zm1&1fAD=!~_Hv2RBBf2&RS#gi<)~)%&T9{9VKK~vU7XDMk9AT@TuVDidW9 z`CkvgU#tp(FJX8JYN0h0e4UjcD^dSjL|P%X!spP9TYTb$R#NJl^jQ9MP>YyPBO%>l0tirx zT{-Q6GWLNl6_41h{aZphnzvysb;O9)`obuSZ$YFbPk+;arUgm)yKn;*A$e4BELw8S zp)CUGt6|~sIql}w>ZuQHu&X)k6j~=<@$Tr|YUd2Iw&HFrP^&?OBY*3RPK_8EYHE_W zG0TJQWmrj*gHEiwVfS6M@vkFU4EYcTtH9ul)dNEt#E+bitqbovlq)@COLv4fitf3n zb=3%=&7(dU%U?>|H`z1vB9yCi|$U4!<%#k3tH40l-wv}k8-n>dP7id3o9*2QgiV}1v(h8b zw}ScMc5Bd6XXd=VFlKX#x7p9}+g-hX-{(<{RFBL+?%B_56>@*}?tC z{7;z)8z*bJMh}~-lPYnJh!P-AMHW|H;1owC%}>P5(4Oj#prQN<_@sN@%h0Jg z!qS?GU>59WIqtolXr0}-@yUDf%DWJF@1ZGutfwulGOIijH*qayVcf#TsQ|-Uw!u~h zmuSq*Dto(kx9P617H3!&wi9I6uR^>maZkwdw*z8Tl&|Vm{>5C{E(qL|w#k;m6Gm(O zVWXIN4Vt-WZi6TH66CjRuk16}f^%8!p1)~Y1}!9Oo>g5@%W4N{x9qOuGyF)*CB~rY zGu}ers`Rt`Il!9gv}jfFN?>WIAbS(g>LV8)Kxm^`v2E^+e*5Ybx{13a&|i8lYU7bX z2o407lN_)s8yooHh?)Z`nx;~`b4|q2i5^Y zIe?x4#cG5(9q8rwLk|kh2t^N6a{#pi>aoH*&RCYcx9h~*-YzrIUS;ahpyhzF15yJt zPQ)qYfw12Nf)ibgJNfE2M2W`DUvw>ymkuCHJdq59SEP=BQb&p84!K%cq!9I?2mQ** zc^(GQk<5UaNu=|wHS(8I7_G(Vs5;r4vR&m-ZJBN1yt{F72uGwlqz_38R!5=I=7cA5 zhgFkng*q1T;R2YHQwyjwElu$@4Ana?B#OWW5GVw!htfe#%Nk=v;2p$p#_5f{| zYFrpO|DOkT5J-56emmQvHWo3lFh`IcBpy|S_8yBPjupJMnr{NuINTq|dlNWa?-#r@ zi7#4Fs48KEbttPKZLdmjUQl_!8c^5!uf!DY0O-7t^?9+4&w#)@yj5Z+20kc58fMoM z|C^bH^RhipVAiy8!KyR*Je5m)Y3dJ#3-B?JZBim=D$grUoLq1+qxr>6y4nK%KdgK7 zTu7t+Tm?tRH~<@KKW#e3fM}3HFIc%?ZU7P^sE$h?f&~Nrdu?lUfJKV61n@7e(UzsV z?`qvj?&9DBtroZksHAe;WzR`$7dipIAV_>ak0*RyfP8@L$yA)M<}bQz#k_|yC9e$K{kPRfN-5I}uKM%M91-q}*4s+mLXGwoZ|w0jhoCi^!*ib^ zS<^mGk>*SvlQP|#WMF8}f3>1t7+Tc-S2O$vyq)|X52F8T{eS#FEq23_H|BO8i8Msf jcm5B@{cnnMaW!^!^>jA3fMsT3=Hh09rJxX3lz{zTjY8vt diff --git a/man/class_var-class.Rd b/man/class_var-class.Rd index 53eedb5..c2c09a1 100644 --- a/man/class_var-class.Rd +++ b/man/class_var-class.Rd @@ -19,10 +19,6 @@ S4 class for variances \item{\code{resid}}{residual} -\item{\code{diagG}}{average diagonal element of the G matrix} - -\item{\code{diagD}}{average diagonal element of the D matrix} - \item{\code{vars}}{variances for reporting} \item{\code{B}}{var-cov matrix of fixed effects for gain} diff --git a/man/dominance.Rd b/man/dominance.Rd index 5fbfbcd..53aeddd 100644 --- a/man/dominance.Rd +++ b/man/dominance.Rd @@ -4,13 +4,11 @@ \alias{dominance} \title{Report dominance parameters} \usage{ -dominance(params, digits = 2, index.coeff = NULL, gamma = 0) +dominance(params, index.coeff = NULL, gamma = 0) } \arguments{ \item{params}{list returned by \code{\link{Stage2}}} -\item{digits}{number of digits for rounding} - \item{index.coeff}{merit index coefficients} \item{gamma}{contribution of non-additive values for genetic merit} @@ -22,8 +20,5 @@ data frame with estimates and SE Report dominance parameters } \details{ -The dominance variance (Vd) and baseline heterosis (b) are quantified relative to additive variance (Va) and std. dev. (SDa), respectively. For single traits, the estimate and SE of the ratios are calculated based on the delta method (Rice 2007, p. 162-166). For a multi-trait/loc model, the result is just the ratio of the estimates, with \code{index.coeff} specifying the coefficients of the standardized true values (see also \code{\link{blup}}) and \code{gamma} specifying the relative weight for non-additive to additive genetic merit (see also \code{\link{gain}}). -} -\references{ -Rice JA (2007) Mathematical Statistics and Data Analysis, 3rd ed. Duxbury, Pacific Grove. +The dominance variance (Vd) and baseline heterosis (b) are quantified relative to additive variance (Va) and std. dev. (SDa), respectively. As of v1.11, the variances are scaled to the population (previously, it was just the variance components). For a multi-trait/loc model, \code{index.coeff} specifies the coefficients of the standardized true values (see also \code{\link{blup}}), with \code{gamma} indicating the relative weight of non-additive to additive genetic merit for the standardization (see also \code{\link{gain}}). } diff --git a/vignettes/Vignette3.Rmd b/vignettes/Vignette3.Rmd index 8189cdb..3f75517 100644 --- a/vignettes/Vignette3.Rmd +++ b/vignettes/Vignette3.Rmd @@ -68,7 +68,9 @@ prep1 <- blup_prep(ans1$blues, vcov=ans1$vcov, geno=geno, vars=ans2$vars) ``` -For multi-trait analyses, the `blup` command was designed to predict total genetic merit, not individual traits, so the user must supply index coefficients. The `gain` command allows breeders to compare the expected response for each trait under different indices. The argument "merit" specifies the contribution of each standardized trait for total genetic merit. For example, if genetic gains for fry color and yield are equally valuable, and maturity is ignored, the code is +For multi-trait analyses, the `blup` command was designed to predict total genetic merit, not individual traits, so the user must supply index coefficients. + +The `gain` command allows breeders to compare the expected response for each trait under different indices. The argument "merit" specifies the contribution of each standardized trait for total genetic merit. For example, if genetic gains for fry color and yield are equally valuable, and maturity is ignored, the code is ```{r} merit.coeff=c("total.yield"=1, "vine.maturity"=0, "fry.color"=1) @@ -123,6 +125,12 @@ GEBV <- blup(prep1, geno, what="BV", index.coeff=index.coeff) head(GEBV) ``` +As explained in Vignette 1, the `dominance` function expresses the inbreeding depression and dominance variance estimates relative to the additive SD and variance estimates, respectively. This also works for multi-trait models once the index coefficients are specified. To be consistent with the above `gain` calculation, gamma=1/3 is needed for tetraploid breeding values: + +```{r} +dominance(ans2$params, index.coeff=index.coeff, gamma=1/3) +``` + One last comment: groups of unrelated traits can be analyzed independently through `blup_prep`, and the outputs from this command can be combined as a list in `blup` to make the predictions. This assumes zero correlation between traits in the different groups and requires the same genetic model to have been used. ### Genomic prediction with secondary traits