From c51375d99b2c5119eeeb12b19afe17d0b27b25c1 Mon Sep 17 00:00:00 2001 From: noah Date: Thu, 13 Oct 2022 18:29:27 -0500 Subject: [PATCH] Finalizing NPDR --- .Rhistory | 512 ++++++++++++++++++++++++++++++++++ .~lock.lab5_expression3.docx# | 2 +- Schrick-Noah_CS-6643_Lab-5.R | 183 +++++++++++- lab5_expression3.docx | Bin 138090 -> 139763 bytes 4 files changed, 694 insertions(+), 3 deletions(-) create mode 100644 .Rhistory diff --git a/.Rhistory b/.Rhistory new file mode 100644 index 0000000..8fe61f1 --- /dev/null +++ b/.Rhistory @@ -0,0 +1,512 @@ +################# Max-Log-Likelihood ################# +n <- length(g.breaks.clean) +kmin <- g.breaks.clean[1] +alpha.ML <- 1 + n/sum(log(g.breaks.clean/kmin)) +alpha.ML +lines(g.seq, g.seq^(-alpha.ML), col="#D35FB7", lty=4) +# Homework 4 for the University of Tulsa' s CS-7863 Network Theory Course +# Degree Distribution +# Professor: Dr. McKinney, Spring 2022 +# Noah Schrick - 1492657 +library(igraph) +library(igraphdata) +data(yeast) +g <- yeast +g.netname <- "Yeast" +################# Set up Work ################# +g.vec <- degree(g) +g.hist <- hist(g.vec, freq=FALSE, main=paste("Histogram of the", g.netname, +" Network")) +legend("topright", c("Guess", "Poisson", "Least-Squares Fit", +"Max Log-Likelihood"), lty=c(1,2,3,4), col=c("#40B0A6", +"#006CD1", "#E66100", "#D35FB7")) +g.mean <- mean(g.vec) +g.seq <- 0:max(g.vec) # x-axis +################# Guessing Alpha ################# +alpha.guess <- 1.5 +lines(g.seq, g.seq^(-alpha.guess), col="#40B0A6", lty=1, lwd=3) +################# Poisson ################# +g.pois <- dpois(g.seq, g.mean, log=F) +lines(g.seq, g.pois, col="#006CD1", lty=2, lwd=3) +################# Linear model: Least-Squares Fit ################# +g.breaks <- g.hist$breaks[-c(1)] # remove 0 +g.probs <- g.hist$density[-1] # make lengths match +# Need to clean up probabilities that are 0 +nz.probs.mask <- g.probs!=0 +g.breaks.clean <- g.breaks[nz.probs.mask] +g.probs.clean <- g.probs[nz.probs.mask] +#plot(log(g.breaks.clean), log(g.probs.clean)) +g.fit <- lm(log(g.probs.clean)~log(g.breaks.clean)) +summary(g.fit) +alpha.LM <- coef(g.fit)[2] +lines(g.seq, g.seq^(-alpha.LM), col="#E66100", lty=3, lwd=3) +################# Max-Log-Likelihood ################# +n <- length(g.breaks.clean) +kmin <- g.breaks.clean[1] +alpha.ML <- 1 + n/sum(log(g.breaks.clean/kmin)) +alpha.ML +lines(g.seq, g.seq^(-alpha.ML), col="#D35FB7", lty=4, lwd=3) +plot(yeast) +hist(yeast) +hist(g.vec) +g.pois +g.mean +alpha.LM +alpha.ML +degree(g) +sort(degree(g)) +sort(degree(g),decreasing=FALSE) +sort(degree(g),decreasing=F) +sort(degree(g),decreasing=false) +sort(degree(g), decreasing = TRUE) +head(sort(degree(g), decreasing = TRUE)) +stddev(degree(g)) +sd(degree(g)) +tail(sort(degree(g), decreasing = TRUE)) +plot(log(g.breaks.clean), log(g.probs.clean)) +# Homework 4 for the University of Tulsa' s CS-7863 Network Theory Course +# Degree Distribution +# Professor: Dr. McKinney, Spring 2022 +# Noah Schrick - 1492657 +library(igraph) +library(igraphdata) +data(yeast) +g <- yeast +g.netname <- "Yeast" +################# Set up Work ################# +g.vec <- degree(g) +g.hist <- hist(g.vec, freq=FALSE, main=paste("Histogram of the", g.netname, +" Network")) +legend("topright", c("Guess", "Poisson", "Least-Squares Fit", +"Max Log-Likelihood"), lty=c(1,2,3,4), col=c("#40B0A6", +"#006CD1", "#E66100", "#D35FB7")) +g.mean <- mean(g.vec) +g.seq <- 0:max(g.vec) # x-axis +################# Guessing Alpha ################# +alpha.guess <- 1.5 +lines(g.seq, g.seq^(-alpha.guess), col="#40B0A6", lty=1, lwd=3) +################# Poisson ################# +g.pois <- dpois(g.seq, g.mean, log=F) +lines(g.seq, g.pois, col="#006CD1", lty=2, lwd=3) +################# Linear model: Least-Squares Fit ################# +g.breaks <- g.hist$breaks[-c(1)] # remove 0 +g.probs <- g.hist$density[-1] # make lengths match +# Need to clean up probabilities that are 0 +nz.probs.mask <- g.probs!=0 +g.breaks.clean <- g.breaks[nz.probs.mask] +g.probs.clean <- g.probs[nz.probs.mask] +plot(log(g.breaks.clean), log(g.probs.clean)) +g.fit <- lm(log(g.probs.clean)~log(g.breaks.clean)) +summary(g.fit) +alpha.LM <- coef(g.fit)[2] +lines(g.seq, g.seq^(-alpha.LM), col="#E66100", lty=3, lwd=3) +################# Max-Log-Likelihood ################# +n <- length(g.breaks.clean) +kmin <- g.breaks.clean[1] +alpha.ML <- 1 + n/sum(log(g.breaks.clean/kmin)) +alpha.ML +lines(g.seq, g.seq^(-alpha.ML), col="#D35FB7", lty=4, lwd=3) +plot(log(g.breaks.clean), log(g.probs.clean)) +g.breaks.clean <- g.breaks[nz.probs.mask] +g.probs.clean <- g.probs[nz.probs.mask] +plot(log(g.breaks.clean), log(g.probs.clean)) +# Lab 5 for the University of Tulsa's CS-6643 Bioinformatics Course +# Gene Expression Statistical Learning +# Professor: Dr. McKinney, Fall 2022 +# Noah L. Schrick - 1492657 +## Set Working Directory to file directory - RStudio approach +setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) +#### 0: Process and filter data +# load gene expression data +load("sense.filtered.cpm.Rdata") # setwd! +# load phenotype (mdd/hc) data +subject.attrs <- read.csv("Demographic_symptom.csv", +stringsAsFactors = FALSE) +if (!require("dplyr")) install.packages("dplyr") +library(dplyr) +# grab intersecting X (subject ids) and Diag (Diagnosis) from columns +phenos.df <- subject.attrs %>% +filter(X %in% colnames(sense.filtered.cpm)) %>% +dplyr::select(X, Diag) +mddPheno <- as.factor(phenos.df$Diag) +# Normalized and transform +if (!require("preprocessCore")) install.packages("preprocessCore") +library(preprocessCore) +mddExprData_quantile <- normalize.quantiles(sense.filtered.cpm) +mddExprData_quantileLog2 <- log2(mddExprData_quantile) +# attach phenotype names and gene names to data +colnames(mddExprData_quantileLog2) <- mddPheno +rownames(mddExprData_quantileLog2) <- rownames(sense.filtered.cpm) +# coefficient of variation filter sd(x)/abs(mean(x)) +CoV_values <- apply(mddExprData_quantileLog2,1, +function(x) {sd(x)/abs(mean(x))}) +# smaller threshold, the higher the experimental effect relative to the +# measurement precision +sum(CoV_values<.045) +# there is one gene that has 0 variation -- remove +sd_values <- apply(mddExprData_quantileLog2,1, function(x) {sd(x)}) +rownames(mddExprData_quantileLog2)[sd_values==0] +# filter the data matrix +GxS.covfilter <- mddExprData_quantileLog2[CoV_values<.045 & sd_values>0,] +dim(GxS.covfilter) +# convert phenotype to factor +pheno.factor <- as.factor(colnames(GxS.covfilter)) +pheno.factor +str(pheno.factor) +levels(pheno.factor) +#### Part A: Logistic Regression +# make sure HC is the reference level +pheno.factor.relevel <- relevel(pheno.factor,"HC") +levels(pheno.factor.relevel) +# also rename levels "0"/"1" from 1/2 +levels(pheno.factor.relevel)[levels(pheno.factor.relevel)=="MDD"] <- 1 +levels(pheno.factor.relevel)[levels(pheno.factor.relevel)=="HC"] <- 0 +# Fit logistic model of first gene to phenotype data +gene.row <- 2 +gene.name <-rownames(GxS.covfilter)[gene.row] +gene.expr <- GxS.covfilter[gene.row,] +gene.fit <- glm(pheno.factor.relevel~gene.expr, family=binomial) +summary(gene.fit) +coeff.mat <- coef(summary(gene.fit)) +b0 <- coeff.mat[1,1] +b1 <- coeff.mat[2,1] +b1.pval <- coeff.mat[2,4] +summary(gene.fit) +coeff.mat <- coef(summary(gene.fit)) +b0 <- coeff.mat[1,1] +b0 +b1 +coeff.mat +b1.pval +## GLM +modelfn <- function(x){1/(1+exp(-(b0+b1*x)))} +g.min <- min(gene.expr) +g.max <- max(gene.expr) +curve(modelfn,g.min,g.max) # plot for domain of actual gene’s expression +abline(h=.5, lty=2) +curve(modelfn,3,8) # replot with extended domain to see the S shape +## Predict Function +Predicted.probs <- predict(gene.fit, gene.expr=gene.expr, +type="response") +if (!require("ggplot2")) install.packages("ggplot2") +library(ggplot2) +phenotype <- pheno.factor # just rename for legend +gene.gg.df <- data.frame(expression=gene.expr, +prediction=predicted.probs) +# plot predicted probabilities versus gene expression +p <- ggplot(data=gene.gg.df) +p <- p + geom_point(aes(x=expression, +y=prediction, +color=phenotype, # shape and color +shape=phenotype), # are true phenotype +size=3) +p <- p + geom_hline(yintercept = .5, linetype="dashed", color="red") +print(p) +abline(h=.5, lty=2) +Predicted.probs <- predict(gene.fit, gene.expr=gene.expr, +type="response") +if (!require("ggplot2")) install.packages("ggplot2") +library(ggplot2) +phenotype <- pheno.factor # just rename for legend +gene.gg.df <- data.frame(expression=gene.expr, +prediction=predicted.probs) +predicted.probs <- predict(gene.fit, gene.expr=gene.expr, +type="response") +if (!require("ggplot2")) install.packages("ggplot2") +library(ggplot2) +phenotype <- pheno.factor # just rename for legend +gene.gg.df <- data.frame(expression=gene.expr, +prediction=predicted.probs) +# plot predicted probabilities versus gene expression +p <- ggplot(data=gene.gg.df) +p <- p + geom_point(aes(x=expression, +y=prediction, +color=phenotype, # shape and color +shape=phenotype), # are true phenotype +size=3) +p <- p + geom_hline(yintercept = .5, linetype="dashed", color="red") +print(p) +abline(h=.5, lty=2) +predicted.probs <- predict(gene.fit, gene.expr=gene.expr, +type="response") +if (!require("ggplot2")) install.packages("ggplot2") +library(ggplot2) +phenotype <- pheno.factor # just rename for legend +gene.gg.df <- data.frame(expression=gene.expr, +prediction=predicted.probs) +# plot predicted probabilities versus gene expression +p <- ggplot(data=gene.gg.df) +p <- p + geom_point(aes(x=expression, +y=prediction, +color=phenotype, # shape and color +shape=phenotype), # are true phenotype +size=3) +p <- p + geom_hline(yintercept = .5, linetype="dashed", color="red") +print(p) +abline(h=.5, lty=2) +## Predict Function +predicted.probs <- predict(gene.fit, gene.expr=gene.expr, +type="response") +if (!require("ggplot2")) install.packages("ggplot2") +library(ggplot2) +phenotype <- pheno.factor # just rename for legend +gene.gg.df <- data.frame(expression=gene.expr, +prediction=predicted.probs) +# plot predicted probabilities versus gene expression +p <- ggplot(data=gene.gg.df) +p <- p + geom_point(aes(x=expression, +y=prediction, +color=phenotype, # shape and color +shape=phenotype), # are true phenotype +size=3) +p <- p + geom_hline(yintercept = .5, linetype="dashed", color="red") +print(p) +abline(h=.5, lty=2) +plot(p) +abline(h=.5, lty=2) +print(p) + abline(h=.5, lty=2) +# Lab 5 for the University of Tulsa's CS-6643 Bioinformatics Course +# Gene Expression Statistical Learning +# Professor: Dr. McKinney, Fall 2022 +# Noah L. Schrick - 1492657 +## Set Working Directory to file directory - RStudio approach +setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) +#### 0: Process and filter data +# load gene expression data +load("sense.filtered.cpm.Rdata") # setwd! +# load phenotype (mdd/hc) data +subject.attrs <- read.csv("Demographic_symptom.csv", +stringsAsFactors = FALSE) +if (!require("dplyr")) install.packages("dplyr") +library(dplyr) +# grab intersecting X (subject ids) and Diag (Diagnosis) from columns +phenos.df <- subject.attrs %>% +filter(X %in% colnames(sense.filtered.cpm)) %>% +dplyr::select(X, Diag) +mddPheno <- as.factor(phenos.df$Diag) +# Normalized and transform +if (!require("preprocessCore")) install.packages("preprocessCore") +library(preprocessCore) +mddExprData_quantile <- normalize.quantiles(sense.filtered.cpm) +mddExprData_quantileLog2 <- log2(mddExprData_quantile) +# attach phenotype names and gene names to data +colnames(mddExprData_quantileLog2) <- mddPheno +rownames(mddExprData_quantileLog2) <- rownames(sense.filtered.cpm) +# coefficient of variation filter sd(x)/abs(mean(x)) +CoV_values <- apply(mddExprData_quantileLog2,1, +function(x) {sd(x)/abs(mean(x))}) +# smaller threshold, the higher the experimental effect relative to the +# measurement precision +sum(CoV_values<.045) +# there is one gene that has 0 variation -- remove +sd_values <- apply(mddExprData_quantileLog2,1, function(x) {sd(x)}) +rownames(mddExprData_quantileLog2)[sd_values==0] +# filter the data matrix +GxS.covfilter <- mddExprData_quantileLog2[CoV_values<.045 & sd_values>0,] +dim(GxS.covfilter) +# convert phenotype to factor +pheno.factor <- as.factor(colnames(GxS.covfilter)) +pheno.factor +str(pheno.factor) +levels(pheno.factor) +#### Part A: Logistic Regression +# make sure HC is the reference level +pheno.factor.relevel <- relevel(pheno.factor,"HC") +levels(pheno.factor.relevel) +# also rename levels "0"/"1" from 1/2 +levels(pheno.factor.relevel)[levels(pheno.factor.relevel)=="MDD"] <- 1 +levels(pheno.factor.relevel)[levels(pheno.factor.relevel)=="HC"] <- 0 +## Fit logistic model of first gene to phenotype data +gene.row <- 2 +gene.name <-rownames(GxS.covfilter)[gene.row] +gene.expr <- GxS.covfilter[gene.row,] +gene.fit <- glm(pheno.factor.relevel~gene.expr, family=binomial) +summary(gene.fit) +coeff.mat <- coef(summary(gene.fit)) +b0 <- coeff.mat[1,1] +b1 <- coeff.mat[2,1] +b1.pval <- coeff.mat[2,4] +## GLM +modelfn <- function(x){1/(1+exp(-(b0+b1*x)))} +g.min <- min(gene.expr) +g.max <- max(gene.expr) +curve(modelfn,g.min,g.max) # plot for domain of actual gene’s expression +abline(h=.5, lty=2) +curve(modelfn,3,8) # replot with extended domain to see the S shape +## Predict Function +predicted.probs <- predict(gene.fit, gene.expr=gene.expr, +type="response") +if (!require("ggplot2")) install.packages("ggplot2") +library(ggplot2) +phenotype <- pheno.factor # just rename for legend +gene.gg.df <- data.frame(expression=gene.expr, +prediction=predicted.probs) +# plot predicted probabilities versus gene expression +p <- ggplot(data=gene.gg.df) +p <- p + geom_point(aes(x=expression, +y=prediction, +color=phenotype, # shape and color +shape=phenotype), # are true phenotype +size=3) +p <- p + geom_hline(yintercept = .5, linetype="dashed", color="red") +print(p) + abline(h=.5, lty=2) +# vector of logistic output probabilities for this model (glm/eq. above) +# prob = .5 is the threshold for prediction class = 1 vs 0 +predicted.class <- as.integer(predicted.probs >=.5) # predicted class +pheno.factor.relevel # true class +# vector of True (correctly predicted) and False (wrongly predicted) +correct.classified <- predicted.class == pheno.factor.relevel +# sum correct.classified +table(pheno.factor.relevel,predicted.class) +predicted.probs +predicted.probs[1] +as.integer(predicted.probs[1] >= 0.5) +# sum correct.classified +table(pheno.factor.relevel,predicted.class) +pheno.factor.relevel +correct.classified +length(correct.classified) +sum(correct.classified == "TRUE") +# accuracy: +correct.acc <- sum(correct.classified == "TRUE") / length(correct.classified) +correct.acc +lr.fn(2) +# logistic regression function for one gene row +lr.fn <- function(i){ +gene=rownames(GxS.covfilter)[i] +gene.expr <- GxS.covfilter[i,] +gene.fit <- glm(pheno.factor.relevel~gene.expr, +family=binomial) +coeff.mat <- coef(summary(gene.fit)) +b1 <- coeff.mat[2,1] +b1.pval <- coeff.mat[2,4] +coefvec <- gene.fit$estimate # intercept, gene +pvec <- gene.fit$p.value # intercept, gene +c(gene, b1, b1.pval) +} +lr.fn(2) +# initialize results matrix +num.genes<-nrow(GxS.covfilter) +lr.results.mat <- matrix(0, nrow=nrow(GxS.covfilter), ncol=3) +# for loop the function to all genes +for (i in 1:num.genes){ +lr.results.mat[i,] <- lr.fn(i) +} +lr.results.df <- data.frame(lr.results.mat) +colnames(lr.results.df) <- c("gene", "b1", "p.val") +# sort results b1 coefficient p-value +library(dplyr) +lr.results.sorted <- lr.results.df %>% +mutate_at("p.val", as.character) %>% +mutate_at("p.val", as.numeric) %>% +arrange(p.val) +lr.results.sorted[1:10,] +lr.results.sorted[1,] +lr.results.sorted[1,1] +gene.row <- which(rownames(GxS.covfilter)==lr.results.sorted[1,1]) +gene.row +gene.expr <- GxS.covfilter[gene.row,] +gene.fit <- glm(pheno.factor.relevel~gene.expr, +family=binomial) +predicted.probs <- predict(gene.fit, gene.expr=gene.expr, type="response") +gene.gg.df +# Plotting Code +gene.gg.df <- data.frame(expression=gene.expr, +prediction=predicted.probs) +# plot predicted probabilities versus gene expression +p <- ggplot(data=gene.gg.df) +p <- p + geom_point(aes(x=expression, +y=prediction, +color=phenotype, # shape and color +shape=phenotype), # are true phenotype +size=3) +p <- p + geom_hline(yintercept = .5, linetype="dashed", color="red") +print(p) + abline(h=.5, lty=2) +p <- p + geom_point(aes(x=expression, +y=prediction, +color=phenotype, # shape and color +shape=phenotype, +main=lr.results.sorted[1,1]), # are true phenotype +size=3) +?aes +?ggplot +p <- p + geom_hline(yintercept = .5, linetype="dashed", color="red") + +ggtitle(lr.results.sorted[1,1]) +print(p) + abline(h=.5, lty=2) +# code for accuracy +# vector of logistic output probabilities for this model (glm/eq. above) +# prob = .5 is the threshold for prediction class = 1 vs 0 +predicted.class <- as.integer(predicted.probs >=.5) # predicted class +pheno.factor.relevel # true class +# vector of True (correctly predicted) and False (wrongly predicted) +correct.classified <- predicted.class == pheno.factor.relevel +# sum correct.classified +table(pheno.factor.relevel,predicted.class) +# accuracy: +correct.acc <- sum(correct.classified == "TRUE") / length(correct.classified) +correct.acc +if (!require("glmnet")) install.packages("glmnet") +library(glmnet) +# alpha=0 means ridge, alpha=1 means lasso +# keep has to do with storing the results from cross-validiation +glmnet.model <- cv.glmnet(t(GxS.covfilter), pheno.factor.relevel, alpha=1, +family="binomial", type.measure="class", keep=T) +plot(glmnet.model) # plot of CV error vs lambda penalty +glmnet.model$lambda.min # lambda that minimizes the CV error +glmnet.model +# Convient way to extract non-zero +if (!require("coefplot")) install.packages("coefplot") +library(coefplot) +extract.coef(glmnet.model) +# get the penalized regression coefficients +glmnet.coeffs <- predict(glmnet.model,type="coefficients", +s=glmnet.model$lambda.min) +top_glmnet <- data.frame(as.matrix(glmnet.coeffs)) %>% +tibble::rownames_to_column(var = "features") %>% +filter(s1!=0, features!="(Intercept)") +top_glmnet_features <- top_glmnet %>% pull(features) +# apply the glmnet model to the data to get class predictions +glmnet.predicted <- predict(glmnet.model, +s=glmnet.model$lambda.min, # lambda to use +type="class", # classify +newx=t(GxS.covfilter)) # apply to original +glmnet.accuracy <- mean(factor(glmnet.predicted)==pheno.factor.relevel) +glmnet.accuracy +# get the coefficients that are not zero (these are the selected variables) +glmnet.nonzero.coeffs <- #glmnet.coeffs@Dimnames[[1]][which(glmnet.coeffs!=0)] +glmnet.nonzero.coeffs +# get the coefficients that are not zero (these are the selected variables) +glmnet.nonzero.coeffs <- #glmnet.coeffs@Dimnames[[1]][which(glmnet.coeffs!=0)] +glmnet.nonzero.coeffs +# get the coefficients that are not zero (these are the selected variables) +glmnet.nonzero.coeffs <- glmnet.coeffs@Dimnames[[1]][which(glmnet.coeffs!=0)] +glmnet.nonzero.coeffs +#### Part E: NPDR +if (!require("devtools")) install.packages("devtools") +library(devtools) +install_github("insilico/npdr") +library(npdr) #https://github.com/insilico/npdr +SxG.dat <- t(GxS.covfilter) +npdr.MDD.results <- npdr(pheno.factor, SxG.dat, +regression.type="binomial", +attr.diff.type="numeric-abs", +nbd.metric = "manhattan", +knn=knnSURF.balanced(pheno.factor, sd.frac = .5), +dopar.nn = F, dopar.reg=F, +padj.method="bonferroni", verbose = T) +colnames(npdr.MDD.results) # column names of the output +library(dplyr) +# gets genes (attributes/att) with FDR-adjusted p-value<.05 +top.p05.npdr <- npdr.MDD.results %>% filter(pval.adj<.05) %>% pull(att) +top.p05.npdr +# grab top 200, remove NA, remove "", get att col +top.npdr <- npdr.MDD.results %>% dplyr::slice(1:200) %>% +filter(att!="") %>% pull(att) +write.table(top.npdr,row.names=F,col.names=F,quote=F) +# grab top 200, remove NA, remove "", get att col +top.npdr <- npdr.MDD.results %>% dplyr::slice(1:200) %>% +filter(att!="") filter(pval.adj<1) %>% pull(att) +# grab top 200, remove NA, remove "", get att col +top.npdr <- npdr.MDD.results %>% dplyr::slice(1:200) %>% +filter(pval.adj<1) %>% pull(att) +write.table(top.npdr,row.names=F,col.names=F,quote=F) diff --git a/.~lock.lab5_expression3.docx# b/.~lock.lab5_expression3.docx# index dbc8386..9390f1d 100644 --- a/.~lock.lab5_expression3.docx# +++ b/.~lock.lab5_expression3.docx# @@ -1 +1 @@ -,noah,NovaArchSys,13.10.2022 18:05,file:///home/noah/.config/libreoffice/4; \ No newline at end of file +,noah,NovaArchSys,13.10.2022 18:27,file:///home/noah/.config/libreoffice/4; \ No newline at end of file diff --git a/Schrick-Noah_CS-6643_Lab-5.R b/Schrick-Noah_CS-6643_Lab-5.R index 48e295c..25e1934 100644 --- a/Schrick-Noah_CS-6643_Lab-5.R +++ b/Schrick-Noah_CS-6643_Lab-5.R @@ -59,7 +59,7 @@ levels(pheno.factor.relevel) levels(pheno.factor.relevel)[levels(pheno.factor.relevel)=="MDD"] <- 1 levels(pheno.factor.relevel)[levels(pheno.factor.relevel)=="HC"] <- 0 -# Fit logistic model of first gene to phenotype data +## Fit logistic model of first gene to phenotype data gene.row <- 2 gene.name <-rownames(GxS.covfilter)[gene.row] gene.expr <- GxS.covfilter[gene.row,] @@ -68,4 +68,183 @@ summary(gene.fit) coeff.mat <- coef(summary(gene.fit)) b0 <- coeff.mat[1,1] b1 <- coeff.mat[2,1] -b1.pval <- coeff.mat[2,4] \ No newline at end of file +b1.pval <- coeff.mat[2,4] + +## GLM +modelfn <- function(x){1/(1+exp(-(b0+b1*x)))} +g.min <- min(gene.expr) +g.max <- max(gene.expr) +curve(modelfn,g.min,g.max) # plot for domain of actual gene’s expression +abline(h=.5, lty=2) +curve(modelfn,3,8) # replot with extended domain to see the S shape + +## Predict Function +predicted.probs <- predict(gene.fit, gene.expr=gene.expr, + type="response") +if (!require("ggplot2")) install.packages("ggplot2") +library(ggplot2) +phenotype <- pheno.factor # just rename for legend +gene.gg.df <- data.frame(expression=gene.expr, + prediction=predicted.probs) +# plot predicted probabilities versus gene expression +p <- ggplot(data=gene.gg.df) +p <- p + geom_point(aes(x=expression, + y=prediction, + color=phenotype, # shape and color + shape=phenotype), # are true phenotype + size=3) +p <- p + geom_hline(yintercept = .5, linetype="dashed", color="red") +print(p) + abline(h=.5, lty=2) + + +#### Part B: Logistic Regression as Classification +# vector of logistic output probabilities for this model (glm/eq. above) +# prob = .5 is the threshold for prediction class = 1 vs 0 +predicted.class <- as.integer(predicted.probs >=.5) # predicted class +pheno.factor.relevel # true class +# vector of True (correctly predicted) and False (wrongly predicted) +correct.classified <- predicted.class == pheno.factor.relevel +# sum correct.classified +table(pheno.factor.relevel,predicted.class) +# accuracy: +correct.acc <- sum(correct.classified == "TRUE") / length(correct.classified) + + +#### Part C: Logistic Regression Ranking of all Genes +# logistic regression function for one gene row +lr.fn <- function(i){ + gene=rownames(GxS.covfilter)[i] + gene.expr <- GxS.covfilter[i,] + gene.fit <- glm(pheno.factor.relevel~gene.expr, + family=binomial) + coeff.mat <- coef(summary(gene.fit)) + b1 <- coeff.mat[2,1] + b1.pval <- coeff.mat[2,4] + coefvec <- gene.fit$estimate # intercept, gene + pvec <- gene.fit$p.value # intercept, gene + c(gene, b1, b1.pval) +} +# Testing with just one row +lr.fn(2) + +# initialize results matrix +num.genes<-nrow(GxS.covfilter) +lr.results.mat <- matrix(0, nrow=nrow(GxS.covfilter), ncol=3) +# for loop the function to all genes +for (i in 1:num.genes){ + lr.results.mat[i,] <- lr.fn(i) +} +lr.results.df <- data.frame(lr.results.mat) +colnames(lr.results.df) <- c("gene", "b1", "p.val") + +# sort results b1 coefficient p-value +library(dplyr) +lr.results.sorted <- lr.results.df %>% + mutate_at("p.val", as.character) %>% + mutate_at("p.val", as.numeric) %>% + arrange(p.val) +lr.results.sorted[1:10,] + +## Scatter plot of model fit and accuracy computation of top gene +# Get top gene data +gene.row <- which(rownames(GxS.covfilter)==lr.results.sorted[1,1]) +gene.expr <- GxS.covfilter[gene.row,] +gene.fit <- glm(pheno.factor.relevel~gene.expr, + family=binomial) +predicted.probs <- predict(gene.fit, gene.expr=gene.expr, type="response") +# Plotting Code +gene.gg.df <- data.frame(expression=gene.expr, + prediction=predicted.probs) +# plot predicted probabilities versus gene expression +p <- ggplot(data=gene.gg.df) +p <- p + geom_point(aes(x=expression, + y=prediction, + color=phenotype, # shape and color + shape=phenotype), # are true phenotype + size=3) +p <- p + geom_hline(yintercept = .5, linetype="dashed", color="red") + + ggtitle(lr.results.sorted[1,1]) +print(p) + abline(h=.5, lty=2) + +# code for accuracy +# vector of logistic output probabilities for this model (glm/eq. above) +# prob = .5 is the threshold for prediction class = 1 vs 0 +predicted.class <- as.integer(predicted.probs >=.5) # predicted class +pheno.factor.relevel # true class +# vector of True (correctly predicted) and False (wrongly predicted) +correct.classified <- predicted.class == pheno.factor.relevel +# sum correct.classified +table(pheno.factor.relevel,predicted.class) +# accuracy: +correct.acc <- sum(correct.classified == "TRUE") / length(correct.classified) +correct.acc + + +#### Part D: Statistical Learning +if (!require("glmnet")) install.packages("glmnet") +library(glmnet) +# alpha=0 means ridge, alpha=1 means lasso +# keep has to do with storing the results from cross-validiation +glmnet.model <- cv.glmnet(t(GxS.covfilter), pheno.factor.relevel, alpha=1, + family="binomial", type.measure="class", keep=T) +plot(glmnet.model) # plot of CV error vs lambda penalty +glmnet.model$lambda.min # lambda that minimizes the CV error + +# Easy way to extract non-zero coeffs +if (!require("coefplot")) install.packages("coefplot") +library(coefplot) +extract.coef(glmnet.model) + +# get the penalized regression coefficients +glmnet.coeffs <- predict(glmnet.model,type="coefficients", + s=glmnet.model$lambda.min) + +# get the coefficients that are not zero (these are the selected variables) +glmnet.nonzero.coeffs <- glmnet.coeffs@Dimnames[[1]][which(glmnet.coeffs!=0)] +glmnet.nonzero.coeffs + +top_glmnet <- data.frame(as.matrix(glmnet.coeffs)) %>% + tibble::rownames_to_column(var = "features") %>% + filter(s1!=0, features!="(Intercept)") + +top_glmnet_features <- top_glmnet %>% pull(features) + +# apply the glmnet model to the data to get class predictions +glmnet.predicted <- predict(glmnet.model, + s=glmnet.model$lambda.min, # lambda to use + type="class", # classify + newx=t(GxS.covfilter)) # apply to original +glmnet.accuracy <- mean(factor(glmnet.predicted)==pheno.factor.relevel) +glmnet.accuracy +table(glmnet.predicted,pheno.factor.relevel) # confusion matrix + +#### Part E: NPDR +if (!require("devtools")) install.packages("devtools") +library(devtools) +install_github("insilico/npdr") + +library(npdr) #https://github.com/insilico/npdr +SxG.dat <- t(GxS.covfilter) +npdr.MDD.results <- npdr(pheno.factor, SxG.dat, + regression.type="binomial", + attr.diff.type="numeric-abs", + nbd.metric = "manhattan", + knn=knnSURF.balanced(pheno.factor, sd.frac = .5), + dopar.nn = F, dopar.reg=F, + padj.method="bonferroni", verbose = T) +colnames(npdr.MDD.results) # column names of the output + +library(dplyr) +# gets genes (attributes/att) with FDR-adjusted p-value<.05 +top.p05.npdr <- npdr.MDD.results %>% filter(pval.adj<.05) %>% pull(att) +top.p05.npdr + + +# grab top 200, remove NA, remove "", get att col +top.npdr <- npdr.MDD.results %>% dplyr::slice(1:200) %>% + filter(pval.adj<1) %>% pull(att) + +write.table(top.npdr,row.names=F,col.names=F,quote=F) + + + diff --git a/lab5_expression3.docx b/lab5_expression3.docx index c8805b6cf5144947bcac1f69dc116a924b897c14..15df232ee60f93fa8d0db42d1caa600aa0e5fde5 100644 GIT binary patch delta 15366 zcmaL8WmH^2vo4IgySux)1P{U81_f8?CF#Bh9Spiec+;#adTPHL~IJkQ%yWy)q;^bk~{&Kn9CGI%4NyK3}L6<;rDyx z4dI{iR(N?O17mon(`&_8A%HS;s?(z!JeAaf29|Y=JO){da2Y?pNkI~k#tUJ-E>JD? z8V0Cxz_O)9-W@q7$p)h8CHIQl23M{2IylT=xw1s zmI3-Y+%L%SG`8a7T;G_BviLeu10+G3Lz_L3+&W4{Sb-6E?hz#}M> zLJD|6LO`5AK|p-m6U6^>Q<8U}Xn;MAp&Td?7k)@#zO0b4j1`jQm2=ATd^Qd1m7^gB zToPEc`Akk-s(k%Eg3E8`Qwi7|zxf)yARSt8TFGJaFM_|;?d@bgx9?cwKIPVGrZl8M zxpqPCR*#SdE~*DNk;EM$^zH>(RO(Zj>%yY)1}BGWD68{>g>>?gL~;Z|t$-sdaA7ho z1jPtdc+!abfcQ^9iYm7Iu4(nUK(Tsup?z-{{OQt6eTwP;3(}NXK7x2Pz-NWi5dKf! zS=@V%cI1j|Crdp}DGlC=w$B&PSwm@}LcE*y3p&Zf>u=NPHfXkw{@kFB7tmyjHg>Qb za>{Dxs%)V9(a%sJthds>3;-{MMMEjoB4)mBhk}bJT~gRMi0nJcN{h$au&`(!Tta4O zxYW?sw8>-Y$XWm*{c|1B2fntx3z}LHG-^v;WC0yX@dX9c!aYv!1D?4eJ~9TxC(d03 z%!HdsFCWo{M!t7`uy z%fI@7ZkW`6*B!FC{?+UAH}}5$ZyRWu5DK=riU|^s|6yrH8HZd!|F56$o8MO<#o+!~ z#djc!{*i;64k7WIUw0u(|7p!7`;deGZDt3a-i9?47D%t}6KYkn=ryUVvQ_St|H+6+ zK#K~a@w{1}*;mL@xZM@ISdA7q49I(kJWD(>uKOJHvozKsHqE{@JAo++9|OdW9S!qs z_+_=ZIZeJQ z_+ein&X5Hkkl$4#7Vc$%_;8;-r{uLyiSXWJO1 z*Y<+-QmD;0_JhFct4*Yp+3MBN^qy$+V%Vn}gsH17=nGXTxvtQzs>k`8{yYhwvlW-C zo43-$K@aUGE3&=qczDF?`mSwwqn89L1}wDIz`(TY1CB>$<;b5-sQIy_vM~DSQ*_aQGIw^ic8c`;w`gL#?ri#TGo*w!gV&OD zJQU#2zX`JAp-RP1SNwBmjVzz*I=W&Q>@QrpBGVhB;9SJ;hm%L?~b`E$VuxI1=VcgfQ%l z?OpHCqcrGrJ=a~>#Et>uyEi<{0j=?j-yYI^{V#Cdan(j3m1=-T(9BR!A~ zJ!nDYx|DWvHzne8qN8t`RQxU|*foqb$wt@VG8p0$H8g0_NcFabm`!YKjv``mbAPPJ z6*uMcJr?`x2m!?R84{m9H`60?9Wrt(45))Pz~Uu7lLVl`tK*!9ImJ z1r|I-l69S7lB_*y#9C6$B=6?e@e7bdXi)3V%g>bH!p zUNm5*jpd+A{*Dh=h5YT1^de3ZL@G{66I7#8L@vof2^Yi-eDgSg0p4YtbIW*=6*Fhr zF$M(b|DM3QM*DdY|3=^ccDzLGy_$r|{mNu^>>F5Mp(bLFzCRyhQ98-UgCLUkvW)Wf zH)S z#vl-68{FqpF5gyp)8m5oWUf%y!FgNi1)4FunUvCO=7Vv;+h!6N`2dkskx{~3bR`VxvS)n z6(SPIFl5roU-a8$zI(}UE-WN-y0qu_F!RT@G?s{wP16zYE(j%K?lLtaSz)%nFKNhJ zm!4LVfT%FyZigW=0ooO;4geMTW}8tY+--?*2^bBCf|i*wd1X|9)u~%%PD$%27wu*R zKq)i19vd&``5a_u%W{E`=d14UUAGx;4wdq3Kn&D@O41q--M>W1xzoB45I{@pnWj4> z?Q1VDzi&S&Z45~h~SzV62fh|?TFGeY^5Mebet-?^<52gl#I zIScuKU4!{<3bXm={%_G;nRV7)No`d=Ue&*b%LJw0miciWYD)qMlpIV5hWx0XK>KKa zot9+0D%Stv5vcjnq+dp24}sl!XI)I==y%F$c1?332eZCD2@>@vXEZEyS^l_}B%+yZ z;$slpebz`sG|*F~pMZcC3iDF1MAl&DCNBR9{Hcz#j}(($>B6on!Cf=zWpC5Tih4m0 zrVkH0cqAK7WI4>Fyf!ek50JWxWl3xo0RM3>4e*9~r7bqC*L%bSv3dS_3X4Sh6oAR> zg@-D#i>J2v{j*RjDMO|kX#olQSf%+Z`&cN~6L;fU*Y_5MAW}KX2WoXIVekwpGq+qm zu(_U(-+!lk>8*Nkd2t*Wix>|w8y|)n`Al^6SQKhy+TCI#g?8g-`h{_+;W>nVzq{T2 zO9oBV@UPW6Ba_AP+veEwvABk@Qu|n;zM-Sj%4bkD`HQw*^Q~xnITHM`Qq$+|oBh6h zm9pWOGlX4q7MDGe11hP(Ln^krM+~SEAkH3CjKxv}Y{_5BfwQm%5fpK*5Npzf47g0Q zaO})+Md)e8lV`!%Nu)%ym3|lhd_06^4z9&mWw}R_^uPWr)aS)aEF{XSZs^EJts%RA z@m|AqKwV+%qI4qA*3N%kiO>#1tsn)e@@RkmN;?&lViFDRhAKIq5``?mK*=;)0=g(k zhv#!VrXNQT*bG3;PMgo3hw2&$QS!@jY-=fz@u$tVtZ#DPi`&g;#P{f>Xv*kEPQs$n ziB2*l0n-{Ly&lI#dH@NcY7o^Wjk;4W>P#-k_&bxy=zl`^2!q%2RPBbRI9xA9Q|<{`PkD~9d_8;`WoQHY zNLP=pv{*et+(Y^^PhC(2GAT4X)RPRPzjs01mim?y!G` z_aHy20Sikt@|U4^r4c?jBCo9jPB)_OGO`5|l$bY`CX)gMX6|$*)t}D12qsHVr_g&# z{MneuqEm^sa~&sUd(cb(?bEH8Jc6!$5+zN#)DSuz?28nMh~MU#VKX5V&|wsRYOBu& zX+^*0JnmX96dyQqHmaBRiKJOw}e8almTnbKU$ud}-ln1JnoG-zWs#`^Z z^=jKOFd1WqohHiFynXhuBK7l!ucG%ZLSB+VK+PvPA!PL%8Jh0N`)O0Q_qsn8^^A{H)!S#LL_F~?^MoXnB?Qsc9!1ABmk?XGST$ygs zHt9?A28XV%`{|_R#om)j!{q~h%uib_Iq~;)R>BQ6X6>GNz9yiT(7!NpH+#>+q+*&> zF5Pp7JjrymzmHpHqdP9sN;vtW5NONjvU3#Ew=o|&d1hYHNMzXdWmflgkC>B8lus^Q z-xplmGZvfM=d@}Tfn5V8HV!H>fV_wb)lk@)tE31Yg&}s#7wgK-v5~w_gujd=+486* z6+_n&1u)kHaX(ciNfTr(dh-~D6xIxw14wUi61$F5MyF69&5|r4j znh*$$PrVhxxdrT+V~$adfpZ3N+nd*hTF!mGvRcbjAwpcJOht3f_}D#)1oM}S0g=K( zOFf9Bs|?PjhPIzzxDedKN&IqmuTEZN_pZ^7eFz1zpR2oGPw@udO>P~fmGVY(nDZtc z8Lx0U7SdmUK#b0MQepo_`9gBK{XvQgvq=KL51F_{iMB%obbZiRouASNu6|S|# zX=C5lQq^Ba_BO6+S~VP(8qn`HfOsL(mx#j>Xr}pCK(geWcYFGDgVT)*fi(8qr_OO? z`ktj~jbm}^Va3CUN9?F`BtK+i;x9YHBuh}PFHaqsVNcQbyamH;(a_dNPkVC0ltRed zk){3EkL}L>tRnGnNuzxtK9mK`7nL)FxqPdsO3opiTw;xYWL|cm93ZSG1D{j9pR?n3 zYho2EFvC{-9r6cAK&z?u^JaiHb%1AL{O@OfxzM+Vdvv4UJ@shb2j=2OHF;LbJKSHN zYTOQ~_3GUL_69>5&ah{MX^P<+LUnz#b74+{K(!mfJ9dAn44*tgyf-e*_h$#btk>ph z8^SwU7jot_%>2n}uB!&p=DMaxj?Q?D;^wGJ;QfPTiT|(RFf1B+0U>(TDXi;+mSXIv z=<+|hDK~B0RGP@eV4UOy@xXU<@ohm^`8p{DOd*&&fed+#vb`w6HQrcnSu_5CA1M`9 z19+ptJDf9+VaXJp#FvclbA?q93L=VZFDOys`G1UB3l~U|A69B;bi-M2jKp)aZz22r zfl}kmTJmAfxZSVAyf^5qtB{+|6?|FliM@mmd%7`8JF-llL7u6mbvol2f@A=L&$;_ zHFk3%4CGNK{es^OuGUFkCS!G~c8NLQfkN?J=e`N^B3pau`u7qk2fNmL^;=R!M(m}| zo(xYMZK+S=`+t2XbuOFv?@NXC##;boxX5K1J&&a&zGsE+#KoRP6DOFo&eVh^E91{{ojqKM#GW4D9_XgxBq*Top`9U?~x z=OK%5b+m4MF}w*>nfsViAL z4AgCAkk=L5jK|!4mx+OCr@f%9R#0O$cVP(z#Y2#7G<4;A{U}4?UBGaW8)V4$VtptMR6%Eza1#oiV;Wzfy;65yzn|v$4t%p~cuB?M(7Tc#; zr`Vu5fVp<@ddIc=T_}GW3jCEA-!u5ttcG)mz#fGUjyUxHM?GNWvlNf#XENa2_)*N+ zBtD!R2pk(ie`|s*jJGYjZdPL9W3p&8stS7^X!I>y+v2aI6|P%im^^>NpE2Fl(zE|n zxDlvlpMuB2m%E6z`#Z3a9_UG6JVFut|IukWzP(}r+mwTH*NSQ66b1 zyc&vvGrTDm+olK2*#@>}JaGF;cSeq~RpzAy^K8c4@(pnBfg-2wBqQ#P7w~mRo!PTw z)%}oHlf?Z#&%c6@&t1r4b9vA1;7Q|{TX5SKX@Bl5`8kwWMiATA1a2_<0BBFsOZFV| z`RIGsCEJm%cx!?QJmKiexP#jQsoUqV&XRr0lU|)4ENi@0Cs15RZ`fxOoTbB^gj4-) zTfwRVcsMh40|w@+o%Fj|W(GY4+)Wb)0;@4wmm55o1!RI=7mQEdlL4=Q_O{C3SHyN{ z&Y*2v3}+uEF3d@ji#P@bGN29CHhr47l0CyQIFG~7=y2}`krQutq6^($qT!34{<(m@ zh85=6NRZ9nz+C6@Ko9zCqJ4E-1SU~a?bv~pgn-$!gUR8P%doi$q3wnQ-)E2BvX?sg z>_5jk)em8F%0m@*F=MJUfcq)mTb%b_c635>rMX+Ai{`O^M|^NJ!9cTJ>*Scf(%g() zWQ3UrGMT?EktoUK07sJ{Ai7abf3CX>HWpqKzkU(6Y*bq;@vS0BQ*k*KSqJA+I~bSi zZc1t~*Y9>GgznZ{#u1XR(G!teL&&A9p*@NYIWEP8O=yCT2?HUBfEh}1#vy_$n z8ix^AJzVg#HYuTB{{qp-3oVjJFp121@1jv2nbIm6!o1v^zf55dZ>6iMAV`>TbKPf` zDR?WDx&!j|oj*k_2A+*RsA+(GDy{T(_0+Y=aEqvZ0skP{)7KITEwT(9kHS`+vxk45 z=x&-=- z;K|`tCE7FnS|n;36P7ZKUD;djzJ-LdHisEd;74Lg&(+A^74?9-=W@$-2xgE~yj+)G2h!Gqy$LnR=VEYz`wCb`$u_I_y;- ztcb*r)SMn3bq3r@S~~~LDi$}D2H@X_$FV&ew5O=MiKqsTw=e(52eu>n?}KLG^xN1t z0wuS%+VV12z;PW9Qx#ZYwkzO5B65PQBGG-JD)?SYTyxIiHW!CaiWS?W8ShMr5U+k5 zjgn+~$uouhOZ6~uJH|IQEr=RdN8z4w=lc_UGei~(X+Tu-_ONp!FwSfeDHY0}Nf~k6 z4yBPLm)o%38`hDAPq5eN)3ZEc8yZ|cM`-p2GctwTV=?nQTO_A;8aU2KE(#;~5I0qh zc&Sr+`MJd)a<^je^MWn%NOZTB2ac+c*SJtrgrl!_)XJN%@-InSEqejkUSZzaz4S0Y z;8^?dw}H{phu@ftEbaAu?klK7XiTD7-E*jK zT>fNgy_bh%tRY1jr~86s+-fEh0Lxc4I6l@Adw4R~W1RIxX3!|Uk?M+b)?9%{KHFZc zr5vfAE65<7nPS!;RlFW@8)S(O7yKL&kgOkc_md>a5Ygd|@gQWLBvEnu^fe>1CE~{~j zh15Y$0@qZC3Xi4D;UaOoS_>{O-W6lMkqEkx`2A^0=C`IFnmDeojfPCYIw_li_q;c~ zTGDTP;96W)IKB6)XW*|X1AA4OELL~5XE-3i`{S$l{Kw}Nd#iC0wdvvU;Fj;06gl4K zgNDv{nRNV95oS%q)mjaX6Fk4o9sNJK{4Q3nZ-yOvUM;X9+FCa=TEHR%*-VrnPECf^ zz+}Azr%%->bOk!EDnm&M-$ghxMA7J&#g+yyU+@s2#wchcaKkupmwQvLd&XXn=Gh5OE%5O{iAyF(mC#TgOX zcPN8sfN)VUqbDP|#gLovdDQ_6(Qv%gFym+YMA+3uS4*l&C`D$I1UgUZ(lzD@?-xI6 z;^mPb5QOPWH@Bmk`V_tfHq$zWU;;2J%0|V)c6$i!_u*YurbN0Qc$8K5*IK)BaC^R! zOY6b<7wYl!o8afL87gOfD@OqBD%x3bN!*#Labou*!7@<7->Y@&?aYye++lQSF911# z^f4oOHMl@<5#X}DB6lBqMP1jQWfskwaM8WTn?mfCeLXPz8Z}7ol&EET`Uf~Ho@!IM z=Ku7NAuhf5J;!xKp`9%dXECnSBHe2u-$Xb~LWxJqJM_Dh+miLq%VQe`p)Ijp6;_2u zakPOiczV02b%Ux(mtk8;>;u^`-_Vo7`s!9{#dpOr1Jm=Isr?~;0 z<(e8Z)Ek>bHtlOc5Z61!`*>V)7m3K zo(!R%6t#%rB^!+@wV38*TeaiWN4sEUDpQz>e9OugnC4B3=*X#+)MxezRSkZ!Q8zFx zC4vsWU?pp9lSY@gZpJRNASG#9P-?di+1@0%y2Jr0@0$UO9a)Qy+vDoy1q4KDlHymS zD!yQ32SWa9T&;l7adY6uE(?YjH3_IRyvFil&(dhRITtbO3Jc-Wfs3-=pSvHAyb~Bk z`D8I$@~3KjeBGiYR|aciQlVCRQ|F)TRaj-B*lV*L{U#?uqA8#{8;d}nMW*bjXxOcg z1fhocb-JfEtrM16C`p^>G%QP~|7dVwoU&#FyZq#a6PVRRHhf<816W6{J{mR*rCWLlPf)zwj$W(Z}Bc?akhtH!X=dnWIB`rebe1UkJz^ zCwe3Z6*%alQi>Y&fD%sfY*d@_L3OLk{+Mh zD#prxs zr;P2;!_za=Iq|{3V5feB#mDsIct{QXi+`uex_6m9Sy89ADjD79ztg<50K}tDS5wQ% z+|Znfv%INWAeS)a#uJIw&XNPINVxQh*cm3VLp;s~1K)tLqfKL7Q(+UT>*D5rB+RXb z{V{w87FjiN!xGz!SzKoXZOTrU_Qb6A2lQ)oe`(s8-$;FIn3j=2W1XJuxQPv{-r4vx+_x{CKX?aRuq~@9=k-Wc+yWKoU;Kqm?8UBeo`yASQclB78_V&C2JX0uw9B;Ym|VvWaLO;Rn@Jj{Vh?ROC(K*a z8879P%YL|8GM1jU6qVyJ?W?28{lLMoJ;QFk%-a}B2L2=r{-VXEUT)jEz_2h{(Gu*1t>C^h2Pby*dvUMdifjp{dzoCFb@`ytnH3#ksK9sJ5Xu+4#qp(i%0O zH{8!5w&SbW8N%5`e@^qhHzAn-hwN(;ymYdiuxArv)V#rbjGrZ@GxAx{e~8vvS-|81 zm4Ko9Wim9UG3@eJD!sp%>roo|tcLWWVvtA_9OhY5?zup;Hsi0u^3*ViDD+;>{qauI z2n0QTO;u5_rop-tcu}x_B!W-n3goJQS>k3=PG}y98_ulXjoTu@3v_@;qPl34Q&KrQ zG90uh$HyKme-9kvoFSb6mwn8Y;KGkHqvTMA<)N8jtDZ>PGAs*hv5c|1`p=V& z(3hD%u10nKEk-=+c{U>SWmv%{Y6?r%{S8~X#i@@lJI(Zg=vEHDw|NO)go??vgk5tP)mJR&=1MY>SJ9q60Z>Xj4QwrO z!6moO+1{>`jK_U6>-Xq3H(#+}{+Nhe*kN0$8M^>1Z@j3#4hL=e?)nd_rDd?=Ap66ufEGd}k!V>> zvg$Fm;?X~P$7x$w2}}j%1^n?1bnYdI?S`%9s!kftHOhN7Sudb5+4tuk41ULF#8OP@ z0Fkg16@M$M<#-j6{ROXc+_H4dU{x8!xR>J6WBw>DK<&9hs`I5EP8yK7XuNVsXoDRIEBWZONp`m>A|wTXh-cy9$~dy4YDC3M&;Bb8uQ-<)>c9+;d9O zn>U`1^`8wMEitBnpVFJvle67qYkbOYan*)#A2YnZqM~^w4on_{t=_P0yL~(p zTxY*xS|+B1IdTp-R*l(XI!GB?>|`w}jlM?5_+;VJLpp9uKkn}1h?iH7)rTFS07e_;y+aWK@|fr|8K7dX=Tk~T3igEE zq6V5$blwHR9 zE5`&56*ZbEle3H$Y@KLTDp*-q^=RIL8;V4BVqLpy}3Ib9+==GoWs(hTig$MtHIe){R>c zyZXv8EZRQ5--)W5-?PD~3WFzHZ686z2OXgYPiSRu8r8*d@EsAiuzHdgCbR-IaLfBg zJnatrsKr&7>sN$xh~t_*G?pPM8QTOCo=7_xu?ZLjU~P@|a#PpNpWpnr!2U4E_MQaX z@Ko{VVage`V2=JT-d`zP$58K2kjMvL%J7#fV{OcC$--Q2knSc~XP}h#Ex)Vd;JhD| zIzDroTFqIy{k7LR@rmyg07s}V$%> zYsFn2JQRI%Dy5J7>7hOShqf@|Gq2T)QA_;Z#}TQPIC`-w$}=OXGtSMk?tJc`SK z74m(-bDxv3*5R_GS?>nuDD?RZ(^8GWzHB55xfGu*P@J{w)VPt>qe#6Pk0`a~W25b7 z)%fy`_Vo9^T6*f$|2EpxbsmPk3(Zt-B+%a7%cUNzll=|ybLY(+VS=rd{Dr~Xk5nrO z{NM5yMg?m%xY9`3xb^*pgeGJMOWNQn0@|xdr z#L(Vi_*~@2!fnc<)9#r^p!jW_yz3pmHQ(zymJ{mrl3F__Rff3F*<#st?erEwI_5o* zRPydwK*R|0@dis+60KL*(5ehDB)A=#RTU!IhZ)5BscuMrXahCb3AN&}l+JYvIjkk{ zC+r@#&UO7|4%Wh9&lX+iu9k$;3U?>U6TMdTcy>>5%itP?0ysvU3Gc+VbAj=VsTF25$qO0_#UP9&&hwE(RjTn&uW+-oqCf zNXJgy3eD40!WI~Hik9Q<*X`XOhA}-_$aj*R;B=QN|HO}x@FS^`gcz{#l<(_{l1=w6 zA~-TR??zqcH>t9wmfry-S)us-<3Df@1#Yu3cK|KM6x*UNmATVyGEgx5q%A`r&L#@H zFGa3p8`d;%JfDOALF*RJ>3a5)jT_e(=)1#l{!xmfwS10mfYmC~sR;klsTR|ddV%>d z&U<~3cz-g2mYP#_qcU3)69szd`jV!8=y#oawA^whWkd ze%HlC{s_u(Hz^=H1uNTti5KkJnudY=4WdU0*0GD%nikKod_cBR{TFRF?tj8z6c7W# z`x>FWAEVo;rV{k8m$L3xcQB< z4>xgJ@`NirO+PVf5vRuJE3>4*zLTF{_-)Ev_9{!~ z>$89@=Zv4o{;~e&(|a6_M=ZWnbH^d~k~bu)nC0IS{_$k%1gGj+-@DGPo-MUlHaSmw zL=}kRru2T{=NSj?^Ws(A4IJZ>Ln=WZ9nA3LwlZOu%I0FoFZqqaZ@ZcNBoIpy-*&Sx z0eMriQ=ir>Bm2<)%KZLBAs3KCq7^TI#8P=cK<~>o9`@a7&%=7Bu3s(RV{q2ZL_3J2 zk+}y$0gi)1S4M7$*!T|Mc*ZuTV=rIzjix4p?}6aBMM(%a_{zDaZ*U@Ks$vam*zY+Q zb7)E%IB=sPm*eC4NysrNVW34k=57!7F}7$Lcgh6QNlQGHS#|7{2Ox@QOo8bF6{+hG zqosLy57RCBI~Uzp`tDg(^~APF7l_%6U!h@oWfQYWzBtE!bo+J}OfMxFuSenE-ChwQ zK$lY@VyFwG*h7~AejoEt6Yxzuj^HDoqk)&ODEP{SxFgMJ(Nx$*vh|Zj zf9*D~CWAn(&t)bj(~|B>_*v}P^Yz(tg6bmgg*lJr$S#ZoOAJ09!;|~Ltx(#hM z0(_Guxh58ZVOhtgcp$aiz@(yGQLUfa%{v>k`v59%GqhS%Qr5RDG62nId#j$BZ)-y_ zzzvz&x{Gf`--@vGBPSdkvm_4oJbzeEHfBEn!|guTO736`aSui)DedK){~#NE7m)7M zWMGj!ACIyP-nGExbEEHOg3~Ha1l{Ahn2xtV4WB-v2%&eaAVdZgku~c zQwm>>1FUWk5~i)Tc>MBGr|7+;At19eNjKk4WE-!+AN(Dza!}W)Bd9v;W5@CM**)?J zXlh2r??zf|rGg za?v|zov9NR;NCAz@sLJiTHO=YwN%~5&&7CH4p_iud# zb18kh7CSUBz7lle4j|4sJ`U~U=vvw)$3rZOT~$1B|083qVSX@{A$xa&>StdEl#w*# z+5o)aL?p=h=HArrjVElo<^Q~z6y$cC0Ne|`+W##N5h$!MPZK=QNcd_hFGa^j z-`b5*VHv4)6KxXf$dQ>@jG-dEhrT{Cok7}d?j?~HF~iLWG)eRZ1|MJ;jv7JpHI@RD zQ)kDwTm6zS#E`u!{Ua3_G?(vy%bH$M#H00?axsP2&QU1wnf?5OO#&6&$hry%GIjhZ z^g%v!)Z6mW;*E$Ha%+#+Y;$HKT(OreM!ItGyM-YTlXPG5uWucjtH{q*OJE{A}WP}{52aFaiX>s zeUv%$*HE($p=FK&DnF1GXI?LzZuN70w!NzQdP-#`KL6xcW{AS(0W$$^65S7IRj}3{ zg7w*P`JW4PnHN7w&C0Zo5>STE_gtXvXuouAtG8aW->B}=m$n$Yk7_W9*kUp^FPY!O zID;HN7~!BZFKf~_2&|CVctI6soFC2N>jG5%vJVVg$jhz6s7VT1xleJZ{<(Wn`zQBD zbrMRvVumVHfu?j>l=@Xga*T`a;pZJu3u*`%(ikBe4~CO_4RECv1dxch$6N?yUi@vE z19f6#wd#+I&{iP3?w0KKi#J@#E{`8SBGjVz_#tfo?nLdi8xTBzx^liH)*`CrPL#v0 z&jTZo;>FM#6rcKWNV7tOXakYRJ%V~of*2h(SabxU_|w9TPyMAA%Gr7;!9tar2TD3n zA`tN@94@%zXwW*v9N*8Db833SEgV-6PVj&;bU{^z1<$!SdMfeSrzgw}zma;(Y>#p; zU-jh0bXDmM2)r!Y)QSLyy*dX-v|f1#%%Gk9)U0Cwe)FQ*m>~uE#yGf59);u(uenr~ zEX&YI&_$>B7yX^$pim+aRSzcpI2W6$jkgwQv)q`zB0MZ)&Yp2Z^qSxdSBW0qVc}y& zBSvJ$mi`$PCz@lkr3~Fqo~Yj?a5I=?Rpk>I3o!7I0Uk!)hZCkd5%t-0zld~=e+lh# zEO`%{`Hi97nIZ{WR>R9q5MxlvKmo#Lj&Vbd&&Q6N9*Z4wgEG7ShOsP2PWVZfX2NLX z4A;-q1XZWQ$C=sb>MHaR6`#Pd$6Qbg+v_AQa4fjM_tmJm4M*Bs%(_2-(zF*t?R6L1 zgE4a85=eP{@I1MxB2JGrR*}4+R?&5PTkg?7bvADLgUQDQ+;&D6qSJ|rHONDP0`}$d z)cM>u=GfTs0orWSEncu__DGH6s}ma5A^3c7RH~!r`?nT2qgN%><)H~sUOn|SDg%AV9S;`MQE(rKEL=2~iQsGBS2CsWJruGWxyr4c)s!G|!E_SzOG;aYdKp=wVdV zpvjHY09-b0^YEwU2gtkvJv9HY&*#sJpLR3`>c-eCwTT=XvX0cw#U_uTin zbMqR?T|O67$3~FA&AI7w1&XpC`XI<@qjZ^rQ_pU1!zMx7iL&s??P^pKDXQO1eepJ9C zowaU`42}KW?~1%TJeUlE0{xU*8D6gI z87>LjFa&}SwUlMAo=#*-t_l2A{-S9<2Ag8^O-)=VEPb#ck)q>K4tp;Z$0|3l0dY?L z!1QaJu-wn7ud%Z`Hp|~W*+dN1Aw7x3RAN|=?O&g-h%)z~3><*3fqKB>*UP)AwO zJ(_V-)mSUeckl7fuj-hH-th4S7{i+_28~vsQydEf549$ec6V_>jI-P(G;yxEZiNVL zsh0dhazT$o>2bX623zGrkoB-g9ZkcN{`U`564?)4-ulaXa($Obg+G8L#POcmAarIH zYOrA(k5j^Q)VyUJLmcY2DRmS4Omhm~vT5i!+%Ps=_uk*|wF=?uC532IAS6K>#)z#V zC>0_wrFU63bi33)M7|4M0SHOdq(5e0&o67sR7_V2Iphtlp5G28C3Ubt@8=u7N z(A+{i_3vlH_1r-LLYsn+gneeAdlZJeohjfxddLyr3;Wleq)PH^)&!qyLh^Ce!$6=f z;-4}1R_ancm>x}<%CoE{HX)wM|9JZ}o7(pnZjD-+jns8p&q10U*Y{dcI0S?w;&Mw~ zs}PyvOcJzJL~KX>jud6g^dhN(6J^A7E15|8A%M4%Q*{4SDs&D6W2VQj3X*;{QqKTq zhOn(7XqCT)lDF#l_i5gpLugdlcwz!@5*q*w?*@TeW0f)Y+psU+ZUmuurcd?w)K~tg7E30P-Sm&xwBbD#$LJqxZ`Y*_+dnQd8PI|u0bw@ zrSLy)NMZN|wF5#DNB9&xU3A(?^#wizgw_%06jXJi&Q=02#BuRWLTj&3fxeyo=IQG) zbJFbcOm%)g;8Oj5bwYv*Om)9}-QBsspC6w;VYk@%TkU4Hi_I_8a@kyJ+*E*bJS&ipQVlR z{&R5%=(LTP=01Iw$i`%?pZ0sTX|?1WE(GN0>9YLy&=XUS9Z2e1^ECXPpP8+$^!6zA zPS6?llhj3Jk{z=>C3*FJ2+!7+v4lX*p)pjMXtqlJb*hK+4`h*HMAu8zCL^_PnJ^n0 z>d&yMbs5?x#$8ynCdNfLeYx@&0|G`JIJnI`!ykEnQruSW=OF_6OJJO)#ov@=D)zt~ zUV#a3>~n0$*+7rj(u1P&rqHSw&Qy(ZCeS*T=ZgE%E2$`DkrB24X*PxXv9EpgeBVHd z`%-{NUz7`DzTE@x)7~tlPS~`@u-EZ)d)rylq%Q%!0Bc}pAUbac*iH?!ZVS6upnm_6 zzZJ?uPNu87 z$IpVov0I2hrlz(N!OY$gsf+dnEnKhmIg`bU9MEN$>nP{$btDxT<>$zCi$jeEgvi+2 zF{qZ^^rpoo-y*2LGRp8q?d!c;O^mhAhl$W4P2W9Pt2a*ye6z9IV|S4o^`(ZO7-fu3 zO@o9tY=ASu;=e^CNCinf4L6)>7mNF<<~3Amri%OiJlF3$o0e${Y2b!&kcQ!J-f>PE ztgrEH?&dY*+>lIdpQE4 zCLBY_NF(Pxhe9_Wv$76?^;yzj-IN|ZiYBZ3V1W4GT@@uhq{rAg#a)T1^!0=iiSCYr z;FnxA(cg$9UTOAw0`=xn^av!({)BiQf{j(?T5@<(Mpl`a&|EuPw<7eq!2#V_U`~+@ zKwoX)L$MccO0|ca@D<9DToAOmaDfJQ(3kXx318VC5?L0DvvgqM7ELo%B}Bh4Gcu;i zO~z>C+(WRb=jHZ0^uM$@8$^!pOneX!w-V4FX>^cKI1vAxN($1zhvG-{Y**6yPhuk| z86S#|{Ga+iGwm?`lU(=T#UkhkABy~61{@**6bYmYh>ie?>Yw`MpDGWsB7kE44@G1+ zSieIJ0)ojH0s{BHTmMI9f@%q%X#NA*EgLJckb;2dQ2u{G8KBz_&_7}#!Vh8z5I^C6 zsfHga6lgQKb`u}%|MBL>KV$%?^8@)Gb`3aV>@fJ)#YAh7>0P!8yw2#V@Ip#OPm{%7mYHHEnVB(WhQ!QvOffSvv-6X4&#n6J)vvE9 zNwd0l?e1PZt=Y5H>$V^6qy-L9MII6g6ATUx4y?4T9+3nbqO>iU37i$^5P@rBLTTfz zijF~wzX|XPI$7j8-#^`4MdDVCa(c{l^SzM1boQgYv4NdB{TTyy$%-DNYFRqzE0Olq zc8w-c+Lk>71S7Qvmm4#eRnQ=oikxd``@H~$1e#(utbwB8b5pQd9k+_iCj2Li5zXUI zQDh!=fRdLKyL@`>|yHSYoCBW%oA`Vm- z^mn;ml3=SRQ8Svd5e&fVUj?muZnlf7-##!Ey0Vv&uq%bRLB!VEP+Is6qio&=-!ng=IUaro%1-D`bRulv zRMU97vYD8xXnTH}Tc9dkH!^IVtZI3Gc!$KywE`rg77Y<9Ib(eAP2VA?$iu?bzWHu? zf`fsbK!Ab$of5eJ&Pnnv1SN3X@dqoa(($lwS_+%Itm=D=l0VN@ZMq6Hh(q}r<&&APG7Eb8N@W>}I{+=x*NroH2a-xH1gX)r~>-6GRG&JlAqH}2?)T#9hS z{ropAy%Jtf?w0C0amzPJz+z#!jH`qmU|bW zRU(|bM9yi`$u8?7w;V_U_@#F8$!9%Tf_oMvdXsypj>T>7I9^EeWVFh$mc);CFKOFE z$0&cvU&l)?*Z)Jl;lJenmwQD`uzz03VNg?l)n_3v8wa}G)% zfq-tD&;g47HB!vh`1}30qH%EZhupvCtyT2krUidMEAj;R3d(;8zsjtwnccyAb~2t39fuql zO6i#lm6dN%Wqlwaa&gts0|A5jlX3SGiLnt8kZ+Cn znK2T64>givgDIwPOXipOf%^5D*LftuMtt7tn#j(3k6q)QGTAfv33j*(A+aX?s5A_U zF4*OZ;EgR}>ui3ZFrVUBrG{doZdAIod3hafM}dqRK^+6|?2Bzj)oP0GAO%fK+hiyh zsh@J{sXEU{L;Zw4C+_WN51f-aumNMX${?{3i%*lKj4wKju@Qi)5(_1ifv{(t;r48Z zY&%lgj&zZ>S```pgOg{oN72q7|6TACT0JW5iklEm?nt;YgXO}dR+Z00%Q(xhh?aRH zVx=xjkN=$&>Q&@};E}8fM(`v>g84SIdpWTe^%GA@?-a%JA$zWOi~&(Buwsd+OZuow z40=-{J)iwO%MaD^l2kn7n-@9#T#B~z_1D%WY5XTjt!eAt@iNX425i&Jiya<)zyZ8CG-@yz66SL2_jE%{V+418|vdk(iq~xoIPV zue*n!Gn}27^rUogIyMAC;1B7ioqxbTEQk|tr;(rxkx~g zp-K3r$Zp=7k){Si!GUjrJr8rbR){S;NI|Gas$SBF>ZNd?@41ZMPRYpVf!hsCn^W&J zFyh)OTJniD;!5Sg0n4divGG-$$PFk39uaS-1D&>~in<%@XAoOWVG`}*#07Z*tsv@& z4ia~&S}*0Gehv2X_>pErNZ2Q6Fzlnq^wot<;QEjUKJY29#en5Pxdbz2fQv?RYx7+F za03pl7=sdgf(<*gDX3j9 zJ@cLcc@}9C!YM)|-r8CLA`^?Qj@sw_GyZP6_BGxlP(3JTrV(L;YAzInN?c7>`6=r>nrC~Gv7WtsddMm2JiX6WrR3c^}O1&^TaN~C48=b@0}2d~@} zr9YWY`1@5ErCnCn*xXc;q3Av~T0LAl`|A#%+1Tvp0JBvE4*Jv$T|M6(ZY*VQN8_Y? zT5)D=)ZE3K&}}umxS@g~PrnPIVp0cTu%}!4DMv?nQ`S<*xeB|U`~;`lxeV+#XqI?E zdo+c(hNp~kNfdKLkM&{^oOSy`2FZ*#M58aSH z`Yr-s_cM~m?tXJ>xA4;mbTsBDQSDLoGRR(L!kO&r%z+O=mE{gp_8z!Pv5~$CmI&}6 zk}BgY75QK<{CoysBgb$3B}(eWtRks|u3kj03%DNw5JN&?{yP5C2I5TJ+(*jfME2a& zaa~V|r}!J_)!>Bt0Pw1YVIDs)a^PF_B1d=3n2}p?Tu*U~d0n;oIge_sU!@*7J4DSF zCF-Y_rnAG;%+sXgX%3z_p5)3f$1cqeIwKH$6Rg674cLgT%8o-x?!>D z5l~TUEV@B+ab)u~e245|9_L+^6lph~C#jdoa|U?C++lpK^<*P~m*@?L60zns_VcI9 z5pwuVMboZX;@6=po9E@rU8RhMMG_8;FKn^YKWn*ZgRTc@&3y7HaBt^S)YTnz|z2b zaZ=z}Cw@MOaiolg)oFMO+wER&JxRmo9(eNAp#9dCuzd;to+zqVdkGXiISfCx+_0v$ z_`8XhHd#bfHW!&97S`2H&grZ2ypTem=MG|`8zA*+>Y_UgCD>n)0fh7gb~zRGIRP`B-|GW47&3)`9{%Kl|w5{efY3J|({FCRk^Rpp+1Es)^tss{8>|KN6VH0NbpK#%* zf+hLbDt%2vdgFwc?0_x>*&8aPHfAhhIMaC^Z=;i3NLBVDlg zQ&6PJ!Y7K2qA?6OW@i|)19<*0lV5~9BLH%Sn6T~P3735t2Al+6Vs%H7y&K?92G{)| z@p8aNTYcrw{dR8s?@^_^6~yU6d(C3tz~-#CH*$5-i}mv2GjbOpww0?O2wGI={n=Lw zTWwUoE5R$#7cQk{Yco*#q)4Nq8kK$i65IveDiy@}JlvjpTz6thOk!eMPmhkle*=0c zABLgnr*%-TC$$@apn3e3(gb9C4&JWou^$;xUTJ%j%AVlb79Kkj7R^|MOnH9&{UU;* za#BSjHhY@E&Y)Z8qLVN1CpzTxIDDR<6Sb1H6hSf`^99%oec@Ew>}vKd zHAx7Yo9AE;8kJ<9#gezKV0}HOy0zs&ce%klD7;BDs7SmxUSuRXLmO&F$RYTwPOS~b zkFbC@osY>D5Wb)cNNIY09370X2)^Q7ANmuRLq1VQ{-1GPGNhL$|&Q? z>hQw7F>qn~bihr&z)S(sPjx@}rhjf+>c-56l-%iiYEUae*a$RVh>0|&s@5~q9VHk( zSL-kqST7hPT5W!vckA3^sTSbPq}eRFxT;5Qhs=erDOE5F8y4S}F!CL5)Wt0eRnthG zQVO4lm3cDAunCatQ8FAGZ>bnQ)@JW-__@9C!E_imc3Sm*Oy992GH7b+puEtwfuvml zBX813fz+55mj_gPYHofIrU9Z`Hu=i7`Gv#bA&rR~MQv?V1Oj5ofArmg=(~RBs8%Vw z%zJgqwjRQ)c~z9*$a?zp%t;%Y(y2D#oyG;G(id*erA%J;v<0?=f01^6oFJS+e^l05 z{@H(AyVZ=nRRi$$f~!OOQALHjB*GAL@9w~ZDE=k&coSIj-C&bfbHHqb-LB4MP%G^V zL8}l?H8j&G7wp!_ooGw}WKEC{1m_`QfI?$J%hb?_$L;&GrK@m|!BE_1BPq1_eBWO5 z<}UH6B68xX_t!e$=J%&f`OW(ZnsndhyXBNTSLIo|SlHjYcZPi-dkfUDJsBwkA#K@P zgFMG8dbU8Ob|<91(qex{-Ehl=8UQ5i3cH>|6ukZ zj?xrKaBbD#S4($(tH$*tP_t6^dgZQW_e#TwhvbD+jlZ?)rO&T51|IvBh^agv19Slc z+}LNZS$oRi;w|I8-Hhgo{y6jJt&BPgBSH9~06ZYw@6SR9*9_$0xT{!?*azGxdx?Ba z&`@utM)dt}#jp6b7VB%^oFRS`KXXtcs4~nocE~<*1bWi4rTS&|Wph;`&(2}(+)I_S z&c~^FhjAr8zt5PjdzF;|r9I0%>*Qgb9k+YpYnjpko4>*HBX~7i24qCKv?=^sen^RQ z^%em?|Nim_%aXUH3+sQd$9m+PZmrb7VUk+xrmd#i=xkuGJ*6Z{=hFo_u38j1~ z-)0Q3kzVq2nl+J8nB`$Y2WWSgC{B~u!;_AKD;st8y65WqG zJ&LZ^?q#$XtU-z36e?*D)d`hYZsN>teResJf32+Vv`e|1D88JArl$AzgL2|+cA;#8 zv^?61cpmSd>5@%K6pnI-roWWAh=0D8BWMqMj02r*8N4`|)RU+`C4|pbyN`;v6gz-9 z3isz#{+jCQT8L3rN@!I=>GRWZQ!$Fo1kZq5;IdURhN>Vi z1RS)DhIB^OR5~2)xRveNsn3egMoA2usO_1m-LTTexshvY`{6Sc)^*$V&fzvsm2;)~ zM$K#5W(|QtoBR62YLzJ6@A#eG;#$mzW3BlZ`vh7yhf&&LL@e__yq87it#H;c3p*D`G$x|pq$$kN9v}c>cv&To;3S-~HgIaHJXne8 zO*I~=r=sR=DstpnQLIuW<>o9?b!x`fTA{@f+Aoo-%8%Z^l=-+&%mPBEw!c}1N^@S_ zQ0#3od2`~kqS43_I_60acu@A)M65G)T^|R2&*=|E(#6(Lv=n=edJlN}q<;)tu|6_J z&HGvO`C?R`v{GjA@z`PU5#!*d+ue2(DWCn#LM8-U0Go^x=TuvjOIQ2DP=3-f>|+q# zeV3`>lc7o+#mr()XjKtw6}G`TU#p4Vg`C_LlBf&9i6vbbNd+VAqKv8?+WB!2?T)s% z6j36&DK}_Ltl7j5dD>vQAN>lLDIq;^X(|%G!K+!FX<=-_nje#RnxjQCJbu$c zu0DkJ^_*ht+g% zB@;Js>fjLn2~pahl*5$SZooOK)_oOZegBgJ4E#Fo|MLLiuBXG=E7USSIw}*s zD1;LhDo6noiIvD6`X=E;N(}jYVe^An&~)huhZqsDEH6QSDE@SpSkdKVdS7MOS@9(2 z2==iMtyIIGC1{n1OP4pwx1_EoOeEZH9BP}wQGtW) zGI6XVO}=hUf(9>sV|^67_M2_b5fz6HWAOB%Ngxf##e|ZZmQ{g}7u3Ldth4|{v%dYCO$}Yn50KFOP6l4m+p8p#ssKA*5fr$RN^pvK>}=bk*%* z_i$z30qdf!T7Rw+F+MQxt9oLAJF|9AKhuGWD*j&K zjx7Nw@*~RnDb&_O1serr?yg^Od;aZAiPz(#nHpWgXTV;osz&IwxMngtvEmp*DgGpM z?$7)#!k4}D=Q``n))8OAmvzT!pJIV>721yS?oT_PtPMG4@k;I`F58H+mXtS2KTA?@ z%kWf`leS22cU-QpV;(cPwM@KEdqs~{xoYk(^E1{a1oTU?@RQb+UpNk6& zrZ%xay;haUr`F)BWdSF|-MJ@s+y7EUcT&CBSboh&%g`FFmN?(*2iv0l>+UH%HK>-< z8L}Q#uH|WqN#SwDX>JCGDu$>SjxO+RgqcR>u~strQb}3f0?b#QR6pvomoEBMdOkYj z8}!>l0+SMxfzz^vL7hhkD!;hn*Od^U;}(1DPgWB6p_0<}b2MhWDs*%@aQ>$~?j_yvtRK4Tpmc7?dLhbOoXUx=C;|_yIElx( z=}myla17@JOQ~O=Ab2#+M4y}&GNbdQnfoU61#Wm}j)cz@t)@Vf%2P~C1}9ELT*9QB ziC7?D#=s9vhEgCx?K&EAy6&4XB9|25;>9>cDY~8-#F;CB#Glw_rwx1$i0`UAiB$H$ zDz_o3_B95YjSJdSJm^?LcDu^6x}4QH8PHc~#tW=s4omnclSLB7!PRVyg_^^=uTiEcNko@ehi&uTad@hai&9ZKp?K?^=}_3Ivy_GIG9?P$#S%LJ1zy7*IDl9#8t+zaly5X< zDUQ7z$8G_R5o^1$d@t3SXD^9l;U`y(vj=`}-WyXMVXIhVnh>wE8it-<=sX<4_BBaU z!aM8e7mQDc?(dIPk)=Z>p7NiQ*!*?=6fu87DVEd6&5xP~hB1i?LBsa^*0s9=+2}fD z&w999lN{?PdbBGOAcEAo>1CGjxX=~#Ma`zjc-mW>XW(!HV58$C9`n0KGB_bd82Vk- zPGLZZZwWivyZ80v_&k+tm^UNtKtC19nVEM$G}L!X)TE;5a2Ez#5nYk35v|sfQq4JB zNA? zOk6)EcuCrU%J@v=UXX65nhooH7m`vUGI#fwXNCGNSwmiV2p5kRQa6jGF@{I_-e;ni z#)iiD+$1485G&Hh$*&<+Q8dHoc7&S+qjuyxR=oSk`X8SO@&tvnMyos3e{1CSy5H{i zmVKadfzDFpiB~+XYuZYc{g<-bzc}{!gM3#sw&CS~$XSg6n}Nr+e8;n4R8A~_F(q{+2kMb=S+9v6$WbKPt&Jm z4hE#c=4TXc?$Rj@evlQgrKNE^6=g_pjNWvtsNFIq=Z5{>pHbIKm$Zdh?y&w4hxG9x zZG~|H8k#$v;U90GW=v+dn}=&S2}X5Ke3#jsqG@Q{ozrceuWzrr7_({p4q(X;v3N{* z6&vjTB+$fA2hqV6v6wDRBhuPr-X6BMTFI{bQ$=ot=R!$INEUInO*LbSnB;d^Nv+kN z|7SSHkK8A$I+p>QnK`~y8KoLP#eTp%ezRu{(6$+W5>Qk^Y?o;-nT(?5pxO~LU+FAj zl*2Hps8!-cDT5^pYSd+tXrY%kZw6L#lh7RwkL&|}yWsjJ>&B8yGnP;B}vt(0lzGsb#tQg z^~q?(r*h+W3D_06Vlm{bJ}_EB9Ir6LAjxS3ag|-jcy@Nf#ZKqUWF1AF``dZ~*jy`r z?T1gZI5%D*)a(g$KldR3G3V0XWp;Odj2K9<^gSd<%@1KEgvxU!+fZj%Czl$g}OX8x|(LPJEke``eRKeUUFD6>1-Sg#0TQ;_YpnpkSb84>vvM9066>@WR z+X6T;C&%@GQ2vuT&NO-s#D)NrLD}1C4@qILEd z@8BVm`bB8GhYsV-#oa=EbDSbji;tnMrE&pKG~1xZRkfm%{QYakY`tErRE`i2K_(kF zh3b6*<(RRuaj$Si=6IVW7HIyoMylDUCUZfmYQtcdDs9$&@l*?`t+=+`qU5QZ2vuZj z`t`Q4F}7sb4Hdpws%*W5KS!!;<(P4FTN3N1&bg;L$G|~x1xeB`Z#Ng9*89ub`9X-n zI%}0fqvaytzGcCt4vU1VXb$|kpxqMT6BlxnejLk`5m6qCrmU8WY=8Bl-GXoDoV|YY zqOazF{Ze(uK@t<~oFFALpFn{1!pd+~ydWvf1olzQq~2jQZM(HS-ZI(T|LpS>X@`> zn`JdkQs-bY!Jg9_W>W08pSNQYf^hg;Fh%Bq@RkkqzF%^X%Gbk?kejP)Eu{dSpk@rv)R*pMD1^u8gKP(&#~D{0Gw6Ldd=%6oO+rnc1ytV6hKW!Ct5kXa{6jA zE68Q~K&s&F+!tcP<7Uq$FB>zO&RKwkcndXQ=ELH4|b*fyBhlPaY*`WjCyY^x%+D%2b z&W9eq-`D%l&rjRgUB1sBR6fAB3o4)Yo1OjcCfm6&;+C-A=DYU}KO~FScb$u*{G>E~ z(F%r?4htV&&&=Yr)IIx+g!g+yzx?9eALL{e8iL1>HC6=XS{0zCl$kBdcD7&OD~qKa zp@>!`zy=d-UA6EQ@G~Yk7&kKb$esRj8R1h1G($x6f#Aq_%?lyM#s&P0{7#WPx=!Y@ zi@^!8xV5gYJr8{X;Stm!YPbcC*%m-%5_N=Sljw{`JljVcjJ#kj&f~1q zzRr{mz0(FH|429V&oS(ttYzLyb^epQIW!m%h+h@qkT(1YPcrBfj1m6Ie--FU-23w1 zR6|0f52xqw8CSnzonwJ!!osL4<&XPXv zF{!9+Q3|0-H|gUJxp-?-zrV$r} z@?jLSVOHxd1T@sXb*5#GpJ4V+hVZvaYLO{?l|~C$+d840sD#$`nqtDu>r6uThE!P# z_SGFbmZnFD=2Dy2llbyN8~wBZOi-0BmPns0X;7)jJq)+5Jv4>)RxG4Ac>Lxj?UW6p zB`Y-iNFJwSk8+LP>cWzL7p;)cT|QwIR;kq1)9zhqYZgd=l^DHjLoLoEi+c*aY(} zms8%b3xY9%7JqKl-GIfFJk?dtO%V3$>{Q9A(WkZdV(WYM^<%>ceeHmb#SiS>_842n zN&O1x4G0RT1i))F?jFTU#$g*wMCY*Up*v-In>dnjw?ED|?^k|cYU+3LDPlU=HM6nx zwFGGPmGH)ifmX9Gn9MI=h|C=o;Y>PUok8cZgf_$-bILJZnt_D)7{tF6R$%49%IR5h zx86KPx=>wjO{e76Wbm11NvX~8+nWxzn>3p7X*z_o!y$IA6Aap>L@l4O>4ytU8A++> zmFE#6Fj*KwiZ{oAW`8z11gU>oN<<pj3G>FJ#2yzWf&RIIXFc;Pp{Y7|WfIiYaoVbE^jn(GbTM4K6N)sYPLIJ@M^rM+(I-vq%gYA6$*>pMS0AL;)=^4#TGLsN9KTshnjTM+xytxzkiyp|5#q|M zyB2kKfcK~40wpWfJF>Y;FI929L;&4x(oKEp?SoyIl7k{bm0xa{OjV)gUf51Bd7Tsw z>+wI<`|_ke;y;=7dqL)Z4)zi(W(%UJz5E!{xVucD&xG*>y@x~3JQ!k&POCAXwC~7# z@k3OnNDQxqVU}PsO}I=DHyq65;W*cx+gbvsDg)-gfhf7Vh1hf*0V74}FfO*9+Xy&E z7KT=7sEuf)%cU6>D2KB^CZ7)gu-2Lv_C9)?Puf>{rXCb-0QEEV+|WC>uX$# zLfjyo4)iJbRz&T@oKM{L)ZbL>2wY5?kCp1Q^~9Hzf3T2UPp=lo1$zW}vAEDiQ-v9uIwWUaQvXxOwwhB5pWz?ayqnswKYlF# z2)qc9&{_Ml|E6)dT<_lM7Ca3NCGH*@-x)Z9*t`O}_a>|maA}D{oJ@kF#1ioP-1yXY z>B9CZDSGs>_dUyCHmYzyJw(uL{^CLQ8U^jua!-npjaSNQpnGUwudQi})rWaKIgGV= zPwCpbLx$De%VZdaoXu{(0a6B8=aje2^rpfzm1tp}@Cl5p*2xy!XFByw>qSZS#0_Zf zoIDO158lXDxM!af(Z3(tAy~s~)x`xxF0cYL8y%Xnh;Ba8_l#_vTOl5w&aX>S9|Y&V zRfp>~<7+}`zmGuSubub3KfS+V56P-+W!C(OXKZ!p+D7e=F0{6Cbrgy3_3xq9!gWvQ z`87GfdHk!suYX#bYh_VZtG`+skO2m!X>v_g=p?{%w&@4G-kkfy?)5on>&IGkvO$AP zJ|SIrZUy|_z9~)%vdeI2G(1({)9zl{a)fro&*ueib@a@<*AgI@Nm{~RtLd90AEQ|j zYt*u9nbi8^c%v86iPd(gWK3l3TWgu#tW`4$IJ`5+5Dn=hYGc{7{#B?`z5qOF_AE;# zm|+cr@nUe^%*)H!{~k74r+GQ=wq^FGb*-&{HmJl*07BZVBlo0Kr}}uG<|Wsku6_x4 z`*^iSH^ob{r^ot4|)*3+fM#aG(MCr_yBkaYq&wZgwQjmU43qiYP z{a!PSmw+?NW_RU2Z>tvYROjB8fd-p3=x_ThaV+s}2t1R3vW4;%*R8QXuSR= z<*u@_@Qw#3fg=S9*+v=!X_)&~F401J_=NGWpfOlTowubmU9rjT2nx7XhA zZM<6*RN?XKF03O-OYzx{Q`$e5a6y6Z6H|>9E{#9)Way34yocFNL!WAt$IuTqTdjOP zfN?NcGNO&ix0tv#o}RPzm0yO%b!$lm zx#Ou%RMWlBFYakXoJ3wgnPI>p$-eIm`rE$RQPN>sWNUR^C06-hKxTC=B?~9oJ(QAp zGo6<0_f8Bj20tMZp~#cRADk;%MJoZQtvKWACp)clAI=8K#;bk*!TpOLV&VbncQnlgpzxyj^xc(4HpMYq4QeFnWC|sGbt@iVW+(N$HF}L`q@XDQnw>0B}Qi2 zD!k%1uZO?8GmEf)L%`E^a!iq7w(NC8t?tWbb*sFT*&_ho$4qs*fqJyhQ<2%SgfA~= zsT9d>2xbMVAHFp3;=a0zHZ_I$9@^0oSb(+sMa=?5%_DZ581YP68m%`eWFnsJw<DH~j03@ju~wyIA+dVIHML=Ta!q3518M0sr< zL*gU8388968fOmc^__;rHWG{DCsC}5{TiM2VuPu$mjm1|9k=9zQUV-!*Hv9!#=Xb- zhCFNA>a!>-x=A&@8VVEbrS-O_lb4_0o^>W)lemf%$hJeJFZz(m%Y3YP%4u#y5LD|h zKjqEl|0?qYP898QO{t!E_LCvB6>Jco-Tv~XiQDCpOjD8Q?b@QQP*d-kd z`H;z+EN?#VPBg8L6&Utqfd=4SK=nH?A>J~+kqzKIto`D7993vxNQDgP66)!C6D6Dq zT}zin%1gREG{-}P7+Q{m5y4nO@xh!UrDn$k{YjJq8IDL*Guqf(ps%Hk&c+;G(^ECA za2W}GW9-W8nnyo+no^Kvh$r!u6-GQlP4*u13)`B?%OJdllB09JS;KRSsJJ9b%w8sZG5gij;V?D$4^gg|} z0;6Ans_E<$`*g5cAZHT~qlU0>z{2`61bp6h^tS;)%8o6Rfr6=*+Q2%$biL9D zH6yUiPvT%s9M(8F#}^^ub?lqlPT27JIIS7;e1=gor6Ut*9CD{LYbIOuK4RV0GayOC z6i?WG;`kAt0ias8&&elLId1o<3KzTU_vd4vr0r{EA#c&RMH){aowVRxDjCLzENZ&! zLU1u4xGm@m<~{LH_FBT?<5!MB1T>qH=b?oyx=`|uulXqN&{^VyR9_6&<_iIQnz9_o z?a~#f0T|*nP&DGvn#j-ueA%HndDRp%3skEXag84 zE-6fl)}sehJT1!(GFH3MUtct!3DGb1;Gz zqmlFI>o@u~LlD-ui>HQc^I(rVT)xYa~Vp#AUB3+4use!jF^FOc{siTwl89lKV#w%i*u+aKaZYKeQ}jV z0(G^!FgSf&$3ER3;7+k$nl;9h38g&BMI;j6&ZuN0v$}@xHtN)7T9IoOvvELvn~4!9 zzP^f93ca%(4IXKZAk~b%WM2P%ytEg*{q{4WsapCQ=aLmwh#v^ZKcv8+?Dh<(30_Tx z(&f#(+EEnM@Ub*Qv1R@9bF=y3K_Ka4-eE8Pd#t0f3?k!={m-N3tJ80-2nA7mQm*O+ zsDME1bu)+}=*7?Qj7YgKPHu+mh20u2AlU&wVXM7&7tadn$74*A9P*&bU<~E;=f~5; zJMPSqY#kyZcoZ1mdZ7mtC8KS8@neIhSHf5KQ3*(Xn-ekD#ijCG0O82i+&JfuxPAJ( z9tDq8xGO{Ary|#ZI>h$HUqTH1a0!MK$@Le@uC?eBamfZWrZ94SwYzToGSwpe->DWa zKh?-pjY+t}Eya4SD`0dLp?zkU35by@pjQlaX)bLIQ6($_MwPz^*W;w}%mmvu2p`T_ z2od+6QTb$G-UsXfH8js(6T18Dc_>UsB7KOj88fw+TQm_B`8|sEY#eSjz6w9g+eNI_ zk>m6m%l3nfv?za0lV8>1F*$U!M=5EqJEIlHH>9sfPNE*H4^H)@zJa(Q{Il8K_zn11 zqQGDCv|NP60tU7p4Dm-G0uF%*_Ft6(XpIPh2fnmTQS%=W35bvwf}7+o{J&)lsQ*zi z{QJ2Gav_Ew`NwE`5iAh&r}+BM8{>Z%If3ekA;|wRQtP>vh5aLtV1Nb#Bm5sm(x5A1 z2!_A91tJm%Lhy7DKM4foUug0df3JLcrpmG yzqO2Hf8Yv8mJEXQFSPp$vp|t#{|@sXlKK?x&?^G6{9_J06a2s