From 3e71be9d2b9438e78705ebfb8c1983f923399559 Mon Sep 17 00:00:00 2001 From: fouslim Date: Mon, 3 Sep 2018 23:15:00 -0400 Subject: [PATCH] adding case/ctrl preprocessing code --- .gitattributes | 1 + README.md | 20 +- code/20151007_RV144.preprocessing.code.md | 683 ++++++++++++++++++ .../20160509_RV144pilot.preprocessing.code.md | 2 +- figure/boxplot-chry-2-1.png | Bin 0 -> 19945 bytes figure/density-plot-1.png | Bin 0 -> 22923 bytes figure/exploratory-heatmap-2-1.png | Bin 0 -> 126072 bytes figure/heatmap-top50-infection-1.png | Bin 0 -> 74584 bytes figure/heatmap-top50-infection-2.png | Bin 0 -> 145426 bytes figure/heatmap-top50-stim-1.png | Bin 0 -> 41988 bytes figure/heatmap-top50-stim-2.png | Bin 0 -> 61872 bytes figure/mds-plot-3-1.png | Bin 0 -> 26602 bytes figure/mds-plot-4-1.png | Bin 0 -> 48414 bytes ...umina_expression.rv144.matrix_non_norm.csv | 3 + .../GA_illumina_expression.rv144.metadata.csv | 441 +++++++++++ output/rv144.eset.RData | 3 + output/rv144.esetBaselined.RData | 3 + output/rv144.esetRaw.RData | 3 + output/rv144.fits.RData | 3 + 19 files changed, 1160 insertions(+), 2 deletions(-) create mode 100644 code/20151007_RV144.preprocessing.code.md create mode 100644 figure/boxplot-chry-2-1.png create mode 100644 figure/density-plot-1.png create mode 100644 figure/exploratory-heatmap-2-1.png create mode 100644 figure/heatmap-top50-infection-1.png create mode 100644 figure/heatmap-top50-infection-2.png create mode 100644 figure/heatmap-top50-stim-1.png create mode 100644 figure/heatmap-top50-stim-2.png create mode 100644 figure/mds-plot-3-1.png create mode 100644 figure/mds-plot-4-1.png create mode 100644 input/GA_illumina_expression.rv144.matrix_non_norm.csv create mode 100644 input/GA_illumina_expression.rv144.metadata.csv create mode 100644 output/rv144.eset.RData create mode 100644 output/rv144.esetBaselined.RData create mode 100644 output/rv144.esetRaw.RData create mode 100644 output/rv144.fits.RData diff --git a/.gitattributes b/.gitattributes index 11c5925..c287fb8 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1,2 +1,3 @@ input/GA_illumina_expression.rv144pilot.matrix_non_norm.csv filter=lfs diff=lfs merge=lfs -text *.RData filter=lfs diff=lfs merge=lfs -text +GA_illumina_expression.rv144.matrix_non_norm.csv filter=lfs diff=lfs merge=lfs -text diff --git a/README.md b/README.md index 0a1696f..1b093e1 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ ### Fig. 2 ![Fig. 2](figure/20150201_RV144pilot.Fig2.png) -Fig. 1: [R code [MD]](code/20160510_RV144pilot.Fig2.code.md), [Input file [RData]](output/rv144pilot.gsSetVehSubstracted.RData) +Fig. 2: [R code [MD]](code/20160510_RV144pilot.Fig2.code.md), [Input file [RData]](output/rv144pilot.gsSetVehSubstracted.RData) ### Fig. 3 ![Fig. 3](figure/20150201_RV144pilot.Fig3.png) @@ -40,3 +40,21 @@ output: - MArrayLM list: [[RDA]](output/rv144pilot.fits.RData) - gsea result table: [[RDA]](output/rv144pilot.gseaOutput.RData) - slea ExpressionSet: [[RDA]](output/rv144pilot.gsSetVehSubstracted.RData) + +### b. transcriptomic case/control study: +code: +- preprocessing: [[MD]](code/20151007_RV144.preprocessing.code.md) +- geneset-analysis: + +input: +- non-normalized matrix: [[CSV]](input/GA_illumina_expression.rv144.matrix_non_norm.csv) +- arrays/samples annotation: [[CSV]](input/GA_illumina_expression.rv144.metadata.csv) +- features annotation: [[TSV]](input/Illumina_HumanHT12_V4.hg19.chip) + +output: +- non-normalized ExpressionSet: [[RDA]](output/rv144.esetRaw.RData) +- quantile normalized ExpressionSet: [[RDA]](output/rv144.eset.RData) +- DMSO-substracted ExpressionSet: [[RDA]](output/rv144.esetBaselined.RData) +- MArrayLM list: [[RDA]](output/rv144.fits.RData) +- gsea result table: +- slea ExpressionSet: diff --git a/code/20151007_RV144.preprocessing.code.md b/code/20151007_RV144.preprocessing.code.md new file mode 100644 index 0000000..da5249d --- /dev/null +++ b/code/20151007_RV144.preprocessing.code.md @@ -0,0 +1,683 @@ +--- +title: RV144 case-control gene-expression analysis +author: Slim Fourati +date: October 7, 2015 +output: github_documents +--- + +loading require packages + +```r +suppressPackageStartupMessages(library(package = "readr")) +suppressPackageStartupMessages(library(package = "Biobase")) +suppressPackageStartupMessages(library(package = "impute")) +suppressPackageStartupMessages(library(package = "org.Hs.eg.db")) +suppressPackageStartupMessages(library(package = "ggplot2")) +suppressPackageStartupMessages(library(package = "limma")) +suppressPackageStartupMessages(library(package = "knitr")) +suppressPackageStartupMessages(library(package = "pheatmap")) +suppressPackageStartupMessages(library(package = "grid")) +suppressPackageStartupMessages(library(package = "gtable")) +suppressPackageStartupMessages(library(package = "RCurl")) +suppressPackageStartupMessages(library(package = "dplyr")) +suppressPackageStartupMessages(library(package = "tidyr")) +suppressPackageStartupMessages(library(package = "tibble")) +``` + +set default options/variables + +```r +workDir <- dirname(getwd()) +opts_chunk$set(tidy = FALSE, fig.path = "../figure/") +options(stringsAsFactors = FALSE, + width = 80, + mc.cores = detectCores() - 1, + readr.num_columns = 0) +``` + +read non-normalized matrix + +```r +rawFile <- file.path(workDir, + "input/GA_illumina_expression.rv144.matrix_non_norm.csv") +rawMat <- read_csv(file = rawFile, progress = FALSE) +``` + +read arrays annotation + +```r +arraysAnnotFile <- file.path(workDir, + "input", + "GA_illumina_expression.rv144.metadata.csv") +arraysAnnotation <- read_csv(file = arraysAnnotFile, progress = FALSE) +# remove unused phenotypic information +arraysAnnotation <- select(arraysAnnotation, + -title, + -`source name`, + -organism, + -molecule, + -label, + -description, + -platform) +# remove prefix 'characteristics` of column names +names(arraysAnnotation) <- gsub(pattern = "^[^:]+: (.+)$", + replacement = "\\1", + names(arraysAnnotation)) +``` + +read features annotation + +```r +featuresAnnotFile <- file.path(workDir, + "input/Illumina_HumanHT12_V4.hg19.chip") +featuresAnnotation <- read_tsv(file = featuresAnnotFile, progress = FALSE) %>% + as.data.frame() +rownames(featuresAnnotation) <- featuresAnnotation$IlmnID +``` + +read clinical annotation + +```r +clinicalAnnotFile <- file.path(workDir, + "input/rv144.master_wk26.csv") +clinicalAnnotation <- read_csv(file = clinicalAnnotFile, + col_types = paste(rep("c", times = 21), + collapse = "")) +clinicalAnnotation <- lapply(clinicalAnnotation, + FUN = type.convert, + as.is = TRUE) %>% + data.frame(check.names = FALSE) +# remove unused columns +clinicalAnnotation <- clinicalAnnotation %>% + select(-v7, -v9) +rownames(clinicalAnnotation) <- clinicalAnnotation$pin +``` + +create non-normalized ExpressionSet + +```r +# format raw matrix +rNames <- rawMat$"ID_REF" +rawMat <- rawMat[, -grep(pattern = "ID_REF|Detection Pval", + colnames(rawMat))] +rawMat <- as.matrix(rawMat) +rownames(rawMat) <- rNames +# format phenodata +arraysAnnotation <- as.data.frame(arraysAnnotation) +rownames(arraysAnnotation) <- arraysAnnotation$"Sample name" +arraysAnnotation <- arraysAnnotation[colnames(rawMat), ] +# format feature annotation +featuresAnnotation <- as.data.frame(featuresAnnotation) +featuresAnnotation <- featuresAnnotation[rownames(rawMat), ] +# create ExpressionSet +esetRaw <- ExpressionSet(assayData = rawMat, + phenoData = AnnotatedDataFrame(arraysAnnotation), + featureData = AnnotatedDataFrame(featuresAnnotation)) +# save raw ExpressionSet +save(esetRaw, file = file.path(workDir, "output/rv144.esetRaw.RData")) +``` + +quality control: kernel densities plot + +```r +intensities <- as.vector(exprs(esetRaw)) +# log2 transform (replace intensities < 1 by 1 to prevent -Inf) +intensities <- log2(pmax(intensities, 1)) +arrays <- rep(sampleNames(esetRaw), each = nrow(esetRaw)) +scanNote <- rep(esetRaw$"scanning note", + each = nrow(esetRaw)) +plotDF <- data.frame(intensities, arrays, scanNote) +plotDens <- ggplot(data = plotDF, + mapping = aes(x = intensities, + group = factor(arrays), + color = scanNote)) + + scale_color_discrete(name = "scanning note") + + stat_density(geom = "path", position = "identity") + + labs(x = "log2 intensity", + title = "density plot [log2 raw intensities]") + + theme_bw() + + theme(legend.key = element_blank()) +print(plotDens) +``` + +![plot of chunk density-plot](../figure/density-plot-1.png) + +```r +# print problematic array name +outlierArrays <- filter(pData(esetRaw), + `scanning note` != "ok") %>% + .$"Sample name" +print(outlierArrays) +``` + +``` +## [1] "RV144_d305968_DMSO" +``` + +```r +# exclude arrays from raw ExpressionSet +esetRaw <- esetRaw[, setdiff(sampleNames(esetRaw), outlierArrays)] +``` + +quality control: multidimentional scaling plot + +```r +# create tempory ExpressionSet and normalize the data +esetTemp <- esetRaw +# exclude problematic arrays +esetTemp <- esetTemp[, setdiff(sampleNames(esetTemp), outlierArrays)] +rawMat <- exprs(esetTemp) +normMat <- normalizeBetweenArrays(rawMat, method = "quantile") +normMat <- log2(normMat) +# generate a multidimensional scaling plot +distMat <- dist(t(normMat)) +mds <- cmdscale(d = distMat, k = ncol(normMat) - 1, eig = TRUE) +mdsStDev <- apply(mds$points, MARGIN = 2, FUN = sd) +mdsPercVar <- round((mdsStDev^2)/sum(mdsStDev^2) * 100) +plotDF <- data.frame(V1 = mds$points[, 1], + V2 = mds$points[, 2]) +# color dots by whether a sample is a hybridization control +plotDF$controlSamples <- grepl(pattern = "control", + rownames(plotDF)) +qcMdS <- ggplot(data = plotDF, + mapping = aes(x = V1, y = V2, color = controlSamples)) + + geom_point(size = 3) + + labs(x = paste0("1st dimension (", mdsPercVar[1], "%)"), + y = paste0("2nd dimension (", mdsPercVar[2], "%)"), + title = "Multidimentional scaling plot") + + theme_bw() + + theme(legend.key = element_blank()) +print(qcMdS) +``` + +![plot of chunk mds-plot-3](../figure/mds-plot-3-1.png) + +```r +# print control array names +controlArrays <- filter(pData(esetRaw), + grepl(pattern = "control", `Sample name`)) %>% + .$"Sample name" +print(controlArrays) +``` + +``` +## [1] "RV144_control" "RV144_control_rep12" "RV144_control_rep8" +## [4] "RV144_control_rep11" "RV144_control_rep1" "RV144_control_rep6" +## [7] "RV144_control_rep7" "RV144_control_rep2" "RV144_control_rep10" +## [10] "RV144_control_rep14" "RV144_control_rep13" "RV144_control_rep9" +## [13] "RV144_control_rep5" "RV144_control_rep4" "RV144_control_rep3" +``` + +quality control: multidimentional scaling plot excluding hybridization controls + +```r +# create tempory ExpressionSet and normalize the data +esetTemp <- esetRaw +# exclude problematic arrays +esetTemp <- esetTemp[, setdiff(sampleNames(esetTemp), + c(outlierArrays, controlArrays))] +rawMat <- exprs(esetTemp) +normMat <- normalizeBetweenArrays(rawMat, method = "quantile") +normMat <- log2(normMat) +# generate a multidimensional scaling plot +distMat <- dist(t(normMat)) +mds <- cmdscale(d = distMat, k = ncol(normMat) - 1, eig = TRUE) +mdsStDev <- apply(mds$points, MARGIN = 2, FUN = sd) +mdsPercVar <- round((mdsStDev^2)/sum(mdsStDev^2) * 100) +plotDF <- data.frame(V1 = mds$points[, 1], + V2 = mds$points[, 2]) +# color dots by protocol used for RNA extraction +plotDF$RNAextractionProtocol <- + pData(esetTemp)[rownames(plotDF), "RNA extraction protocol"] +qcMdS <- ggplot(data = plotDF, + mapping = aes(x = V1, y = V2, color = RNAextractionProtocol)) + + geom_point(size = 3) + + scale_color_manual(name = + "RNA extraction protocol", + values = c("orange", "red")) + + labs(x = paste0("1st dimension (", mdsPercVar[1], "%)"), + y = paste0("2nd dimension (", mdsPercVar[2], "%)"), + title = "Multidimentional scaling plot") + + theme_bw() + + theme(legend.key = element_blank()) +print(qcMdS) +``` + +![plot of chunk mds-plot-4](../figure/mds-plot-4-1.png) + +```r +# print control array generated with different RNA extraction protocol +rnaExtractionArrays <- filter(pData(esetRaw), + grepl( + pattern = "Rneasy", + `RNA extraction protocol`)) %>% + .$"Sample name" +print(rnaExtractionArrays) +``` + +``` +## [1] "RV144_d825482_92TH023ENV" "RV144_d824374_92TH023ENV" +## [3] "RV144_d849655_92TH023ENV" "RV144_d723863_92TH023ENV" +## [5] "RV144_d203143_92TH023ENV" "RV144_d827158_92TH023ENV" +## [7] "RV144_d736104_92TH023ENV" "RV144_d603208_92TH023ENV" +## [9] "RV144_d627156_92TH023ENV" "RV144_d431442_92TH023ENV" +## [11] "RV144_d127852_92TH023ENV" "RV144_d513999_92TH023ENV" +``` + +normalizing raw expression + +```r +eset <- esetRaw +# order esetRaw by idat file name and features by ProbeID +eset <- eset[order(as.numeric(fData(eset)$ProbeID)), + order(eset$"idat file")] +# impute missing intensities (intensities = 0) +rawMat <- exprs(eset) +rawMat[rawMat == 0] <- NA +suppressWarnings(capture.output(rawMat <- impute.knn(data = rawMat)$data, + file = "/dev/null")) +exprs(eset) <- rawMat +# remove outlier +outlierArrays <- sampleNames(eset)[eset$outlierFlag != ""] +eset <- eset[, setdiff(sampleNames(eset), + c(outlierArrays, controlArrays, rnaExtractionArrays))] +# quantile normalized and log2 transform expression +normMat <- normalizeBetweenArrays(exprs(eset), method = "quantile") +# surrogate remplacement +normMat[normMat < 2^0.1] <- 2^0.1 +normMat <- log2(normMat) +exprs(eset) <- normMat +# save normalized ExpressionSet +save(eset, file = file.path(workDir, "output/rv144.eset.RData")) +``` + +quality control: gender check + +```r +symbol2chr <- merge(as.data.frame(org.Hs.egSYMBOL), + as.data.frame(org.Hs.egCHR), + by = "gene_id") +pb2symbol <- strsplit(fData(eset)$SYMBOL, split = " /// ") %>% + setNames(fData(eset)$IlmnID) %>% + stack() %>% + mutate(ind = as.vector(ind)) +pb2chr <- merge(symbol2chr, pb2symbol, by.x = "symbol", by.y = "values") + +pbY <- filter(pb2chr, chromosome %in% "Y") %>% + .$ind + +plotDF <- data.frame(mu = colMeans(exprs(eset)[pbY, ])) %>% + rownames_to_column() %>% + merge(pData(eset), by.x = "rowname", by.y = "Sample name") %>% + merge(clinicalAnnotation, by.x = "donor", by.y = "pin") +ggplot(data = plotDF, mapping = aes(x = dem_sex, y = mu)) + + geom_jitter() + + labs(y = "Average expression of chr.Y probes") + + theme_bw() +``` + +![plot of chunk boxplot-chry-2](../figure/boxplot-chry-2-1.png) + +exploratory analysis: heatmap based on top 100 most varying transcripts + +```r +bluered <- colorRampPalette(colors = c("blue", "white", "red")) +varList <- apply(exprs(eset), MARGIN = 1, FUN = var) +esetTemp <- eset[order(varList, decreasing = TRUE)[1:100], ] +esetTemp$donor <- factor(esetTemp$donor) +exploHeat <- pheatmap(mat = exprs(esetTemp), + color = bluered(100), + scale = "row", + treeheight_row = 0, + annotation_col = pData(esetTemp)[, + c("donor", + "stimulation")], + show_colnames = FALSE, + show_rownames = FALSE, + main = paste0("top 100 varying probes", + "[Euclidean distance, complete linkage]"), + silent = TRUE) +colorName <- textGrob("z-score", x = 0.5, y = 1.05, gp = gpar(fontface = "bold")) +exploHeat$gtable <- gtable_add_grob(exploHeat$gtable, + colorName, + t = 3, + l = 5, + b = 5, + clip = "off", + name = "colorName") +grid.draw(exploHeat$gtable) +``` + +![plot of chunk exploratory-heatmap-2](../figure/exploratory-heatmap-2-1.png) + +```r +# test association with stimulation +tab <- table(cutree(exploHeat$tree_col, k = 2), + esetTemp$stimulation) +fit <- fisher.test(tab) +print(tab) +``` + +``` +## +## 92TH023ENV DMSO +## 1 29 122 +## 2 171 90 +``` + +```r +print(fit) +``` + +``` +## +## Fisher's Exact Test for Count Data +## +## data: tab +## p-value < 2.2e-16 +## alternative hypothesis: true odds ratio is not equal to 1 +## 95 percent confidence interval: +## 0.07483947 0.20646421 +## sample estimates: +## odds ratio +## 0.1258127 +``` + +```r +# test for clustering of samples by participants +participantList <- esetTemp$donor[exploHeat$tree_col$order] +y <- sum(participantList[-length(participantList)] == participantList[-1]) +# derive p-value using a 1000 fold permutation test +B <- 1000 +yhat <- mclapply(1:B, function(seed) { + set.seed(seed = seed) + participantList <- sample(esetTemp$donor) + return(value = sum(participantList[-length(participantList)] == + participantList[-1])) +}) +print(paste0("permutation test: p<=", max(mean(unlist(yhat) >= y), 1/B))) +``` + +``` +## [1] "permutation test: p<=0.001" +``` + +create a ENV-DMSO ExpressionSet + +```r +# identify complete pair of ENV-DMSO stimulated samples +flag <- dplyr::select(pData(eset), + `Sample name`, + donor, + stimulation) %>% + spread(stimulation, `Sample name`) %>% + filter(!is.na(`92TH023ENV`) & !is.na(DMSO)) +esetBaselined <- eset[, flag$"92TH023ENV"] +exprs(esetBaselined) <- exprs(esetBaselined) - exprs(eset[, flag$DMSO]) +# save ENV-DMSO expression +save(esetBaselined, file = file.path(workDir, "output/rv144.esetBaselined.RData")) +``` + +differential expression analysis: stimulation effect + +```r +# create list to save linear models +fits <- list() +# identify samples stimulated both by ENV and DMSO +flag <- pData(eset) %>% + select(`Sample name`, donor, stimulation) %>% + spread(stimulation, `Sample name`) %>% + filter(!is.na(`92TH023ENV`) & !is.na(DMSO)) %>% + gather(stimulation, `Sample name`, -donor) +esetTemp <- eset[, flag$"Sample name"] +# identify genes differently expressed separetly for the VACCINE and PLACEBO +# treatment groups +TreatmentStim <- interaction(esetTemp$treatment, + esetTemp$stimulation, + sep = "_") +design <- model.matrix(~0+TreatmentStim) +colnames(design) <- gsub(pattern = "TreatmentStim", + replacement = "", + colnames(design)) +rownames(design) <- sampleNames(esetTemp) +# fit a fixed effect to the treatment/strim and a random effect to the donor in +# the linear regression model +donor <- factor(esetTemp$donor) +corfit <- duplicateCorrelation(exprs(esetTemp), + design = design, + block = esetTemp$donor) +fit <- lmFit(esetTemp, + design = design, + block = donor, + correlation = corfit$consensus) +contrastLS <- c("PLACEBO_92TH023ENV-PLACEBO_DMSO", + "VACCINE_92TH023ENV-VACCINE_DMSO") +contrastMat <- makeContrasts(contrasts = contrastLS, levels = fit$design) +fit2 <- contrasts.fit(fit = fit, contrasts = contrastMat) +fit2 <- eBayes(fit = fit2) +fits[["stimulation"]] <- list(fit = fit, fit2 = fit2) + +# differential expression analysis: HIV case vs. control +TreatmentInfect <- interaction(esetBaselined$treatment, + esetBaselined$"HIV infection", + sep = "_", + lex.order = TRUE) +design <- model.matrix(~0+TreatmentInfect) +colnames(design) <- gsub(pattern = "TreatmentInfect", + replacement = "", + colnames(design)) +rownames(design) <- sampleNames(esetBaselined) +fit <- lmFit(esetBaselined, design = design) +contrastLS <- c("PLACEBO_No-PLACEBO_Yes", + "VACCINE_No-VACCINE_Yes") +contrastMat <- makeContrasts(contrasts = contrastLS, levels = fit$design) +fit2 <- contrasts.fit(fit = fit, contrasts = contrastMat) +fit2 <- eBayes(fit = fit2) +# top 10 differently expressed transcript in the placebo group +topTable(fit = fit2, coef = 1) %>% + select(SYMBOL, logFC, t, P.Value, adj.P.Val) %>% + print() +``` + +``` +## SYMBOL logFC t P.Value adj.P.Val +## ILMN_1896145 --- 0.3017263 4.652130 5.976632e-06 0.2828321 +## ILMN_1711883 HTR3D -0.2017627 -4.159871 4.730030e-05 0.9988950 +## ILMN_1677374 --- -0.1844982 -4.068273 6.819273e-05 0.9988950 +## ILMN_1717324 --- 0.2424789 3.875968 1.440693e-04 0.9988950 +## ILMN_3239480 UCA1 -0.1652911 -3.821776 1.769920e-04 0.9988950 +## ILMN_3299441 --- -0.2514595 -3.790279 1.992803e-04 0.9988950 +## ILMN_1691361 --- -0.1809396 -3.727657 2.517173e-04 0.9988950 +## ILMN_1685045 CHI3L2 -0.1708355 -3.724162 2.549989e-04 0.9988950 +## ILMN_1725139 CA9 -0.1671642 -3.704351 2.743727e-04 0.9988950 +## ILMN_1684897 SPDYE1 /// SPDYE5 -0.2025732 -3.684547 2.951235e-04 0.9988950 +``` + +```r +# top 10 differently expressed transcript in the vaccine group +topTable(fit = fit2, coef = 2) %>% + select(SYMBOL, logFC, t, P.Value, adj.P.Val) %>% + print() +``` + +``` +## SYMBOL logFC t P.Value adj.P.Val +## ILMN_2140990 CAMK1 0.4370108 4.388951 1.845103e-05 0.3459782 +## ILMN_2376133 KIAA1191 -0.1687398 -4.310844 2.554186e-05 0.3459782 +## ILMN_1759787 THBD -0.2611175 -4.301246 2.657524e-05 0.3459782 +## ILMN_2093027 MYO1B -0.1768234 -4.278022 2.924398e-05 0.3459782 +## ILMN_1664602 --- -0.1857514 -3.939602 1.128235e-04 0.5696287 +## ILMN_1655230 --- -0.0997451 -3.935295 1.147168e-04 0.5696287 +## ILMN_1787680 VIMP -0.1504375 -3.926250 1.187913e-04 0.5696287 +## ILMN_3238859 FAM120AOS -0.1350209 -3.852309 1.576563e-04 0.5696287 +## ILMN_1657234 CCL20 -0.7692402 -3.850290 1.588705e-04 0.5696287 +## ILMN_2105919 FGF2 -0.1026889 -3.846069 1.614375e-04 0.5696287 +``` + +```r +# save MArrayLM in list +fits[["infection"]] <- list(fit = fit, fit2 = fit2) + +#save MArrayLM list +save(fits, file = file.path(workDir, "output/rv144.fits.RData")) +``` + + +```r +modelName <- "stimulation" +fit2 <- fits[[modelName]][["fit2"]] +for (coefName in colnames(fit2)) { + coefLS <- unlist(strsplit(coefName, split = "-")) %>% + gsub(pattern = "\\(|\\)", replacement = "") + sampleLS <- fit2$design %>% + as.data.frame() %>% + rownames_to_column() %>% + select_(.dots = c("rowname", coefLS)) %>% + filter(apply(., MARGIN = 1, FUN = function(x) any(x %in% 1))) + top <- topTable(fit = fit2, + coef = coefName, + number = Inf) %>% + mutate(SYMBOL = gsub(pattern = "^([^ ]+) ///.+", + replacement = "\\1", + SYMBOL)) %>% + filter(SYMBOL != "---") %>% + top_n(n = 50, wt = abs(t)) + exprsMat <- exprs(eset)[top$IlmnID, sampleLS$rowname] %>% + (function(x) return(value = t(scale(t(x))))) + breakLS <- c(-1 * max(abs(exprsMat)), + seq(from = -1 * min(abs(range(exprsMat))), + to = min(abs(range(exprsMat))), + length.out = 99), + max(abs(exprsMat))) + colAnnotDF <- pData(eset) %>% + .[, c("stimulation", "treatment")] + rowAnnotDF <- top %>% + mutate(`fold-change sign` = ifelse(test = sign(logFC) %in% 1, + yes = "UP", + no = "DN"), + `-log10 P.Value` = -1 * log10(P.Value), + `adj. p <= 0.05` = factor(adj.P.Val <= 0.05)) %>% + select(IlmnID, `fold-change sign`, `-log10 P.Value`, `adj. p <= 0.05`) + rownames(rowAnnotDF) <- rowAnnotDF$IlmnID + rowAnnotDF$IlmnID <- NULL + colorAnnotLS <- list(treatment = c(PLACEBO = "blue", VACCINE = "green"), + "fold-change sign" = c(UP = "red", DN = "blue")) + pheat <- pheatmap(mat = exprsMat, + color = bluered(100), + breaks = breakLS, + cellwidth = 2, + cellheight = 2, + treeheight_row = 0, + annotation_col = colAnnotDF, + annotation_row = rowAnnotDF, + annotation_colors = colorAnnotLS, + show_colnames = FALSE, + main = coefName, + fontsize_row = 2, + labels_row = top$SYMBOL) +} +``` + +![plot of chunk heatmap-top50-stim](../figure/heatmap-top50-stim-1.png)![plot of chunk heatmap-top50-stim](../figure/heatmap-top50-stim-2.png) + + +```r +modelName <- "infection" +fit2 <- fits[[modelName]][["fit2"]] +for (coefName in colnames(fit2)) { + coefLS <- unlist(strsplit(coefName, split = "-")) %>% + gsub(pattern = "\\(|\\)", replacement = "") + sampleLS <- fit2$design %>% + as.data.frame() %>% + rownames_to_column() %>% + select_(.dots = c("rowname", coefLS)) %>% + filter(apply(., MARGIN = 1, FUN = function(x) any(x %in% 1))) + top <- topTable(fit = fit2, + coef = coefName, + number = Inf) %>% + mutate(SYMBOL = gsub(pattern = "^([^ ]+) ///.+", + replacement = "\\1", + SYMBOL)) %>% + filter(SYMBOL != "---") %>% + top_n(n = 50, wt = abs(t)) + exprsMat <- exprs(eset)[top$IlmnID, sampleLS$rowname] %>% + (function(x) return(value = t(scale(t(x))))) + breakLS <- c(-1 * max(abs(exprsMat)), + seq(from = -1 * min(abs(range(exprsMat))), + to = min(abs(range(exprsMat))), + length.out = 99), + max(abs(exprsMat))) + colAnnotDF <- pData(eset) %>% + .[, c("treatment", "HIV infection")] + rowAnnotDF <- top %>% + mutate(`fold-change sign` = ifelse(test = sign(logFC) %in% 1, + yes = "UP", + no = "DN"), + `-log10 P.Value` = -1 * log10(P.Value), + `adj. p <= 0.05` = factor(adj.P.Val <= 0.05)) %>% + select(IlmnID, `fold-change sign`, `-log10 P.Value`, `adj. p <= 0.05`) + rownames(rowAnnotDF) <- rowAnnotDF$IlmnID + rowAnnotDF$IlmnID <- NULL + colorAnnotLS <- list(treatment = c(PLACEBO = "blue", VACCINE = "green"), + "HIV infection" = c(No = "black", Yes = "red"), + "fold-change sign" = c(UP = "red", DN = "blue")) + pheat <- pheatmap(mat = exprsMat, + color = bluered(100), + breaks = breakLS, + cellwidth = 4, + cellheight = 4, + treeheight_row = 0, + annotation_col = colAnnotDF, + annotation_row = rowAnnotDF, + annotation_colors = colorAnnotLS, + show_colnames = FALSE, + main = coefName, + fontsize_row = 4, + labels_row = top$SYMBOL) +} +``` + +![plot of chunk heatmap-top50-infection](../figure/heatmap-top50-infection-1.png)![plot of chunk heatmap-top50-infection](../figure/heatmap-top50-infection-2.png) + +print session info + +```r +sessionInfo() +``` + +``` +## R version 3.5.1 (2018-07-02) +## Platform: x86_64-apple-darwin17.6.0 (64-bit) +## Running under: macOS High Sierra 10.13.6 +## +## Matrix products: default +## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib +## LAPACK: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libLAPACK.dylib +## +## locale: +## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 +## +## attached base packages: +## [1] grid stats4 parallel stats graphics grDevices utils +## [8] datasets methods base +## +## other attached packages: +## [1] bindrcpp_0.2.2 tibble_1.4.2 tidyr_0.8.1 +## [4] dplyr_0.7.6 RCurl_1.95-4.11 bitops_1.0-6 +## [7] gtable_0.2.0 pheatmap_1.0.10 limma_3.36.2 +## [10] ggplot2_3.0.0 org.Hs.eg.db_3.6.0 AnnotationDbi_1.42.1 +## [13] IRanges_2.14.10 S4Vectors_0.18.3 impute_1.54.0 +## [16] Biobase_2.40.0 BiocGenerics_0.26.0 readr_1.1.1 +## [19] knitr_1.20 +## +## loaded via a namespace (and not attached): +## [1] Rcpp_0.12.18 highr_0.7 RColorBrewer_1.1-2 bindr_0.1.1 +## [5] pillar_1.3.0 compiler_3.5.1 plyr_1.8.4 tools_3.5.1 +## [9] statmod_1.4.30 digest_0.6.15 bit_1.1-14 RSQLite_2.1.1 +## [13] evaluate_0.11 memoise_1.1.0 pkgconfig_2.0.1 rlang_0.2.1 +## [17] DBI_1.0.0 withr_2.1.2 stringr_1.3.1 hms_0.4.2 +## [21] tidyselect_0.2.4 bit64_0.9-7 glue_1.3.0 R6_2.2.2 +## [25] purrr_0.2.5 blob_1.1.1 magrittr_1.5 scales_1.0.0 +## [29] assertthat_0.2.0 colorspace_1.3-2 labeling_0.3 stringi_1.2.4 +## [33] lazyeval_0.2.1 munsell_0.5.0 crayon_1.3.4 +``` diff --git a/code/20160509_RV144pilot.preprocessing.code.md b/code/20160509_RV144pilot.preprocessing.code.md index 54c088c..0451e5e 100644 --- a/code/20160509_RV144pilot.preprocessing.code.md +++ b/code/20160509_RV144pilot.preprocessing.code.md @@ -1,7 +1,7 @@ --- title: RV144 pilot gene-expression analysis author: Slim Fourati -date: May, 9 2016 +date: May 9, 2016 output: github_documents --- diff --git a/figure/boxplot-chry-2-1.png b/figure/boxplot-chry-2-1.png new file mode 100644 index 0000000000000000000000000000000000000000..0f7cf5da2799c3ed787de42b6747454cd2d9a75d GIT binary patch literal 19945 zcmX6^WmFqo*Cn{SOL5nr0g4Bg;!s@s-~^ZA8e9vdxVx4XcWWu`#ogVVFMYpX$y#^T z%$zIxoU`|yD0MY?EOc^oI5;>gMFp@X>~{eU4nTv74Er~i*1-t}he&2CBcpCFBQN9Z z;OwdeF}Jjqw|262wYAigmxhA_#l`6u+feEdinqD_N6#>(<-i~AUM7LXh?SztpWoOT zxUwSXl%^=SNVjiwdph>@_ZWz;?Ir(%thyiYFq5#qt0{_FF3S(~<}Q#NTQ-Nq9HiA`_ru?1^EK zc_hgBdKy&h#^oyTf*&xvei*LF^nNz|JG=6lTiql_MMUF5ZSj4#3PwWDWdm9@4$`#4Ea?hxn_KyBP zDc6_MGL!HmkQ0|zztxuDoq>(2mQX4>bBn5fMAHWBg$Z79Y;r=!h4zNEdKM}xQU}yD|~X`w|^Cu#)eFQu!7D?na z;su-XP;6vg!fHaYXw&>pBouqy?=M9fGe&zldsAP3{W1-GXdY{1vx*cCA#Nalj;|t~ zw!WwgH?ky<2;dNcmtZ2J!T%dNZo;d0o8-a^p!11-uX>w8*+5`*_PC(rAZzkoXzryw zmQS-JQ+7qTu`SQ+jgL*mNNc`|1N&zL{U6pRB;$_fLDPDp0}dgLpPOXat$Je*rVmHY zytnMnc17!f+jfGRh zP`jB9-r8M}8?1dM_91Qh?;y>-N{N=vjhB@ygGbXB~w8lbkTJ_NraV~d` z7!|f9^=|2tQ74x!AnK_e{#x}_M@RbrlFCY@>#Pdulp=^_kg@_UJAo+JFzGu3RlQ-kK=l;x^F{Ksh#rz5OUA}{5s6btfNLAJDS zhx^bn(uiLw@itnKi;e~qNN20p3D=ptZ8!vmT=H2@nBAQc7spRroAHfsX_P_1L-bCh zvAbpevzM1s_fzNU@x43I zy&Usg|2g5KDxd0y0!au26NLO268v`br%b%evKyF33W2{O=>bm2V}gQj@LI$hExXZ1 z)&N72)WGv1@%WXL{IkxdQ^=U`-nCVg6!vBA+p|3PUEtf^TAkWOd`DcD$fm>1tR7Bf z^-qK8>iGA;D<1R4<|jYqvHx&)`Au)j`}F0$UNdK}^5s5^XyoiFjGj0LJ{jI0nZA(| zf3!+U`li?R`}t|@4SR^n%Sa3@HghXR@a_4ksck%I(fR7fj)JNam>Hqw*cn{9aMNfM za6oPy5Mr6<276zRcy@^?;Di8jR(GphM(#T4J|%*hS@za&f;*;r9Dup5Mq`&5ge)=? zM4-Rdd`VA*U~Fl&{bFXkrr`?TSvLJaLhSwP{g!##$xnjSA(S2A=nz^!eYdkfa$pwd6duzFs5*?o+wu*FNf-k-mm1$KF=jW97hMi`jqmC|i~{4bkm?ETdU?jw;{7x1R@jkLp89rHK7xae6spHaZW zjuzG2xEH2-cyNvjkI5&KL(2p)UDM+wr?zTyI=8(sER=WNX5NE*Ie53rZsTzO(GEh^ z{O_A2f@5-D&mSe=94F@g!#=D4Z4`)qCQ1EgP^j zgU0u11odw9bM32?6)ws7^PYqR4VSbTF2WGp2Uvg}m&;%wIE~I7PG_{W*W0EW@><{7 z%3j~2>QbD#Oc$olRkZC(J|Av4t0@0Yyp;m6ilKyGaT0rY-A@z$wg8$t2oZ_ywhr1uO{+!F+ zTXByQCNGD2UT%k}2{&eIry4drohd$ZE2h<6|DCHf4a4n|1>*ra&U-PAs(@zCuw`^^ z4i!CZra!J9dS~~0$+2mr*CK|q#ONp(O9w+5gB;=X%MmvU=vc zhs4anLb@pq+s=ZzC|vvbjTCV&X=j=thVgj4>m!W0vj<^wg@5V_Rj^s*6t#~fg z9XI?z#{5XyMw!zsIgWN5mGAE;VrB;0fEzPfAF!HMeJk&r^Rj>b(803g#=2=I)Mdb2x#3%J$T3Xd@*28@hTm6m~6*m_#=>Pf$I6;iSa@(w5 z(|TA@H_TQUCW9?pbRs-`0LNT+aNYO>sU`gZW%ei=42prydS6t>2j`e|qg656TuV_R z%);=faxli(1M?aYc9zAY3UM|`VV4r<3z+y-TIMms#*fg*{c(pR{~h~^i%ml+1hF32 z_=k!8P!Mtja2-)Mj${pXj?rM_z9Dkd{Xq49SYX$FuRA390rv891pJRQ@xgX51ku?V zXc9s84S|j4IbQ2qS z{l_m_9}%;YIgU@d!hF{C_8Q?9a8qd-p5DIAh_&?i5Zd~=F}t`=DZASwxpnmOS5QeC z`e?vm#CJuZ^(Ctm)ka6$%?)F(m16%6Z#ZSnA2U3;p{7efbqY)?VN%C%o1SQ&*J$<{ z)T5dFTMw9SG~S=seyH52eoUo82py7S46g3DT_hJ6cNc%&z!4ffd>9W5 zm=;|7x2F7-dCA@X(nvS;`Fbc$#P=}~Wl_@3%P+N<)X(cEeQr5dS-VMoG5=mVEVv3g z;sv@SMlEkK3|js`~5GZSlHT{rlu<63D>3lp)vVi-DeS7?$glOWo%c=BE{=C4Gx<2fU zWMwWgr1nS3K{$l|rhlO7Wi?^3MoO#p0jz#oU81f4#8mGK#M`8XXyzsD> zOa1zK*A@A}|4>V;ITdGv)Egeb>Q^8vA!4Q?gHP6F;se$=m;J6rvCz7U(8SKW5zM0) z)Al{fIw^B^%u;RXGylS>7<)|Ikh@5U=lpd+(kUhc0s;bBEWhTxoX`J2TQDR4(N!%i zG~a79OPFsYQ;55LUzQu4R18l4veM+HqEqD;O7=TaeZPwVie#TZQ61S98wk|*LD65= zo_B?2y*@9y48>{xluzqcs48-hB$GS9(8$Srt3g3#k*e=OrU?6&0Y%KCm(_@QT+fA~ zVS6-|V&4pxtjgIp!hP}AkX}2F z5VG;)bBTWU8|YNq135(z2#v*p+8fdS>EWPuh~fM>@wDc-xW13@0)Uf24FJ}Cc^1u~ zZ2vAMG}vuSSl+|8d0blwo{pHSt{WanvwY+C8c)6^xacqZ%{mn&!pk5W(wBV%wP_UT z3mYESSsgu9wq@HmfjQL1d0$upc3Sg#RYCFPB#(x>*d7RLq*fP6IKNt3C2H7M2cHQk|LH!xNx2mn5#E5)0z(uB`v!D%_|CfW> zVHb|)v`;Cqhr_V51{A^eA1(&91Jr6#mP7|$bU!aNM3$Ba<;@wX25oEIzqo0>7lWJ9 z304U!mdh>oP0CYO-`Ee&|4!xzcLHUk2Jv^3RLHNY25;H8xx-fCP*_C}X{7kfdq}i* zn7x+r=ig5hE3%z`VTqAE&Dis=HfM@!he6bVzhd?Sip+nXY$H{U!+f# zwSy#0f6LE4valrlg~9~!3Tls2KyM4IYP|6jf1FIP%i21utX*~g?05WtYT*p0KNXtg z)o8fB;yh4lfKOU7mf*KdZ;*h%CH?S|QW8}K-w7q84z)Pwvh`|Av=U^-%ZKQBRNZ7p z2H%gUFV;Q(t9cEU3SFyh2GUE^%n%2awFag$U0qme5kMaW`Y-DP6#}M}dL5k3BZ`kRBwGFk!b{{(GWt zWqf~wdJz`htLXQk7OEh6f+ zpiOQ_LfCPIwMaP9zYi7vqXONphkA%cADc%~m6)@K)u|gMfki9efvXjjH3qmgmG)1! z2XTO5bahgcn}MJTZ9?l1S~Juwsx5!EJ`B@)X%E1{Eo>^=aF4vb{zC(({p5nieQc@P zw$F$wAdoGeoV-K@gE!P2PvakLz zrHXY_cVGmu&baDHa~l%EJe80`N4RDMJ#44^_%CaJ-W=Am9(xDy6G;*EG*7T<@6Ze5!P5WzCg&daMAH7!ocR<`;ByjkN%4y7@KM>y9t5M zTM^2%R8c?b%`;z`VZ8%CAGI2C{T)UF<;b6(Hy9^dKrZh00V+`RzR+eyawkl+&{+JhgNtsmshwrEIu8JoIpHTbfVzlKl~1k3E*di4Dq7Z>>$Nw$cGV_Wm`?B0n~= z6l_W59|0e+2S?B4-Fbzq0|yJ9fV9+Z?mR*J$a?D!flo_obIWW5Xc>%;{8>Tilo%$% zqsLqd5Qp3<6TZs52u!QP4iSb(Fx2KAY}V>ZGHMSR2@jN(5-?=4_O zJ&C}3p`k;70wyQB{d6E$K`yOPd@G;NIs|+}KjA>VNCjB>5zjY<;==Yb5a`Z<-lX{i z6iY%j*1lG4h>vph_&eAKPvSz;cT@waB#DrG7Q|Vz!8DS(drsMe0%!#J4=r&rt+spN zLP$D`Zfo;;2!`AS@ql@=3(o<@hHO1tMT?OIa-rWZdJU3nu{$;{txS^?InqCT4sPVm?Gj=1e%&z=#aR%b8>=C zdn2n$Jr^arQvC$}XFzjmoAr>L2*ir_9-=Fy|$2i?eiIFgq zU_&iz6Of#N7jAfXf8WKR)Md`$SKAv^*gqlYkUNr4gBIJ&2{b;W78)fP3iuv&j*J9* zV6Xd@$2uFxN(-K)wlBZGd&1iy)6}BgUbqArgBJ;5rYW!?bopV zF!{SKMTAM-?28Dgi;-?}T;dzAHf#Kj<+E1@P;wFPlear#Rk)5XmPi6v!`-8>Web>x zJH3dOx0FvpzD~fmWhe{^v7ErFXTyG8ag=ZPH9qL++_IDA8F}Z|Ck>pKL&9*r!zzbe9S$ z%1_=w(~dymiG@xRkWZ6$k1nkIk7ZVf54}qpNin znb@Zr5gft63Rijy0xlw+F(A?T^qbzG9fOTFHUzR`3r0-S$Eoz;w&=B@G~60(B*bma z=7>32zwRh9hGN9T?Pa=?+_il@TSc^=#X%}D41aj21;5g#Ae{?=v6BZ9z;#4 z_-|@E$U1W;#mB?3g69fLnSYxVx4MHJ)AK=;KnP+qjwj=wih#AK%{Lf;F`dl9L&})6-%}+ zxvvV9){$iBke#N6bOYE9lVbm36#-*T#^_&h14%Z?M94*m3m7=#|2q8;B9%v#jobpo ze*n~fA7)8@$M@w%37-D?&7<~nwi-EEv(oSEVcD1E$w9cD<$XnDe$#D6@Saqu9m>AM z3cg+&@Gm!|TIH(+W>1}YD-F2M&+D%>iK>j0d|{rZ9b=I`*PZN1JYh3)!RZsy=q_!` z)J7zA4J)erJ0B~V3aKupjjecWWFCd0errDZldxs(*LB1oe{VN#b0Cd&62V3sE_R=1zk3e(NulOB-xMWCvXSQt_BZouXmf3QGPpH zb>8^+NIjPR5oO2X71A9b>Lvmw_zL}$;1&Ba2A*dEqeO51PSjsNSifp>V}TSJ5?3yGY2ag1+7#r9wNGJTP&V2 z?{tu_>)13%dYdHNMG#!1u8S!?DZ_8|;u4s%)agm9sB^H+4p>nR5+Ljg6AH5(fgs$& zSw9H&Kx1t0DI!$<`@Yt`-T~*%?)z6KD$x2g9AUAUt$#e=(2Eq}oGo{`B%ICh9Z_yp z`~8Xbds(Z06~I z={mc~t^CD1_0ueCbW-6@<2ltHSoMo|D5#D2IKtFi@c9O(X>llcS`%~4N)HB}My=r9 zf&-IBpu*jf2$fbDdx$yLPL5|T#J;x@wTmqf(;@aZ*Z1~sAwvDG`M>ZL6oCYERZ&!K zq#rP!VN9C+rnWu~<6Xr;L=+&5g$5f0r)Cp9*$9uCPjqq7;O2;$c0~J06I@DkAsK}0 zodM7`&a}V0gVA&vE$jmtT6*wcS%0}l)xqx6R(7|1KIyKe&IfJ^dQk@U;%2FI0; zWnST#PJ3HM$kI3lb}VaJXAGe#b}t0!3Jpt*0{0 z)hBpdi=SIOYClj0g!aLxYX#w}{{Ar6%8IQ8WT3Ro%6p;H{MjLz%FrL*!knowk-aA` zq>?+O5D_d?B!W}pbh(GMYhuP^hcSG~g+Dg?m$eAXWX)L93}NyG&b}u_7)yhvk#RZ3dAk$#{}2XS&_MY()Mhk*;5ki^ghV)%Qlj{O3^ahh5;^l5~G!7A@vTNe=LLQBsHO~O?EM3yllM-DYIfsM{Z4Na4ow7ms!7gGBNp=4n)gv21%(IRFogIEVFOK3b{ z_&$b4noAjj_4zOQ1n^+8D>BxeOBO_iNFt-Ke_0LrgEpzbW!P-kr_4?%Z)#;nUh-7-h~Ywxx09xg ziyS6n`&V!nAvt4YA&VJLc3>ekCP415Q;qznTs%pNzY9BVxjxp9SgFbc_$>vh#GDzW zEdF4{p;#8I@y6ZY^f9g+Q+akeeA(AS^wiv4Mm02sL)HKWMl%_$f**1Pbv^6k{Z@OWsQ;qe&=!Ky6C=~(6J&*bdccC zR|UxxF>SVoXDOM?0J-7<6f1=Mi+uMX22}@?V0d75=_*JoV{*!5jIk1V&5S9uuUEpH zA4uhQ-gCGv`_|iR*)td$N|xhehOtiRfPhT=Qh)Mu$*vKtCTPZ6+JWJK<$kUr7QFAF zt=v~LLGl%S2Yb_q1P>*6=g|J#dKIP}##DT`gl8k^*3NaMS#nZ#XNW?zT_=UwS)e&A zd>UFx9@dz$<&Bn-xe;n6uwV^$n>ZjP@gI$EDpve>+t!{cx#PW?@t1Mup+j{`OfxdVB=!loe)K?#a(cN$g}{ z4^5!+gLtZ{Ek{-Y(qzy+b~6Pq7(zjdKDyp59Cz?eEkfWsCPc=*S8sL5*#hV98;lWF z&y-lN8f3%v+WvM;Te_jhj%_UPsdT9Uad6aIT)(u8GIy&8EAg?KroiWF*EU9X^WJ7w ziCKAM@oSqQz=C4ewc^dfq7;jhPM>^Kjpsyg$)RvQI5BvD*idCfx{AB)xM8sepZaT5 zvI92)7{13XD;T0@ILkQs&i#lXy6-R=_de#=@B%c-il1*hM<$+gs`9tp@?V!tW+*5(+G6h^>d z{F6$nN?4s=Lh{HzqcW=i;(J78byzU*U^HVH@)Le>kF#bWDLqnuKU=rH^vqN2+W=qk zc)0g(0dEOB)(rzz9ziu5_yw&Rs7Clqm5!dVfLOcGk+9}~d1Nwcm60LD38vZDF^y(V zANh~<_{L2>dne}d5ldfj}V2#znG^6(SHF%)={N!ZA6gr zowi?%MQAZ(cJ^?L3pZng<&np4MHVfAeHVHri&76KhcAg!WFy;~eg|Qe!ca47g0Lpe=AvD*ejh3AZZi&dys2Y{rvI3*kqf+5TWa*n_uGSKcBnCTKkL ztM~pTLb#t~b)EU+ zr1CUx2`neH#5%e8w{0n`Md&BAio1rJ%yn@_&>H${;1-%w39E`$9+e(tST1gz1M^goR7n zu@OOG7g&-N^^+v!5*QH!Xb8OU*9;+Fuf%Br2Dv$M3V`5XI)Q)Drj`$Sig1>~BiqsI zvIHttA(>H(J%mfr^8HT)+?(YVH;?Q)Z0W)Q?XvIS>Kd%azoYEVd? zj{W2Gfq(Jdt%F$UGFr00tFo1Yj@1qZuHs|gIS&hjekhOUW$nmsmoJw!nJTJ?#zq&+ z8&{Q{Uc;jf6(`;qm?M>WHRUDXP@-%wmeGDzh%BiE^%riJsMs9DeX-BHv%W zxf7nZdJ*rnk8R>GCECA_Ya73c!6LLo-%M&T=&RaY%OEayP)$j@3u7;A&TICgzNICf zS{gSS^@*XTA0CNj($^_bjV#ay0IW>G@Tm%^YP=@paqF-uv|B0N;BAB@s4V<+R#~!$ z!SEdrSDBGk7mirwyRNQrLhRN*(Znv@iJl|+sqUN=tN{LF4wGW*#?1DEx)tV`N-!SY z3;1DLP6dp5GkRc&;g9g5+46_z>y%mB`c#pCkMNe# zn?_u_;RnQ(v3>UzLAv7lD2yrkTZmDv`y}TV5qVf5kLw_XyP>R_H>%&7HtaFtz_r0I5kKNVf!}F3-@c_qC z?>qEGWMS*tmIx^g{T?BKs|hc}`_5!a<@{Bm_Ezys;X1g&hHOAMWo7Rm> zRA$9s1b!jEr~c(CYwHfd#3V;(|6IO?4kXTW`=vJCRQ z=$Etz!n^9&_4-oyj3QR*IseCiV>x;O+rKk#f|SOpGAn|2U)bHPg9F8qh9j~7YzHNz zLSmIZ_@{)rN5B-TjW=l(hJ_T1m!hDnr4Vt^wlL})(yyJQzwk+SqJIx1T|)+N#eu+S z`|&zmmK2_J2#kT%J{|OudR0gFK4z`_A_ppc@?Qa}|A%RX77N9SvHO`>bT3|xu0!)B z=m$RT+c$vBnyG{@B3UJ`gl1_iXtOD-)5MeKM0TTkv5M7;*kTQGb{IPn(ud)v*ca@N(eZn5jm0sbcKLF+RA} zJkNmdY0WrsEN%ejB?``zBU*dQ<&XZ#yYJ&(`lP>Hc^c#uk4lo+Bvbef1WAkR^S&KN z3fOklSmts48>JH*9+1T)s(%&!C#0VITlA3wZT$f+SxmY)3F(}5J5b2M%l7|_8-tR% zV1`gK;3V7y?n{h>lKq#;y9Jj{cbLK8H3qirHA>JT!Jj+2SJgK*VVR-5@bW& zq~~A3QsM;JZT3YfZcxp20-fuQ@>OQI1Z$x{)=`s%c}_c zHWn;_j22 zL9+jBKuf?O$}PmXKPZCX^$mz~s*{cgWT$f$Gvj>Z_)1G zSiS_+51~i>1T71RBC`P}(n9fmO93=vw!L@QP0LGmh9~V>FJLVirs><`{13UUopjS0 zV=J2sdlCF4XlmYd6Ch)?=q!;>zCR|^G7lZ$Blf2x%B|&2)Ls#oJz`cRo~v=y_lPzB zqoQ*na=ymcQc#ss0pm}YC?Yv8+q6pd$?xO8Suc@cK&ehCt8@-xfnKHiG7xGpQpZnu zkTPMyDi`@`6$b&D58@A2OFPO#k_gBRbabW{{D_f!ZYMG1ay7pb1m^VV)iKukrJ}33 znwDLI?FRqp4zsFK#7(myMgF$Jmm$)Uu>8OV17rsgu&erA(&Z{xZ(<+^y_(n)?2=2k zZz+dIV}A_iRVesdTbNXXbbG0~VVKRoUccmS{V4R`MlVrjKymzXZ;t{IhK{2;yguyV z=u=cvY59t0f1;sOinWY$K-!*%xlT<)m>B!J?mEtFUFLO%o^E~3X_U>hZwI;o94*_) zi{xpNQ)GOfY^uEpG7mlejDFt{KN$Os68`#W-$X}}Gnp|>Us}rxsy`$zS6t02$wJ{! zCndQI{$L1G0}axeLzSnb2F=E@gpKh?ZCOb_?w2Ex?EVGCVk1|q6}^q(?Bu?V)cTCt z!|(cBN;U%l*O(>! zE`D=Y?rn6F)WiD_WBOcRv*8W{_{uM?J74sE=0`Jr3NR2NTlQM>|GDV-JInNGah`aH z)`W@pJy z6mK5YIk;dITnqCTeX||`*|IaSF{DsNfm&dQ_S1pK)I+u3bC3^SW|PLI*|fRhn+$i; zWG6HMyt0cE{NtE~L1+fQpOUJ0!{>QlUV?m}l~`!erm%U@u2NR98fq~dPxsdJzDZBC zXkU%cRHrt=o?M&G%1iJ=3D z_LP0XS+G_MAsN2}z{jyM+@0QpAyAkN&D3R8TNiv5Z7~pg3At%CsvV<|E4u&Z<2xk| zJCd6WW#7y8;3wGjg}Hx`eSkCF!P0Lq)zh?V0_Ku2*LN8p`uqSmTG9^_Fop^D8>w1q zXlI|-o_gmDosIo&@ah#j`!^$N(Bx&(d*(f!FNP|A_4eaZPl`?^K}>Nkg1e3rWw0lG z`V`l7k^SN(cW~#SSwsCF#z@i>#{N2dAC0a*p)~A$nYh{xB^G}<_;94Ra$Ef8uIugP zo8I6B+<+6Rts5!TKM!C@RJd}pejXY1k$#1D%Q@LB0QYQyz-ME$J^o8#`U^~|nM**5 zpy>bk&sDOKR@RZQfv^b=X$ZV)Tf>4mk2y`OpeP-wAX%`~V9>&^+$|wEfZ0d>pkH4} zSlO7f-IE=;(g^E>U#_Qb$c-?1iyL@F-d)w77VLc57{R%Smb=y`m^X48wNNo7!$)R& zSbP(tTctUKzTytR6A!3R4TO*m9DWh4E`2Itqf>flkai#-BssbahHkq48Y4=Qg>on@ zXm5Mh)5!Vw6z=N}fTp^|qRsP3QCxImWi&#HZW5@vV9uDBt!ZH%Pa;V5HH;t+Xl=fp zQkI{TFHD`Bp2jsVzLo2%&tg7}VCknAKvEZY{^gYENk@k+E>7Gz1yzL#`yEWz*1w-f zLXHLMNG^)YWl?sU*A>wbYt4S0A!BID&e-g!4 z!-K~Oag|*t%(O!=EPbqAVGgxnjXx1AIXIM0jniscYT5BS{~%8nrO!9gyB8s{8Lj~W zl$l4mU9iN-W>d3*bHUA<6zmVl7Htw1 zAd3|!uC}jnts&FSXQndP+hdFb-B^DIq@*io@igc|NYaW0__A3XM~Yw5V*O=H2v=~v zGhE&4LJ14bYNz3rB$O|Uet(%)*sjf1@_gIWHJUv)?XN|Imof6st}g3ri(ia5fQN^# z!Hxnor2IooED9Fwm)_X5v$k&AS7rP5*TvJ}QPep76M2iW;4wj|2ce~ki+?|uBZV~! z(u>9>=ZfWmQ)a_Wop7C)WJ)cMho`}Wu~``PH553iylkABWBzyRbH{fQ#KEiy>S>bY zIiQ;|ueIl*bvD#(l;3F!6~;BK-Igk=-0W`N?+SI139YzKvrvKQx%G)8RW{EPENr5u zq$`S94f}gp(#BY@FIfh^R^-e1E;Z8gM9n}g9;psx4v2!IzLp>cb^VDY#6|05pR{kd zW69U1BESg?I~#fakJbFI2&G2f^$L7B;Ej_r-@nFVjElZ7z&*Wj4Dy*Yl%vc`*f5nzDCl*1F!Eua8JEo=8-fQT;HKELSg+@J%PpgA zZ+~n=9?~#t#O!C-q{~ysNqXLQ+(J!t#f_(iS^YoU$m9PMno4{^G%HEB`)*f|?kyZ= zMe_||#X+4%2b+nw<5wRE$ep~3%Joo6{c z4hxM!YT6KwN+c892IIa51{MtH zEArS0rHn_JoL7`VEhui?`46rfeK0IJ*R_tsg#AbDmQVuO_B4H zko$5E^zldYe`i77(__X*s-Fmfv0S1^i^9E+jpF%eq+ve7$P)+;D$XA7Ptmy3bbeSQU|D;lr8_!AookiIE z5WiOBEqg4=<_6F5AYWC|yBCDUZTkO?gt$&#`>F0ACGXuM(p=HQ$BEz; z2dg1ODga&(Uw#F>*pABlw0!yVk$5r`gp(azp$!e++@Jl0JDmZ_W%0+SscpX;F8GQy zr;L#qLXgi~fI#P2c$Cy3>5ohIzVo|xXk z9RlUhxQ<~_492TW1w)TaPqR$UGC5Gb`<(w6vS9H5&Lp!;9*}zke-}= zB~X*Dd~64N$p_vv$$G|sp!DL!S`#{XGr>rA#W`vA5LRaru@4fya^9fa=>)hTWpq^$ zKhC*Qs^E+89gNOdT|?XU)x4;fCJQ!1xo?V(uAWDHtdR;yqgaUUS6IOb?)8#5fs21KS$XLRkqOx~LeJHc_8?KKcms}Qb~{q}~J>R##A z@_oOob-3^!cT|0m!$=TlJq*@iDl@QrU0-!KG9;lk0saJM@>3#srSKP81XAe^CUxgb z-0qa?01nvhjeptpx~sB%KHRSq`H%3Yc`WPQ`@!y{Pn?#-%oJpyfH619Ay)Mr65JPs zHwD#)nNsmtZ7J0Mc=s&oa5ZX4;e}AFXAA&4a6HND_5|=O*nI*)Nyz^-Gw9jXQ=KBo ze+NX6=EMO;Sv1YsHN1LDuBeu7meW;-Gi2-gZ3Sz=B4EXJd$NmcivQT2nEv!kR5R>V zY*cdLXtX3))AUU)di_%gl3l54rwzySHGo#dWrtenhV2vVf4!D(jC z8tw1C_MXz@3F^jk6rJA=m|$%|E}+e~kG1Tcml$=D+aqKzIjI?Nr4)vN9|w2Sw~+F# zAw9}G{;xFBU5~Bi;J2PHJlR_e=}nc8`(BvZ%(xxFQ$`3WA#hMGC>iHbQ652>bw{tZXF-ihJ1nUe^y zX8a?pM5>kQ3Ukyfafe$JKpi}Nr}2k)nF6nwLYx?5s}zLmf%=>3g7gSk>>j>-NxP$& zr3Z}VcrxV5Fosq8cq%o;{r!~2Vh_yEVN$fY9)j!Axxz~SBXH4$$Kn9>ws?_aL#Swt zK%9ZQ>iv^%cPtqnzf}~sgi*{I3YB}}Azwecrf@yG%c5v$;#P`7}CealwUw#a0Pjeww&%M*nNrHe14|ANaprOWdX2jSm# z`0$H34JXVM-<{>zw|teZr(460R*p4rt_a6GT;{C*+fV=qE0C2rd9Rpq;#uQ1lhKxG zuxH&D@}5n-r=zZ)JZ~)uRH=4$?RU-BOfQMgR0Fscb>~k5y$O@7U$_;Rf|okiI-rF)tA2)SPYVXDM@)U||A z23rrTk!!fHPuN%2`lw3o+yX;dTx#wE5%u>$Gib0=#x{DbEt3qPA=1M&2LI$H^zrGY z$0>?eN5TN2r8Koa8{BeagIcs4gd_>G(AB*Of!K8ma+Nj%f*XDY*)fEqQ%XF0mOKS^ zu0&@0l2%2Wcuq;~HJ%dGqN8o%)^MotAl|UlVFgxNBIH+INAIw80kqZbiog|zE0oFRWA8`|f zfUoUwsGqIevKPyiLL0eeXefe<^}g(rq7MK>8Boh+X4G1MdxX}sN9pI@uPyy7Oan}@ z{2G`Y{2%6O$z@ff%iG6yjG~|6FQvh$-XoKjn{wriGyGhg?=kHTPpUPLUPJcQzp$aNv|^l;T@NYo@+ok$P>mj(Kp^!6&EB^U=o zM^)qpg-+4k>{_hjHn7epI%w&iIFcJo?fA?vqG0eqFJx)srAx+j9jt`R0saK{l<=#9 z$+HD{!YT0Ox0pS?LU4yN(4~afKm%d|EiE2SMI2r~0P{ITJ&9aL1WAG4G$}U_Bx%+F zA}snDA}rz;nFNyWrKb2)-(mjcFD~LSVq|(P*5NWVjVKr7q1gkG8{+pWbUo&S_=;Xe z%p%+!4w!Ocj^Um#qkA$o)kBO*}|Z9rj|6fl*;x!?-q zQdL1OIGBw~K9``q;ZI)JhW=jMgAh!^F_auIfE!{dgd<|1v=RZfUxI}j&MjXLQ$G#&Wa+ZCBIi z3g$;G!xXR#BI+WZm7!goK7x1C@u@p0gw!w$;x_3lWsw=<48|CuJFtvmKyest4D_79 zM)7tDA>!pUpllkE?B$GS_CxGg?y%-j5}Sy7@~}Du^e6*9Hn<=F+z=2p*2ok%9yfsz z8B5NiSvyVJ(Z8yI`w& zl;#2Fz^#C9kB%v=fW@mo%m4rhe@R3^R2qjID|(jJKYZ_d->W`BQcx(dZ3JiMLr^qh z#}7tXXwqOHXtv|(Akh*NKAewZFa`m{GD;C28^)KTih)s=<9nflNj{hARp?+;f0CEI z;kI!V(4B;ED0xnrQbKgPfjEZMpn*|?{hh@9VGQA~7a;yabtMpPMtY4{9Gi5FB=DIyjZkrn)IqMsF}eHuOVN)%16GmLc? z5xb$>5jd7w8dwDh#smGDsOUZ3cmL zk2)qqXAq~@Yv5ke2oVCOiSe{#&t6~8$D0}k##%%qMFne|0fjhMUU_Bv*|d}@Q>yEv z#AOFz9VE-Uinuri@e|?#Anbnms33Op8)6AXe8uemkPd|INP$#`(ZIkHZ{KLlr>C>P zFhJ`doVJ*{X)}a*yB~l)gc%)xKvOQ$%j=OILnSNxX>e2S4SL*5hEmA!abP>q9hf;9EOVl=U~Idf6uOt!jLg@JSIqQOnl($S0V5hhWDe6015RLldgdf; z(nP}3Ncx1}nhfhAVoe^DPqUsFa0096bMsl(1OrZBHF;1z&3axqGpw<1!g21l*6IjReqtpab|<~>aO68XYk-Ec-29Peorj68u+r=7Zmbu)L{ zOH#LaSG^>I+D2Uw>WS11XAl_gX%&n-fe~^<93PFo>jXs71MRx4)ajo3Jy$+TgqhXR zDNnbaD(G1g^s2$9^R_v2=IHIy>C&My3{Y;4EQFJ#LH)QuHD4Q-yYu^S3=A0qBQGKj z89$T0yr=GWk2;FR~!al{eY6vOWK#v5-8d+)t>L+>@)$8sF-<2qRSL|d%Wuy5wfnGJhsv~b~t7iy7F z^euWRUvR+%0q%y0o^CCdUV5qS-+1GVwZR`2vPU0%G|*_BO27sXEcFZ_;0-p|AlstW z&Z*;)o@;O`{9UB;kTbX9iYw}^&;W59@#7q#1i3SAN!a#H>WkHRGTw%sAr5Iamj8^TG=+9JbhE3q6u-Ww;Q9ZB9GwG~GsPS_p=hh0@v7 zgy+jlLvxy!L!eZcI{*ChmBSJGWd(5<;Hs;x(wZI83tGpfLyGj1*0IO!A01$hJn~4j z7Sdjd4-hb}_;CX3;0>m+!s3fBUL48SmzK2q?z^wXpj~#^WlB_O0|u+dVJ97MzyV?Q z?AdzW=9_QcP{xh1ci3SEB~a<;S@ip|bH{0kxJ{fUJ+61!X{Um~#uYzKV3T-zAs)WC z`|rQMa-x-1TB+?Qa5Oj`oX)tQokUZ;qmMpXV=fsx)z9phaatlK{^&C?|UE^a4X_K(CinPt{89v zo32mZ?|UE^Z~_|$n*Acv6$4ISHG1+VoNz+edFP#LG-IEt6PUjPj)6X6zzM96$nC|@ z#5*p({PG$ha`D9%7eoHw2S3p5GiJ;Pn{BpPQ78tV8a45kUV3Q_+nP6Tp7tnPW|?I) z)~^_#V&x}4`H418f5j_a5!P8}9gU6Tu&Qej8NcePt7@<^@rC&AuD|~JiczQeCh_Aq zj3@>kdg!6J$RdlxTW+}}KK$^*apjd)juQQ#vt40@72-A5ToWZKBJR5Du2JF}b-OhD z)8CS26D963N==g}p<3~#n{J8{Cz-v|jW*gSo_F4PdW<-}qr_ciuO&Uu;%>X`7I)ln z$GH3MyJyFqt~))FV}w7ja0xgM0mXn2K)C<@`@>rp@!RUzn$GJVePfo4niz}a4j+bUjO>nXRk#7G*K)IEU-Wj zVsAO_bl%|vHo{LI!Xd5003w0FgtKY?v89$;N_VkbZ@u;Gup;fAh!HD+)~c*nagwYM z8%RX^*=L^}h?gYghV&VFomE}ypa1-)hKiBXhL}h?ibEp?TnB5!$4LExJHQ!4Jk*e^ zxpU|0{LpNTW}9uoR$FbQ#5v^& zjNt&z!AAI@F1qNV8UzfX!HM_{z?n$-AcY3bMC!{ar<^k7JArG7KBpXz}BQ6 z{(FiTZ~~hm0{wKyKn??LMa*Hx_c;cphyf?CDI(BMcMRk(-~^V#j_-2}Oc4W4U{ge( fpY9mQVc`D(sl!*~m2wcy00000NkvXXu0mjf*ix(m literal 0 HcmV?d00001 diff --git a/figure/density-plot-1.png b/figure/density-plot-1.png new file mode 100644 index 0000000000000000000000000000000000000000..a3492b65cdb6b14ea384f75b9c6399af6bf5ecc7 GIT binary patch literal 22923 zcmYg&byQUE_w~#$z|cc?N~nNz3k(g?5)L8VjkI(KNGOPOiXtH02m*tofPe@HNOwv| z=X(b~-`{%In&lr}?>x^v_dNUTbM_uzYH27D;L_lNKp+AYW%f>PVd$!R&tDapAyxjxfzGqvUah4W^egeNfrbWjgQxRYDcF>Ce`d=dxv#U z$4T&oXOT1k8$pu3;JezU;F%d=mlPG@DVBAs<3EFG-v&hmny=nzDrg1pZDx{nwl#z> z-8~}uEmg~tjLSAP%5+z%_B32;J@(^Px`aTweN)w}^`C&h*}FOZe?HspVy=JwgiJ{! z15;#_b5rs!dPKdBmgG!`h{~hHdiO+|p7Uw38bnwV&!^64PVTF}>EzhVb#7UMB0V{i zJEKL-Y7t)aYw%v7UcWc^`FSFv^rj)hvBi1=AEJ}zdEJD3<)>rZ^GQDQ`$$^AIj=tn1j|L zg#U6aaz_t|to7-%Pn3L0>dI>O?rab;JUEbmqFy!`W$gDnyKvUYIapiS6w2uL_mjUo z|Kg{79QZ6(Y-b~2{cl7HH}+*g>f3Pg&hY5qcZxE{Ym3)a(~sa7)DAQ&D!~ulyeBl` zS$gXHaeYV?pHgI$QI%rrJoItijeJe)vych{tP<_Gc+Dy9&zoCk+FgT`0l5J$giAt$ zTM24#Fz=@EMz!FoBQc1r?jo?}C?@QR!wIqY$;!y6V@%)KP|~iU9v+C-rhjj5{rdL0 zx7RfMtZ}gRzEzY|7)3SB<*QPPpVs>&FN`forGvOcAkysAOeB-xLnizx$BFLzU>3ia zhw8^kbk(F*yXO;+oD@tRii}@_*_enz7EWQZB;!$UFV_E-;_A zTy~pQ8ETF15%(!)oO#u;mg zV(&+K*UK*YzwZ0XPlyzQQog0HRUAHjXyE#QP_Z->^W5P~pIReczVtBsFshcTg(L9p zRd>txz`nrmfrEkVPa%mC{n<;c3lH2FJxsfg9sbf7&iPaLQ8n0Zq&Sv7V!q=<&Oq_( ztMG!hIsX9%KFm^9dm=QYZ0mv|mnTP(p3stUweZ5Y^?^S4Y!BCw&>=0Z8kHf6s`==h=M1^Hq#rAGseZgU zVmTWKG59I&vPJNn+(&gJ$%3Xr_%|xU| z50tVl@KLSO``Z@}7GKqZQs_cLd+xYU4fSEqVv7Ean5zEWPM@e8fHyySvVGik&hHy^ zHM%j~=~24I9a1MEk=-`8>p3okxG4eI#s`D*as%`*jPZ6sAOuK7UPjj!ypxR`LZ)!s zJ*Vhx@J#;=3yj4Ap`~(v6&fQJ>r8{dZtZ2_8jG0@J-8?v(?Nh zn;M0O)28c6Us^BzG;^2#h2kWDa1szUIE}I(^)z0Y>)o4c9&1)6_CL+Lts$VP;F~LY za!7v)cT(`xmMKYpL{k)F5VHL9YDT3qX4K8)nsl|zcbej|vgGDJO{ziHXKR)2OJaUq z&%{`ruz4gCc91mvv*0NWUpO_@-Xoz)H zJ+t_J7uZEX$!opZ+A0(Oi=74x?eN{$njtBpy+?wAP7}k2cNzuIp#kAEzD6i=Au1tR z0SDhkKC07h{Hnwz zgU`z(mfSMuCv(B3-qY?%_S213hks@ycq6Lpez0u(cps&!tJ|x>k3 znI+1nby;dqW7RO@E#Qv9aRlm6y1P~M{uTcUgBuD=TqaOio8>VmHd60wr1M8dQ_{Lx zwpnY{V@Q(Yj$#~auyxiFN$S#aq;J_W=}=O(oO^RwNX=eDYEgtaz{Km)yd%8vv+#A^ zOeJf}`DO`tKm63lwZ&_qNsg6};UTiT_1xj72v6H(b;ucC51r+RpyT+!Vb4Ii0A|+n zqfNGzm6eaUrh!Mv_FG^4{tnD}tPUj3_{=?ZaBvuyJC<$_BrcQgsQRi!{s()Sa($@E zvS(^CLp+Jx8Om{9SzR^kVGRBv zWZ8|s)z8yb`{>alt`9%3%L<6>R2h#`aOqVQA`6IQlU)8xe@1Tqpi50jF)TBzOKEGcWq-L9I;i=Haima?#gx^$N*OnziA-LQqHT8?-k?uj|`i=%Z@H~+cF*R3zrtC zre<6oy|a6zx0m*O-7Tj0*WAsuk197_%6UEXpq|Jx?EwnB=2_G?5>@p%;P8i%k^gQ> zL&0L-$O5%0<6(-S+dv(qX~2)i+LnPpa3cFdu^H?(!@Z-h>YgDTLix94D?;r}OmYhH zQRnmFM5UN`6v+FBT;z@^F^V0*!+A2}#f3QzZ}Ls0&e!jF1dQoDNe|EyO>5h2y9vVc zA|$wdZ$F~iOLYY^g%^XgN>Vm#50}<*gI!{@s;s?srpC^Rm&T2KVB(zCGSS~#4=I8lAG3oO(_Ea#qX=~ z*DB2LejR)% zcRb@`DkEKd%CCRl80w1edbJyTBap&{j5Pa(NBHhUVXw~a#7bK`asS9~x8+1$;jO!| zQXB7oAeM=k!i0C?GE`GNuMWR9Y#{o+7kd4P;>y~^#BEkDHi5dn`pyQA?}gUh{8u+D zD4e0mV{I_G<#~}qz@%n{nI<-jE zA@7ZHV}{;~^kKTsu9-teH68|fjF(Vq{P{8O;fpV0ES8b0O!#AA7Ph2ElE%T2?5IkT z2;PaQ6Rf`KsCl@wDC4xiaj7)S2ERvU;?r0!OsH=31qJcG>s$Jb9Nm(vqB`aIN`iJ+ zcRegFjhW|PEUdE7u_%>-W4G62e$>om0D<}P<%|_!3+>^Q40k!Q++>+F);s7a_bS78 zu^DvJk})~4vR>|ol=IBDM^7SrH}+nw^BPUCy7kpz$WxZ${&sC^3GDl_#2jaT)eaVj zkzK+xt=OgVlih|2bO<^}^X^5=4FXk;wioLce2%bk5D|ycgi(~9hjGs`iJG;DU+#6p zSIVF56;%C=Bv*|37RN&v0V>sjHxNmk=2gva4BG@Mq;D;p3kWj*!9%`dI^Mom%yex{ zRALMigbN8WHs@64$>4WT!cu7HH)HMC@T@32;X64I(zK`ZCLvU_)~Qe9q~CD2hE*>x z-xmCkSXe~xNMU6d)j2i4RrJGy-a9tEpz-+mjqZ*pA`w~Wt4^4&?{3ppdHyv*+wq0h z6EpI53_^`+Du#9~^2pD!7=s;K4Y4|>8G^DzYac(cNp`C$st z7)MWVty+BdEyRGsuMn7m?*pUwVm4XKgKC1j!{(&!WIG5WZIFDy;KZlv-Ot@XRYDDg zmD*M|q8#2j!E7)WFqtt_@j5ZbO~MLzQjXq_N~mQ3b=^ME=2RzW(v&8&B)$ z*J+F*qDV$Q{Fh(ky3(>K=l=1zX&I&+2<&^e99a+*l|#D|+e5RRyjK<7JlJwp^b+2J zL=Q>}!j|m7M5;)J2hH33mgF+CZ~?f{JFq@MG{j3oD3#X(7lP4Yv*-SjtTJfu^IDE$ zF#sXG}yxaOnLG1qb(AR_k+DtHBOI3Pzb-~McPq^1U^pxE^4Tp{wgK!aAlMH8) z<+Y|w;i`;a!vY5h!ZYB6dj+x=+rKa=jm;M2FA5c6DQv#8oZ2&MM#F0{D9C<)4K?S5 zi9b9g-Xko5jj+;&8JPRxP>2si$OW5Y(VP>Xs0M<^5isZF(DbOrldl9Ak+YQEtv~*# zAcKeX2c#l|2)DS3b&hQ0^1P)2IBwal8&N1-oz10?{+re3?YQ86*+wsSGe1;?!^6r& zf=UXRw1|sj*R>Z5^!Ut^9g@DM=1nI6{@W`}Zqv)%-6Q_m#OFKz>j{t9i!SJR<}kBR zG>H0{tn+va)v|fG{gkz+zJ~%_+3lAZ+I?YfC-cQaJLo4W^C)Hitpc2 zNFM}>>QnW}K9&6(smd5YSn}wOG!k~_#qb4-s7b9q$rV8f%ejOo>Ul(K2zgBCEh9%s zLJtHC$dObmn9-^7jEzXfqE@q^uV3Na!8aGCHIMis6D@G6M)BmBkggoYA@g_~=Wt3Q zAFju+{;;P#^AEqCjUm>}uqYyPPz1Y3+g_ag@M`T}$SFEw&B9k5Rk-BTRVn0LGaDs*O%`4S>;C^_5rHka36hw^m;+ z=_>!=w+e38#avG552StlG(CYF-39+#(vi&M{fYl0^)${xTDUFBXS@^uGbW~n` z0AL5-?{B=LgPd*+cvV`-90tW%XOr{G%k|E2iHY2VQ0TN` zh!nzQi+DWhe?1w~3ebc-ewBv#e{pGy4|y+VHc^OsVajaw8|N6i0g`#@kN%wd3(FQeJ0Msw94%PyDyuP5{*b<+iX&yh6B z@mr-1{dJQk6K*w8Y5q|%4eSPBRN7Jsw^1L^l)Q1^u)#OS_H?@3SPR+cguPUgJR(WW z$e?M*FSgpnYiOW$md~}=tj!Id(P%pX^77 z=$*$;brI@1X^T!+Lv7*_Q)o0)*!cNN!7o9z%clYEU!y_3w@1ACZHOfVcBvtCX7v|O zB+(SdpfC6PATQYQ=zi6BSplfeO}f$-o$JaX;2uiZ52q?Vl+|B8z}9)+>zU-}QmCo@ zW1fQ?NPr6#6;*?;;4roDLd;&`oPh}RS@1q;JWikWj6D8D`rPJQFwYeb7go}%4AvJQaKLr&xqL^hM;CA^sV0q);LZ!$)oC%F# za`JmTG5?*S##rE<=JW766%MdJtn}0U?1E~pf0#lE%s`&Wa8%+_vuJUId}tDvBR8x1 zhg<0^2=Lo`2V5dcP)H5xlRdc4*S7NZ1ti!&&g>rykv-ATNz$!Ku~sR&ey%B;`{90K zWT%~_hQ_Pi-!HE_IcKg6j5q?g=pm{~lj#n&G%jP`V^D#y34|Xr$+JSG`Q7Ji^3p%{ z%C~tW5@rGEZ4Kn25MdT6tK5<4>%wv@yjWIr6mU0&Ibr z-H(MCAYsH`ek=YPWGQ|F4%S%tVaxEisTu`BA`)rcXxk$)+p(2V35ANvZ7*c{_b~4hVGSL=!SP|$AR%WT=9)`t#QxXOM`zL?? zap@ zz>Pjf*p5jW5HHBxLPS)%9rzK4Uy1^HWoix}1NI-7Ve%_7^+K%JDp5w=o8QC!F|If4 zz(qk>n|TAs{g7vy0U?!CEdQOjV<4&ea$0Rv<#3>psCkS|uMpY&Jy-OCdji0M(kBmo zQzVB>5=iml#ZlxHRr)5PBNU(rKOd>TG2D?6$ah=+h;fdp24q;XfjM5 zXF|>liJ$m$p0uFX6P>B^DnRUVwuSYlh#@i*Ob;T~VcIsb=>6v~Bf!h)Wd}NXJOo7cQC_`bVf)&KH`RGQkygqq)JYhq!}TP}FEQ;_2htw{U0Iph zq|0kKYyk2zIgk4s7=VHe=umO$$d$YM^~s= zJU^RLzce8$E35IZjQFLf4LW<>)Ii_U_~brnWZw@wKjixg4MI5ot{2*w0d`|O-RQo0 z3NXb){W3$gE^vpVf`YVm@G#HuT*s;ULr?-hZqdS)eaK6&)> zqxz902~f)BEI%XNvqPcAicHtdVvjJqC(Ny4dOCU#5-44;UeLNq5vv z0$+Kve8{kE`b}Q`1z0%#eTVler!}O0DYH}6;66w(9K9xDgSRs=UAR)d$;KA!q~vVLF)?j%zO|b zFDZnzSF@i7jEnD;+ARt(umVA8*RV{cXVKE(7Ms8>)i`y74KSAvMBZ)-DB`^s*`4=4 zMkC7p7?y8=CYw_|XaZtH1&nS_ATxkCl-{@{Zy6vq6jb_Qkvv)&uaI;QRU+S%TbOtN z?0q>*!;G9?9%X{v`&bn|j5n zd8`h=^H@ZBn?DZHs|E!;-iI(e=B9=kq`>X+DTvmZ$+28oIt|XG!#V`R*w~oLnoh8@ z04;`dX;;VO+PjY*yK4m|FIY^}>iQZC&kWNDE2I!YD=YU}3f<=N%QjqWzZk!Xe|<6j za>%dxNLlF|NYqSAJ-sngQmpMHzAqdZ4A$&PV=bWn7ERPiPhh<^LJQ$Uqht~M)BzG zs$mO{M=e8fDaJXev${3)KEJX;z~!RPqDM|llUO!KsbA2pCa+xZqCae|c|Q{CoVve@ z6$6XJjT$_g{9p^N&k+yh{ZeE)kjaHGe2|^@6PBxxi0+JRfb_ObkyagI1*&5yG~S0E zzvg}kzX5H^zn&^V5y&9>i0}Evil`h+EDpQIyX)gr|Fv|JK;GV!I8QIbR{|NCqeQ?I zOwsS|mH zwsPpe65PXPYqgAyo{-j1!7^_rs*`fSKK;)#xCz+cX_&{Ue?4|C z&)!v?iu@b>5TTEFrOxCBDjXCJ%@(MgM3=Jkfd3+47k!2h5GkM@FE0&R|J-CYvkBc$ zVp0aF82CxP^i4-Uq%w*KeFTSNo#ZT_x*+F+%QMc;1+@A6|N5*`G^@P1NaK!^2MM%# zUggf4CF_7N|IYKF7e7q~uA2I%?0iZMNylR33gNv|nx8-G8TYRfn-4?-$F5vGlQv`m z_l`WsV47G;J3~$Bwg$7HUce>Co|cxrhdOv2Q-1tvDS;0ifR&U^dkzHuhfL}ycf8@l zji1)Q(7rXYyZbBxl$FOE0yo3S2TgT`-JmG|f&%U~3i0fr_fft$-mwPCbM-G1|8KxR zJq;5!DL+Jl5)D)oJb0MPWD^WV8D4scvJbh-c+gGOkDtiK41 z55Ao(FtP!v#_3h70tT?CFp~azp_sX_Pn_t7)9%8GUJA9rufcQNLAfC8N&Z#ygQy&BD1Rh zS$G#EQGdPj@1|0l!6&veO zv4`RpWU4d#<_$11Rm%`md4BSTd#2eZm02#bH#hj&-r-~=)f88e0H9q7=9NX#~?Y?*+ulXq@vZ(d(xV6^7&TT3R*5&mQ@nRP7aUA_bGr1J z5_Sg2F(&0-4jz%NO<2B8rgi78 z+74z(j@_c-3%%$K0b?u!G7B$y;&~5`@`2#lVy7IDBr95uyk|LDq|HuvqmmP1Ns$+D zl#AkN1!l+wzuuH20HP?{*%Ti(i849Z2(*!!MT|J5wSbKaR zqh5t>P~vkVun2CL!B2=MeC7hc^ty1p{~17E}p_{?+k<>3Y!QV_Kk3_s@AstMNY%gZ})wd9vltqk_39nAt_@# z%I(@ln?j|Qy=shu>KDm@gu!yr?F|$a0O~zu!lHh|0oMbWOEOLaqt1B=$OIOKOyDow zDp&FGVD79!5boYQ7^m?C3O)E{C!Y5Fsou?jB_~VuAPj-5cSP0W1XPpcZcVXSyBkZXNz2sfZp11kW3Jgks4eAS!Si$oGH{#P@HD_O zg@*3q5>BN?`yxga1s(4P+uRN<`te8M9YWy(&`qH=vr-Zr3}Uy-X^$alqUY1PqLCiJ z?92BRhzpT+ZBk~lIzIEgjZ5?PGds0raeZJ~~pl}Y4J2;sniav>acaOy$!mZe;iQ>dh??ng->csQu@S|Ifr zD{V?23qrIp;oSlYZ^^lW=G(%V`TiU-5O7v}&H;^s!$Nx0Vzt1XlB>pZfu&4n>cXnZ zR0_%vcFr(!m;;*hfklnu#yeS}3nAWvVr#5M4}K;AiU#C!9{rSex97pDhydq1h)yqGzuAp{lQ!< zD7SiwCPGZi=-UDUc(*Y*Y`N}|(@ZtEGCq$9P6}y&O+fCBLO#fu>07^OnEE<}CPW~d zPIchnszoB?T+$$}*j(deE@hAl@i`Tg0JL%i6JSu@VE_pwoMapjaJEqR$F(fj;a&_E zie6Yxz8CkQ51>A%}E2R+L@l~+JorXCV;kOb1nXI+=+)s*3_S%v zV!3b;TY0N9FtCooWH=22V&gf$*8@n^+eQAH#f5!(~5;MkF-rjUw?+F2ma1fT82fvHs3> zn~(5)z}%7HxjFhRE=*o+0569fs7`^$leW(;_j@t}+qQRjhYd^dvx`6ng53JE($G;> zSA}yu8XM$*G=Z(EGO8K+EG7oZ!y%U`>}<8tm&~pT%hfB~% zq?)1ztbd>3gcd>8lQg9gu^D_(YA<|L0u2a&DGtacRDlF>&IV*KYXUvgtT>_>`!XPf z%YfyVu}8eUXyhMWoOOVnBoIkoY7t7E&ToU>D3rzB5Apu(ECP(7NHDc1X5a=tM;)Ht zL_zL4Hf0j9h1x;qxBddgj}Cl$i(bM7zym0<%hPiB#&1pf`VR1vLVz<1z~!Jb>q(%D zLdx0zc>GL>48bCLsDg+>@e%|jH9S(rGnTSs`%r(G1hj5oNKmfRB8X+qnph9i4i?Wk zUokOjZ6>NNm+3wedjqtUc9rdD;7JIC=0Avt4n$Q_z`GTy?pgp0+Um&>s)vggi)2p_ zEGoXXIXdmeQshFBRQS86FX`V0|9-I6{ykU7F~xxH{wWK&W=ggAD}Bq0h(1q047Ch2#0~%FYwrGbWr9nuszgrjik`o}3Ek{J zFTIjj23~Oj(3a`G5?iH82bj<|#>n`;gA$qmUtey=W)vJR%l4l776&qikbMPwWoE{v zs5M|-KABIqU6AtLF*c&47u>Ar2YVo&xdErNsr%C1dKK`HX?)0Rw&!!ypXw0A43PW#G)O*o3ky=8_J#( z2zt+9nU5G%Y>FR3d|gDw6h27|*%E{CLk$vbL#auutqq_N>TKW>()Hm=<;Yrl@)O6B z{N5q1+H=4ylq7jFt)ua9U;YjPm?W5Yexl|3R~w&0v-&5><%pS&&+?5hH{jsqFz!)mLnAq8D~UrGECr53bH;Da&IGjwz020TMe9RSZr z6W-pNcqd*D`zr<*D*_mY;ErUueZ}d{AbI$kW<1HH<@q-tqglC)>dE5j1AubaZK$KI zZGO(yy{SiG<_uC`xN)7%mJKtFA+@QU1uyZe7kTW`E9t{43jYe^(2bdcJW%pHh6r! zf~2HA2(0_rtaaICW|nW(&QbjL2C{3Zx3|^xYF2~kTjX_ZGXv2_`PSDY>b9wwuY8zR!`;8FJzmOIS* z4uc0;_Pi1z&DJ_WyM-Ti>>DV;Qf(VO*qO%qMIjvXgHX^r*w zE1&Oc8lgU9v6+3?7z`|qq0e1v1Y4le1a!59?$A>i-=%nwxM%g85X@)|g^+LeCk?Dl zN5CKWbBVMkOofjLT>l>amgppstv*qBXpGw@{K!SGzzkn>1VOc|T^#gzOqlB>$5xU1 zPwvBfeRc&(_qoO7!D}^boz3qQSN1m+Wg^v4ck;4NnKOHRh=YoM=ejqrpgKtPv0xB< zA53hdBoHo3F3AbpP2h(I!l}#|*T)TQ)FqgauDrV8Za}8d&QT*>Q;Hdxk4NK#%vOaW zNV!<36xwAwE5AhG?WbTk5^?GIV!fB1U&7eJi!7}CLJ;6@_AU$#7}6`^=1Aprs=*dS zSj5N{^!OzbKm6~sZEoPe9%ZWl(H=q>y&d`4P-l^rMne^V(=gzjM5qeZ2pXH46i6ebbevsQ<>QqEC|pP$5Uy9wQ0JcS{O7h z=W8|Sqe`CCeqr{W-*&>T1rJeE@xhrO`#z0p4cY9oJ}_}nDa~A|nNvyg^z8YNMCOM> zmc$e9j-q6sgR{F|{%TBn%|3ny{dIrP@~)C^w`JUel*rC>$@fVime5i|QJG&AS7imU z^qfi5fL73$1hUZ>J3$eVOWUci?uR=rEG~Wd$vpcktcz%qWGZt09Y^72%tOA>7a``2 z!*FL&Zb28fD`w<^(F+6jqa#y0sng%U+-zUWw`{j@=~cH89YL4@;I=YWL0kn}@CA(z zWH0q?);tEN!yqy2>Cc|eS?Q#(taq;#eZ}a_?AqU8yu|8(!UtxKJzL~CJ3Ltvgdq;X z@!ufJC@SA6U>=Xw6GEUC?3MPzJge(*$`C2$5avqpwXGuVPnm4QLWDg3`cq&8ftj$@ z6i4nI0qWG`g7hYSO~bjRb|JrEfG43$37LBcFM9bDWFW2@_MZ#T6s=Q?$Hb87J+2kG{6N6A)e==IBZB=vJfIVeoF}Y<@Lx9`oz&N$) z>nk`4js*JojFrrkFFqfP|5*Z9=W&7i%h$2=S$3eM6PQ|hBwK19TAp$x+^Qr@?7&q% za|_lvVy&&$YG)Gf1xV`gdo4B3*=b@&brIq zWAbLAg!S@{lo+OeQY<(aa9PyIlEq7u#|6( zpWQ5XroqCk#{_KIy><=7dXOt|3LR^tX!j17>cmC3mh|c~pv-*h?~lD` zl@$+Y{XbFC!vk#*P;p14Y5=|rGD=F0LqP7aqqWpy)Tcxz1cSzMl9S})Y$>BDoLBZcLx88n>PlQ`Gi}`a<2e>V zQrnEwg`^spX?EX7ZyQykK(Fv>K4P3!xK!*QQz(=had08{ZLXwEH&_Y$wtd#}2fFu4 z5Cr4qeIVw?&?G5RaWgC*L?*3%FR&R`6Ts+Aa<(n_HyVn7M2pCHlPG(_)a>;D^-9?# z2N_Fg%^^{rf!<31T4=}4?qz31ptdnn5aQ|K&+k-ktkc7_PaPAeZ_oJ<4||S@t@iIA=)TpN#DIGjzHny4PRnf=;_52T-2?s zw|f33`p&YTn0Z}~C>^&z;C>vit^u={IDM=41_`!9C)dB&j-ZA6nN9x*fY>@v(;S&*aHe%o3DP9nevHJ&rH=q<5`0S%&kw1M%q>x3As_s%;*?hzx} zZmydfiYt$=*|na&U(H$!U@wZI-bg0lLQzrKMcO zt7Dt@za%FotE?=%0IYg0F6A^dqZSKGT6%hH>p!h_W0AxL6F+8(^G^IHFVEimd}rAr zwulXx$|(09H+n~=fVzsQKX1uKJt@1knjF?Rv{g`eVek5-fna3e@Av=gxHM7$jrMau z;jcH4NL|-m0Dhv=k;QrUZbCeuS$a#(O7CH2u1E3BR`XLcn$p0Nf++ouG0^^p@zy?V zux?(WNa{v#LH9he)-kwCzSFpPShC&0sCIhW^Mmy8xM!WfWBSK|3@6+SX@Ev6Y9kk% zpIx8VMO5u8Ck*TKP3LWufQjA@VG%0!*`4LL8KCkw1DJS13b&;{Qohd>7!{2GZNs~2{ep7?Rx$Sr{Uj22)I*QT9!aG&IQZH)uYC-tvS55_(hi2o+e)w~Re z7eT-{Sh5zmqxMypOd$LHXTGXt0q8NS0Kx3LM;nF+p;!vBx70QNy+**1m*<&iBUa9A zCbj`BlY@+hot4|ua{eoSaN5wfuslHaV67jaj)n5epj7Y@2pz{hB}xA?-vP7Y5Aq;C zJg{S?7qpTvK}QBN!D%{AqT`8J@}2*s3I?!-Ai{?D56k?-qhaqFUi|tvBZ~>q(_15g z3EB_y!o-3;n{ECB5!@ltyt2!fcxqC^!sRyY4S`cpywX2Cqm25-%v-Ce)b3z2@OO=P)1Nd_ty`us%M7I8rC8;&i-JFH9j@+^Zk`w;H>?Qj(?x_%-xGj zCZF~|oWb*^#Be76#ebB87JaLtbXAHu?6TUX%vjv-)5E!n*N2WjraRYveo6Oq9Ytzn zABG1!ENZJIvcT6Lwjt7l5Lsg$$)!X||+ef3@$e zA4(Z;odG`6Np9dbfshitiHD*d_0zRWgl{mmRMWJI!j@A^Bv-~af|c#pXH7LU!+X4` zHepZu$MQN3p_*t650L=s$>$nVWXWA9gZ9?Y{5ZGunvMGVP7>;gPX*iKrRkTp`n@!g z;sN@u(M+W*ALVjZRxI>%<>S8M{d>cI(7wx@Rx@=Arb-PF`ap+20o0Raj`&py<(pM4 z52b9&sR?86h0^2bXG3MbzcnRQ?z(lL0l$ktn@U|4xp?kDh6auDG={!^uubabTd4`} z{$1$B>8D$!c#lAQnf(m&@zvFbw-g-%|4LUDdPRJt)i+pxD|kTA=Fda1!v1>qWN(#U zJ8@&M@%L%r;>aFyUaOSMJ3<1soyl>D|8C@$W7^pQjcv*~0d@qHx_ZNG!+y`Pxlg_f zb}Gqf>78~QSpofH4V^>aE2{p${n={r5&6Q*U1KEtm*l8+oeErFd8c$SF`dtbtD=Vv=a6i9u z6b~1dr-x8dIc+Nw^YLagPtJ{3MpoU^)aOetN`~7;DeXE>q`xUhc&2W;5g!^O;451_ zvW!mq*yO@>DT2mLx(7{Ozgff>SnM2^x~mB0oD)NpG!V;rMFPHoR~w0|@@4FungA7> zQvfK@V|azN?nc?uW?|+9GV+H!Z$dDji>ZGJ$ugWM#`Mbu#)y`sZZhxQJZ)u}+)kUt zfU<9;t`Ce5@PE0|#~G`yC?}gK`0(yat(m#g26NcrRUKYQq7>7o>k~gGYhDVRgy{;= zNDzpe8U23&0O^P5hWS1QIgjOqPM0{M9yY4Fp7jt(E6ttUv$Qj)OoHZQ37BeZPRPWJ zv833EwmDMz`5x01`+#se3@{zM94B~N^O-Vqq3-O+fC;KhG;skyG6~X)mtT&xDHP%Un8xU2-D0M zA2Y0b_{f!~Ty%XN44E!tJs?Xp}VzyWW2vocTwT&SsquGPE&zBGXi4Ug*os z&pZ|-ch*HtD?8|NW~(7VZv*_4NIa08OwjpAiIMf~+n4__5joUf2rSKIva0Y)iycnX z*Ab7Q4|QZ}fx*vb1I~Ktv+3_O%q_k+A`D%rx99wPKkUGJRw#%QPB-o)nfD3DRP?dw zf89+24Wt#0aE_tr?cHfeaxA>gp?NeQWC@-9F&Q**NcP)Ml!Su2TqKh&yCCgnc{bT0 zBhzDJW7o^?ENSMIxK)5pcPB$A^XVev{|@pz9vGTs>+;h;zi1ANO9ua%d9O&EyQg{2 ze003Rg!+g-D~-F7B6z8ghz=d_Q`(-O_6EP^(1br zW;U$q;$%#+ssDSMX6k)tyP#jMEzQ<2B*?$!8n&XVFL1${yk$5+6rwQrVl(O?; zhuV0Zgtaxl_8Y#?IubE`sbITaRP@g?&aZ$HfxtjfE;AbcuZZ&wr~3W>xIIH8vR4@y zk?a{+3EA_YY&rIbI425~?3Hn{v&XR=BYUrmV`cR*Lb9^O?>=;WfBrZw=i(ga+~>aE z_x*Z39#4}33Ic!STMb$Zb95C;X{bpXY{nP_iPw)?x{I5Ju4rF`q}l-3)Jz~_B@h!< zggu*0imq`dkwW3j39OLs$REqJ8<9h|u3otV*TO>VolP$a(*-Zo$W|k~pGq@~BCJc8 zi>qLNHF7`4-DGv?OOuMD>$cuYA~vZd@0qzLfuR&pOeBYL^4Z+Tn_Y!naCQ;(^+T!h zz8_adH+eszBTcN$sq)|juNrV(bD?TZpsYDN4>cQs)@Rem7M1R|uy;qKjF#std5g(7y=r{I3;sr4`ObC))1y^OpeoZL*D*$-w~f)3Lc1L zZEp$EhA1qA2Zz~&4}Ni4Qb?32b$z*WO0{^QA}Pn~v+q}$k` z$_7LR+DlV4-%m1z3LiTYk|`SQ`9}66?A})$5-sT|Cm*De@;+ji?!!%=*^)4W@PvXb zuJ`zIZ2zd?kCfXvSU3Zu_0C+pUIc@m>h^>Q67jpV7J7L`tzVNXZsyJyh%?=t^ObT; z(E2l@>$uvV@SRzPy^~k^=qZo%@tDS|V5N^Yms>NHbaoa3$`|`bv;UMw>zt2}>c||^ z=8q4U)BT%bD)K%qf9&j|5_vDCyzC^_H(YWo!fZ$C9dzNMV=Bdd=ArywO{+{gDbV;t zx|FtOl%j@4fn~0(R%9uUPZ6ixNR7;S-e5w;+xVwUtT%_G?~GgcH#u|~PNr^1)H&tr z!M8vA;y-_Wf&47>_~E=vzvI}pQg6?MAqONCX)N6HO^NN_e-gvBV2X|RKuQD(99@~E zKWcKl^|%Ld2RkL}Dbf}z*Es@-3n#ks)`XYGWlkSgv{K6~WnGh6UdMSCDUi$y-x|?K z7JzV0+KZXI7UZGd5zD7*sEM0}g(dAeA5gs6t1N*~^Einv(FjnjA8lCYj~JR|Xj@T2 zYK7nnmOS>ekBd4rZfp4f&3z?;S5EnIF-yC#U4^AEk0=<;qy72Ar|s*xrukG>X|vqB z76{%P1Z7roZ1qJNZAUwmW9+%1o+(nwKXXfF?La1bDBERu=LLCu1uU7*@L#K_KXP9tH{->UGDjA{>ka`MHyvRI&jqag)8vV)jk14b`&x2ztWS0 z?UKDcH>X8zY!aKU=${DkqS=sdxM2J>n`62m(Bk?4Ve8ap_o_Hc)M=M`LQj@VeQx_! zm-AAmsdG6+Thr=lyj9%#i zjOAacf7PeyH6BJ7O^d4k^()ZV>yg~w`bMVXq4#|;H7}k^g2(R zkb5-~i><~o(6=Rkq0}B*PFU}JoTaStwX-bmjJMF>>T$1JsjqtAsDzUK-=kXqBxiL4 zY-TzL889pGJlT!8qi!83a}d==F2%Q#ERnh!;|#H(-N!|pO*rg!7fCqMt=}fh%i2#UGyj=G5@eeNHPxR%q5_W_{k!9qu&w)G<*YA=6_> z%uUeqr0#_6*`4(_T>kv8)=)Ys1@_F71;Letb`KxD%V=#21)_e0)WbdVyAqDRMNdOl zr^ykwef}fJM@Xr-Px}bYI9px_ark*MM5OvxF+Nk&H@y{)?jdS#*{1h@9hU(c#<8j0 zVO^PPH(-Mol#LjMo0P8KFyMLdB9H{hDfw}-Jf4Wp^rSa9zR=fWv!1F4nrk8 z_j&E!JM%xGiW&5b=SPaYU;momRopZzQkAwtJYvNqC=iBg_2>sCf z;;kQ3fs`3Wh^N33m;H};Nf64j3<})ntyRam?N-tOA@>w9m(LC0_{{mAafNM9o_1Eh zTjoVv8a(25wLqz?Ar*hznvMb%#p8Od=P-LoNG!B^eAe?qVB~op)2d1O0^7gqXDAm|tBt*tXp2%q(ONKY8vc%M#+|4VFPA<@3&x=Qw)|Q*~ zdg7ux{feKK7LSJdzBh@Ye{o{wwbdDSs!T7ZH8lR*!!TfiUbV?2dD}VEb-$`~{E(Ta zc(z8Bg-0{HAzF0qLG*6<5Lb$w>#JYV2!9_`=VKo9F=9O!6WQh)O z9;9B-;(XoQ-ft=_Om}GRMa!QMi`!zGH2PINC{A&RFSYo2lv$&g0~0FEJ5|DU(_^rC zC@G~;FgLn`A6tMr6`6ms|30UpxJYsA@sw|hhJD-jL=F9-)Gx347lSptSL5UR(j*(I z%(m9N#zo4txAq>(Sd^5S1UTk>7u-s@e+~B0Nm4g!XV+5Vwd#n)O6P8tZxYtV8>Z>M zB2Vi3Wob8=XxR2hE%r*#R^re9OHJT%T)dCwU{rj`#vC%~_$McbxXoj7J?p7lijpIu zcT`;D=vrY_zACL-k+z$Dv6nLhrAbm*@?xA`r*Psy)sC2Bp3O~1<754sBdU>oqztM# zr_Ck^gU%yUVO2bly(l zNZ{Lv=&UEN}2 zXw5?YYD6@tF4ap#Kj=ZHt*Non;6xC!o)wTaePavyHC${eR`(mV(VvD1?2zZ4n)olk z!gmNarqp7+LpE(({|g$+mH&w~X~@eATR~Vp`JVhOu{W%)zV_p0*;+wf?YMYQXwf-8 zjQkdXMU`WGCA6?fMwVg&T?}dw$w#R8TGPaJL1&KrU|v;AIs(C7FW1`PiL$BMN~$ur z_#Fs$6CN7Pzf%1WE^-1%y-oI!?h2y7)T%!=TiDVQ7Px+Q>!MSq2x%LXJ#`{#%xSsaDDQ=-Q^(|^IM7F9PIAURc7+PC(P#nM?}f${wu8kbMnu@WAq6s zDMd8;?Z*H#x~|DyyY$|$1|4#}a-Np8tu@9rYXDVzGZQ$7N_@$o-?Y7j$~LTAyOvMy`CrYdp(st@vM&R%i-_$&vKJ_0n~4eYUFs~qB} z)4m1(Xg3CW31#-PM}4rC;JLkx@vA@52``~dT{NEY^h-kUh7;{R`+j-yNstYqRGvzG zqj{He@)zpW;cp&+_fxU*Cq%A9d;vwqRF_p;%)wioXeoSTAqLfM7LI90Td>1eaKh-` zb#Y{})O&Bj_`T;y?F}x5+8yxTh#??;7$=b+Tz3a_QZa}mXF7@48$;l*IN#ewDg?$#b z(@cK=Gr;ARezgr^kOuN9Pt-kKTQHDR}6B|fYma_Q1`xvvp&n2VFDjya63Jj0x`=shy&x< zW^kA+%h*`&BI!_Gz%(;qa+oYCY`H?Do5zbH0I$=*S@_zX5yA|hRQ`vg*jmKt9bvSD zj5agXP+s*!fqiv3Yp7INQ@8MtG>6fQ1BhC>?=)FfWV5&&V&S)z(;0I(AU6NWcDrJI zd5)S6Q_^e(B~2Q8BDq3rYUyoa@QAE@tf5L(*8#P}9_D_WJ?}sI5d-W190c#29mAYl zW0D-Z=63QW(0KYf%;4|D&RNk(sa8d!(r5P|V7hoT_qc;^BVY-3+LF_%VaX&6V`ZNV z5}jQ7z1C^tyA0l(39{y)g0CmB#Z!-(b!pO})l_8Y+82^iJ{;1XQG-r1x% z_6Gz90QLC=@DN8J7>TzftPdaZ_{bIG>wrL~-m~MqW#C4rinDOu1`9%<6UIu{XdI~K ze6j$*q#M72@v!&6XEMn!FRJA+8qhO(WB1Md*6tYBx;`5d-!X5m#+TZxjLjbfX;4fq z?wyxxZxL_`Ga;}u zj$m}}W(yi%1}3VFs7;6MXN(TS=+Q4^n-UShBtuX)0&kdQ0lCMx-1@x) zKChOW8}dogeI)>aGu54+Y=j4|0PHFTq~mpi75<}qU|tensrysJhT?*Y_!b_n^78EH3vu9R3vVB&C0NEFor+cm$EsVK(^bNWQK60RU0luLLissR-kXl*WwKS! z8T2FM=no1hnbcDF`LP`bx{RzJj5TfQtZkPO4-=ZxC`i!L2cW{v-U#^RG^+*~5xoae ztm-oq8-w7bI-7iTo|Oz#$^%@1EzLyOy(W86yI>{zogW>Xq+YqnxTb6+I^)KB(f;3(m?}VtmSCZ`LbM5*I{J>{);JT5T+WThY z%c*i3_u+?xY8V=xG-bDw0n7o4x*s2S@hVyb`S4%d@DCDGn$@! z+g%mGY$Cb|aH2Oir-uEig*JdUDJkB};|mgA+!q1~lamZSJ^$?~PuGGCQ`vb7BsG=gh2 literal 0 HcmV?d00001 diff --git a/figure/exploratory-heatmap-2-1.png b/figure/exploratory-heatmap-2-1.png new file mode 100644 index 0000000000000000000000000000000000000000..6c44dd1a656dd8655cbbebf41a0492f7d198993f GIT binary patch literal 126072 zcmYJaWmKC_)b5?2!5xA-1P^X4PJrT2++7M3D;gj`a4Rmw9ZHL~XdytMI2340fnvog z1Su|W{?9q*eb%ge$=vr|lbJm`vwzo~Bm+GS5<&(-002OurKx5F001xl`+)E;-_$}L zq5uFKdJk1qgJ-H5s?WWi`y2Z@IJs!J__+9cI2mcE001&+X(l%AOeWL{y#a2VTr^qb4gcHic7(g%-DsYdNe!ee zdMiXHwiS?&k@#5PB{HF$4zI$@kXgW{K?h4ppD6g9)w=3XxZU?JKkm!kbi$Zv*!)-> z+jr}U5;K74I+MvDVAypgtKy+0%emuzr>N3Lk+7~!r*Fdn3Xy3ZuG@lJ3yHWU#-54! zVmoVZR(Uz}v<+`4~8uxc99bdipQQqo%OEXHTUe%mJaVwU5s>8G z+cf;B-I!Gdy9I4M^Nx7Ww6T)+Dw9GH@8#-RQFATDk(zbW3VvgMS@oND$?TDPNU2l7 zEuq9&>AeZ#!KK%Y*(1e+K_5FMEl*D5aZpJ;Rz=63{BOLCOHcN84ke2w!$Z|>ulqtZ z&m;WHWd9t%_fL^|Lio2e`8lyPA7c}vDqzazdt3Ldt53+W={*^Cw4#Kwsz|Lwe%g4~ z?$2ry(@8I~YSaF@jy5WPP-{=|M>JdDy*5lu(_a?qKKylMI5I;AFN4QPy^fCRCut|Z zg%m&&J`n1fV^cUol@6qLqNOW*Vdduw}#QAlsjlVSslD?1UoaDQ~wfS&6YS(Ig zASiA4ev7`O$86@x?&{aA*tx*1XZ1$Zw%Ku~|0m0--Zq`#-(Kbc!6ppect2^;p9@El zMwn0ijMZ-etKP0tj=pc5& zOAK;SFu%kRb-%DFFzKch>rD6F#N_FLZw?=`>2qC~er2ju zE>-8$9I2-lU1p~h9J%^CMjqtG^2SsRt&7TJ3GKB(<=$8~x6EG}eA{mc1<*agWx{GR|I!^^#Fp32PMJ5Tp7P zkAL3?Z>4tt@|Y0OpE!NMvlIAhxH6~6<+hbU=1fgE(Z=G%@AH9cv5?5S#e>z40Zn^C zhVVPhNE1< zf>v1xM_t)#N$I!a@LJl3Z(zw?$=3X(s#d)4;rQ0z*bfyP{FO$6^#~9|eLM2BtT$j_ zcsO8EIHkpRw==jb=&)E;Jit|tZxW)3OpN~DRzryWBCm1zvOyw6cg}XE3(7)Y#&hNgMX7A)#37bt3Q1OhRH-~GQnsh8X+%v;xr&H<0x1%0z z%HR)o7xT5cuheZvuf+62FaLC|Zy$eR%`c0*)-ttu@g;{``|E4F$Sd{IutUC0!h26j z+o-!2p_sRJ+GiEA@kkiPIy@G(ck=DtyRF435^D=ty|=Fa{V9!+x4G{mS;BKslR?_Y zx`R2Mgpo(_79AV7Kj%kAzX zZZvuaWnHf2zqOw#5FPzxYVLsRvE1Qd()K23sk#4pA_Q{hH<80P-yd<~Sx3Vu>`H83 z+04J{_ldcXIcTY(@pB`!ExmWGO>a;anb^-y8M7nY&uPxjZE@kdF+{>+5jR^YB4RW& zH0Nm{jL-Z0V7MlKZ??KVlREBvcYWnet8s2Jb9y?f*m?4|b3M`|-AYj9;r=ebJK)c| zZ>~Q5SG%j?e`xOSufGk@YQNpTVD#x4obSU-w^&UF_Eq!1{k8gn%_4BfANc~3cQmW` zP%L{;oh_Dc6J*g^e%$M5KPl8Hab2Jg5oQR7Uv4H8PD`&ps>14YUvfr3%Z&D+sZqxx zbXN4NVfz*6e;gkMc8i|*NBE%-z!Zt*g3HPgrFZcg_^j|2`02FtL;@`dF|*v$>wmvG zl~}~xW^%Ii6)fa}enXD>Z{F4e@*i#n9_ z|Ka}88=U_mc>PDV7!2S8iZwdfQ+%)p?C@Ma>fU2Xi8}glk(>^~S$O4nDeHn~e7_cc zBs5(lnXM3UJ%+KvN))x?-MKS=@d`8RGyQjWM^X9Mq%ulH($9Uom6!qV;7h^GB#iyO zFLa6|#&v!su!Pmuzf-mB^bx{m2`T?@5!gLHc(LR}_Dp~IA7+3v5xbGYB+?Su?0?Wx z;@$7g0z_};%MOGhQ;BkKWErm-UDXnoEFp`VHln zncb?}nccM~(lR}|kpUE#Q#>f~mRY{4$QBkC0h2E-3^cWi_?MbKD%M!zdibrmj9%8i znmglryU&jv`hUu=2n-{7X>B7uG<+XG4-uW)#?{6JP<YQYM&NsM4Nwj>h2aOtTHTe}<|VGt!(>*EPHI{B$NqfRJ>xI_L*G{tx?4!U{QCKNf8JZm z)Fzv!`-QUBxlBZ_X+byXd}e!fU?1V0{mN8+=()`v^8N}w28gmvL7*4cD@73DIV5W z`3F`K)9{w-9% z|8T-A2S`3HnG3W{d|L)T8P}=PD-rbe{Wr?^+?nvR1|`0O!vA>Xs#+Ar;R>~*&soHr z6J0j%G-xJ>;T92RUq!AGdttd!KaF{deO2NWfNLnSK*Rce+sJqCT`H5-S`Fo!-J?zL zz;aX5${IcN%$bD-ndLoGenROz95IM0aR$pA;;N}p+=FY*vwwg8> z-k+A!8Tl1?^nqWX+V3)Pxvb73$kUC|6YCdWOwMiI@ShX6eb4k{o}bD1SfEXR0yvCu$0eUW??ITQn573<8#=5dryB_g9S!%XcB_n4vSnuojCEj7$a~OWn-wr0RVRqm z)3)h|E(ZcLut?$PL3~qEY8EAplk?Pae9TfCplK_g&HHnY9<+fY$Um`9@oFUqZ!BB_ zk8Mz~R=#h)D(Qc9w|>JwpCjt7ER;kW|K3Qt?njjub_c~1IS(bOkwfib$i|?`pLqUP z`?*TH*sjb@{iBbntHf{OlW4(@!#Zs4SRDY&-@ku}*gH2EQYhvQBC+ zpl`dBUkhkfoBo64i#bWWN+J^HcXehS2L%uT;DZJ6Nht>`2D_eS4|O^t&v83(`z3WU zx+tL0o30QX1>+Sn;LHA?$N*XUu$o5zp7)MQ;6qY~ar?7GE1i(iQbL=Xzuq5X3Aw`S zhJL)gB|7~VR8>;=6R=;8I&o&fL4{lM7mSJCImqE@^bIpvgFutXN$5Jx+N(-;>Gr)j z;wF_Ri~gJif*&FaTtvS$Qp(xgZpY5DANO2+88&1+`3NWT`GK4bcgd@)avlD>{H$51 zvCFc93v(uO+QQFT4>fCILi*WrC+U@ocg|kzl|BDcn_T@y;uyO5F2{e1&~KT*!%zs4 zNBs}GK)%A1=e77-R+L#=OA%h!&kvp1xSCh_I&^X?Q4e?1G^YICJvVq4v<8NT4rZ&w zN|du}@8d4@*n8fqQ9WR`CWpILhMs@}Ba)&}mkberAJS}DHie6{n%6yAjN1B`L_2jw z?VTz=UD5dOU&yg#FVsK?NY)yn5wzLZt`l${@}+}>XpMrI`<}&)0rAkuk0_&=Gwie} zv3N-~MHc;=7#v?!C8w1)$Z@E6TX$)89KC65oZTfFAgP*??yl?0hBp6-D?`3JiIdayN^_V&dtxD7+h6Yd`rLeQ$hV{f1+aj=iV-2K zTN4YLz3bVv;sneV1Ddk@af9;A?-+G+lYj#qX+W^g0Ra}6BMSw?N`__K|u#`lVixUK|PD zv8jxqXgL3q4pC;hAv z6$o>=*E=O5$>Y@tj9yptoi9Hdx@A_f;T`JZlw4m~<$72AJ$!%eDB>?(^0C$P4=(Fh zW=;Lwh~rvQr;P{geEz4{wMuoa)WfMZJZz#s9HFcv-kztQVvC;XTDkcb7rhr0;W-SE1DQye@w^iq!YP2~;RTxWoU#>@)556sk4@aX6O z7k|5$J7UZ+I^jhKWZAypj zJ z_yPv%z}Cqmuz|>!dlbmvf4C8RvHa@!!~Mm8kY^Tb=yiWHk~0z)i&LNI$rhB!=Rr5%e$zF(-fdHazl7=o5lcGdDUZ(9-+W98yt(UJEAQ zsH?X1SV83vs5|fn-6#@5LwauPT6%i7gT>>Mr4XXt^q`h5B>2ItXi!x|=+{=japWA%4szMoK zWH=y7(GR@<)JbgL8QXkaZZa0muYtCGd}^zg#4!!bQzoby=cH%<1`PyKZ@E0yJq%&ux|s2g7KsUEI&{J1s$ z)8CB!?sa6oY6_mZEzfn{KCPMmyKvu-`)T+0OW^>M`YKiS+AD*Ei&{vnyHya+oh;EUb!YdV> zWoJk6SskpiP%OPV9X#i^)Q_UI5)XUYT6D&62*2bXH7j)+VK2aLQhp=cd-VF2v2m_Q+Og*Or~@gzzwk1`O~YXL*Qb zj!V04z6y4i9N~_2wMJEefE`KtQdVrf6GY1u<`xv?SMzGh{+lH_YlWJZ3Va*G&qw|? zSh1H2&zYJ;hs5nuMCX~(okUyoqF#IUkE2H>GKrff7mZG9-xPqs?E&4CN!QXqmIJvFb+8>*dB6AdDzYo0~r^r(vt|uSYEwCzn zXSIB#$G$2?hovZg`h4jTsI*d;_LaV2s$#3^F7@#5r#gGI!w5v!}z&D6FF4O#0N1@*a?8s{P#qTqY`hPP>p2D z7fV`E@`g|guYcxy|BtNIuwY|W;h!#(z<<{MYb}C~Y=oCCtzVoqtZHdt`<39h_6gM$8Q+1`)g7CRftH+l(I$KCj?6#v2(W_}+F_)G1-gM7o z9d|Jp#1QrE?|Dj12H2b*ib@KsDK74M`t&{WiLclyCUJk(JDW@X2rmh&Iii${`>2Rp zL+iSGrlH~MRzB0aG0}=1voN8z>i3nuJZviKkRxgElPXB01u+EDw3HYQ)t}{p#DD~QSb#gG5dF3 z4^J9RGC1-7Hu5Aa1taDEe|=C!iP4aa{2^GH$^Sl6-$tfT*(5=n$}Dec?%!TL$Xf*ZkLBa0o-q`Ta_rVCF6j4ubqS7qV`$LZaohe#Dg2_zkGqAKg@fJ z)F0>DZbS($)fdluf48Ol#!Q7#r_T@d9MykWHcn}hO*&;!_o4qN&DC0MBmS9$VqnBh zalnj)n(RksPCQygJiW*Np30~9`28HkB^p-SH3HxwlG>%$;Q@SuP6aN+|1EGOo`{G@ z{xF%MwW+n=hjmVWuRXe2S-G=l(K2jXZ?FoH4?+4v94=G%rfcy)Y+ zX}eP`TDM|ICZp`#*}^m7K}9xhsUP&Bj>zKRwNS}Y#sD^2+!m^hfZ4~q7#(V_QnTN* zx*LNQo~Wyu#UR>Cz)ld!{u21-jU4&rC$J5#v=x^ zG0O#-m)lPD!&X>id^fPx=k_u7_hbw`fm0vp;-Cnc^5CPjxEMjh^~zaeuR!=06EiIg zjLfTlUjfx&Q*X$??8qi3j6sb3Z*3T6^!@i#Bb)nTz2Ql7)ZJ@-L&v4WX!ZE$&s7?e z-Fu}1l7C7oq69lv{kesL8=vgR>;dzQ9w|5Ro2NMZx1CQMFsAm1nASz)!Kf(<_q^L` zt!(G*vCJ3|DUt&Cw|2ML4*jU5UOGGY=dWM?qe@Zcr9|1Et1=Y$1?g5JE`HGQ6oP~K zPyhj6-|jE3%?i}TgIzFaT57${{GDIi#IL{H0p{iJ?AkOU%AwOM#v!P*|8)I z48_b=TvF%=;=P(5|H$tjwYrQ4n9;U!J^L5z9*V&<`MIW+4Nv|d{?_RfIwm-paBpE` z0Sso9qHquJ!_L-!-nEHmfl1h(H9@+bmbXGSo)On*!&Sb}q6?@fwlUrox%R}T?|Ae* z>Ft8PqNHmFhI>wU)Rp8(!j~rn5cKzYr?9U~r48xTBCUk1$T;Aw0 z`KG)tld^(0e>iW{%c8{x`hR!HH<4HPH5+w(P}g@Hb-9V(y_GDmOJ(F%>=&PAnM@Fp z1K6)zBEMRGI(A}+GWAE`+2CHVl2 zR-e>q$cQMh69t<4+F)&9DX#+EedmNO2_;k&`NR6do!zoqxJCh$o{aL+o}(5vo-ssC z_mz~c5W^is9Fw#GLkr)x|Bz*e{;`b^2E60K`ToUqD`4!?o?+t0+pqM$nw&5Pz-$)` z;TGVY%o^U%5OPU78+f>TE)d)s^oOu^8sg+Py;$c|n^t`)+7?68fXnJ&0FfNYQy%>N z1gKMt325sa6zEi*R^;u~Y?N>fggfC_C_y{jKpT5Vwej^T?`!z9*|o$fvA4{!a#S15Wb$hzi76Y=kh**K)zZMnn1 z2E6U{{<~{8`q*{+6N-{^gIzdp_nqe#PlYCSbN1!xLQ=f534!J*y7&BK1?6P&zgrkq zCaBE-RqsX!F4n@E;pgAox;gA?P#-qu`a%O_F@gX|+5Owvw*|s}gI(OwX_3$x&OV1k z!MNVrlS$+D@a;^EsiT9Tt#RHm31bi4x8iP#FJ>QKRDUgLtNKVYaisWgrEi#`?72g` z#qEtxx{zZ4dxp>Ct8^%hLw5Y+cxO+5McbuNbx>i?0~P;-w9Cr3xoe|(R?LB1SyfRT zr1C0?i*3v|1Bt4ZT=lLXcvXYKL{1;DS9Ev?)!VU_861@>*{wuQT2R?}3=rT_8P%r1 z`fH)kt+enNXOrg9z(b@k{~dmGA4Z#4lAjpw$NqwQ*n;_`lNN3zxzb6PV6>d4v`2*a zLagqFnZ{WNEW;gO+uftu9*q?Dta;h;vs{~`b%V8j~BM7PsCq(F?&h%b+D`PapIR3{fkt}wSE z`b$ghYekrzwc|kZG0)1U5dHxr9yqWbdGkFTIiBIA*!PJgoXf6G0i!1lpO%!yXGJUZ zVJ3H%Xx2K(Q3dZFkt^IJjj!~X?%*{k)BVYRPzl zF>IS;%QK^q&0jQ5mj>-MOu&vpA5a5GCUv~aK2K=MH+;rG%I{d4Jitw5-n+}689Qmx zYtSvt_+vLSwp-ag|79>$9ROChcFcZ`k^NTJ_HGjFozymm$l}M`ztRIzx4BDb{{GX} zPIx0hZ&HznM5#RjJ18uGONipLr6`?_5dTR>u-l7*H-$?Dwax~%Fp4Y0c}=lx6pmOG z2=STilCS!XDxHGUG=>W*=hw$Oz+GlYZE(J%layJ``Zc!2iGzu_@cl<4l%UfPH0ceP zoW+^7b=FEa1)81(hwPM z3zf|}x}#T5h&;TbfQz5h&YaO*pBQQ|?+nt7=SFTvBoXhQnlF)~iAhT%uRSL7X;?k} zDLYh6rBX6mzW3t0{5ZJs>UqD11Tr^OVPD{dc?n85QLmFF|Fdx5AZchEBJg`qs3G=~ z7h_dfRD6TXwAC1$&D;=0zZsFldsK)^hA=}$7>%jsGpQ4msVmnLl-q0&{Z1JoD<>&n z`{ivwb@fq#mITiT7mhhe(%LUVpUdeklLpyJz!*IcbDsv1?MQ1=&+(zhq*yUmoW8A~ z%O?0qEMGTnG+O>|my=Hp!h@Ix%XgOey9y_>963W^KWRoc)o7_Rb)Ir&6!d9T1XEZ&i?g`=5uGXuKkuQ^VKHQRm0fI| z!JuaC(JV->GZ0P6f^|fYON!CavhwOIekZ?1JR(0l6Lkk@(@h9=tkr1e<7SmxgZ_2= z3G#U6&Q72Ja1e`E9B#1L@q(t{WQ zdz(&68jFxE&*bRy>HQ-m9#zp`>UL!~#>JFSozyI}UcWrN(`zKMCDMdQ7-wsuZ!Nz{ zNH*^8lNBXp*MEiJ%n@mVXsx^UU-oH44@aei$DEHG8jbBGoge_n^yDjP2Wq4?HjXG^ zHFlLXYyn~gxTrv%`9`|XKXsqM$L|u<#J9HR8R0j6(#s)eKhf0!z!A>wPyy6|LRNuS zfM^4E0y8j~>nu#|)_xS^xr}(RcC}I9lrHL;u7dx=p+RLVgH>@2ewaf)sh5}!_oyg1 zKU$CK2LfeusdTR9(HA-R>g?T5a4*<~fWhqR3n8&}gEUMFovIFTlriPU-JFoNYAXIL z2ZkoXArq~TslAQSANmkW9UUFv3RO@WC5Q+&dx`uEA_RKWM5P7P?>Q6C!CEHc2B18X z{85%c*EDysfx{J`cx-~4r+822SZPIVh=NfrKLB;G!O#|8lCzwqT! z>9?$_;@gA!X$q#Gaz!X*T2jR-3g_mx)+;n%xmrb7QK>T8CDc>fu*299>dwWxi}mLd?8!J0OEym->~sAB z09(?57^lmQQ!e-(q<@aKtGrnZ9Q)j!PgBqLJn+c}BU`;WBr7TbTvLT3zDbFA-?}yT zcUxN?iWtUCY=hlHxBxK>r3Y=;uSGik5b`H892pc>NH-~-u^z^(qLr(--xEIwU`{Ut z<^L)ki-zJa$VM7WnG3pY(U<#?7tFYe#>EC+e8^39!db@D3uJ325Ipigyo~Ao?5Bd_ zQp2tz7=#wDI|ZG+Yc|KKXJp__Wh++N-^@0X16OfMILfiAB^tOH(ILK)7^8}H+)fRx z#Y(otcTSk=5>wt`$H$-Rl(It!QEKY+1?D)aZaF*R<%Tm5Mk=sJ)!f6UM3z%JEFc>x zVT%6T37mWl#3AwgFrEe#f#|?4Az}|3V^09D^>K#X!p$aN;hAss+a~Fx;AABRW?)cfy`O;5rOltcB3oeT#;oTi%+ zLWh+~zmswVvcQp2vNEK*E?W>is*o(OlTlIqexhZC^iztiE8}=cQ2?08=xv8EdQi~e zF=PE9ovsS(+KRT1ehwWPvW%rn3ZMD4&DrgURr=!l2F?!sU@$ z{uBY)paI6k0d{|0E!?M;e!m1c5MtN91sRpyQr8;qv-H&YfUw?z0tH>>Kv+|}((xX1 zqf@X?m6ZHR9dwG+f94)Bj;dEwP0h(x@icj{qf6keQx2}&f`K-KD1<{J)AX+pgtYx* zB?}e9MqxGy6u}>4xe|KQ$pvWzGB5|Oci{?f45i5R?4aUo3rO$ElnI7jQ>P$4aVrWY zrNXZK$o?Re%9QpiAaY9jHCx|*`qqeSz?3Rt^I8&fuyQ7IsgN?>$RXJG@__VTBsc75 zu(Hv#q&5nyK5FaGhw_q5vNM;hOEpkQk1p{?MAAbKSgaR~HGaYBg! zZAi_qOh`n2LHd`Jn<-i53Kd(QA9dO?dGd}6cmx^>{6warB#(1oLSG^Nolu;;h^dO) zmBEQ2#>2r?3u_952QK}U#ugtk$3j8+*EnWFGrt2$AwGd+pXYunp^%&uY<=;K9pnH2 z*`eO8akW?X1?ru!t#h%pj8m1#4;9WNT+znFb@7S&0EPw)B3`J7F@t57g>YuQ94;)v zA}AhRRvEBlZ3vU)P@$=3g3&DysSPfYbrzmK$ycTxpd<+eN)dHkox9@!!hKs-_uwL{1vQIS$oAj4=KVdOs1B`I@@CYnW^A{zFCzS|VP zTvkv65(@}bF*dKS&%Sx}*t_&|Ix0IFmh?y8tT+rDSog(&)43I$?sIjASuu;b?%>rIh^-H%^1tg>nE7WQPj zZm@2}z{A(^FI`>8;KX6I&z0;t$U-ll@LcmH3wZ0ST*5&jdu@?#vJk<}6$F*KS31E_ zcE1T#19zjPG=8GI=w&y;^|6?dqp=M#bM#5tOa*LHlW}a8axq}MKg`sQD6?T_YB{MZ zrB#q%;^0C|t1L7OX+Aa#90jP3B1y0T2Bh{|1*ZpR?xo~mpz6te#?ADYXh0ojHr_63 z+a39+RnR_-O}Huoajmq0IWNT;9KjJD6gvnPad!JBYR^V>X2(E%4mRceIVsnuKG7XI zN@0*8Y{8sY#ojsERWHJ!1w@a}6P8JLIUzlSzt_bG{sMhc^P^tGC038TSnGC(=?a*` z3?p(H{H-ZlZU|NnL4Z=z4Cy+sIJKObKp25PKpFg5^a+Tb;i{XARh@Kq?+J%Wa;Ace zea>jLZ~^y(kW8I*YEQW;{xTk`u&dra!8IFmlNh`GW8ekU8O(91;%jeFS%zYBC^iQh zh-DvByv@Po<{$ybiF0Ta#Q(vG$QF)B&-4{#6P{p`shnBUXkd+@ivy_BD9LJU){0yv zLWLr7UUqz@cll{XOpXO%nv?#M!b-=Wj*9c9sEuV~z*H~s(RLky5-J(>Hgw|g^&yqE ziv6kyXX7uwUZ&3<{*|5&{-#~PbFqEiH-*L6|o*s zw9Ey7Ko5=PgpY;0G-I6~IEWw3UU^$GMI~r|4ecLDG0SBfd5IV$QN;v-I6Znj)RTC= zb7H5_#b2&o1vJRO^z=1*6$!zT>{M zQ=encu|E%lB5|jHKqBMXc1{EO5pJOJ6olW|@^W^N(SyTTC2j3zRVws^OCb~q@Ies~ zD?C*ifrL?B(za@H3x6z{mDAK$J*ooPflJDTrF@5Zk!6H0RZ4%Vqau3O;XwjVPn9mT`346P^-i%!n!sVcp)hZ;i8m^T%hGXxf)lcScTIZV)RFo8`9i2_09 z%!h_3IKB1cE6L|R!Jq2EVHVbY)1oRGUEO>%ao^m24i^}1U&?03_ZwIay*SQ^qhR}^ z-Lp}kjH!?;kB&?MJ`q3Dwb<{4SV4N{J`8IgrhgoY8NyL?c&b;z$GGFTNc;yI4W>4t zHV33o=dsUEcAn9|<2uM5bFLnoRya;SPR6+dMq)GbhKdq0cBsfB0qgOz2;Lm4{USMI z+s?pD9*Lg-eGam^=TRA~tTklZ{M!v^`+fPLtd!^&EEXyc35`$66&ng95W1K+8PLLiKqd#PS)7yeUlc{%5PbVG+>N0-DTg=TQyNu~4C!E9X7$2i_&{-)*V0=B^eK(67%fB2w>igLkw-^{A>~(GbCDoB_pBby zCBvs59u_psbhU&3VEp;znxMH>pbe8P>#+5wdA2jW=L2&@3mrmV2R` zL=YVr_YTlMl3}sjp@@}vT|uqV3G&20K?Rob^cm>!Nhd4b5w4+(Y$vE#zNc;xGRK+X zE*G3T;VJZQlb4fR0=-A;m+?8~gb10GeiY>6kRhZtEj9qxZ}@u;6vEOfIL&ETpI|Hj9zXh$9m~-{W zoYX;_x0+NB{5lXm%||wJvQ>}KuW$~rvV-iMnZ;m84R^df`6+A%JSvVTR_<^H>jolk z##)+4{ocA_EhQi_9*tkcCAtNQBlf}djX{tAS4EAK-OFXfimG)D#agP_h(#5L7Ar(x zzaT~%KRmKpKqP3@lEGH-C8!9`y6AkMver#DozBZXlsc$G{-E*_xFK!riOVu{KN4&$ ziO(GP@-vGcgsN*qXR4mn1rcO#)wjHW;i$lRp-x8u*SZ{8XyKTvW8#&pSDSY|W7hgI2*1ytinC;x)uYbrW#Wfr};k zF{O^i_CVxqV%IZ*Xsw5uW{UM(|3e*C*zstQa>;?Q7HVR!o`s!+JVC0)>+CVL6gAI2 zrdCs7n*>dQr%@g4h{wj%%2f??fZlN~ZL%(o`(+bnWJULWTq(Y;0KFU#up}PM`pZ3` zAKv?nj5|&lYr=-opkPC*A{RgI(+`7j^l|9;N2L#qv?(yRj|s4ZU7|S;HYg~AG=7|4 zHkm(%tfbIT8YdYw)Pz%os6avH%m|X!r z9^SAAq~1gD*8x(kY<7oAb-Y3ofS-=eB;vlC)Vqa2WPedqRS5>{$^__>VYK_YFBiC0 zol(>x^}Sw9hTI832pVsASmVGU+XI~ny(sJjVJ$QNQtm?HSKwj?)+XDlS2&#El7_~! zUnPR!4w)5^=8OhkA=oMl7mT3E7y$A>E`W3(lQ)JCh)kj`pPA!(f-6jiqXXK;{qRyT zH81>AcY>2M2wNS@S$sewasWwr&mjU3z=?S&FoC(E7jd$=y0FVYC>NzKOHK$WKg%HB z&n_4%Oo>6M461=DxI~BQ`q<}=$n!c_Wv~+1=Y%)nqXC4noI0?jv*4U3htWebV*H&G z>|G6tx2U_dMb>*NtV_r8=qnafmRh$5{x#{bO+^e3B=#cMw*UvN`-oZ#(Z@#>k%%@g zF6aA1aYtB@*an*8t)O^|7De1140uSm0QnyVgZEeFWct_%<|Mt7%10|{eMCx#U%^NJ zDSgMTFed;}RBj8yeqTHCaA+9+eytB;N#nDk2kD^(coh2UyUnYP2z=c}&ymn?S(v zh@9{$M}SVS!sO#Z-CCe4%w4BCCds}i-p>TTWz3#|G+K;fHAS=wlb+?^2CMc{kT}@) z%xHQzl~Q!gk;SC+sj}%R7RMPvt8R=Yabvd=$^mAp*q0(Dw!?ix;nW7y5+ctoNp|(& z@JSHY1a@p)8U<1kI43kgsnl)QRR93#Lh%7R-*D{;STBd;P&Au_5%c|L-1zR)F6;-d zxJie|ELEbT6Ux1!6`Lr|qz}z2u+pAs{xaV)fJcx5LpMz;){%favE@1&=3=rzGV-xC z5pq{^jr|mJm3c%M+pWP9W&@U?v%AM;*DO!;6Xh>;jm8On`jY92kD(m50D7V$;ff?e zVJz_ie)Mc1U1+Ni{jiSQjMXc=57--7Cfw$J2Gxa(NJ-jaWuhs3-*2Vri4UGMc#T{e z6f6!2!luTq@5Qi`ZuFhpN4DxBGz_*x;*|LGDeV=66SLRnOS75D@rcDMhj+1I)0dt^ z7%9ES7zSos?#tjw!p8>MLocfw4A`MYY>^F0EFch=p6AaT&XEU7=G5aSl)!iM{slEZ zn8{U7POhrC5rCEVxifaQ6TMf|?>O#P9Pd%pCp!P6=yVubFk8V0j%q;s;};rmJf^SV z@eodpHCV^ozT0c}g%(Pb*%SQC^Uol!@djLrWm{I5>fp2|!<(p~x_$B>m%bsY~>-oB(ZH zoJSvlstulQe0lggDB@SMSq6fU{YMZE_t|j@JIj!c6?pc70l7Ut$P>(W|I0wuv3hF0 z6`oA?M9qYnX5!PRGSnKTVXejb4C9V~qv7Y!A_I5!kC&ti!JQh3_&HfdXf>v?=MBGJ%>2c+30=IFH|s2)N3sB-z+%$*Rmd0u}RiNc6yuy_3Fz2 zYdW$uBkEM+mt`-eV5IVVmaH*S9+Ve4n%LQ=m4Jz8MGrsE{G5uHM^se4{pl&)#rc*wWyPy=exPx(fE~Mb(iYiME7;irUmPa;J-txs^DGV zwlQ1dE1y(4Q3EScgiE@t+7;g8ckgUlSKm=#hvLE;z6{oZ*&W&bPm`yY?g{7M7UTn) zy8{QRpK|4K{!a&h7#c?UBZWc$Q@zNKhbb2`wNm&h?`Y5W5%BsMG}~;CJg0`7x%J7U z%zvGu1Z<{?j`kiyfb11>Tw|~Mx%sDi=J|=yaBw1)0;}uR&HpP?OpLZnjeZK8t!14e zcf(?S)@gCL`>Nd68TSwkJK{!jYJdx~?9o+gdQOX%tBUm@ zN?d;ZcnSuw)+%FO@7mJ{89J^aEoFH%8ltd089t#bGofxhp%u0o%lp@YZZ_A~RzGhl z8r!udn>3{27FogEbr^G-?ItBwo#UP~{=^i*(rs<8vE@bh984Dg6{r4L zBMf>kDEN{*y$$W@$3f{dI0>IS>L3=1jSdg0J62<;19-pMKd*QigbBO0O)-%cu2(?j zG$=Zk%)9@fQ7_NCCS5o1bMWMaOQnqJ0>R~GM8!B8!r-zKD40^q(ciC>^324SlM>JA@vt1`N{C4T}f;et<8 zmennNW14!;dDp!c*9bpv&bQs0BFoo4WU4W=eCFckt6vZ2`1L<{^E*8^F$qO_WoC39 zVnRv4UQG8=%Aqv|>O+GOrPVKLu=#Yn(snJBx+3PI%+kmcYQ!}$xiE3z0d2$hzX)()jg{#? z2rw-PsaV{Q=DoG=>&0aaAEOgKlZ-pV@LpF&HaGrt0rZEm*J-~<{8GpF>i%KBjzf2Y z!as!VYsyb~xNl6(c8YIZDtkF0MnCRI@_PIY*D?@bG`v9F(ly>2WNs{k6z%I7Y8Hte z0{JIEz`5=Y_4*w5aydP?rby!MC>-S0a!#i#MxZtKZXU$*8_(UT4O{FDc~8Xrt6$~e z@ln;DV8hm9L#LWa>u}+zC_)r_#7ib;_krKL{lCH7VqpSp+iVn=CwKo)S>?y@9g=^> zWLwLmzPVHOT7Cv-PJ_)RX^8;waqz@hPECpUx7lK_4%M-d}*n*H!vQvDZ2Xy&*PFF&2!M{?8X^LEc7`?TKdJ znk-HT4>JTKvffgl0^Mz#0`DH?qS9mdu_xPB%aY#vsZ37y(xu6X&@M{XumVfQ2*Tyw zsFDf_QmjK!Imq4I4GDg4K}d%SJCK~SpRkZpGjp%6uR|h z-2Fm4L?lLFgPg1*EbwQrE2mufkr&B;*wpJH=GPLhZ=SJ!QwREV+rHP$*N7t^r`Y`Y z(^)`99DzX6C}%8Jf8L^$CZ&CIrhxm4=d|gx2;{L8)5D7AlbNqV++LMGA?g6{5h!*i zmsW$zf2MYh>m-<(VHvDr3P&Ul3qq`gI|(a!&{U^Y(BC2O9UyY|R4J0&;A zxq|WE?>g`ysn@Ix-Eydfjid>g!faT}EGfOkZ>AR7UTQ%oCJWK$RJ#!0*i9D+OD%Pm zZ<*L{QA#EE=zZq6-v6EX51aE`pS29`w8^U9)Z%ZRUy#z(L9(OtqXUql>saxBVu~D# zhsj&IQFqP%LCpX9DhSJTPV!cutVd2$y%k%raRQG`<E#=U3fXPXzgMNoJVgJgiW0V}HCqMlaWFrP z8wi!pZU3RzU&rebjgK(We+o)@(g=HIETMMK5ZfY}+v$>TggYy%SB_wvQBOC(zGncQ zcA@vnOsnoa8gj(`qy+i|%-YwYcJbwr!?j2w!GYMHs{(Nq`lY;)`_;|!5<-&otD|Qc z2|oiwV_b*RZI&e{n>xy9WL{_q7@Jzx#j_+H_QBaPjF(m3_t~}R^>E@& zp-$mg@5Vi`@kIZ2t=)I93@!BvDuHK3S6}%7?Q4&;&M%G}4XrO2(akG$;>S}dp*