From 704bfb3b224153e5bd4364741d97907c4b8747c6 Mon Sep 17 00:00:00 2001 From: Zhiyu <zhiyu.zhao@utsouthwestern.edu> Date: Wed, 26 Feb 2020 15:34:44 -0600 Subject: [PATCH] V1.7.1: first commit --- V1.7.1/ODA_1.7.1.R | 1075 +++++++++++++++++++++++++++++++ V1.7.1/ReadMe.txt | 25 + V1.7.1/data_template_1.7.1.xlsm | Bin 0 -> 96603 bytes V1.7.1/run_analysis.sh | 17 + 4 files changed, 1117 insertions(+) create mode 100755 V1.7.1/ODA_1.7.1.R create mode 100755 V1.7.1/ReadMe.txt create mode 100755 V1.7.1/data_template_1.7.1.xlsm create mode 100755 V1.7.1/run_analysis.sh diff --git a/V1.7.1/ODA_1.7.1.R b/V1.7.1/ODA_1.7.1.R new file mode 100755 index 0000000..60955fe --- /dev/null +++ b/V1.7.1/ODA_1.7.1.R @@ -0,0 +1,1075 @@ +verNum='V1.7.1' +#This script was developed by zhiyu.zhao@utsouthwestern.edu . +#v1.7.1, 12/04/2019: Figures were saved as both .ps and .png files. Euler plot was disabled for >6 samples. Samples / features can be clustered or not on the heatmaps. +#v1.7.0, 10/22/2019: Figures were saved as .ps files in a Figure folder. Correlation plot was generated with the old method for better speed. Euler plot was disabled for >10 samples. Y axis of volcano plot could be raw or adjusted p-values. DESeq was renamed RLE as I would implement the DESeq differential expression analysis in the future. RLE is the same method used by DESeq but I slightly modified it to allow for zeros, missing values and / or values below the value filter. Due to this reason I couldn't use the DESeq normalization directly. +#v1.6.9, 10/4/2019: A bug regarding cpglm was fixed. The bug caused an error in cpglm fitting but the class(fit) was still 'cpglm' and fit$converged was still TRUE. Visualization option was added. Version number checkup was added for the data template. Multiple other minor improvements. +#v1.6.8, 09/19/2019: A bug regarding featureDescription was fixed. The bug was present only when data format was "table with feature description" and samples were excluded which happened to cause some features to be all-empty so that those features needed to be excluded as well. +#v1.6.7, 08/29/2019: Correlation plot was enchanced. Euler diagram was added. Software version number was reported in the result file. +#v1.6.6, 08/28/2019: Data table with feature description column was enabled, minor bugs fixed, and volcano plot colors customized. +#v1.6.5, 08/27/2019: A bug regarding number of features per sample when data imputation = "NA" was fixed. +#v1.6.4, 08/08/2019: A bug regarding groupMeanLowerBound was fixed. Definition and default value of groupCountFilter were modified. Data imputation was implemented. Missing values were allowed. +#v1.6.3, 08/05/2019: A bug regarding nSamples and formula string generation (only when technical replicates are present) was fixed. +#v1.6.2, 07/24/2019: A bug regarding normalization size factors was fixed. The bug caused size factors to be applied to incorrect samples. If you have used versions before this fix, please rerun your analysis using this fixed version. +#To run this program, first install R (version 3.5.1 was tested by the author), then run the following from the command line of your operating system. + +#Rscript path1/this_program.R path2/your_input_file.xlsm path3/your_output_file.xlsx #Depending on your operating system, / can be \. + +#Use the data_template_1.6.xlsm file accompanying this program to generate your_input_file.xlsm. All the parameters can be changed inside that file. Please read the instructions there. Excel macros need to be allowed. + +#############################Begin Reading Arguments############################# +if (1) { #Set to 1 for actual files or 0 for debugging a test.xlsm in the R console. + args=commandArgs(trailingOnly=T) + if (length(args)!=2) + stop('Please specify input and output file names.') + inputFile=args[1] + outputFile=args[2] +} else { + inputFile='test.xlsm' + outputFile='Result_of_test.xlsx' +} +if (!file.exists(inputFile)) + stop('Input file is not found.') +if (file.exists(outputFile)) + stop('Output file exists.') +figuresPath=paste0(tools::file_path_sans_ext(outputFile),'_Figures') +if (file.exists(figuresPath)&&!rlang::is_empty(dir(figuresPath))) { + stop('Figures folder is not empty.') +} else { + dir.create(figuresPath) +} +#############################End Reading Arguments############################# + +#############################Begin Installing/Loading Libraries############################# +myLibPath='~/R/3.5.1-gccmkl' #Specify your own library path here. +dir.create(myLibPath,recursive=T,showWarnings=F) +.libPaths(c(myLibPath,.libPaths())) +#This program was tested under R version 3.5.1 with the following required libraries. +if (!require(openxlsx)) { + install.packages('openxlsx',lib=myLibPath,repos='https://cran.r-project.org/') + library(openxlsx) +} +if (!require(data.table)) { + install.packages('data.table',lib=myLibPath,repos='https://cran.r-project.org/') + library(data.table) +} +if (!require(gtools)) { + install.packages('gtools',lib=myLibPath,repos='https://cran.r-project.org/') + library(gtools) +} +if (!require(matrixStats)) { + install.packages('matrixStats',lib=myLibPath,repos='https://cran.r-project.org/') + library(matrixStats) +} +if (!require(cplm)) { + install.packages('cplm',lib=myLibPath,repos='https://cran.r-project.org/') + library(cplm) +} +if (!require(ggplot2)) { + install.packages('ggplot2',lib=myLibPath,repos='https://cran.r-project.org/') + library(ggplot2) +} +if (!require(cowplot)) { + install.packages('cowplot',lib=myLibPath,repos='https://cran.r-project.org/') + library(cowplot) +} +if (!require(pheatmap)) { + install.packages('pheatmap',lib=myLibPath,repos='https://cran.r-project.org/') + library(pheatmap) +} +if (!require(ggcorrplot)) { + install.packages('ggcorrplot',lib=myLibPath,repos='https://cran.r-project.org/') + library(ggcorrplot) +} +if (!require(eulerr)) { + install.packages('eulerr',lib=myLibPath,repos='https://cran.r-project.org/') + library(eulerr) +} +if (!require(GGally)) { + install.packages('GGally',lib=myLibPath,repos='https://cran.r-project.org/') + library(GGally) +} +#############################End Installing/Loading Libraries############################# + +#############################Begin Main############################# + +#############################Begin Default Settings############################# +#Below are default settings. DO NOT change them here. +#If you use the data_template.xlsm file accompanying this program to generate your input, all the parameters can be changed inside your_input_file.xlsm . Please read the instructions there. +dataFormat='list' #Available options: list; table. +valueFilter=0 #A cutoff-value used to count the number of features in a sample. If feature value <= valueFilter, the feature is taken as absent from the sample. +groupCountFilter=3 #A cutoff-value used to count the number of features in a group. If the number of samples < groupCountFilter for a feature, the feature is taken as absent from the group. +groupMeanLowerBound=0 #A cutoff-value used as group mean when it's lower than the cutoff. This is for the adjustment of fold change when a group has a very low mean value. +pCutoffs=0.05 #P-value cutoff for the differential expression (DE) analysis. Multiple values are allowed. +fcCutoffs=2 #Fold change cutoff for the DE analysis. Multiple values are allowed. +batchEffectReduction=F #Whether to model batch effect. Only available for linear models such as cpglm and glm. +averagingTechReps=F #Whether to average technical replicates before any analysis steps. +normalization='RLE' #Available options: No; RLE; BySize; ByMedian; BySum; ByUpperQuantile. +uqProb=0.75 #If normalization=ByUpperQuantile, this value is used as the quantile cutoff. +logTransform=T #Whether to log2-transform the data for statistical tests. +logOffset=1 #An offset-value added to the data before the log-transformation. +statTest='glm' #Statistical test method for DE. Available options: t-test; cpglm-glm; glm. +controlGroups=NULL #Control group name(s) for DE. A control group will be tested against every other group. +dePairs=NULL #Group name pair(s) for DE. Each pair of group names is a DE comparison. +interestedFeatures=NULL #Interested features to be visualized separately. +excludedFeatures=NULL #Features to be completed excluded from the analysis. +dataImputation='half minimum' #Method to impute missing values and those less than or equal to the value filter. It affects statistical tests but not visualizations. See data template for details. +visual=T #Whether to visualize data using graphs. +volcanoP='adjusted p-value' #Which p-value as the y-axis of volcano plots. +hclustSamples=T #Whether to cluster samples on the heatmap. +hclustFeatures=T #Whether to cluster features on the heatmap. +colors=c('Red','Yellow','Blue','Orange','Green','Purple','Magenta','Cyan','Gray','Black') +heatmapColors=colorRampPalette(c('green','black','red'))(256) +#############################End Default Settings############################# +message('###################################################################') +message('Analysis started ...') +#############################Begin Settings############################# +message('Reading input settings ...') +sheetNames=getSheetNames(inputFile) +if ('Samples' %in% sheetNames) { + sampleIDs=read.xlsx(inputFile,sheet='Samples',colNames=T,rowNames=T,skipEmptyRows=T,skipEmptyCols=T,cols=c(1:5)) + if (!is.null(sampleIDs)&&dim(sampleIDs)[1]>0) { + for (i in 1:dim(sampleIDs)[1]) + if (grepl('\\s',rownames(sampleIDs)[i])) + stop(paste0('Sample ',i,': ',rownames(sampleIDs)[i],' has white space(s). Please remove them.')) + sampleIDs=sampleIDs[!is.na(sampleIDs[,1]),1:4] #Excluding samples without group ID. + } + nSamples=dim(sampleIDs)[1] + if (nSamples<2) + stop('At least 2 samples should have group IDs.') + if (sum(!is.na(sampleIDs[,2]))==0) #No biological replicate IDs. Automatic IDs will be generated. + sampleIDs[,2]=c(1:nSamples) + else if (sum(!is.na(sampleIDs[,2]))!=nSamples) + stop('If you provide biological replicate IDs, they should be present for all samples.') + if (sum(!is.na(sampleIDs[,3]))==0) #No batch IDs. Automatic IDs will be generated. + sampleIDs[,3]='Batch1' + else if (sum(!is.na(sampleIDs[,3]))!=nSamples) + stop('If you provide batch IDs, they should be present for all samples.') + sampleIDs[,1]=as.factor(sampleIDs[,1]) + sampleIDs[,2]=as.factor(sampleIDs[,2]) + sampleIDs[,3]=as.factor(sampleIDs[,3]) + sampleIDs[,4]=as.numeric(sampleIDs[,4]) +} else + stop('Sample sheet is missing.') +if ('Parameters' %in% sheetNames) { + if (read.xlsx(inputFile,sheet='Parameters',colNames=F,rows=1,cols=1)!=verNum) + stop(paste0('Please use data template ',verNum,' for this analysis.')) + parameters=read.xlsx(inputFile,sheet='Parameters',colNames=T,rowNames=T,na.strings=NULL) + if (!is.null(parameters)&&dim(parameters)[1]>0) { + if (!is.na(parameters[1,1])) + dataFormat=tolower(as.character(parameters[1,1])) + if (!is.na(parameters[2,1])) + valueFilter=as.numeric(parameters[2,1]) + if (!is.na(parameters[3,1])) + groupCountFilter=as.integer(parameters[3,1]) + if (!is.na(parameters[4,1])) + groupMeanLowerBound=as.numeric(parameters[4,1]) + if (sum(!is.na(parameters[5,]))) + pCutoffs=as.numeric(parameters[5,!is.na(parameters[5,])]) + if (sum(!is.na(parameters[6,]))) + fcCutoffs=as.numeric(parameters[6,!is.na(parameters[6,])]) + if (!is.na(parameters[7,1])) + normalization=parameters[7,1] + if (normalization=='BySize'&&sum(!is.na(sampleIDs[,4]))!=nSamples) + stop('If you provide size factors, they should be present for all samples.') + if (normalization=='ByUpperQuantile'&&!is.na(parameters[7,2])) + uqProb=as.numeric(parameters[7,2]) + if (!is.na(parameters[8,1])&¶meters[8,1]=='No') + logTransform=F + if (!is.na(parameters[9,1])&&as.numeric(parameters[9,1])>=0) + logOffset=as.numeric(parameters[9,1]) + if (!is.na(parameters[10,1])) + statTest=parameters[10,1] + if (!is.na(parameters[11,1])) + dataImputation=tolower(as.character(parameters[11,1])) + if (!is.na(parameters[12,1])&¶meters[12,1]=='No') + visual=F + if (!is.na(parameters[13,1])) + volcanoP=tolower(as.character(parameters[13,1])) + if (!is.na(parameters[14,1])&¶meters[14,1]=='No') + hclustSamples=F + if (!is.na(parameters[15,1])&¶meters[15,1]=='No') + hclustFeatures=F + } +} else + stop('Parameter sheet is missing.') +if ('Comparisons' %in% sheetNames) { + dePairs=read.xlsx(inputFile,sheet='Comparisons',colNames=T,rowNames=F,cols=c(1:2)) + tempStr=unlist(dePairs) + if (sum(!(tempStr[!is.na(tempStr)] %in% sampleIDs[,1]))>0) + stop('Unknown group ID(s) in the Comparisons sheet.') + if (!is.null(dePairs)&&dim(dePairs)[1]>0) { + controlGroups=dePairs[!is.na(dePairs[,1])&is.na(dePairs[,2]),1] + dePairs=dePairs[!is.na(dePairs[,1])&!is.na(dePairs[,2]),1:2] + } +} +if ('Features' %in% sheetNames) { + tempFeatures=read.xlsx(inputFile,sheet='Features',colNames=T,rowNames=F,cols=c(1:2)) + if (!is.null(tempFeatures)&&dim(tempFeatures)[1]>0) { + interestedFeatures=tempFeatures[!is.na(tempFeatures[,1]),1] + excludedFeatures=tempFeatures[!is.na(tempFeatures[,2]),2] + } +} +#############################End Settings############################# + +#############################Begin Reading Data############################# +#Please read instructions in the data_template.xlsm and use the RawData sheet to hold your data. +#If the data are generated by the metabolomics core facility at CRI, UTSW, usually they are in the list format. +#If the data are in the list format, in the excel sheet three columns will be read: feature ID in column 1, sample ID in column 2, and value in column 5. +#If the data are in the table format, in the excel sheet all columns will be read: feature ID in column 1 and sample values in the next columns. +message('Reading data ...') +if (dataFormat=='list') { + rawData=read.xlsx(inputFile,sheet='RawData',colNames=T)[,c(1,2,5)] +} else if (dataFormat=='table') { + data=read.xlsx(inputFile,sheet='RawData',colNames=T,rowNames=F) + featureDescription=NULL + rawData=melt(data,id.vars=1) + rawData[,2]=as.character(rawData[,2]) + rawData[,3]=as.numeric(rawData[,3]) +} else if (dataFormat=='table with feature description') { + data=read.xlsx(inputFile,sheet='RawData',colNames=T,rowNames=F) + featureDescription=data[,1:2] + data=data[,-2] + rawData=melt(data,id.vars=1) + rawData[,2]=as.character(rawData[,2]) + rawData[,3]=as.numeric(rawData[,3]) +} else { + stop('Only list and table are supported data formats. Check your Parameters sheet.') +} +colnames(rawData)=c('Feature','Sample','Value') +if (!is.null(excludedFeatures)) + rawData=rawData[!(rawData$Feature %in% excludedFeatures),] #If speficied, features are excluded. +rawData=rawData[!is.na(sampleIDs[rawData$Sample,1]),] #Samples with empty group IDs are excluded. +if (sum(is.na(rawData$Value))) + message('Missing values are detected.') +#rawData[is.na(rawData$Value),'Value']=0 +if(sum(rawData$Value<0,na.rm=T)&&(normalization!='No'||logTransform)) + stop('Negative value(s) are not allowed when you normalize or log-transform your data.') +#The following code is to generate automatic feature IDs in case they are not unique. +if (length(unique(rawData$Feature))*length(unique(rawData$Sample))!=dim(rawData)[1]) { + rowNames=rawData$Feature + ind=data.frame() + for (i in 1:dim(rawData)[1]) { + r=which(rownames(ind)==rawData$Feature[i]) + c=which(colnames(ind)==rawData$Sample[i]) + if (is.na(r==0||c==0)||is.na(ind[r,c])) { + + ind[rawData$Feature[i],rawData$Sample[i]]=1 + } else { + #Currently the auto-ID function is disabled so that problems caused by mistakes are not ignored. If preferred, comment the following line to enable this function. + stop(paste0('Duplicated feature or sample name(s) found at ',r,' : ',c,',. Please check your data.')) + rowNames[i]=paste0(rowNames[i],'.',ind[r,c]) + ind[r,c]=ind[r,c]+1 + } + } + rawData$Feature=rowNames #In case there are duplicated feature names, they are replaced by name.1, name.2, and so on. +} +rawData[,colnames(sampleIDs)]=sampleIDs[rawData$Sample,] +rawData=aggregate(Value~Feature+BatchID+GroupID+BiologicalReplicateID,FUN=mean,data=rawData) #Data are stored in a list. Technical replicates are averaged if present. +rawData$Sample=paste(rawData$BatchID,rawData$GroupID,rawData$BiologicalReplicateID,sep='_') +data=dcast(rawData,Feature~BatchID+GroupID+BiologicalReplicateID,value.var='Value') #Data are also stored in a table with features in rows and samples in columns. +sampleIDs=unique(sampleIDs) +sampleIDs$OriginalID=rownames(sampleIDs) +rownames(sampleIDs)=paste(sampleIDs$BatchID,sampleIDs$GroupID,sampleIDs$BiologicalReplicateID,sep='_') +rownames(data)=data[,1] +data=data[,-1] +data=as.matrix(data) +class(data)='numeric' +sampleIDs=sampleIDs[colnames(data),] +#############################End Reading Data############################# + +wb = createWorkbook() +style1=createStyle(textDecoration='Bold') +style2=createStyle(textDecoration='Bold',wrapText=T) +addWorksheet(wb,'Visual') +addWorksheet(wb, 'Data') +addWorksheet(wb, 'Overall') +addWorksheet(wb, 'Group_Features') +addWorksheet(wb, 'Sample_Features') + +#############################Begin Sheet "Data"############################# +message('Writing the data sheet ...') +if (dataFormat=='table with feature description') { + rownames(featureDescription)=featureDescription[,1] + featureDescription=featureDescription[rownames(featureDescription)%in%rownames(data),] + featureDescription=featureDescription[rownames(data),-1] + writeData(wb,'Data',t(rbind(NA,sampleIDs)),rowNames=T,colNames=T,keepNA=F) + writeData(wb,'Data','',startRow=1,startCol=2) + writeData(wb,'Data',data.frame(Description=featureDescription,data,stringsAsFactors=F,check.names=F),startRow=7,rowNames=T,colNames=F) + writeData(wb,'Data','Description',startRow=6,startCol=2) + freezePane(wb,'Data',firstActiveRow=7,firstActiveCol=3) +} else { + writeData(wb,'Data',t(sampleIDs),rowNames=T,colNames=T) + writeData(wb,'Data',data,startRow=7,rowNames=T,colNames=F) + freezePane(wb,'Data',firstActiveRow=7,firstActiveCol=2) +} +nSamples=dim(data)[2] #Bug fixed on 8/5/19. +nFeatures=dim(data)[1] +groupID=as.factor(sampleIDs[,1]) +names(groupID)=rownames(sampleIDs) +groups=levels(groupID) +nGroups=length(groups) +repID=as.factor(sampleIDs[,2]) +names(repID)=rownames(sampleIDs) +batchID=as.factor(sampleIDs[,3]) +names(batchID)=rownames(sampleIDs) +formulaStr='value' +if (length(levels(groupID))>1) { + if (formulaStr=='value') + formulaStr=paste(formulaStr,'group',sep='~') + else + formulaStr=paste(formulaStr,'group',sep='+') +} +if (length(levels(repID))>1&&length(levels(repID))<nSamples) { + if (formulaStr=='value') + formulaStr=paste(formulaStr,'rep',sep='~') + else + formulaStr=paste(formulaStr,'rep',sep='+') +} +if (length(levels(batchID))>1) { + if (formulaStr=='value') + formulaStr=paste(formulaStr,'batch',sep='~') + else + formulaStr=paste(formulaStr,'batch',sep='+') +} +formulaForFitting=NULL +if (formulaStr!='value') + formulaForFitting=as.formula(formulaStr) +#############################End Sheet "Data"############################# + +#############################Begin Sheet "Normalized_Log2_Data"############################# +message ('Normalizing and transforming data ...') +normData=rawData #Long-list +sizeRatios=data #Wide-table +sizeRatios=sizeRatios[rowSums(sizeRatios,na.rm=T)!=0,] #All-zero rows are not used for normalization. +if (normalization!='No'||logTransform) { + sheetName='Data' + if (normalization!='No') { + if (normalization=='RLE') { + sizeRatios[sizeRatios==0]=NA #Zero values are ignored. +# sizeRatios=sizeRatios/exp(rowMeans(log(sizeRatios),na.rm=T)) #Buggy: Each row is divided by its geometric mean. + sizeRatios=sweep(sizeRatios,1,exp(rowMeans(log(sizeRatios),na.rm=T)),'/') #Each row is divided by its geometric mean. + sizeFactors=colMedians(sizeRatios,na.rm=T) #Median of ratios within each sample is taken as its size factor. + sizeRatios=melt(log2(sizeRatios),varnames=c('Feature','Sample'), value.name ='Value') + sizeRatios$Sample=as.character(sizeRatios$Sample) + sizeRatios[,colnames(sampleIDs)]=sampleIDs[sizeRatios$Sample,] + } + else { + if (normalization=='BySize') + sizeFactors=sampleIDs[,4] #If normalize by size, size factors are taken from the input file. Size factors can be e.g. cell counts, spike-in amount etc. + if (normalization=='ByMedian') #If normalize by median, the median feature value of each sample is taken as its size factor. + sizeFactors=colMedians(sizeRatios) + if (normalization=='ByUpperQuantile') #If normalize by upper quantile, the (uqProb*100)-th quantile value of each sample is taken as its size factor. uqProb is taken from the input file or defaults to 0.75 i.e. the upper quartile. + sizeFactors=colQuantiles(sizeRatios,probs=uqProb,na.rm=T) + if (normalization=='BySum') #If normalize by sum, the sum of feature values of each sample is taken as its size factor. + sizeFactors=colSums(sizeRatios,na.rm=T) + } + if (sum(sizeFactors==0)==0) { + sizeFactors=sizeFactors/median(sizeFactors) + data=sweep(data,2,sizeFactors,'/') #Raw data are normalized by the size factor per sample. A bug regarding this line is fixed. + names(sizeFactors)=colnames(data) + normData$Value=normData$Value/sizeFactors[normData$Sample] + sheetName=paste0('Norm_',sheetName) + } + else { + message('Size factor(s)=0. No normalization is applied.') + normalization='No' + } + } + if (logTransform) { + if (sum(data+logOffset<=0,na.rm=T)) + stop('Offset of log2-transformation is too small to make all values positive.') + data=log2(data+logOffset) + sheetName=paste0('Log2(x+',logOffset,')_',sheetName) + } + addWorksheet(wb, sheetName) + if (normalization!='No') { + sampleIDs=data.frame(sampleIDs[,1:4],SizeFactor=sizeFactors,OriginalID=sampleIDs[,5]) + if (dataFormat=='table with feature description') { + writeData(wb,sheetName,t(rbind(NA,sampleIDs)),rowNames=T,colNames=T,keepNA=F) + writeData(wb,sheetName,'',startRow=1,startCol=2) + writeData(wb,sheetName,'Description',startRow=7,startCol=2) + writeData(wb,sheetName,data.frame(Description=featureDescription,data,stringsAsFactors=F,check.names=F),startRow=8,rowNames=T,colNames=F) + freezePane(wb,sheetName,firstActiveRow=8,firstActiveCol=3) + } else { + writeData(wb,sheetName,t(sampleIDs),rowNames=T,colNames=T) + writeData(wb,sheetName,data,startRow=8,rowNames=T,colNames=F) + freezePane(wb,sheetName,firstActiveRow=8,firstActiveCol=2) + } + } else { + if (dataFormat=='table with feature description') { + writeData(wb,sheetName,t(rbind(NA,sampleIDs)),rowNames=T,colNames=T,keepNA=F) + writeData(wb,sheetName,'',startRow=1,startCol=2) + writeData(wb,sheetName,data.frame(Description=featureDescription,data,stringsAsFactors=F,check.names=F),startRow=7,rowNames=T,colNames=F) + writeData(wb,sheetName,'Description',startRow=6,startCol=2) + freezePane(wb,sheetName,firstActiveRow=7,firstActiveCol=3) + } else { + writeData(wb,sheetName,t(sampleIDs),rowNames=T,colNames=T) + writeData(wb,sheetName,data,startRow=7,rowNames=T,colNames=F) + freezePane(wb,sheetName,firstActiveRow=7,firstActiveCol=2) + } + } +} +#############################End Sheet "Normalized_Log2_Data"############################# + +#############################Begin Sheet "Selected_Features"############################# +message ('Writing the selected features sheet ...') +nSelectedFeatures=0 +if (length(interestedFeatures)>0) { + selectedFeatures=rownames(data) %in% interestedFeatures + nSelectedFeatures=sum(selectedFeatures) +} +if (nSelectedFeatures>0) { + if (dataFormat=='table with feature description') { + report=data.frame(Description=featureDescription,data,stringsAsFactors=F,check.names=F)[selectedFeatures,] + } else { + report=data.frame(data,stringsAsFactors=F,check.names=F)[selectedFeatures,] + } + addWorksheet(wb, 'Selected_Features') + writeData(wb,'Selected_Features',report,rowNames=T) + writeData(wb,'Selected_Features','Feature') + addStyle(wb,'Selected_Features',style2,rows=1,cols=1:200) +} +#############################End Sheet "Selected_Features"############################# + +#############################Begin Sheet "Visual"############################# +if (visual) { +message ('Generating figures ...') +options(warn=-1) +plotID=0 +plotTypes=c('raw') +if (logTransform) + plotTypes=c(plotTypes,'log2') +if (normalization!='No') { + plotTypes=c(plotTypes,'norm') + if (logTransform) + plotTypes=c(plotTypes,'log2_norm') +} +if (normalization=='RLE') + plotTypes=c(plotTypes,'sizes') +if (nSelectedFeatures>0) + plotTypes=c(plotTypes,paste0(plotTypes,'_features')) +for (plotType in plotTypes) { + if (grepl('sizes',plotType)) { + plotData=sizeRatios + titleText='Log2-Transformed Size Ratios for RLE Normalization' + } else if (grepl('norm',plotType)) { + plotData=normData + titleText='Normalized Data' + } else { + plotData=rawData + titleText='Raw Data' + } + if (grepl('log2',plotType)) { + plotData$Value=log2(plotData$Value+logOffset) + titleText=paste0('Log2(x+',logOffset,')-Transformed ',titleText) + } + if (grepl('features',plotType)) { + plotData=plotData[plotData$Feature %in% interestedFeatures,] + titleText=paste0(titleText,' (Selected Features Only)') + } + if (sum(!is.na(plotData$Feature))==0) + next + tempData=dcast(plotData,Feature~Sample,value.var='Value') #Table format of the plotData. + rownames(tempData)=tempData[,1] + tempData=as.matrix(tempData[,-1]) + nonZeroFeatures=rownames(tempData)[rowSums(tempData,na.rm=T)!=0] + if (length(nonZeroFeatures)>0) { + plotData=plotData[plotData$Feature %in% nonZeroFeatures,] #Zero features are excluded from visualization. + plotFile=file.path(figuresPath,paste0(plotType,'.ps')) + plotFile2=file.path(figuresPath,paste0(plotType,'.png')) + p1=ggplot(plotData, aes(x=Sample, y=Value, color=GroupID))+geom_violin(color='black')+geom_jitter(position=position_jitter(0.2),size=0.5)+geom_boxplot(color='black',outlier.size=-1,width=0.1)+facet_grid(.~BatchID,scales='free_x',space='free_x')+theme(axis.text.x = element_text(angle = 90, hjust = 1),legend.position = "none")+ggtitle(paste0('Violin-Boxplot of ',titleText)) + p2=ggplot(plotData,aes(x=Value,color=GroupID))+geom_histogram(bins=20)+facet_grid(.~BatchID+Sample)+ggtitle(paste0('Histogram of ',titleText))+theme(strip.text.x = element_blank(),axis.text.x = element_text(angle = 90, hjust = 1),legend.position = "none") + statsData=melt(do.call(data.frame,aggregate(Value~Sample,FUN=function(x) c(Median=median(x),Sum=sum(x)),data=plotData)),id.vars='Sample') + colnames(statsData)=c('Sample','Stat','Value') + statsData$Sample=as.character(statsData$Sample) + statsData$GroupID=sampleIDs[statsData$Sample,1] + statsData$BatchID=sampleIDs[statsData$Sample,3] + p3=ggplot(statsData,aes(x=Sample,y=Value,color=GroupID))+geom_bar(stat='identity')+facet_grid(Stat~BatchID,scales='free',space='free_x')+theme(axis.text.x = element_text(angle = 90, hjust = 1),legend.position = "none")+ggtitle(paste0('Statistics of ',titleText)) + p=plot_grid(p1,p2,align='v',axis='tblr',nrow=2,rel_heights = c(3/5, 2/5)) + p=plot_grid(p,p3,ncol=2,rel_widths=c(3/4,1/4)) + ggsave(plotFile,plot=p,width=min(50,max(dim(data)[2],12)),height=8,units='in',limitsize=F) + ggsave(plotFile2,plot=p,width=min(50,max(dim(data)[2],12)),height=8,units='in',limitsize=F) + insertImage(wb,'Visual',plotFile2, width=dim(data)[2], height=8, unit='in',startRow=plotID*45+1,startCol=1) + plotID=plotID+1 + tempData=tempData[rowSums(is.na(tempData))==0&apply(tempData,1,sd)>0,] #Clustering doesn't work with NAs. Scaling doesn't work with zero-sd rows. + if (!is.null(dim(tempData))) + if (dim(tempData)[1]>1) { + plotFile=file.path(figuresPath,paste0(plotType,'_heatmap.ps')) + plotFile2=file.path(figuresPath,paste0(plotType,'_heatmap.png')) + p=pheatmap(tempData,main=paste0('Heatmap of ',titleText),border_color=NA,show_rownames=grepl('features',plotType),treeheight_row=0,annotation_col=sampleIDs[colnames(tempData),c('BatchID','GroupID')],annotation_legend=T,na_col='gray',scale='none',col=heatmapColors,cluster_rows=hclustFeatures,cluster_cols=hclustSamples) + ggsave(plotFile,plot=p$gtable,width=min(50,max(dim(data)[2],12)),height=8,units='in',limitsize=F) + ggsave(plotFile2,plot=p$gtable,width=min(50,max(dim(data)[2],12)),height=8,units='in',limitsize=F) + insertImage(wb,'Visual',plotFile2, width=min(50,max(dim(data)[2]/3,4)), height=8, unit='in',startRow=plotID*45+1,startCol=1) + plotFile=file.path(figuresPath,paste0(plotType,'_scaled_heatmap.ps')) + plotFile2=file.path(figuresPath,paste0(plotType,'_scaled_heatmap.png')) + p=pheatmap(tempData,main=paste0('Scaled Heatmap of ',titleText),border_color=NA,show_rownames=grepl('features',plotType),treeheight_row=0,annotation_col=sampleIDs[colnames(tempData),c('BatchID','GroupID')],annotation_legend=T,na_col='gray',scale='row',col=heatmapColors,cluster_rows=hclustFeatures,cluster_cols=hclustSamples) + ggsave(plotFile,plot=p$gtable,width=min(50,max(dim(data)[2],12)),height=8,units='in',limitsize=F) + ggsave(plotFile2,plot=p$gtable,width=min(50,max(dim(data)[2],12)),height=8,units='in',limitsize=F) + insertImage(wb,'Visual',plotFile2, width=min(50,max(dim(data)[2]/3,4)), height=8, unit='in',startRow=plotID*45+1,startCol=ceiling(max(dim(data)[2]/3,4)*2)) + plotFile=file.path(figuresPath,paste0(plotType,'_pca_cor.ps')) + plotFile2=file.path(figuresPath,paste0(plotType,'_pca_cor.png')) + plotData=data.frame(prcomp(t(tempData))$x,sampleIDs) + p1=ggplot(plotData,aes(x=PC1,y=PC2,color=GroupID))+geom_point(aes(shape=BatchID))+ggtitle(paste0('PCA of ',titleText)) + if (dim(tempData)[1]>2) + pMat=cor_pmat(tempData) + else + pMat=NULL + if (sum(is.na(cor(tempData)))>0) + hClust=F + else + hClust=T + p2=ggcorrplot(cor(tempData), hc.order = hClust, outline.col = "white",p.mat=pMat)+ggtitle(paste0('Pearson Correlation of ',titleText)) + p=plot_grid(p1,p2,ncol=2,rel_widths=c(1/3,2/3)) +# p2=ggpairs(as.data.frame(tempData),lower=list(continuous='smooth'),title=paste0('Scatter Plot of ',titleText)) +# p=plot_grid(p1,ggmatrix_gtable(p2),ncol=2,rel_widths=c(1/3,2/3)) + ggsave(plotFile,plot=p,width=min(50,max(dim(data)[2]*2,16)),height=8,units='in',limitsize=F) + ggsave(plotFile2,plot=p,width=min(50,max(dim(data)[2]*2,16)),height=8,units='in',limitsize=F) + insertImage(wb,'Visual',plotFile2, width=min(50,max(dim(data)[2]*2,16)), height=8, unit='in',startRow=plotID*45+1,startCol=ceiling(max(dim(data)[2]/3,4)*4)) + plotID=plotID+1 + } + } +} +} #End if (visual) +options(warn=0) +#############################End Sheet "Visual"############################# + +#############################Begin Sheet "Overall"############################# +message('Writing the overall summary sheet ...') +featureExpressed=(data>valueFilter) +nValueFilteredFeatures=sum(rowSums(featureExpressed,na.rm=T)>0) +if (sum(data>valueFilter,na.rm=T)>0 && groupMeanLowerBound==0) { + if (logTransform) + groupMeanLowerBound=min(data[data>valueFilter],na.rm=T)-1 #Bug fix: this was added on 8/8/19. + else + groupMeanLowerBound=min(data[data>valueFilter],na.rm=T)/2 +} +report=matrix(nrow=19,ncol=1) +colnames(report)='Value' +rownames(report)=c('Samples','Groups','Features','Value-filtered features','Value filter','Count filter','Group mean lower-bound','Adjusted p-value cutoff','Fold change cutoff','Normalization','Log2 transformation','Offset of log2 transformation','Statistical test method','Data imputation method','Formula for linear model','Visualization','Y axis for volcano plot','Clustered samples','Clustered features') +report[1]=nSamples +report[2]=nGroups +report[3]=nFeatures +report[4]=nValueFilteredFeatures +report[5]=valueFilter +report[6]=groupCountFilter +report[7]=groupMeanLowerBound +report[8]=paste(pCutoffs,collapse=',') +report[9]=paste(fcCutoffs,collapse=',') +report[10]=normalization +if (normalization=='ByUpperQuantile') + report[10]=paste0(report[10],'-',uqProb*100,'%') +if (logTransform) { + report[11]='Yes' +} else { + report[11]='No' +} +report[12]=logOffset +report[13]=statTest +report[14]=dataImputation +report[15]=formulaStr +if (visual) { + report[16]='Yes' +} else { + report[16]='No' +} +report[17]=volcanoP +if (hclustSamples) { + report[18]='Yes' +} else { + report[18]='No' +} +if (hclustFeatures) { + report[19]='Yes' +} else { + report[19]='No' +} +writeData(wb,'Overall',report,rowNames=T,borders='surrounding') +writeData(wb,'Overall',verNum) +writeComment(wb,'Overall',2,2,createComment('Total number of samples',visible=F)) +writeComment(wb,'Overall',2,3,createComment('Total number of groups',visible=F)) +writeComment(wb,'Overall',2,4,createComment('Total number of features e.g. metabolites, proteins, etc.',visible=F)) +writeComment(wb,'Overall',2,5,createComment('Total number of features with at least one sample that\'s greater than the cutoff; these features are considered \"expressed\"',visible=F)) +writeComment(wb,'Overall',2,6,createComment('Feature value cutoff; only feature values greater than this number are considered \"expressed\"',visible=F)) +writeComment(wb,'Overall',2,7,createComment('Within-group count cutoff; only features expressed by more than this number of replicates in a group will be included for comparisons',visible=F)) +writeComment(wb,'Overall',2,8,createComment('Group mean values less than this number will be replaced by it',visible=F)) +writeComment(wb,'Overall',2,9,createComment('Adjusted p-values less than this number will be considered significant',visible=F)) +writeComment(wb,'Overall',2,10,createComment('Fold changes greater than this number will be considered significant',visible=F)) +writeComment(wb,'Overall',2,11,createComment('If Yes, non-zero values are used to calculate a geometric mean per feature, feature values are divided by the geometric mean to obtain ratios, and the median of such ratios per sample is used as an estimation of the size factor per sample. Feature values in a sample are divided by this sample-specific size factor.',visible=F)) +writeComment(wb,'Overall',2,12,createComment('If yes, data are log2-transformed with or without an offset.',visible=F)) +writeComment(wb,'Overall',2,13,createComment('If >0, data=data+offset before the log2-transformation, otherwise only the positive values are log2-transformed but zeros are unchanged.',visible=F)) +writeComment(wb,'Overall',2,14,createComment('Statistical test method.',visible=F)) +writeComment(wb,'Overall',2,15,createComment('Data imputation method.',visible=F)) +writeComment(wb,'Overall',2,16,createComment('Formula for the linear model such as cpglm, glm etc..',visible=F)) +writeComment(wb,'Overall',2,17,createComment('If yes, figures are saved in the result file as well as in a Figure folder.',visible=F)) +writeComment(wb,'Overall',2,18,createComment('Y axis used for the volcano plot.',visible=F)) +writeComment(wb,'Overall',2,19,createComment('If yes, samples on the heatmaps are hierarchical-clustered.',visible=F)) +writeComment(wb,'Overall',2,20,createComment('If yes, features on the heatmaps are hierarchical-clustered.',visible=F)) + +report=matrix(nrow=nGroups,ncol=2) +colnames(report)=c('Samples','Value-count-filtered features') +rownames(report)=groups +sumGroupCount=matrix(nrow=nFeatures,ncol=nGroups) +colnames(sumGroupCount)=groups +for (i in 1:nGroups) { + groupCount=featureExpressed[,groupID==groups[i],drop=F] + sumGroupCount[,i]=rowSums(groupCount,na.rm=T)>=groupCountFilter + report[groups[i],1]=dim(groupCount)[2] + report[groups[i],2]=sum(sumGroupCount[,i]) + groupCount[is.na(groupCount)]=F + if (visual&sum(groupID==groups[i])<=6) { + ggsave(file.path(figuresPath,paste0('euler_',i,'.ps')),plot=plot(euler(groupCount),quantities=T,col=colors[1:dim(groupCount)[2]],fill=F,fill_opacity=0.7,legend=T)) + ggsave(file.path(figuresPath,paste0('euler_',i,'.png')),plot=plot(euler(groupCount),quantities=T,col=colors[1:dim(groupCount)[2]],fill=F,fill_opacity=0.7,legend=T)) + } +} +if (visual&nGroups<=6) { + ggsave(file.path(figuresPath,'euler_all.ps'),plot=plot(euler(sumGroupCount),quantities=T,col=colors[1:nGroups],fill=F,fill_opacity=0.7,legend=T)) + ggsave(file.path(figuresPath,'euler_all.png'),plot=plot(euler(sumGroupCount),quantities=T,col=colors[1:nGroups],fill=F,fill_opacity=0.7,legend=T)) +} +writeData(wb,'Overall',report,rowNames=T,startCol=5,borders='surrounding') +writeData(wb,'Overall','Group',startCol=5) +writeComment(wb,'Overall',5,1,createComment('Group name; a group is a set of replicates to be used in a statistical comparison',visible=F)) +writeComment(wb,'Overall',6,1,createComment('Number of replicates in a group',visible=F)) +writeComment(wb,'Overall',7,1,createComment('Number of features that pass both the value and the within-group count filters',visible=F)) + +report=matrix(nrow=nSamples,ncol=1) +colnames(report)='Value-filtered features' +rownames(report)=colnames(data) +for (i in 1:nSamples) + report[i]=sum(featureExpressed[,i],na.rm=T) +writeData(wb,'Overall',report,rowNames=T,startCol=10,borders='surrounding') +writeData(wb,'Overall','Sample',startCol=10) +writeComment(wb,'Overall',10,1,createComment('Sample name',visible=F)) +writeComment(wb,'Overall',11,1,createComment('Number of features that pass the value filter',visible=F)) +addStyle(wb,'Overall',style1,rows=1,cols=1:100) +setColWidths(wb,'Overall',1:100,widths='auto') +freezePane(wb,'Overall',firstActiveRow=2) + +if (visual) { + if (nGroups<=6) + insertImage(wb,'Overall',file.path(figuresPath,'euler_all.png'),width=4,height=4,unit='in',startRow=nSamples+4,startCol=7) + for (i in 1:nGroups) { + if (sum(groupID==groups[i])<=6) + insertImage(wb,'Overall',file.path(figuresPath,paste0('euler_',i,'.png')),width=4,height=4,unit='in',startRow=nSamples+4+i*20,startCol=7) + } +} + +#############################End Sheet "Overall"############################# + +#############################Begin Sheet "Group_Features"############################# +message('Writing the group features sheet ...') +report=matrix(nrow=nFeatures,ncol=nGroups) +colnames(report)=groups +rownames(report)=rownames(data) +for (i in groups) + report[,i]=rowSums(featureExpressed[,groupID==i,drop=F],na.rm=T) +groupExpressed=report>=groupCountFilter +if (dataFormat=='table with feature description') { + report=data.frame(Description=featureDescription,Groups=rowSums(groupExpressed),report,stringsAsFactors=F,check.names=F) + freezePane(wb,'Group_Features',firstActiveRow=2,firstActiveCol=3) +} else { + report=data.frame(Groups=rowSums(groupExpressed),report,stringsAsFactors=F,check.names=F) + freezePane(wb,'Group_Features',firstActiveRow=2,firstActiveCol=2) +} +writeData(wb,'Group_Features',report,rowNames=T) +writeData(wb,'Group_Features','Feature') +writeComment(wb,'Group_Features',1,1,createComment('Feature name',visible=F)) +writeComment(wb,'Group_Features',2,1,createComment('Each cell in this column tells the number of groups that pass the count filter for a specific feature',visible=F)) +writeComment(wb,'Group_Features',3,1,createComment('Each cell in this or the succeeding columns tells the number of replicates in a group that pass the value filter for a specific feature',visible=F)) +addStyle(wb,'Group_Features',style2,rows=1,cols=1:100) +if (nSelectedFeatures>0) { + report=report[selectedFeatures,] + writeData(wb,'Selected_Features',report,rowNames=T,startRow=nSelectedFeatures+4) + writeData(wb,'Selected_Features','Feature',startRow=nSelectedFeatures+4) + addStyle(wb,'Selected_Features',style2,rows=nSelectedFeatures+4,cols=1:200) +} +#############################End Sheet "Group_Features"############################# + +#############################Begin Sheet "Sample_Features"############################# +message('Writing the sample features sheet ...') +report=featureExpressed +class(report)='integer' +if (dataFormat=='table with feature description') { + report=data.frame(Description=featureDescription,Samples=rowSums(featureExpressed),report,stringsAsFactors=F,check.names=F) + freezePane(wb,'Sample_Features',firstActiveRow=2,firstActiveCol=3) +} else { + report=data.frame(Samples=rowSums(featureExpressed),report,stringsAsFactors=F,check.names=F) + freezePane(wb,'Sample_Features',firstActiveRow=2,firstActiveCol=2) +} +writeData(wb,'Sample_Features',report,rowNames=T) +writeData(wb,'Sample_Features','Feature') +writeComment(wb,'Sample_Features',1,1,createComment('Feature name',visible=F)) +writeComment(wb,'Sample_Features',2,1,createComment('Each cell in this column tells the number of samples that pass the value filter for a specific feature',visible=F)) +writeComment(wb,'Sample_Features',3,1,createComment('Each cell in this or the succeeding columns tells whether the value passes the value filter for a specific feature; 1 means yes and 0 means no',visible=F)) +addStyle(wb,'Sample_Features',style2,rows=1,cols=1:200) +if (nSelectedFeatures>0) { + report=report[selectedFeatures,] + writeData(wb,'Selected_Features',report,rowNames=T,startRow=2*nSelectedFeatures+7) + writeData(wb,'Selected_Features','Feature',startRow=2*nSelectedFeatures+7) + addStyle(wb,'Selected_Features',style2,rows=2*nSelectedFeatures+7,cols=1:200) +} +#############################End Sheet "Sample_Features"############################# + +comparisonsSummary=matrix(nrow=0,ncol=11) +colnames(comparisonsSummary)=c('Control','Treatment','Control_Only','Treatment_Only','Absent_in_Both','Present_in_Both','Up_in_Treatment','Down_in_Treatment','Significant','Adjusted_p_Cutoff','FC_Cutoff') +options(warn=2) +if (dataImputation!='no') + data[data<=valueFilter]=NA +if (dataImputation=='value filter') + data[is.na(data)]=valueFilter +if (dataImputation=='mean') + tempData=rowMeans(data,na.rm=T) +if (dataImputation=='median') + tempData=rowMedians(data,na.rm=T) +if (dataImputation=='half minimum') { + tempData=rowMins(data,na.rm=T) + if (logTransform) + tempData=tempData-1 + else + tempData=tempData/2 +} +if (dataImputation=='mean' || dataImputation=='median' || dataImputation=='half minimum') { + tempData[tempData==Inf]=groupMeanLowerBound + tempData[is.nan(tempData)]=groupMeanLowerBound + tempData[is.na(tempData)]=groupMeanLowerBound + tempData=replicate(dim(data)[2],tempData) + data[is.na(data)]=tempData[is.na(data)] +} +#############################Begin Sheet "Control_vs_Others"############################# +message ('Writing the differential expression result for control vs. other groups ...') +nControls=length(controlGroups) +if (nControls>0) { + for (k in 1:nControls) { + report=matrix(NA,nrow=nFeatures,ncol=(nGroups-1)*6+2) + controlGroup=controlGroups[k] + treatGroups=groups[groups!=controlGroup] + group=C(groupID,base=which(rownames(contrasts(groupID))==controlGroup)) + colnames(contrasts(group))=rownames(contrasts(group))[rownames(contrasts(group))!=controlGroup] + if (logTransform) + tempStr=c('_log2_Mean','_log2_Adjusted_Mean') + else + tempStr=c('_Mean','_Adjusted_Mean') + colnames(report)=c(paste0(controlGroup,tempStr),as.vector(t(outer(treatGroups,c(tempStr,'_log2_FC','_FC','_p_Value','_Adjusted_p'),paste0)))) + rownames(report)=rownames(data) + for (i in 1:nFeatures) { + group1=data[i,groupID==controlGroup] + if (sum(!is.na(group1))==0) + group1=rep(groupMeanLowerBound,length(group1)) + if(sum(!is.na(group1))==1) + sd1=0 + else + sd1=sd(group1,na.rm=T) + report[i,1]=mean(group1,na.rm=T) + report[i,2]=report[i,1] + if (report[i,2]<groupMeanLowerBound) + report[i,2]=groupMeanLowerBound + for (j in 1:(nGroups-1)) { + group2=data[i,groupID==treatGroups[j]] + if (sum(!is.na(group2))==0) + group2=rep(groupMeanLowerBound,length(group2)) + if(sum(!is.na(group2))==1) + sd2=0 + else + sd2=sd(group2,na.rm=T) + report[i,(j-1)*6+3]=mean(group2,na.rm=T) + report[i,(j-1)*6+4]=report[i,(j-1)*6+3] + if (report[i,(j-1)*6+4]<groupMeanLowerBound) + report[i,(j-1)*6+4]=groupMeanLowerBound + if (logTransform) + report[i,(j-1)*6+5]=report[i,(j-1)*6+4]-report[i,2] #Log2 FC + else + report[i,(j-1)*6+5]=log2(report[i,(j-1)*6+4]/report[i,2]) #Log2 FC + report[i,(j-1)*6+6]=2^report[i,(j-1)*6+5] + if ((groupExpressed[i,controlGroup] || groupExpressed[i,treatGroups[j]])) + if (statTest=='t-test') { + if (sd1==0&&sd2!=0) + report[i,(j-1)*6+7]=t.test(group2,mu=mean(group1,na.rm=T))$p.value + else if (sd1!=0&&sd2==0) + report[i,(j-1)*6+7]=t.test(group1,mu=mean(group2,na.rm=T))$p.value + else if (sd1!=0&&sd2!=0) + report[i,(j-1)*6+7]=t.test(group1,group2)$p.value + } else if (statTest=='cpglm-glm') { + testData=list(value=data[i,],group=group,rep=repID,batch=batchID) + fit=try(cpglm(formulaForFitting,data=testData,link='identity'),silent=T) + if (class(fit)=='cpglm'&&fit$converged) + coefTable=try(coef(summary(fit),silent=T)) + if (class(fit)=='cpglm'&&fit$converged&&class(coefTable)!='try-error') { + rowName=paste0("group",treatGroups[j]) + if (rowName %in% rownames(coefTable)) + report[i,(j-1)*6+7]=coefTable[rowName,4] + } else { + fit=try(glm(formulaForFitting,data=testData,family=gaussian(link='identity')),silent=T) + if (class(fit)=='glm'&&fit$converged) { + rowName=paste0("group",treatGroups[j]) + coefTable=coef(summary(fit)) + if (rowName %in% rownames(coefTable)) + report[i,(j-1)*6+7]=coefTable[rowName,4] + } + } + } else if (statTest=='glm') { + testData=list(value=data[i,],group=group,rep=repID,batch=batchID) + fit=try(glm(formulaForFitting,data=testData,family=gaussian(link='identity')),silent=T) + if (class(fit)=='glm'&&fit$converged) { + rowName=paste0("group",treatGroups[j]) + coefTable=coef(summary(fit)) + if (rowName %in% rownames(coefTable)) + report[i,(j-1)*6+7]=coefTable[rowName,4] + } + } + + } + } + for (j in 1:(nGroups-1)) { + report[,(j-1)*6+8]=p.adjust(report[,(j-1)*6+7],method='fdr') + } + report=data.frame(report,data[,groupID==controlGroup,drop=F],stringsAsFactors=F,check.names=F) + for (j in 1:(nGroups-1)) + report=data.frame(report,data[,groupID==treatGroups[j],drop=F],stringsAsFactors=F,check.names=F) + sheet = addWorksheet(wb, paste0(controlGroup,'_vs_Others')) + if (dataFormat=='table with feature description') { + writeData(wb,sheet,data.frame(Description=featureDescription,report),rowNames=T,keepNA=F,stringsAsFactors=F,check.names=F) + writeData(wb,sheet,'Feature') + freezePane(wb,sheet,firstActiveRow=2,firstActiveCol=3) + } else { + writeData(wb,sheet,report,rowNames=T,keepNA=F) + writeData(wb,sheet,'Feature') + freezePane(wb,sheet,firstActiveRow=2,firstActiveCol=2) + } + addStyle(wb,sheet,style2,rows=1,cols=1:200) + for (j in 1:(nGroups-1)) { + expressed1=groupExpressed[,controlGroup] + expressed2=groupExpressed[,treatGroups[j]] + fc=report[,(j-1)*6+6] + rfc=1/fc + pRaw=report[,(j-1)*6+7] + p=report[,(j-1)*6+8] + for (pc in pCutoffs) + for (fcc in fcCutoffs) + comparisonsSummary=rbind(comparisonsSummary,c(controlGroup,treatGroups[j],sum(expressed1&!expressed2),sum(!expressed1&expressed2),sum(!expressed1&!expressed2),sum(expressed1&expressed2),sum(!is.na(p)&fc>fcc&p<pc),sum(!is.na(p)&rfc>fcc&p<pc),sum(!is.na(p)&(fc>fcc|rfc>fcc)&p<pc),pc,fcc)) + if (grepl('raw',volcanoP)) + tempData=data.frame(log2(fc),-log(pRaw,10),F) + else + tempData=data.frame(log2(fc),-log(p,10),F) + rownames(tempData)=rownames(report) + tempData=tempData[rowSums(is.na(tempData))==0,] + colnames(tempData)=c('x','y','Selected') + options(warn=-1) + if (visual) { + plotType='volcano' + plotTitle=paste0('Volcano Plot of ',controlGroup,' vs. ',treatGroups[j]) + if (nSelectedFeatures>0) + tempData[,'Selected']=rownames(tempData) %in% interestedFeatures + else + tempData$Selected=F + p1=ggplot(data.frame(tempData),aes(x=x,y=y,color=Selected))+scale_color_manual(values=c('black','green'))+geom_point()+geom_text(aes(label=ifelse(Selected,rownames(tempData),'')),hjust=0, vjust=0)+xlab('Log2(FC)')+ylab(paste0('-Log10(',volcanoP,')'))+ggtitle(plotTitle)+theme(legend.position="none") + for (pc in pCutoffs) { + p1=p1+geom_hline(yintercept=-log(pc,10),linetype="dashed",color = "green",size=0.2) + for (fcc in fcCutoffs) { + p1=p1+geom_vline(xintercept=log2(fcc),linetype="dashed",color = "green",size=0.2) + p1=p1+geom_vline(xintercept=-log2(fcc),linetype="dashed",color = "green",size=0.2) + } + } + plotFile=file.path(figuresPath,paste0(controlGroup,'_vs_',treatGroups[j],'_',plotType,'.ps')) + plotFile2=file.path(figuresPath,paste0(controlGroup,'_vs_',treatGroups[j],'_',plotType,'.png')) + ggsave(plotFile,plot=p1,width=7,height=7,units='in') + ggsave(plotFile2,plot=p1,width=7,height=7,units='in') + insertImage(wb,'Visual',plotFile2, width=7, height=7, unit='in',startRow=plotID*45+1,startCol=(j-1)*13+1) + } + } + if (visual) + plotID=plotID+1 + options(warn=0) + if (nSelectedFeatures>0) { + if (dataFormat=='table with feature description') { + report=data.frame(Description=featureDescription,report,stringsAsFactors=F,check.names=F)[selectedFeatures,] + writeData(wb,'Selected_Features',report,rowNames=T,startRow=(k+2)*(nSelectedFeatures+3)+1,keepNA=F) + writeData(wb,'Selected_Features','Feature',startRow=(k+2)*(nSelectedFeatures+3)+1) + addStyle(wb,'Selected_Features',style2,rows=(k+2)*(nSelectedFeatures+3)+1,cols=1:200) + } else { + report=report[selectedFeatures,] + writeData(wb,'Selected_Features',report,rowNames=T,startRow=(k+2)*(nSelectedFeatures+3)+1,keepNA=F) + writeData(wb,'Selected_Features','Feature',startRow=(k+2)*(nSelectedFeatures+3)+1) + addStyle(wb,'Selected_Features',style2,rows=(k+2)*(nSelectedFeatures+3)+1,cols=1:200) + } + } + } +} +#############################End Sheet "Control_vs_Others"############################# + +#############################Begin Sheet "DE_Comparisons"############################# +message ('Writing the differential expression results for groups ...') +nComparisons=dim(dePairs)[1] +if (nComparisons>0) { + for (k in 1:nComparisons) { + report=matrix(NA,nrow=nFeatures,ncol=8) + controlGroup=dePairs[k,1] + treatGroup=dePairs[k,2] + group=C(groupID,base=which(rownames(contrasts(groupID))==controlGroup)) + colnames(contrasts(group))=rownames(contrasts(group))[rownames(contrasts(group))!=controlGroup] + if (logTransform) + tempStr=c('_log2_Mean','_log2_Adjusted_Mean') + else + tempStr=c('_Mean','_Adjusted_Mean') + colnames(report)=c(paste0(controlGroup,tempStr),paste0(treatGroup,tempStr),paste0(paste0(controlGroup,'_vs_',treatGroup),c('_log2_FC','_FC','_p_Value','_Adjusted_p'))) + rownames(report)=rownames(data) + for (i in 1:nFeatures) { + group1=data[i,groupID==controlGroup] + if (sum(!is.na(group1))==0) + group1=rep(groupMeanLowerBound,length(group1)) + if(sum(!is.na(group1))==1) + sd1=0 + else + sd1=sd(group1,na.rm=T) + report[i,1]=mean(group1,na.rm=T) + report[i,2]=report[i,1] + if (report[i,2]<groupMeanLowerBound) + report[i,2]=groupMeanLowerBound + group2=data[i,groupID==treatGroup] + if (sum(!is.na(group2))==0) + group2=rep(groupMeanLowerBound,length(group2)) + if(sum(!is.na(group2))==1) + sd2=0 + else + sd2=sd(group2,na.rm=T) + report[i,3]=mean(group2,na.rm=T) + report[i,4]=report[i,3] + if (report[i,4]<groupMeanLowerBound) + report[i,4]=groupMeanLowerBound + if (logTransform) + report[i,5]=report[i,4]-report[i,2] + else + report[i,5]=log2(report[i,4]/report[i,2]) + report[i,6]=2^report[i,5] + if ((groupExpressed[i,controlGroup] || groupExpressed[i,treatGroup])) + if (statTest=='t-test') { + if (sd1==0&&sd2!=0) + report[i,7]=t.test(group2,mu=mean(group1,na.rm=T))$p.value + else if (sd1!=0&&sd2==0) + report[i,7]=t.test(group1,mu=mean(group2,na.rm=T))$p.value + else if (sd1!=0&&sd2!=0) + report[i,7]=t.test(group1,group2)$p.value + } else if (statTest=='cpglm-glm') { + testData=list(value=data[i,],group=group,rep=repID,batch=batchID) + fit=try(cpglm(formulaForFitting,data=testData,link='identity'),silent=T) + if (class(fit)=='cpglm'&&fit$converged) + coefTable=try(coef(summary(fit)),silent=T) + if (class(fit)=='cpglm'&&fit$converged&&class(coefTable)!='try-error') { + rowName=paste0("group",treatGroup) + if (rowName %in% rownames(coefTable)) + report[i,7]=coefTable[rowName,4] + } else { + fit=try(glm(formulaForFitting,data=testData,family=gaussian(link='identity')),silent=T) + if (class(fit)=='glm'&&fit$converged) { + rowName=paste0("group",treatGroup) + coefTable=coef(summary(fit)) + if (exists('coefTable')&&rowName %in% rownames(coefTable)) + report[i,7]=coefTable[rowName,4] + } + } + } else if (statTest=='glm') { + testData=list(value=data[i,],group=group,rep=repID,batch=batchID) + fit=try(glm(formulaForFitting,data=testData,family=gaussian(link='identity')),silent=T) + if (class(fit)=='glm'&&fit$converged) { + rowName=paste0("group",treatGroup) + coefTable=coef(summary(fit)) + if (rowName %in% rownames(coefTable)) + report[i,7]=coefTable[rowName,4] + } + } + } + report[,8]=p.adjust(report[,7],method='fdr') + report=data.frame(report,data[,groupID==controlGroup,drop=F],data[,groupID==treatGroup,drop=F],stringsAsFactors=F,check.names=F) + sheet = addWorksheet(wb, paste0(controlGroup,'_vs_',treatGroup)) + if (dataFormat=='table with feature description') { + writeData(wb,sheet,data.frame(Description=featureDescription,report,stringsAsFactors=F,check.names=F),rowNames=T,keepNA=F) + writeData(wb,sheet,'Feature') + freezePane(wb,sheet,firstActiveRow=2,firstActiveCol=3) + } else { + writeData(wb,sheet,report,rowNames=T,keepNA=F) + writeData(wb,sheet,'Feature') + freezePane(wb,sheet,firstActiveRow=2,firstActiveCol=2) + } + addStyle(wb,sheet,style2,rows=1,cols=1:100) + expressed1=groupExpressed[,controlGroup] + expressed2=groupExpressed[,treatGroup] + fc=report[,6] + rfc=1/fc + pRaw=report[,7] + p=report[,8] + for (pc in pCutoffs) + for (fcc in fcCutoffs) + comparisonsSummary=rbind(comparisonsSummary,c(controlGroup,treatGroup,sum(expressed1&!expressed2),sum(!expressed1&expressed2),sum(!expressed1&!expressed2),sum(expressed1&expressed2),sum(!is.na(p)&fc>fcc&p<pc),sum(!is.na(p)&rfc>fcc&p<pc),sum(!is.na(p)&(fc>fcc|rfc>fcc)&p<pc),pc,fcc)) + if (grepl('raw',volcanoP)) + tempData=data.frame(log2(fc),-log(pRaw,10),F) + else + tempData=data.frame(log2(fc),-log(p,10),F) + rownames(tempData)=rownames(report) + tempData=tempData[rowSums(is.na(tempData))==0,] + colnames(tempData)=c('x','y','Selected') + options(warn=-1) + if (visual) { + plotType='volcano' + plotTitle=paste0('Volcano Plot of ',controlGroup,' vs. ',treatGroup) + if (nSelectedFeatures>0) + tempData[,'Selected']=rownames(tempData) %in% interestedFeatures + else + tempData$Selected=F + p1=ggplot(data.frame(tempData),aes(x=x,y=y,color=Selected))+scale_color_manual(values=c('black','green'))+geom_point()+geom_text(aes(label=ifelse(Selected,rownames(tempData),'')),hjust=0, vjust=0)+xlab('Log2(FC)')+ylab(paste0('-Log10(',volcanoP,')'))+ggtitle(plotTitle)+theme(legend.position="none") + for (pc in pCutoffs) { + p1=p1+geom_hline(yintercept=-log(pc,10),linetype="dashed",color = "green",size=0.2) + for (fcc in fcCutoffs) { + p1=p1+geom_vline(xintercept=log2(fcc),linetype="dashed",color = "green",size=0.2) + p1=p1+geom_vline(xintercept=-log2(fcc),linetype="dashed",color = "green",size=0.2) + } + } + plotFile=file.path(figuresPath,paste0(controlGroup,'_vs_',treatGroup,'_',plotType,'.ps')) + plotFile2=file.path(figuresPath,paste0(controlGroup,'_vs_',treatGroup,'_',plotType,'.png')) + ggsave(plotFile,plot=p1,width=7,height=7,units='in') + ggsave(plotFile2,plot=p1,width=7,height=7,units='in') + insertImage(wb,'Visual',plotFile2, width=7, height=7, unit='in',startRow=plotID*45+1,startCol=(k-1)*13+1) + } + if (nSelectedFeatures>0) { + if (dataFormat=='table with feature description') { + report=data.frame(Description=featureDescription,report,stringsAsFactors=F,check.names=F)[selectedFeatures,] + writeData(wb,'Selected_Features',report,rowNames=T,startRow=(k+2+nControls)*(nSelectedFeatures+3)+1,keepNA=F) + writeData(wb,'Selected_Features','Feature',startRow=(k+2+nControls)*(nSelectedFeatures+3)+1) + addStyle(wb,'Selected_Features',style2,rows=(k+2+nControls)*(nSelectedFeatures+3)+1,cols=1:200) + } else { + report=report[selectedFeatures,] + writeData(wb,'Selected_Features',report,rowNames=T,startRow=(k+2+nControls)*(nSelectedFeatures+3)+1,keepNA=F) + writeData(wb,'Selected_Features','Feature',startRow=(k+2+nControls)*(nSelectedFeatures+3)+1) + addStyle(wb,'Selected_Features',style2,rows=(k+2+nControls)*(nSelectedFeatures+3)+1,cols=1:200) + } + } + } + if (visual) + plotID=plotID-1 + options(warn=0) +} +#############################End Sheet "DE_Comparisons"############################# + +#############################Begin Sheet "Overall" (Part 2)############################# +if (dim(comparisonsSummary)[1]>1) { + comparisonsSummary=as.data.frame(comparisonsSummary[order(comparisonsSummary[,10],comparisonsSummary[,11],comparisonsSummary[,1],comparisonsSummary[,2]),]) + comparisonsSummary[,3:11]=apply(comparisonsSummary[,3:11],2,as.numeric) +} +if (dim(comparisonsSummary)[1]>0) { + writeData(wb,'Overall',comparisonsSummary,startCol=14,rowNames=F,borders='surrounding') + writeComment(wb,'Overall',14,1,createComment('Control group name',visible=F)) + writeComment(wb,'Overall',15,1,createComment('Treatment group name',visible=F)) + writeComment(wb,'Overall',16,1,createComment('Number of features that are expressed only by control',visible=F)) + writeComment(wb,'Overall',17,1,createComment('Number of features that are expressed only by treatment',visible=F)) + writeComment(wb,'Overall',18,1,createComment('Number of features that are expressed by neither groups',visible=F)) + writeComment(wb,'Overall',19,1,createComment('Number of features that are expressed by both groups',visible=F)) + writeComment(wb,'Overall',20,1,createComment('Number of features that are significantly up-regulated by treatment',visible=F)) + writeComment(wb,'Overall',21,1,createComment('Number of features that are significantly down-regulated by treatment',visible=F)) + writeComment(wb,'Overall',22,1,createComment('Number of features that are significantly up- or down-regulated by treatment',visible=F)) +} +#############################End Sheet "Overall" (Part 2)############################# + +saveWorkbook(wb,outputFile,overwrite=F) +unlink('Rplots.pdf') + +#############################End Main############################# +message('Analysis completed.') +message('###################################################################') + diff --git a/V1.7.1/ReadMe.txt b/V1.7.1/ReadMe.txt new file mode 100755 index 0000000..9aa68e0 --- /dev/null +++ b/V1.7.1/ReadMe.txt @@ -0,0 +1,25 @@ +To run this analysis tool, first make a copy of the data template, read the instructions in the template and make necessary changes. + +To run the tool on BioHPC @ UTSW: Open a BioHPC web visualization node and run the following from command line: + +sh /project/CRI/shared/Metabolomics/run_analysis.sh input_path/input_file_name.xlsm output_path/output_file_name.xlsx optional_BioHPC_queue_name #input_file is the copy of the data template with your data. output_file is the analysis result. + +######################################################################### +Change Log: +#v1.7.1, 12/04/2019: Figures were saved as both .ps and .png files. Euler plot was disabled for >6 samples. Samples / features can be clustered or not on the heatmaps. +#v1.7.0, 10/22/2019: Figures were saved as .ps files in a Figure folder. Correlation plot was genrated with the old method for better speed. Euler plot was disabled for >10 samples. Y axis of volcano plot could be raw or adjusted p-values. DESeq was renamed RLE as I would implement the DESeq differential expression analysis in the future. RLE is the same method used by DESeq but I slightly modified it to allow for zeros, missing values and / or values below the value filter. Due to this reason I couldn't use the DESeq normalization directly. +#v1.6.9, 10/4/2019: A bug regarding cpglm was fixed. The bug caused an error in cpglm fitting but the class(fit) was still 'cpglm' and fit$converged was still TRUE. Visualization option was added. Version number checkup was added for the data template. Multiple other minor improvements. +#v1.6.8, 09/19/2019: A bug regarding featureDescription was fixed. The bug was present only when data format was "table with feature description" and samples were excluded which happened to cause some features #v1.6.7, 08/29/2019: Correlation plot was enchanced. Euler diagram was added. Software version number was reported in the result file. +#v1.6.6, 08/28/2019: Data table with feature description column was enabled, minor bugs fixed, and volcano plot colors customized. +#v1.6.5, 08/27/2019: A bug regarding number of features per sample when data imputation = "NA" was fixed. +#v1.6.4, 08/08/2019: A bug regarding groupMeanLowerBound was fixed. Definition and default value of groupCountFilter were modified. Data imputation was implemented. Missing values were allowed. +#v1.6.3, 08/05/2019: a bug regarding nSamples and formula string generation (only when technical replicates are present) was fixed. +#v1.6.2, 07/24/2019: a bug regarding normalization size factors was fixed. The bug caused size factors to be applied to incorrect samples. If you have used versions before this fix, please rerun your analysis using this fixed version. +#v1.6.1, 07/03/2019: a bug regarding numerical batch IDs was fixed. +#v1.6, 06/24/2019: Multiple normalization methods; visualization plots; sample exclusion; feature exclusion; cpglm and glm for statistical tests; +#v1.5, 03/13/2019: If logTransform enabled, data=log2(data+logOffset) before statistical tests; p-values are FDR-adjusted; Bugs fixed; normalization enabled; +#v1.4, 01/18/2019: Added functionality. (1) Averaging technical replicates. (2) Auto-IDs. +#v1.1-1.3: Bugs fixed. +#v1.0: Basic analysis including counting features and t-tests. +######################################################################### + diff --git a/V1.7.1/data_template_1.7.1.xlsm b/V1.7.1/data_template_1.7.1.xlsm new file mode 100755 index 0000000000000000000000000000000000000000..67dd2c1d744578a2f2ad102812218d7fbe42a7af GIT binary patch literal 96603 zcmeFa2V9fc(mx)gNmr?YfP#S15fG)S6qRnFN>@RpN^c1wO(~)xO`3>^f<O=vsiBB~ zfGCJW5F%Y5bVz`Z<o`Smmc4g(fA8MAdwutRU;LOcKIN1-bLPy<`JRcP4kZ;Ej0UzD z27~Q_-NWc2Mk!#hR(cqW6}EZP8Fg0|&+{&x)~Eg4&U;u&`Z_!9yic`h&pp^C;Q5U| z|A#ryn$YLMEwd$TMr?`KxH|uX@I@XYb{)yPi%(<jEK)@K*jTQK#O|8K@y>)l%wIlJ zbh3;hIC&v$Zg6fprXYW>@s<+-CPiYZcGL$(9H&Ov33H1jYcJ__JGwf2Y1w;WPkWnt z6^3&3VB(K@Y#Z+0sdLz6C~qU=&R_HHb?Z?vR(Q0>)(7TI2d_=O`=YgHX4{L?51-Tu z4=E)zXCH4M-mRWa)h#*OkT~tL-AHEp)py#$`$UZkNVF$cKkAA-sxmjfJ<27}#LaW% z923h4RvFbJ+V8pc6~%3@wNd+`W0ER*<}CWlm6=CNDxvK8sh#Xy7|j)36)eTMMMO<h z(ad6HRzarkO&NGx!hVw{F}}FVGE0bBmZ;&XQqwx~;DIR`4>6-fyuLou+q>^wDyV6X zISWPNrY5L{G33Yfu}VH2^!UZg@Lj6?*4JWLI0uJ@!j@mG7^rTGG1@XUxpFi+{-C?G zv@_%9nm+IQjZHDphHApx=!IIQwa4(65%c>of(I5a^h~5*<M$U>ddZ@4&3`Z=T}@=r z<kZR*qJiOc{@7if?B4G}PzfU~e%)P1$3sIYPr@SKl0BU}wF)FD=??31BgSGvG_(yI z_#K+&o`jc6osp373)o*B_mP6rnImltnBQbFCCu<2)B~$k_4!SJdcXj8836TQ?S9_L zLrM~)8*1VIurByT*Pq<h)1j2v5<YX4ysNcJsaI5>!%0f>i{3dRg-kNhMj`p;gK{EX z>#%3$b3XsnLZ5WM<R$y3S-x3Y0=%~+XwqFxwePGjuG1mDctN{0&Rr`8d#o;4z$te! zZ%H;rbR-kKq#F5T@P4SI5cvw1*HPn!t14v{7Y~@myK|W@YO~}DS;UX<hGV!5+vFPF z9NVnqL3v+wxpJj)zk*)j#TI|p=4?T+n-4YK+%#Lh`SHlag8Y5s%<SE@NW|;$Yp)3c zXDRta%!wIaW*pz0JgFrVc9$}3bOn2W`dQ@zO%K+fPsv(=W73uTSaM8Vq?+&JWr|eJ zzA<T6*!HIX=_<g5Ur!m?YVO2OlrY!{Iv9)vc*567%E#5+(azP?5nN3l7(2RV9bh9? zcaUf&(hAvH#ppzN!?wu9+k5Bg=ENJvi_+>Ay>@iY$;=%^&B%w|c~G3^B65jhC?Ixb z=u=!O+^OYAzyk|pIsFz@`CD3f5A@Bv`0ve_6sc<7z4<`<_z4raLk6es#g3uK-p}(_ z>dOoG<zSCAuDI{z!ANV|?!0}ivpka>drn^JdZ-XeL8E@}aZUZl;tcis^xF5(<r&6) zUMRJBZ}w^%cUSYn2l3d$2PCeA&T78lmbn+-eCD%u9smBA^!Mi}w>TP0ue=Wy$k$zq zBIJ&qvpm!z;a!qYpFvAThepY#Qk^NhA!=6Vc*0?x?Ou%P7n7B<oi)q%sHP^lE18;Q z4#3alMi3soz`uD^sgi)-6y9d7Revb|hL|ebOxlHM)%^@<8Rwte>sIlS?~n_+ebOuU zkon<Ii-)=$a)($=Kb-NOz3O_`t)3oTf^?>+51xJI6O?}7RKT`3dVL;(1{ohe9EqSc zDXo8eQL_GmOwkkmrUMw38X^8N>k-{gSYh?;6j7sL4=M$2OHOESy6B_1Y;@Hd9aR~A z`K|Hd4AT(X!!1pl!-B<*2R`LU)EMmv=cHI|K^a{>Z^&?&z-}reAkbVX&=7ZT&0E8M z-rt=>n|Azu$Nrr&s?SL)FKD(1;Avl5X&2NkXfAH=!Nn*aVGR}?cX=gsPX}F{w{3c; zVyxj)^qgn;d&?c=<!2jAGX}%5kB+{K99(HFNnK$%79M~&85h&)vu5~~!hP1@dG$?i zbnQcw?YItC`{;C59^u`2^kmb$JF|D+<@*kVnO~Y>%Iomw*~xq6=KIb8_g?3I<mu5v znjUr`$H~*I=kW(MqGfQOvv*I%cf5d=35(0V$~>uktaj7I-M6r}tZvC)c|;`<uX%;} ze(dU(v&&y5+_zN8(4#sXj%8Ii8Y!F-dpjj1XPi@48mwsIMLir$nAnN9W!JSt$@$cX ziYHn9UVz6WdhNE`0XiCLs=y%F?qWPLR)s^m_#o>&f*@zW@ZgS$mv@?&qMMugw|}f? z#9@woWTPFOZu&HNP4$jg_V&CBxf)cF3g%PVzTF%~e*T=Y*)+NP90eZTei+JTDWa4v zDB*BXyGubv(RxWpd1cEgf1T{v(3Kfl@iWcA%#AfW=_}6ZlT0Os;sl?`?AKt~C8XW; z`Gcl>{X_XvA3KZ@98W7ka;5K?e$iL&OgAfM9WwQqZGB(=*5lB=v%P7bFBEy6m5-fJ zxQY>diB0QJlUJO{47<0pug9c2vxtyr=K4I37#z5gcP#tq0!gZVxqNcgpk3x3+fHZo z!_foXZ&%cQQxUW^N{=Z6YbTBx2IGcprT`Sd`ttcpb?|K|r2s5!;MxD%e=6_jxr9h> z3B$+x)3y-ocG3|~=<YYl#^xsQ)yPo3py2uFV(wxc^(gGv`@Ws~kI8b#W@Ie&d*i#_ zf1&ExEVTE=)h&WIEk=$cu%XTyDEM$~Z#nnol2Gs<Gh=^tw@EMOqhv3&WbJ1(WqUge zEzEjVx7hhoGsy;DI$R!iWIWu=Wa9R-Xnr%g-OTl&U*;H7ty(n*sdiQEeWWLjeWPTD zy!%^r@UfFG9O}8w{t{i_&3R~hr;ur6yVqwfnPe-}T@PYT$I(+1=-|(#=Z^(Fov;1$ zhVE=1Pc5t?@c~&&HQd+supPhdGwEoP<a4%3FLPS$ZIfc~$rib~cVBF^5X++HaqGHS z9F!9nponzH?+=bmm0i3Vax-YfrwYmJn8A!Mwl`~#8!j}ic*J+K=LL~cU9v-wNb_5Z z5R?6yb2}9bma_#0;{t^0>muZF;rw||4=Ly$8E`d9HMR7}+Rp}`(Ir#X2QY35QG12c z-+zNuq%SSar4%k~`dpsk<>rI>J>x+$@l;Kl-W~rCXR45xmnAs<De6u`PPej$m2}=_ zi3CaRm4c9{JJZK5o0hAQ&=z`@IfKWGboQBl&L)J5;X5a~tvMWCr{@X2`a}_W-@N?j zo0hUMqw%d{$K~$}Nq%llf2uq%E289U)a>>&E^WS*Ay5rBAIes6B&is0^Ol(Yz*n#H z)D4D$LtJmRF|Xwb?#-i-Fu+hIieOhxTV=~hQ4C%`7m7M=*xlT@f-^KsP(Cj`V<M-2 zaxA{q0EV)ajb6lSy*~c>f+jwYswsW&QwFIu#qW5Jz|t65l1mQ9h1g|BDSvv_<Bs8b zJ8Ew6wd2Sa)T}KHqashwl#`d2#^}v=SFeewZL{Xg2<q!PGBM2KaVUD(@pe65Ag6nj zrE3X#D;h3yHe!I-QKhvoDHS_ls(ibTw0qgY5`DvPk~O=1u$IGgS>XL-gUbVR+qb;^ zB015Rw^*sCLyg%=9z2+t#lYC+fu>!S0u*r_jIXY`iiikA98zq~b863>@sx|Xe_1(A zLm;s~G<D}56{FNqdF9+-HmWC8p;w-VZ+<8!ha37_wHf(&8$}N{#Y&b9cjgx&*AO>l z4-I9^tHW~10wEFRn;Ld3H($<lUZu*qFc@K<NKOf~P}eoT_VH-4dfRRBd(Nhz{O2eL z`QGn1J&2#R{I$28w3-YR$)nD4pq=M=UM0$ULgL+yU0eJ_E4%BeJSp!n+|R8cK8>_z zn4#W%w@%8_W$U4B`>6&JHXr$rIJXmarqWG)t>Ma9m3zCF%ILQCikbF@p3JJVa9l7> z;c@Q2xL}oX##!~k!il>y7Y$Qmp1Gww5JXuDX_nJnJl=A5XGVUBUtU|x{+VZfg>Cu! z`Jb)q(F=KYOs}P%-p3(3hqog)vPu5KsCi0h&!-Cu`YHIHPnQ-Br^I&1eK~UX;}`nU z;|f~5<$A|VF5N`#*BF@INyGb6_(9a3H+%^~?9<G4CbW9TIJb86gzx!uWwTz&mARMc zS2pV(`>=z9l-<iIJ6u+!p0<>DHasot7QEqf#lj&<uS2UNz1k+XFLOSV<bb~owhLbE z5}0tY<>@@Dc(b+Xc3a=e3upcBQ8v>^ofC=?CaYQdGZM3cscw?=<Hwy_4C!|Y(t5&Y z#tuCwe*NJ51uH>Qxk7xGhO}J4)=M$guO+yB8@aZ*ZrzkU+mfiYgJmbz-R{lW^T-AB zH1!>f!H!hfeTXF{#;FNkAxo}OUKuX`V(uedo3(Rj#xl)#a^G*(zLU7)h-G1GX@#X) zu%;U2$nhNlo_-xW9#%PVcD>;^IH?)#t#YhXx$~JL5A!3%XOoA^l{S|nlg79z_>N@L zo^|$mes{BXSSUO!lrou`ld(7QxQ_#erBPm9xFW*`PM+77MP2gs!(@?&tZ~`?zR%*U zIqF=;{Dat%$0O>Dg=yD>Y-vxP)(^)cmr2JgB$CSSJW@L#A=Jt`;?{C=bgLwWOW&2W z=Lm_bVakDK#K|DIy{6v3btWl8+rhxJ0Ywp3n6q{q=f?Ybv&E|s-oWvMJC_r-Oe^Nl zyOzI`E*1<W>+MoKlStK{q}s$lqO3jcKd_Y@MdYlX`0WMQa1%V6wHSzUc~JkJ65IzV z!L&23ybj*Z7|LO_&Hk=YMw%yX+bZ4jLSYJBFD1J$IA0InU0Shq80GLpv?JH@tmMO6 zAs&hc)uY7S5gu-?`J%ZfLSYbMdFhkd4OKWIXl;(LP_Ga$*3JtL^jYg9#0C*aYYPqv zB=?oIvUcl!#M<~&VIzNRwHmoFt9^_-89`o)ogp0{E#lizJ9yy-8y%LHR|pHJz?J#l zQclT01is~6L}>=<AZhB0V9@-+Tn>vXYL-k~EM>?DB6u$1a#liA*Zgfi^(wuX9hf4- zjwQCwPUQrRdF&vsEf9Q;cHh7dNdsSo5=Y_!6<3$0IIBtRs7mtG!ngxrb{ezd$={Dy zCCs3RBYkrbm}RygIC+vFMZT*xzT&~(4-X)?^2e%G5hn`+G2|UMOtlie-CzxKi>#1l zWj|Wp5rLmk<qc9Lt+Xfvjxw}#V_nIw!^AlIs|gFTD+F?>Y8x_Wp!1zte>Vk!czTFD zqr)vpS}A)~nBlY<gm@aFre@Z#G~SEBkQa$7<bI_UwV;6@<|Y0BYx2-yJ85i<1=FVD zw>ZAzyt4l+j?!WF*~I=-)NH+vTGiO%q>?Ir{L9!(J3ljqPqor><sRt*3CTGVOIog6 z#k3opA^Qie;a5X?1K9)Q{SOx_j6QdNmvIAuAP@?z(|3?BPktd{5Gz-?IVaWodC9NG zM%3EjzSh--ai!!`@)-Y2(CRz`xtF9+PwKDrEqb<MBJc3Ro%E>!qZV|KG&_H8c$ifp z<7(=dK~ejdB5?;Yj)gBXNM&_x<P~9c+nY5FLQcj-bY5OYprj;X{;@D=ExIryR`CGQ zh52C%!CI{z$1})k<g3Hx5!&64@FUdwWZWPhI2y4csz$~=M+h)tJ}8ug;@h&LQpYQ} z4?MX(({5mK)pbI4nP1JKhdf=UHW8QJy5_sDIMj;Y%V_v+lh-XC0R>H$>z_Ev1l{qW zt+uV0Zk1aE0TvCuNabyQ;`18AkDFv~jR{C;SVdND^OKm@8lGKDLeY5o`sa#y`ac0a z1;7Uhe4YZIXTaw<@Oc4z%Eg9<R2EnGR<E$Y>mrnn22$veJhL^_dgJe+%iS*zy>4!O z&3?;Kpj*@B+9$R$A@?)GH=6obpSo)fS2t1KiV%p^;EPh;>L)s{J6zVpa4S)uP{S&! za_d25xz%@Ha^H=Y6<E%_A8&@}M~IGtR7H?~Hq&6)ZS0zl#em97@_@FUljTA53-^7p zD5IpHOV;dV!tTUSi@VlklU-lJ?Oh^2Z7Y*-KRP^3XI3WQZao~_q<-tDK&)o^HDz`` z;dzVUuqK^bCIW?;{nvKzaT6~9-Ja}bH%Q~du6i|`$~jtr0NXK8U~Qs2GNe*Dw)YqH z=XdqQcPdzgq14nQT@l@d>Ps$OONyR~*1m%8M+OWsy$hNG8W*?IH;P%kUI;4_2pV5o z{V;p7+K$*D7uXgy#6FqkmL#=fjFu#ts%kU3A~>OE_~Iq9@*sT~%hK6>@R##LvT%)f zA^6La_ZPS7Z4ZCWv}Pn?eRf>ET+3rkJ%ChS4y*4*af?$QOEk^EPf+$O6rg;F$?t)0 zlw<eDDeZT(SCHH}3Z#eDC?<zUt3okT7h@+ot>A|O&_k<86aKg9OKY56ovQh&_f#fU z&)BU#?2Y4FINp9B=Z=BDlYxJ$h{}*}*xroA=mpy0xODTPNvYml+1tm&fdBNw6`)je z#GZ+$<mZ$8J7nKjpHv-xtTGNg%|3s;JC?jA>-QQeffxwf=f7!C&ar((C$~iBKm=4q zzpH`28dx~+_)Mw(%vfCd<?17<<3wBz9-%t^^_d`l>-Ba-tvql@Roj`d*Z|yb0E;wf z5`h;5VG;~h?E*SD!G;GNUwI()Tm&uy6kAreHtjo>HUvI;!@A=zSRxDX##=@BO#}Z| zEA%T5Sj=-+=7HC&cV#Rd4(Ol-Ybd;FuzF8rJmgyiz~f+_hS()XfcBBD1a>5v=SYFy z!M5t;cZ$F6WOxmA!g_){9dEFDGN2<7EEQ@@i3t4Q`s-l#px1ljs#Cu<VA5C{G3f?A z05ZX2D6ax2vt{bI^&5b3YeFM&)t}NA<G~J?ytoa%Lpguk`l7+A_E+{)cdXxrSdpzV zF7~Yl+9<VV5%@?nwYKVnZK3bl{qU+lG=r<^Zhi&4EFo{0k{`D!^x8^w-gdN&Gl2z; zp6uQeu!4?ARxr89u{VGye7?{uyS1`RyIA+*q(P8F(DNQt>U8^>zFQLKc`Ge(<se%n zo@Yx&tDXFJJ2U(|;zt+T(V0`??><<Y2CXhU7l#vw9QkLm+A3G9?soF8>Y0~nr}$L^ z#qq2<_q=`5$<c3-eHn#$#MGwkJkVwY6q4xN79ls?KKXo9of|(Ec4lE$41wt4oK$)t z^DNQZ`N9WBhIVVyVjG@Sre5ipjylc{Gw5*p=I67c433{>@J{N*7tZ*C4Ix<eO1pM- zlwv9{PkO@<nEkCOC7$v9MEjBbi%JUVL9bKrh{^q}gp%2?rGU~*vq}l&jF!r?Iws@R z_>_^}x2~9QzxITJ4)U-ohnqSom}k`TOu06>93wsuj$lw|B|P_xCwZaFUNRAdT`m}3 zW@tZp(6tX}kLyslJ#XpksA!u32_dM^T2Z<fwv>Yr9|=c@0k4$;uPvrDd~z#KUa%a< zG;5Pk7Wp7P7LHI-XbpHVTeV0i%QSl{QMpogH)D->!{tJPCxdH2$J+1=-c!B!5!5-r z9O+_ed+GwR48KOKcDcae&CqUg(b!uZRm?LwWwA#WRl;+(q^~YT>IuN#PfTqS(8E4( zBF-fV?>XK()91?JtB$G#x@8G;D-LBg$VB9Dz3{=8q5XOrPX)l0l*&xA5eZ=QJywbP zfjfR+J5mv-cRZtgma8~Ir)z@0I;x9jln9g%wAyCp`b8l?y}0L$uf9|NV3^bs%rsM* zzbiDZGYWy~#e-i*MgsHp$qG31D+B~;oM$x6ioQ$*j!4e&^L^O26w_AVIxrWoeJm{< zg$z^Rzio|gM06|cuc!`FFljC%oJCFYq*!9=Qu2noiA9+8=i7|tCrRj;)mQf2rL%$Y z_T36NmKLZaP^xBp;6s(247`&*@Uz8}XM7jJ&f^<N*UO4NE`j}9*kkmy%mt%z+Oapj z?gas!(@z}6<_lbvhsMjM=gN{5Fgm1WHC#n?Q9x(`P-&mKC$Yd4%a^4zbg&$#E`4sd zC4P@3S-cO8UFAquckzb*XsFyY9(eI<WpSwQZ3aiUz}3kc{&j@EHV>E}lRN}V3|JDt zrGrG$!HTbM14B{l4fO$%@v$kW!+@-f<J&(${aq0S_=nA?_{t->9E@6k!W*Tk*>cQR z_5d@*3jd>czdnFHP?7&sFK{E{70+0j)q&UF_-}1Gg!s<>U-s(z<pra=K*PV)BnnJk zS$~lCVdLL<FaHZ}{6F6N&X~%vZw*|Z*FX63wOM0%asYdOa1|WyA3XX6AHQ|O<f8t! zL6{Rg1xduNHYvk4rJ7Oo`EdiH2fVYgb0KyP6u2C<rS>Jg1$`;Y?5cVtc~A<iO%Z(= ztCpDe(D(?UrVW5myabFQFa^Egr%;nA01wvX^D$;g50qy>J}as(_EYTVzx~TsE7X=M zsviS_p=_P~aY^rMF~Br?0SwkdP?kqQ7Ot1RKKs5E4h-zqM*`!gIUU-U+dcwfm^N2Y z+$YtZzSLVBLvUClIpG<e$+}h(UlC;y5zm1t9M(9T@s7_N*7%%@CsW=m>-9{VTm*yx z>`OYU<1d^4qXtO6)*A{=7f_PEHm+?FkO!!{(7?a{abp}1d9Y{e`eUX?zs&;@^Hc<x z@ch(*1;`+v2J)LJsOrKF?CL{OS@mN-h+EK{g+}J<n6L91m<zxdLW8hAh7fZhk)QB^ z|K{@~ANL%;98S`0Kq%K6&EVL-^gC>aMD1HKDnZ{zQf_ADE03W0__Ym?-2E8gf2qxL zq6A1TBIL-7JqQ%G7vqPFMeSHyk@~RITTM)J#=L11N4=@X5GSmi3W5se8qX&MnN5pT zaHC6V3TV9rktxA_^jL2jpG_4!=&G7&THH=#M6e1y_Ogw^rT}hqaZNfcP5_x0Tt$!d zvT@y{q^v6~KOk*1Ani3E9X}vlG9cYGAiXpo%|0kCKPYW9DD5>U9X}{tGU#?HsxL^( zee5`zw;>^D!p1~K#+_ZzSmq?V;AxqC?1Cp{wy_Hu%FMF~>dSPn3F^u`WmA`tVPpMm zQgpI)O9SI6c5NAhhx_i>rR#d_X;4n4KN0D-&*(ne!Z}gRy`#-~68j_b$u31%Xm_QH zKdU}W?986$YwpAXJG)^~ttWO5M718-ofFl1V7FUT>z>{4UaeHS2Ya=W?9S~Cx3k+L z{5vE&X;<X-hoWHD&VyZ>Fa)~B4t6b91L)f6$dj#cDdtXWXIv$H*v=p%F8{3l;Z>9( zCCI)lKdXP3?Y$?QRvN}E8(uaUMVmGpGm18BXgVEztl`+{=tB)nr=kxu96J>)-OzNB zvEikm?(Yz78Fe;UT_cc{uR&IJWr9)yvGO{=%ExTkNg@g*Xrm_tSEG+=D-OGS*kmiX zM+>OERggTveT-ORn^MXOR&-y@G_|)Hk~Ua{5qrWWmNH-)y0a#o8mElp46b6t8rc+5 z2C#BPuDoRl(v>!%XsC<G*>}%0N!RO8k#sVBu$}EbqfHHU=R~O^!yx(e!ejLLO8bqR zR*ug4Oh@_MFqJ==YZof2_0Ud3RO`MSm8e#_-OIgNckM#=Y9;)>yzCJ+j@*795*kT2 zYjNF&`yr|dQ7wq-L)6&K+R4E^3%z<VY>B5~rJKXsMV1viZ?JVnRQS(`tL>#trRnwO zCXVbKZBOn2S-b?Y_!7wC_YjNOv(f%w#m&{c4RxKoBB;UQyrj`<M$xAmHXB9jHe5R$ zt<kXgbhL8AwNue@4VzCzi#Po4cpB<HVUyK>2GtTA)c;3~rhU1pa4oN+GJW;JrhiiD zNXuth0kxHb{VM_cH=z#p?;0pVtbhoe64NeP8rCW*H{xaQRB^h>E@&cS$1Z3jbA(;+ zl#Bqo;0c*c?1Ba|gKUC&GIeZ%I=`#ofE!|fhGPXaoFu3QoWX9cf)&eyYG49lY~)EN zm9i60T#Sx3c(SC}H2k?KQvW^_><9zc5gxE&Rj}e1uwq3}wU*pxJ7eI&mUpJ>#9zcr zup@xPG`|lAD`o*J76dES2P@9s)9@jVr(wd_>$EN&NoO2w(eP)aMOP2($lq(|^3ZQk z%q9JM($^LYl}gjw3Z_#m8L<8~3Y#jnqw{NAXuU;{6~TQBScDDBrT{K<Rt+C5ZZ|R~ zScL)WXCt*KU^_am#)=jvEX}q+vEX_1%42%vzOBNiw)>U~>ONGQ^*oxS8}@nphIF!M zZ6~kozWA#-N_E2zIRLKS2=r%HS}D{=8|jwLH2>#S>mFmD>@2(TB0c1vTo?gO<H`nT z8i-z^guX-c5*73vqL-*$*?fc*BD*ZdA0@dsu|DA{i`<&W^*C~CE?08o*2i4eBey={ z@{8P>!*wxo>q9QH$jl_J$m_pNO!wL+-34c!7o^glDIWm}>{d`<ML>ZyzGFAp9AoEY zqA2>1r<Y1htHAE1sMf!&$sY=X6X1W=Z1FcD$6xKQjl7k|RTOzEmn$>!)?=<)k+&Xk z1xMb>;kq1o>miqI<gM)A?HX_+-Y)zx9oP|0up|4yis9f);RdI?7pQA4Y3_ZV{^>6q zKxoQ`fgEH4E9M6))&eU|11mN^5t%QtF?4_7&>1<6p@<iWo{h~x#LhSQ4#jx{`|eE@ zeCU`O7aDIdWN&aE9oEGLzo~*B9a6(bgWH2_3|67TUbJc76u^g$s<EQMi6Xm$tLPXT z0&4&bW0MQ5_M(*Fnx6|%Q7Y&=L`A8A@9S#D*b7oKI=x1^_>!NZX7JO1nt5!np=SPP zgK<O6Y^a&PRLukxo^T?Cw4pA#dV688<~H4w6<g75H7?ZNN=T*PJ|?V{4W2S!3%b6B zj~aIXDH5#0gf+8irwrJNZmF@N#wj8df-{p?l5&(b=^B0h;E|+z9->|l4TESrM6*Hq z>AhVRVBYH&+9~ULO)nV#oR~sa6#=j*CjJ7OqK2nII=LL8JrG@mD0eM*T@j+@08RV< zg=GN*hYWQ+ei_jKTvABGR3_Mns#l#@jL1=Z&!cJZ@CZb;A$k&|E|-c<Jek%8oS3u! zk`n{S&Yp({e;FkD;WNF@11t-k2_w+2SP})DG=Fx`ZK{GC8LM^IqQLXCuEO)5=30mW zUSW*$)t}@Q{u>Pab@W7A=8t-X6&-+Ai2f^F?El#-tcbM4rcWo4Z9FDFVkfFU{(I3& zbs2tX10SlS+?xC`>EzcnS~7soggn6*_jrpF;9=xafzwdC7Pt%vsEbjKNl{dJ#9!gT z6atYTuV51(DAZGv2gHZqpN)W^jn=Fvhwv#8Uv`GL$Ly7nQ@Vh&9e@Ye`d;KOh+O>7 zYz5rHpFAUeZk*)5ct-sEq5JPuX49hHZ&mta$hD*6%zuTeWZ)Cpe3~YGdUM=kW+*(f z3kuI@L*bchFg#QCxB&<n`eax&FngZL`15f*Dk~t!Jh2JXG9{3fS%tJr45VfBAuS{3 z%if#3;ne@J$oht(wBacIm5vhjB;K=ejW#7msr4s@eEZ@4>=h4ER{w(VB&0$(Awz|! z64=XnPYY_COH`0ua}~5}%xi(*Hn<h?-@8r#*4fWAdN-E%KeW;Ncl1|&^3WC7C067p zjjExQ`)vP;Kn%En^tuMzK<WwxHjvJVQb8sS9c0q5fhG;GcfQ9MKi!#!KL3|EE0EDr z4jC=ZkkO(FX&gaN<B)Z7;AeZioIh^MKN%UpCVZRzp7Fihpu-4V#Mq}k+P)X8`ZgN9 zuM{0w;SH}w%r&ix6cS@!P8<0PzpD4Pv-JzhW=}9OJ!yL|ESoLigy|96mC$U~1U=Jz zwq>F5Tu)xw2+@eim|m~dW;C`#ZVUN(sp4g*BVA%rq9X%Hm_Wh;lC2<N1qpj1dU8at zC?G!l`W!ON)YQuMV%T~1gwv+_wtK_Qvn3dss@jf+o@Y(aHI=me%{riJC1nQ8>x4!} zdO%`c4}JrQ`AhH{NX#3AaQ-XDuvXa)cFIm790?|-=C)>GLhK2rOto#dhY7JI7?>*B zehd|2P0;yOKYnygNB2d7+jYOP=OFn25*$c4E(2sgNQ^IEud-6zR$(%p%*E1sJ#2(M z!PwN;_ITI`Tf#|GEnB*<5!M8K)BU#1p(ERV$*UiAfb7F=V=WM}u(K+V0Es+ER6wE$ z5<QR@S=beB4AFN{`@z1y*$}z15->rlUX7R_-?ege+(*>83>I92jWsi_%Z8tZZ~)yM z0+&UXhQJlj_z<ErIyQtThi(rc?n4)bz`3!dwFrJ}Y%OUAwy;*amHZjeC-eF45!5tS zq9Y@*{H#U6!NtiNm&<706)L?NxvK%)U;}Tc+Rq)s_D^@4V1xYL8;FY+BA?YgrPmP^ z&q3a)+efdnM?3@RS@)DqM_4@Z$Bg~K#;=@BOmdUREJVMsIP*npzhpP}YGh1HNkp~X zkh{cv@Nx=tS#&kBsw{A~cVUOD4BvBBh88jLeB_h5+w?7a#j}yob=&A$M8wmOPIb5G zT7-U2TZqY^g@i41Uc9WdI0F(#kOY7v3M45Y$$RnAN{at_U*4i_bi~xlx_NpXG4VX) zgF1J5oxS2&$k4iZIvo-56y&)&ce-D7=DPy%zE0P@O;Sdxn*$`=AmIav07!&EA}U$b z-E!0pFRH>iWUEw!tgP#xZxI!Lj7+UNN#C+ZJQEpM*Fo1JES`+CsQX8)_seQ6?dnQG zby{0W!ax!Zk_eDo1IcxeMCxFd=5L{^wA)|$qV@mOI1OxDaCZhtK1ilOA~6VF14;hi zwt%XQA^P=N^>-VhICMo2A;Q{cjjUEl^lG%pV8X4~STo{AY^<4Zw#e8JxC&Y!1g?bk z34yDjSwe^h&;}tyWwcZXQ4yU}i{QcHYe`(#_F57*mPnu`yVlBw;T;)fy1Dub<twc{ zd{MYr1k)hV{3SvDI^L0Wrh89+p;2W}>mVS7eX)f=^+PAr{+z_)9qAL3$bm7S2)2QQ z2^2v+ki-CFUMA|VoLYoPa4oh5C-Jlxxaz|jKy^_F+822C4XOWTj%_HW|Hie$*^W#y zKus?bSYoYxnp6Ws2!Ti`8iG5rkcJ?QZmg~5!=}_$^I|J%tNF1JwTK;9KuPdn6KfH? z*iIW16OLdL#E2WS31Tv~BhjMBDC%xDe}6t;_GadcT`n7a+qv+<BIL`spGCp(#YJVw z?mfWW#8@-sZ525aOY8)ne&nB9@F5F=12QbaLBnDSvNxij%fis*$Z9)R$&FS2*De-l zoq_7r86^22nF5IfsA!!*lE1O)|F&`acU$#;mhFrY#_^K@mjOkdUG|E=7W+JHrzO}S zrDzG4kSrna9oWg*>h0Lx+G-vwp%%e~ZLCFbW5;R{+p*oXs!X_E8x#w!(I$u)*KH#- z8srA}8*X};);`fK^5KAMBH6tIaCJa+(ak&Lu6G~05eeSlFZSK-?TVii!|^Lm{5k*o zt%;LGfy5#^4-irSqZK{s)<u`gINpJm^`Oh_@P>Yi4JZ5`VlmL-zY!<d2xWbXZvC-1 zNwpqUjXd7X8$=$Mn`HjzF!R&H7it^fi+_%CZG<lvHX=-4gHC@pLnZ}{BK@!G3+oO0 zVq+g>V;|<1!Jc34(EZI!K)}A}QzMoS&8R6O`gRepY1oC@YA$R`Z8bNx-iDu+;E!~r zA$TKaXb1=-uwQZ+nL$JFLt4`iypTZHb3b|}geZf~2qDUwUSHw=RH#H_aSC!iV%EKq z2*@kh33(-P>o*{;q*0vjdS5RQ*thu8!A>w{;S9wr9HE#69~85=23>~sNb<gj;YK%N z{J+>3{wkdRK|3Y`2$oEL-<SChZl|gOVY5%)C)mag)_-~jD+LXl$17AKk?M!Kd7rlT zH>>hvqiU;nVDoBK8F819g|q~Fq$@4K1vx`YutfqI<RUVImT(?vO-qnNyM_=|(flFA zgJ|oJ_#|>Iuv1zi9|q|NE=WgsKsthHJ#uj!y0Hhm5xBlRQdo!%ntPP-r+^Aaz;MNL zD2gEug)3Mf%!LMcc|R1zFj`~}@Yp~rYye6&043kI?0<!o`Gbbc{vt9d2iQ#|Ml=>G zoy(|-gR7!BL*NI|CLwSx?2wH$3$D!u!;I^*!LZ=!ZPb`?oi=JLxE31}Gwvc1LqmWg z)o2JFNEEQ3kV6PwjaMA$A>f4!1iTo13j!~=*MS$?pc^N^8!Jf-0IaFP+Pa85`e%1j zHw=sooX59ajEx<i|1cB>*iK#e>$X$zlQ$^3*&xf}EM!@<t-B@HAj@JmbR!GAu?~sw zMpfI9>$4Dl4si)>ia=Jz2gu4$fvk+J&}GQVFn(u(ZQFod{u;IVAM_Kx`yvqh?f(=q zc7^=qS)xSlYl$;;g?CK$`P@2vVirNUMC2cSIeVI8VbVCnR_Mz+`}&~#|7K*XpB;G> zGDK)WoJRGs>&L_pxI7wwWq<&Ueii`R^*e5hGZI4!07Y}Y!MdKjTw|E-=A|sy3z;CV zAQR-SI)we)0+}Gkp&MEdDD*4`0OX!Rz#B*nKSi+=Z@?!1IsUqVY~4V%{x~7=2cur? z+X5C`$+MF;$Yk0Ka{1YBpsoOgb;vjYjCV}gV3=@38#PAUkc}D>Zovk{h?`sobpeqN zxo@G4DtN_bvc!vrv);}FQ{45vCE=xx%&y&h{RM|BEu(xZgy+1S2dB6jdP^cpA!Ju6 zVU8WB{?Bfu%zEnqcqBFJvmkcc10*3Ji33R%K-Qxwe;M?-8aNjr|H%XdS%a8Rj>{k> zR2ag9HozPH8^M1Y!0~U6(f2^XADmu)B>)-3@&B)5``46<g2=xxFnk3e4Cue3F5;7d zYF#1VczSwcQ1vSa0_p)lKo3I@(A^*e)R9d^CAJ0HNIA8%8c9@gMgO_;PLP|z00l!K zNEZ(BUNoW0-_NDBZva8Ri~5h4*hX03KO7eLuT!*Y23WOqphhjYks{UK1nr^pGvStO zFpM|^hz9-bGbYd;N+XB{{SE}R{`&$!RVy|QJ^lR;J^d{B^Tr{9fA<gpU|x{d{yOvG z-`RMb{}~&v1ndS<<EJ7u{?VHGyTKVZ6w|McYyZI^T6jsC%s}+b;rStC*I>*#sH<;% z#}xo)NdL~McL>x4oJ8sT#kn-<e`h0QBi{Kd5a5QNu;C|sozQ<bs9HNK0r**0Fi{3{ zGA&BssT6o94LIcvo_DX}2jHN;1Ji&GrCkQ$px=?8oc{(Al&BE9x)H$GSdKSRjQncJ zhw!a`v$X>VV6?3x$DzbJ4Lm=4Kec~$!yoylBc^}1HRW#}<X2_>fgS$`M92&Wuv<pX z`Vv%4bD4bfU5os6f?A|i!HG}g6N`Uw@)T(A0^BXbK^as8z}>PQD1%Dgi^BiL<ZOs^ zG<gj|kcWVXQ0vND5E!(IXZ`3q5zNZ<+viOnFes3Jhx<oDF4S)&<Z=ao+%)e{z|r)( z8e|UOm^vHrjMFa@a{2y@gj{g%4FviB9}wi}xej)fwtp@q!bVEM|6tZg0Jue1M}+dB zrf<UMwANE`xvU>G!GkdJ-_Dl-rQ%WrVdOs$<X8V41o_BML6CpkK!k1}LN{`yem~3i zEm6kz><j;3qWyk4+p+7!|K-5v;D)~V=Pl3Q%XTtH-cymtr``UnvMwmQD<oQ=eAg~5 z`km~m8SRb9tc<3ams^~#-m%VDV!Qs-(Yln@7S@s57-Qxo1wYv@G4jhaTa|&@KPA1e zij}I+Pf4`OCzbJY@)*Y^8vLA$u8l1@KPRX4(8#HulR9>8KK`W@fmpa!K7=8E?y~k_ z#I4%6GU7&UT$yk<n;Ayjw9O0?j%35nh#RosXTr_dWH92!|1qg31o(iEq=OAtgGxW& zt&1uN<8)+{`oL}V<Y1*~fNv4aoSnpA!`+urCD%C}*`z**S}C^;?gb9G@mm*lBijEg zuw3c`o1_eRP96l&$%6@U<-r8G@?e5o`AC4|m?{62G7n}3Wx}(A+44Zhs3rs%y$(9T z#+R*&{u`wP7^<I<V%>M;&d<qYyrFD=BS!m&PA|nK77F9WP`Jq(i9xfAKuWL4!sdQ1 zY+G$VH+HJFe>=9Xwx0(}tX1X0*4L_XV~1*0w_`hNRe64qnB#knK=33AV$O#FIEk{T zOakUFfEOLPt8UOf<j!;-%pNNI4$K~!1wuzvA%&4^WnH8;#`GsMYX6}Vv9Eug65g0z zzs@oIdluI(e=1?amieauI^=WwAq4wRT&41^6B`PF96cla%0%EN3|1YVD$oRYAWLZi z0+2p50X|4PO+XNmg@)jUjHMy?A`NH=p2&8ZU*sYFp2P<@L;o*CK;)tTdks3ecJ(On z&n5eZa`J?MdLti7)e``v>Jfl0UxzLOsd~gV4AI~4e*X~OXKsE&2X5%V_1~p|OpFY~ zR5z~_`7@9U@M*VoT$KzuC8SCo-5XLRi>?T%Qa}?z;L_-b5V#z=F$BI3of87zj~)ww z%b=lSvEAQ(6$MD^^)G-qAV&d8Go0~LHd(6NH!R+N^R(+mc;VMjY$2LI6l45};f0^3 z)i2qRSGobX=sZh0kO(k<DF2Zsr0O90Oh}aqS}CMT2^|nprG};rfgeC04S_48MMB_; zXorwrWb^s{vnb&4wLgw8zuN%BYM{g$$s0z-KhNd=v*NfB$M~n;^815141S_k2jUos ze_<Tsk7v~A6i?q)QAz-Q`$UIt=_ma5Nuk9?#Q)zslJXP1<ewJt|7lC|f1M7jzKaf8 zDY5pcTnhnCvNc%+>ET}4xH98D*toLb-q_4A<N9r8Sa7vA{LHv_HvBBOW}6IV+((;? z@4wh+{XJWt(CEK#=rzNB4q#><m_2mhr|ivc#6ExITlp7H$!}~e{QFxA|9$(odDp9O zj;Lwh6<PoScC(6f7aLRCzXO|J+s}thtL^8-R@L_NV?%0HcVM$>Rr#=SwW_?>;@Tfc ziJ|OX(9f@FxBezt^ml%~b$xpk$RFCh5y<#eP;z5$;Xj%Y^-ow$g|y%DuK!;hBfu>& z`CoQR#5UZLe>}$a2Y*m~1BmqNEsSj;e{lZ?!%FDOAdN@xBm9O4*45-*6ln-gUapV8 zkk?|>NG%zNg?6>Ida0m=UJUu8V^9}a1j|D1YUycLBMqq$doqago&4m{(jbxpirj+w zl13n}nUI@;2;+&?h}HVlxS-Vm0*X8^X1#)0WKkn7wiFTvW=KQ~et{pc(yc}s<wvZo z7@z_>P^6>e)d(_fK4W>x6|uU~$WIzcoC#W)%L!UZ!I0*1GJ;lC2>d}SQ?8^X1GUPf zz&`Ssk7`~d6k>vufmt4ktvaz*eKeZ9Qs0hRTgLDoBP)@Cdh!#=c0q)0H6js%T9_wa zK@nLn<bhc7(o!OsSZckt9KlaoOo>G-mtq2Oo&440MDp@S7V>f-d5zePg0JKRkvOHu zjZr94ADIv_qc%ofJ4y;#8~a4ALizcTn_I|8GC^vo!i3Blgm+aVuBlm(B*;@_8W@aB zri2;lP*SnMXkeRRFxW1bSy*`LW?C4mjtd3@{=$5nq<mc69X&3bKkw-w1^pu{>Fey& zYG&rqyMGH}Mwd(((SHA>&U=~L`|j^iOn>FVEi%5i_2JIz^2eojUidh^i1ML$tGx5z ziBwPX7@@;aBZIxH48FLDU3$*~bfc2adGt-Bl!q5za~gC~=vJ#_i76^zXFShy=O$zO z(ou^qg&b;RIi)TBCRiDJ?qf}nWz`Pi8O$`A&l>lNuI_H!-s(iBw>2cgQ(=#L<U1|p zJHnK`a#e>T2Xzy-JzU(e{ai(i;EtP*bT7YpuqBf9B`l~Yc46>%PyML4A5SFny-C(| zj+cktyq_1xZ?@(>dWY>?%j(wSwA*@$65?tfWI1rBWnGB;Vwudd^8oFL%g}E7*dvyf zI$7f*l7eh~!k=%^wrJkM9U5p4FYVkEeZxeNe_!OK<BbW*I?*G{d#YpoF38`<CYrtD z{_qs5m9%ig)=DBeB}%`=SGJqHFE5b-!$Fp2>Kx?H@5Zd!2@5PsbnVoemDjGNRZ!mQ z1V<)0+g;NsDAJnqE|AI}pP%4t89Fss6f2|3mS6Ph?G%%~+e_&~NY9wDE{R<}sdL+o z30oUVyJhaNn^UnG*quD7daWzcz>ZfnT!?0gdup!JKeI^_R$Y7Mw6K%rZOjl;Z~M}9 zl19K$=EBRBy|OUbVFAJ0H4C$kvt<gmY(acEED+!BUtc|PH(g~oZ1C36-O=}-+U(AM za3+4oUW!j5yKC@=wUJY<aWK_5F5T)2Lgzm5vqTgq@WH7Y4~gz(WZxOiA*;>It#ROj z@$fa1eU6HovW}*{7=B+ieI>X?fcoYP=ULbNFK7u`?6FdQhcfOqVBMcQpng&FIq9QU z`KQTnC8uysy8|$!yq)PBm7PAj>Z5n2^j6KbhaVT3X`Q7m+H7xo_0ZGOu2U2Wt+uk2 zTk4r@+&<COq|$5A*1vaiqAohOJ?g-qvBkk9>6;$Tv}KVhJfq0%R@ZE4(*w5g<%rHI zH<(usUmWn@33M-TuFkw#Xv>^35XR*_JJ+Oo!%#nLvvZ2;g9A;8$`XP5O3L<jw@HI7 zLYt*=rVMo*WPj3WY@in#*|4?wEKX08Xe8gx&_}GpioN6Zd>G$clFC+MK;`&oo1}T@ z4aKau3Afz)Z;oAWh}6u@qMFFPd#4-?*X+6{n{R(B)oJ`b^q6GjQ6+DtoHt%Pe3!Hn z@0#9jbWXflbk$MF?}_m-oybna>y-EQ6Nf@7l__FQ)wbT`mX|wP;?d9<{E+RjIU^=V z!^`=U;sM_hC8Pd>W}5fyOMN+#9+<>$u~2Dyj_us#Wvg7=a=PM3uZduZ(*k{Sk@p1d zo&G7Al6TLJTgv75Y)NmsG#vfD>2Xr|g><@2lb!>d=j4yDhpk0EF4iBqSD=`~aDb&a zOT!Si8D%#i+4;cLC#^G#Yc+ew>;0pGq&FwWJR8_8GcJ{6-OSIk7h8#EjBmS-4bDwS zkiMYdm$6srjl%BodY*pg`4_BH5}vyJC$_l`2Q!#|-ZzwgpSpqPLjJ)RmC?5h!k243 znS|&LF4Mi9K9^A<644?q^^skY{u$+dzWawPPHA6=FfA!%W!gm1U;p@vzrk}Xj?>&5 z!6?CcPGV8b?OlcNs*_=yeTd||gWSi@d-PA-cXWH2eQhgR$;U1Id{Q+g3u{tNIz9h! z#@I`4#EYcNG~iihy=HhakMjJ`8KXG;LFVQCg|m^Ke1-?G{1jcuwQjZkC#rQ$QueF( zqH1K+_Pks!nty+McF&X4Q{|Px2WG0qMBG>C8n->~K5{Mgww*+crp8ce&Ij^=`JJ`l zH_swbsSca*p}Xz{q7D>QD|uD!b;ZA-R9oX-)@BNLxBCuLPId9623vz|YlSBC=<YY2 zU4&cb2P>+m`pz=-D!?0yjni!loO|7>Y%WAsl{h{L$(-yBE_*5xr|;ESczmWW#j8E2 z^19mD<?^}|uJh8%H`Ogk`5N1!cG1~Z2OgJsoVSM0ZSEb_YU^q;3lF|19+7_jg^q)> zc;XxM=E{K!Q-y18gV9Qt>SXAq@AZ_Oj>!6$QInUhs3aU@m^=)(&QFxLh;8k1zesh1 zpQ`|6v7NO&jVqJorH$NW-k^_O7{@zKC3y#;I`>Uoe5-xOlwpX<`(3(sSnaLZy*qM) zvq~mMf=c;kZy=4|GkmC3^t2G1<-OMVc0UuwI$R4W*H(k?vzZc&lvD_+<K^nVP<FY! z1bfjfV{LJ$#eA<kO##OJ!tH`BeGO)rLcR>Vi66%}El>BqZ|094zJIG+h;h3Ar5g)= zii?}QSnZi^_xwaR)x>bZp$dHK?j4Vo?<yq4mpzZUSleIMKDn=o-Qdv0ESXQJC`q@1 zvCw^yb{7_3Z;B9AkfRH*A?sdIqnqcF6w^hY{hYxpexrv?aikl$=q_m9Mm1x&Xxdlj zxRhFHluITcmw`3%*Q+Pni=L2RN*Ju3?sr#DIdJt%P3XM=te)XB48y#e`wN9u*?N&X zcFS;B-ebBde6jFi<n!csjj#tSDfegkx1anFc}qds4%^dVFrPvhAt-@KZMnK5yGYw) z`}vSN+sx?F@-8L5lggNxvA<*5n3u{~etSFoc0#pT)vB{!$LPE1Y)`7q@oiuF%lvQc zHpZTBo;uY1^y;nr{X#x&j^5YvV}{GDPjfU&>*XA8HqxZDwRN%mGG$Cx_F>TQ^=s;D zx2SzSqx$lCtS*eL30t`=)7e!mcV52Ct7Xo|gI^Pim+qO!=h<v2Tw8*;vv7{D!7TF3 zCW$gX<@Y|Yf?Y3SMDGO551c$1*-&555$;Nlz9q}!arg}lDSTFK;y^1Q!Nn&zD&c+m zupZZb{b*%T@qCB(v)xadge4r88F`2TPCa$fmiCto{j|t>Deb&#`_{B)PChyr<h7j; zTaZwFRgHU7qXy!$5VlrCd9&c%CkbU&%U@o*eycaS`h1l90wv-dMmyb!vc~ah&V{Tc z!G7f{EPDsf@1p)7J22~UXXg9~7+u|S%qGuk3AYo_Rn>22h-H)_9%UmDbyeLBYN9@b z=r`k<1f2{-;-yj_w^@;EY~@#$X7#kg`e52I8M5>V2DVd9(kZ&HIzLX*d54LJqxK37 z+N4SHC_5gapcr4;oWN|`+``7SMn`QNk_@|h!+t1%sz}bommIg{ib3mLrD}6-xPW2R zmh1Tkz4t!nQjesLbUo4b#Zh}F1;_iLOVzqJZtA^j_3h56?eKlAs2j{yvgK9JUOKdE zS`8CcDOmJ{{K7@`(*C78UWPl-A6_k&(>1lSvp%|deZ2U3!M2ragdY6XF0Qvm0x)bR zJ1ke`UY-UFeM1a(`s^W->#*gTiHIwg7{7>ag+=AE9FnK;jNe@Id6$3wC&L?C5u=)s z+okR}?J1qMdRo9HTVS%ok^adkmRoGHpQ)^)s7Lti#<hdDnx8&!U9a<f!MHQ2KSQ@_ zg`T_P%6I^V-(30i4nn#~rfgTB2zl5(M^uNtEX?C3Gm<X(ZFOBtXGw(9vk=6v6n?gD z*MMPm2%S!&Qx9U)@X8sc{z_>ydTn=Krtq2bPY9bO0}tHcL0#SE%vO128Gj@CL(hTw z_^e@@5n^odt^yTZa)B+PsK7|;=-&66D!mR0l}DO=_GnTbuD;IU-q>(zf%`~L(c87& z?!qoIVAB42+9qM$<4!WbU{-p+J8k8mX`9;XGStMbiJ=_kr>ZJ#$nEGtswDAc-J-Z% zyECpu_RM|fOrcYmU%Y*;t_4JnHI-`w(7lqg5BzA?r8+p8*)}5JidcSbwN!1~=53~i z!POhRn=xu%@GlBPR8EnJt{JMjgevDVs#33wh!%Y~lu0j%#V2U3zRj^O@T01cQE!W^ zva2h>z`3G~kMa3J&LvBi!U?^lRBOLSDPA8(T~iusy*p2j)Y0UdgqX4o)s|bR;^E<w zqbeTnSy)>9ysX=;KXiFa#>sW%I(ZO#b4Kh5`+G|pDn`%9jW{9-U6*#m`_2?%#~A91 zS3k7)opPOwz+JK`q_yf<sYqqA7n$);32NuDOl(U`Y>QBy6GoM_w{$%vy115N&a@W> zAp#3h<`+ao<|P<g{ASL?k)u|~o~PaXrWZ-KTUHPc7gSaXF4nodJ#C%Xm+6%xN~*5n z(jt@%+CE)W65*K<ne;aCAU?JDgd?BAlHg4l>iEoM!w1v#?@(mij0<)%sbqX)|6|17 z?t>wiaxE?#t&c*~c~!{@pHGgt5xu0?fmPL*hkauEC-_??yaE=z3LUfBy@=TX17j8% z@7k>MbNze{;a{(fJ1(6nF&}w6Pv|GLIK=dYh7hA0tx;{HW^O5yvPXQc6O9O`s#eF` z->LTD9Eb{2g-eqexg$fTQj~(;>m4jK@e3rO&q!5t9#;ME+7-9b)fn!{9^P=^MtN}^ zetMP%AFz<<e5TN+IHFHkRlVS&eNzU;p36$_h>T7gUjG*Byg=j|FNK!HON!h~jh$I$ z^HH4ibegU4s~SUcyweV^1J2_0omdzHInf1GQoZxj4<B0xlu^uAqg`H4Yj74s9DLv> zLP9?-C9X}+OMbFnn(UMQ#8^FgXtY94HBxVITo^x=^57#P?y#hdwWbeas?l80DG&B$ zZ@5V}QNO5u*f+30*)A}rI!3LirQ*0#|6}(B$&Q5yMI)DV%B?9rxt{_B6lPAB`Qbw* zE0gr&Q+yh0S4PuB6<Ahiv(f$MR1WbS;;)wsXqi5Tp-n4k?l~O6Bi(mzNZw58%aw() z5Zvt>2W(yqY;mu=;p&Wgch0h7*JvC?lE(om_&!+Njfe-av`OMRkr59JXly#>ZZ0-o z8G5d&!m#7E>h>$g8`V2RD3NhfR5Md<z|+%j3Mvo>fa^!oo;CZIPCO74=Jo0j)`@r! zx8IEA?Qs!Z#qkGvAD9)#LywC9kMs<@?QE<S`f=4^8hlde%L7*2qM?ZMQ%|l5;CV{V z4t4B!rhd=4se5sJb_Og!nU1t+i(3NxAFV$gf$mfX%ImQGT2U-A=cwYvxRMZk|4<c` zj^im$cs)A=vE#SCmSJl!=D2fdp%yC8)*u#YWD#&v<&cQVt8~7}^;XK@7EOUxc)vX! zTa?rIwZLc<l@G^L#Ce5*|1ZJP7O4XK@2|H)1l`FB_Uf>RGiQze6@CA@ZWG0EI^3cL zP>AGl5tZ8Xbf7;8%%?#10QVzRRJzvN#(VV&gGt&VTYEovuY7mf;+?{5uZ~mL@y9@` z%sO;|68P@LrP~^`d#-n$F~I-fdbJKoi{NdbmF?2w$?H8#x7;bBIAQ2nzPLGzxz~&9 z{?MBpH!VFsXWCpn`K2jj+%Sn*g-!PaJ#T<CuLbM@TT<X2SauS#B<vKuK6!I$+H9iE zFk=U*=FN@IA0IR1xTa!v_=Pg8VK?7PEB{=k!X<V}GPU~MKKanS+_ujcc1RpreK@&D zmz`V?0n=sB?cDTPuDh1nA9IL@F}T0oFil?@)^&CK^QJ1-7fj)YW_a8}cD*IEez<x< z?sBw@CBv0&i|#NA8wOqU)grkQT(Zw$H`~1SFm^b+;eb6WW0t(BTKpV_Z1ZvpUQ8~f zhPf1kFr=;|myMu!w5dLqdsARJk}9YLod&Pi&U3)3%9z7y+u;vgDP~h9fhyanBw<%Q zcr8A+wzEmX!aXc2VAmLQ2d~N~Y!xbg0h4|M)a0MxX93$erugE}jBFeOY{Y_=^9y{` zouVccD#1w-mI}U?&d@q`_11yQHPV?3SH4(8)ZN1;cTm@O0_|B;jVp&edIL1E-@%*{ z)?HR$2}4s0d)J~<uiZ|zgk2a@WTUBS|0D@}0kj2pua7#*b7E>?X6@DZ`z_Ujlvg;D zK2x9ddL34{3>&0==u07;KIt#QQ^Xo3)68&X9%|7;zuFgRqS3$;U#NxsYtd=$1B)^J zk9Qqg+GL~By~Ia}<V?a+pY^>jZ684KIepTQYf*Kf3U;L(YSt@G*gU|@S&xi^*08rg zXKUKw_He(LLMx8CCsLaxJuEpmV#d?>r3>n}q#@3I;0#^@9#<qbK{t6VJXdd-SyO!g zUi`8sSX{mNNxRpSovwjK^Vj=3fzq~`W&j);U+^yjZ-e8tq-g^O`17FsiBv~j8H3o^ zh4Iw!a%l_Y#=NoU@#|-!&8#0#d$z6$NMxS<z<F^iZ*ldpmM@Y|vM4`rN<DdymF-!X zAPF|^DBTmOeUni&s}L)`RaFI)3*3cR`n9S>)vt})^y7tN8QmW>&-xE~)~-KJ-`P$0 z-kZs<T?9sWG74x3QnkBQlN>KH|GMJBkN(kz^gDdz2_U=cPf`14Kr8`d@B`<!t}dz; zHb%W#d`VINwXR>kSi*p;{>IjOc)u_!xrZY6zj8e6@Qm#DZd&E7kHT9*r4=+pU=u6t zzjfKK4}6oEYn>b5aIbeAkis$cKyWa6s7<ZE4JA;>Z7J?f))VV&gv9rscm|5^hs2pX z3Bc%$-@XC}fV$`E_3kC1v;Je_70i-uJnV&&I<FS*Fu0SbfM@pT=42!-34uc#rLQJ4 zsRIm8w`9%gw&H5=w%a=g(#t1O6yN}P&WPB<>DJ{#7J87D^BBOsO&`WU_9cET29gx$ z1yG*;`FV03yx4OH8>ZwCyp+v=f78)coH(5xwkY*UO-JFa)r*i1N$Cd`tY1VVxZJOD z<II;)8MzuXpyw$#6{d7CW_&N_8UX=TsQT*4lF+Y;{{pWr+`qk(RVf*2mLl83pfQ8i z)op-4fFt&;C!b4IzrBiDnoE1$oRog;_SZLP91!;KKJ|$AOr-@G<;x+VYZnzN9C*5> zrs*_rRcgdn>Dr59{-kVkPlMeP>7kt}*7m1Vjf3{4G>$C%e|)`TbSBZ#1{!B#+qN}v zCf3BZZQGt`V%xTDzp?SgwteS2=l(qR{^;(#YW3>9R`;%|r=DF^-O^e?co*ytRT=BX z{$)V&T~f`(Iey$R9S7UMlzZ3MSMgDmv1gTJp9-`8*?X07l0Spd_lFm-Q}JZF;H|;A z@wB2dz{GS8;4Ijw%ry1mOyftrXz4A_Y@tcD?R{>j7Wn*Lew|4u?RxD!@A^FRYwWt* zIF$Ay*6Pe6qMLSojwm}i4(&cxebkugz0N^i-snAwGhG+NxxDiKpE%s|_PQnP9}kT} z@&Ad#Isb#hwYHo$#t{7fk<ZYGL~o0d9|;tlDfR^?>_&T5tnJqCvdT%*;;5<#4cX@p z{!DII<e)~w6~|KKcptmY^&lw+pYqq?-S-*Es(SdxHG+3KRJ<j_oPq*?t=`X$H)iMT z3c1(~6Mg)~_II+%&0NB-m+<%2f)BlJ6V1zj5MvH^{j)b0E<*KF88^uaMaPL)2UGQ= z@S-rm$zw-X+D4Pgbi{W@peNV9NCJ^A*=|jjxOn}>p=AsG{5=WpvbbdWWg_lHt=^sT zlBD;J9)-O&r9IHv7>UKvvZm?oLvq!9Y@wknN|#^y&){krWp@|B4yz#2DMu_}!rlAp z11LmpUwf%f?uuYjXN|Mu_1>|uBBi{y#$R_qX{F=t5P#4Qq{#lrMVIbzsdeduzg^PN zL@|>W9O~suk+6RH=$}$_bSEL9y~?Y<>10DIyLgW`+TZPw`~~`2YX!QH^oj(=0REfd z-r4>7GSDqjEmDYeP`za_%#Nc-&l~u6N}d-<Yj#(D57(S%`nt2S{(#2TsE?~M-3Jx) zouyE0o13H;o7`gB+-!&*y+-*pKS6gapxYx81G*n%;JzB^XZPQp`FF(VO|#B}u3ABi zcmIQpH`A-_4unhX+m@NT<5XPE9!>#<r#__I84~p+5G4lBl^y#qu8e3(Wu}MV$n&D( z3Pwa~mCjYxpz;iqp1*T(;>tUn2&H`*OM)uf+wf5Vu5+ZoP83JpJF;u>@*Zj75LA#j zRI}_`btw?IW`vuqc9bc~uLCc-5^C@X>WP3C{bB+K$~QeB_?u9q>51J8&tHY}&p~qM zvP1qJFoEqHK+3FN%%QJZx2p%yGsC7zR7lUJw?{(#0SD!Bkyh|u*+>_(lR*xPrgE8H z@pt8WFJ&4a2RtCv0dH@}+naIvB;Cp(^yEWG+y{w}sE%XvVCfEH08~^@d{pfsm!nt` zcw73#Wr}(s^OPv~-SmpBDzcl#d<a7sO*q+TY}Dei`?&JZZe!)jF`eps2*px01EZmn z2@&7m;$|8}lbRs~gHJsB6)-e1sQH%Pp#JM&+1#ZQ0Td-sQ;)c?u5$nGR%iy*zXV`P zGTrxR5O+^`Hwk&Ax3VC=*oYucU62Gc7{CRR-XmkfA_Vpgg$XacwRR&*O3rSe;Vp_{ zF`wa^!2K|tFzY6DK1T@EVq&h(ugiPHngu)6j>|}U#A<26{ljotI=1O05{SB8(%m-{ zk?$4E*sigt@?3YKYErYJP^K3ciNc1KxgKsD@H?i&nP0MCl1$Mt*oaY2625nZC&Q`y zl>bm8Sn4pL3|p$H{^j$c&hYJhBka`P0k5QY5M}UbJ(S@}-+~}x`wt{16fWEvf4rkZ zt~T>$%vbSn8RzA_@a<0*N8`k87=!UE{S$+X?R&i^=E$hO>>5-yX?poD^F~$Mc*>Ms zL8qZ+!C6+U*5}?X`NQS);~j9O&Pxhu*(ZCGPf^rIeoW(n8_WZ~Wb0?U#&kl%B^zcw zhD<K@Kf@l*u@lnF8v@CqJeM~za&579H5YuV;bfTf?k<FtpG?)oQSYM*Jcc^laLPP8 z7@`{V2ZYLxEi&L5j5(>*{xY%i_MW@xwA;AjuFHt7-QX_H_L&YwL`7#OPM~0{1^qh= zYAVKbp^xh08CFNXq_KoK4Mh>uO$JBiwTtY6>x-^N3dq|OT!L$!*`PVP-E|ZdG;CNG zWK^Py`U|WF9|&w%Hh=lepP4~heULH-gKZxvlunz96pS*l$W@W=Yq5?<;?h@fMmx-5 z8^E52!lJp5fx4q)`nF=rAY3w}r=hALMAZ}=u3Ab4V_;Gslvw`3u|RSFvBM>DM~N$5 z6!Xkj+N^@8ju%bw-->^@yT?CUXIzOzNVZDN#-8x3eXEHdYMMOSWNrxZF=U<P>Yri> z$c)JTn?bEzadZVK$y|dXW+y*{u&=PzV)$nI@7t!@;+5#P;7xH1CUr!H4Wc9AwCLcn zu2pilGnoWA#!zE6(9JB3aW>f!YELQLHUfUYxk@NxBbNVhvJpbZkx)x;+=}n0z#+42 z>-oN3cs8apFJoP$`Wu>2`p!793!ykA*577lqDdRv{misjJ>TQTjYu42i4}b=p0Y@| zj(0D-@OsBdKEDw!+AK4|VSK=acA4#pV`N@jaZsL>-=^cus~vWCMV)Zc)~}$$cuCUe zriaMjE_iz)`9>t-Sl`5rT;G5aau-IGAnOY>QvI}r<w2qQ<8C2I{{ww^c)nmIN;jn8 zh=aT%daq?ynzaUpOkHjgMUk;mk-8EkL8Alq|0bQV%qN~O&HiDUErhDfcn<M@lI(wy z^dF|_|0kJ*FgBSW)1u@#%$TYa6(O~bC9lx=@3Sf)6d6aaC6^J%Sq_67r?K0ew4@`P zc(9$JtYb7pU8WA;m$!N&PE5vds@}a&<{ba9t-9NYy3?E2T{lG7Va^!r;Gq0=m=S8^ zwxB&FEsz?=rep`_h5BvB8qYRcMFQ_~ixBL2VoS)rXZ%tyQhNAiAAij=ydW4@!o#Yv z#Li1;st446TYEG9>j!8#Ywsu-$bMD6%lhP?$-Or|yl5oY%(xd{SYz*?5Jb+K=VNHq ztNepTjqjOgu>J)ysjN6??6i@L<$*&PoLHv1UIWt0s3Yp_KYn&UzIJG;>JNF(bv&4O z1oAzJuNOOaJv6nYthj8Nuu|h()1XR(fZ_%EVh;&?o>Ufa7ePuv+Wt{-pQJRw`xq;V zWB*v_XZ$XW_v)Z;(`tP=HJF<)c*>OfQ2I{S1*^!Nc`^N}hM%yX+iMCq_EqY5s~|e5 zjeFjY)g1rz_d6W_u(R-1zqxLtjHT5vm!Iks$)|kOvYVw=sfme={{!oPafoM>fZRz) z5D-or5D>KgZk)SVnA(~${pZZ`f3W33Q$7xl1HG5*mO$LORLbqXRFgyd@<<-5V?%AL zsZ7FrB;9czC(g#awLSkSA|!-50uGf6@kiOBw9Z7SCUp?$Bi+qK`qV03xb37GZP_@* z{%iN~@pFzp0MPF8swdn{HU(*D39#2Gdnc3w$4p>9du<&lL}Gl)iKj+DF{QbBg`@wN z+_>&_hp`>UrdkDmRk6v>A(BbAU9ut#FeGQJgD<_YoA3xC=esMzbaH|BLettPTiIxa z_`CS5loYUjzCZ(0*W0baTa?F}EQ|n=L_}Ncfq^%J3rm7oAm8pDVgU%p5q9Gnc`mwE zBq#4Wd3F7~wcbdzvZzEZaEhrS;hN7Cu0j4`j^`yu_<KSsb6TyjTls7SgjU^gaHfvm zCiRUS1QN80y>)nTZT<2JQ44Yz<WP{H13J6kY;XYYjZ58T+p@eibR|#6-x*R6_p$#E zre3g<&*X06jCil*E(i5!tNKm$+f3eC{S7X<a>9(-?4z1@W3P@aOgd}FafrM2w4!;T zf?hr`-o%BHvaBZU%Ar=!6{`SP0GJZF)Y-ULSSfYD$0eCGWtPV~c~o%*l&zsj7a$z0 zp0&vAL9|bb()&vSCx}4N()acu97~?an<+<tr~Uij;0r*7M3}YtO|xPObCu7|3YqKw z{Hnj>v-kP>vT=;m>#`pj`F%Mm@BeuS$A$c)Xgb$eOZfS9@b4n{zW*MM_<|^bh9*D> zgSpM&bFpOtQ}vs7qy=C8rBW)|Sd<<HJ=iY%EI*n7*bG^~xXIX(?Vbr9+7_UG<j$at zWGs+6#B#s%vpxVidQNvF8?mk8h?S2or}VG4lK__^7YdKl^p-?9#`-QF_!T*0@+Y>o z*>)avGW7sWN{-j488)w)dpPoBEi^e$wyyl)3fUeDrBr6b#E3$PqDcZZVsgNgIAr~9 zy5H75RdJ;P`ejUfAWUe{@B>?UGDdw?r(meqDH_#^$A@Tq)%0LD2%6H4np~7;AjxwO z8E!YBea5g`e6|+Mw{6f*2rDw?AYn>}=NZMA(y7nv6PqG^!7O2f@=1gx7+q`|fl~Jd z-0T|p5#CYgblqWxauxP@P1FG{ZGpX98u`&|hM>MHD7~8q%O$GL42Rh5Pc2$Pj@0(G zg}3Ik+fPTkoNi^4Zp3Gzq7%@JQZD9;s3p|L*!&%qrte<!6Ibuo56|`jBZJekp=!X+ z`1<~Hl7bgC7Wy@uFBq0WzyJ#T*|N}cxG8VS@J-{-3+=WDGKPIeN<m6ifoMgj(&sQc zJ~G9!*pO5fP~AepWdsI+V(n*;`-6Z`&y_WcGDJ2iq(Bo6J^SSSS(Y*Utt*vmm2$o3 zlD}6TsZi8w>MR3kQ=CQ%hN9-&KNMf6$Q~wUhNEf*qy9*LvZ_o{r__KfOqBJ(j0Zcb z0_-{N?FiQ?MIFoblDM8GPbz>})<rAai5WNw;W9C#;1!bY$<@}WrE00A|M0{*>wqNJ z#t$+CWAr^z-s$tW`)Joz*mBbxr(q{96*Cvco=6QkRPJE}c8>I*l5IzTy?qPNj!sAu zkSqL*D3}QEBvF>gP=KQ22bpk;_vM;Ku`i4mAK@HsJwLuiOoHkhmr!(Ed6>E;q&!Ws z-ZBD(JfLl@gJp4O%Z3<a2aO|N`ZL~aH3r3(pb)h0hFksTj@9i7G@TS@xC7D1Zar7N zvme^8sAe97!@l?*ZF`mn+wd@((kY^9hdbp1Dzi173}w#LL=x<O6crU^SpmBlqB718 z$~_jkhkb!5a2Lk~$a18r2H&nr{GBRv0GuMtCD}i!DoZYjr<rtjVXCSq%~ccS+N}r$ zdmUE^=fsNHs(aBD@2HUd%+(oFtiLt8qz9uK|0r%WIrAZ7=UCL<YHN~Z<#=pOAUTyb zs$oPNwkz3#e~yN}B7tXaLf&Eb>0V_?|J~W%4%Wv1`Ec0UPho2@ci9X2zifYlb7eQn zzYSqJCJ+#W|Jr_M7f%~g=l@uJ+>O{h&WCQLZ;}ab&Sx7wFf|dzdFhQMX}07y+@L@U z=X@!RUzMLh0#AZ((4U!4f?gNMsC|g>pk{=)<Po+r!^cO)@9+Xui5f3q_oNf^`DH!I zqM)UQoH6sS3Lmd)sWyt#6U5n2jwyTjZT+Dyh7bEl!3CQ`FP4FhG$2t;*inBPbOgqn zP#7Ih(XSgQ<P1$x4y6=?*>Wd`Sx5Vuh=ru+t)q;_rI3q_LU&1E0KXC-URRcsU>i+v zlIF;x3Y|n8BZW}~e6mV431y88T&r<^QCF$bt%SJhq2+f=xFk?dHIil~UxU<UO`9o3 z#I1?3CIr{@L!`YMD*e`+dphcw@it;T!1w<X@isJPkWy?)60%0i5MEGX?=2p@$F-}C z9PFp6NQf{%KjTD_D<4*HP4_^uWf>1VvzJ4bNaQQ{)`jo{<*Dg&4SpC8krj)M>41O( zfmAS(t!gBW-1{A>LPlLaXH!mvQg-9dLvQ;KzY!wY+WHpv<LdMa!ft0}cQ-?g>k2W4 zuJAJQnH(@gq&LQ9)kQB|>Jd{aGAD2o=KuZn4-9<sO%8j|?@Uh5{6TkDJ}$1`TSqy; zzZ8J}Lk|<%?2W#>-qDS-M|+V=1#RZg8{K#xJi>9S$1)b`{b&2{Ou<$B+cs71y~voQ z_ck1eOLvzEjNk%GZCI+70_f{6^wHm}!cHhB*a*bdbdEW1v&b#ifJ}-aI2!V_<f_we zxexCgXgpOmZ4g$9pSZFOKCm#h*-+&|?Ps_MMtpR?JU!X<<9ArK2cF64WF+6jgDugg z8rdVHqLf8!^bL4OV)v;M9==#NcEtm{`O)T8_oxmfBJSb^5)aD!LL+VT>3md;*#kbT zawva>cEKK-q9<&^r{yPV0$zS12mwZ5`0lai-a#yV0~WauL%gJdPLdd0NHiGtr;k-x z``hK8xfj7KnCh=fI7&Fznux5yEjoG*G$RKSt_Cf@NLv|d8hTSsafWJ`TKNt6y~|PD z#tOytU`C+ji8})-mU4)8yK9lsd|Wq!yB<3Zc6u3K8(sYwY<GfyPREwQS)uw;ZhG80 z>L-_$ZVVE=Lny);cn-}lS4I`3bUkV1dZ;0Z<a+qs;rmZXsu-izTaQVbZnoDHEi3c& zKH!6$QOx$|c<N@iIkdec<=evVuBZ>f=2j2>^1FJGui6qp`s(pwe&e_FTF=tmUXkcC zSat;K@3*d>X*=}*ufm}dfsj+DaA!p+J8!4wj*VY+Zh{=ypNqX@%(@~u=JXTYQ{lai z+Y=z1<1V%CX%|myYF{bu0Xok;6-2j<$6j!I{DY%g%bmPKtuv1)Je0>j_7Lpbce?!a ze4UrS(GM<$&Fo_d=ANAJ&H-I+^Sg$%8rcPx1SFT!hmZ&I;_K&$r6eYuEml;?>I*x9 zS|*HppTiEfejg0EYt3riV(i&QN>fgEKfCR9dUt6+_#Z-zx0ihfw`~WnWHpFaJVoZT znf8ehCg^v506cVs9bx%_kcV<=iaO#n+*cf@3?DQapdST6Vag+}Uv*I+5AM@W{Q7h8 zML1CmI8pmjEk^v(CEq6<Io5(5J<r+OF3VWt+RFloXt!<RFq}I^5jm{{cgmpbIsx6+ zZDkWG=wd-wi9U0!XEGa@^)3!176@8xqHFohs`85Qg;JPzBMF9a4C8LH<oOcA%qm;X ze1J<~I7?*;7SawaBi*jEjefZzgj$EQHjzpzOx*Jty3;6Ysj#E=gzCx)lZ=YPpDp3> z6>j+y_a!aA(ggmfHygE}Sq{=&5V~3L*WvowWF%a2OR~Z*B*fcWB?$O6E#0op$-2wV zai-!M+%32|H;XeAktgds%-bGA7)U7d!KpdEcH!YiO~1=~;8rIbv7*h|Su^7+Jwg1Y zHqWb7g;c4gH|h^36YyOnb4?#(;Z({Q2qHli=*BXeGaN+LvB91P5Z1e6=#F;%Wg)cZ zbRJUkU@po<5P4f8ixn22NBjYS2h~5ALnGW_r^z?Xjc7gO=;72zPc3<1uzxE~?BVjG zc+zBgg>yH(w4IR`c$#G2{gHmh@Rv)ur)vC%ww(!Q#BwOU7Lza^4@6KYWX>h>n+?xU z$!f-eD>Pw>9E@)c^t}Y~UjJgKb2w`qa946v_{MMp4<O@0K+-Zy7llZN0f)C6&A;bC zOJ8bpdc#f3cXOgL96YLcA#hJJ7s{{fw@EV3%|Y~Pz@K2rf<!MH-{VW0SBL2rs!I6l z>M&-9Y^;ZUtZgLT`*bxap4*LL!C|E4oJ>qP>$-a42y2mJLD;5@S|KK!H0<oyr}M%N zYooNpp1K7$SCZy;vL*P2?e<?Ynz)h*$)1|B1VhLc6YUO~*OH3%6qYC19h+Obrz(Fn z6V+<ka+v!g5gdx2vGY9EWzL~)Ohs_B7il3n*1~obn-bsu)3+YRX@hP2OQ-t$7gwSG zzrNMN(8<(9*~Q7y&ip?vwo2319*+z4%dhV<JkZG^Gq@}2vS_+E7H-t72I_9!ISUH( zUl^q;J?W9UsqI8T*GA+IfOn@5Q%S&VL;7_=!JUOO!*@URa**2DQ4ng^?f7(h!tBjV zYU`|y&FjzO<$eDFyp8ZtXi-RRd^@!${ya%t=>0kzdKda>{kfosBYbbddh1zkeaN<A z@h3Um+8X1=nm6~QWU1n&csaj@QUZza<7attJ=uh=IrtU70KTI6v_Fq8VK3j;=*h!S z=0c*wTW1B2FDpzb71|fFpk7%Xwcam5bgz;QTE6^=7YBd$%B<4#@*um`Cjm)LJxof4 zLM>J6?GB(fKZA_z_?eD(?YgQsQ<IhjY+Z6mazm?$S(6h1Loxz7zN?;Exq1aw;-%H& z{4dy<W))t6AboA6r1J_RnhPe@evxQ^_5wtgJn8-(N4J-+h-di)pWfLzn9*kRWAJP< zivafqeuQt5d;hfAhTk|LpR~^Xrs}YRywe`_+LzwDt~VeXp}NvdaHnrp!5S|iblxK& zg5)D06_aX^t!J)1y&9Q^G#)t<9g2qCUUy~av%S)=_jbmM4;Cx~oSRr!wgzZ?xvRT# zNI@gNE958N_0YlLY^>_)T3@*%_*(FA&$ds?b&OZzs=WK1_iUp?!n6S3+F^Y?1Iysa zutIDQ&T0q<jUck>XmKRm<<STsw#9JC9h~oKGIHiEI~ZMbTfVa@l&yp-B2x*Hg)~;* zHUR(N30$dPJ?wIv;FScJ@kpp$ckJKEE&VP&y&1onz|*{Op=NRwl?nb44x{)>m%Lyv zaIhOuz6WtWPV8+^koMD}H)Zr)Q1Bz6St7DIm?*Zl?T%T~^U%w@Z#t5z20~Yh7UBk8 zC)$NM1r5J_qeYT{_9n#B;Zi!Q4ne$Ul-5btf(NKvuH8mh29<@g&m|LzrpczRY6pXW z`e^5nn8B(;QfMXDiV9XGauPCc{t#PyQykj311_Z9$_yF1ntWk!imAW}QV?>ZcZs&K z<Hdhov8$D95muDL{l=7_$7KTha<mP1z(oS%|9P1u&_R=5@FOBKfKQ)Vxtyj&6I$ap ziKFIhH&?OtmAhR-@7<@{BBlwUH<TU(i2rLKDvsx`?WH&H_k^o#`*=QIZ4*Cg4D~Vv zKN3qNS<l@g36}-hbK`nZoV5Wjpi9tJtvOpQ_)2_jF2AgtA6l<M{UhqP$_a8PrYk4u zW0f^GGs-etUPC%~aXw|A#Lgi`p11s|x3K_47I31#MkIi7-X}DR^z-X?3>T!zY|5Y- zs$;<T0YMS{>=NYRgccU|38~Bpy%EGOVF56ThM&Uit3S3ZmOHda?HRI2JAUCKI&BxI zftTS-(lG|a;KXKkR#Y<Z7ZUY}Yc)q}?IB`tqwwD22Ku&{ivW|_b`rxp6wd67@Eoi< ze+<pnrKX8?$Kg~!-J@Fmqp!wsTdc{=APWt|=+`{iYOR(DzbQz!*?be816tRXa4i|x zVNF~aQ|=WJd12avx1Wt1#RJs*W~CvAradyn#p*D?I>mmd<^=15lG4rgaq+O$>_GLb zHYSUX4@T=%X!8*YuF`g-6X_S~1rhwPfzMKOBc?xEFYq!}{83OJa1~vXj>-xFe#bv) zSVOH)h(YUoAT65|{g{90<d=?q5X%Zy6}7u!MZ?aMGt}&)*p<j1pLX}d%=s5ebxKWl zHjg_Vf(0#ANOCaGfW5fKKmNyXXp_u?9R?$zr$U93qsHcLipiQ8L?Jp-;h5SqNV=hk zR3Y9~o*3#ggpnN=!4DBkBSQLw56EHYMb{_9UlK_jF!u7?vj>(3m6M~dA9H#+o$22B z82X+`)uYJ$M##>r*9yna(;e%hWANe-qz2kt@6<g08BXFT_<*45#?To{S+xaydj;J> z-C60wAT#PnqjvnbPH+j9#$*{Z@P=aQ2j)xF13`$Mn5x|>-*>W%Yx1+tbq}Eb`Q7GO z+kWS`GFv1=?dnI@zlN5gp#c2_a9J!|*tcUeo8eIsVkR~6;hhIzRI?VcLq$9G-+Mc_ zaSp?c&mj>T1@7H*W3eCXrUm{mF?m7I(Ii1=k5u%`lP}hET`N`ME@;F%3kv@{O!%W| z%`JEmM`7G>?vZ!}uBuo1J}V^!7!so_TnWk*>2{oVE6~gc;(0>uyb9?*cDt9~Apu3$ zN1F*SOmE)`WuTPed4lDa5r(>A;zA17!&{2YqNc=~n_gSJ!+epdI5SUby+-0H#+&U~ zf^C38^NZ=BE<%R!fncSVfi8e+GX5C~TPD4NSEQ+hvOX@GFYnGXLueX~4Z+;OhH%Oa zK1<lnJF2rQAnQt?WqELRT;;9&37w#e1DEE8*z(NK$Jj;T^CH1VybV!&l$|qwJ>p)$ znoF7~Y(FMSyUlhWxv!Ph0v*je=^QfWdDsC~T%(IUB`7ZC@)-hWJJTG4!;$J8s(qw} zHE$3VIXHY1RS;?kp0o#SLqII*mbjvL$ZoNN0Z^s#>1jinuWZKHXEl@A%Nj1-{JH$8 zo`Fd`<gWVidbc!v_~?Frb6+Af+yF;gXBV_wjhMTcQ7bu#mM%c1P*k^xa5<W&E(fcM zzu@pUH7OVrd5!d%dnPQ*m$jrf?$o|WB&V8TRcAi<kI5<AbTnj|QKuKJU|q(Y86kS8 zH%gsYtr@<Mjvt#>Pch>YB2Er?2x3>uHbmnr4dER#BSAI+)nae1C=MN8=9laC!e)%8 zr0yb_K|v)ZNQ9&OpJ6SWdPhX<5}hAwIYt*0B*Z~dsXc5fBr1pV7}_z2h{Zh461G+j zm1!iIEJo!OI3)A3Tg9WKgS8wH2)Je@<LGkhJ?E9B{oKajfzW4zE<^n1WkhDIYB0&+ zokfX!yltCj6hgS*W9IVDfU>LED?31N*byvJ7v`g^L)mdp1F<5(b$gOn@N#6Z0Y^05 z=^f(=!Zx8J7RQ|nXNfs0R>9;gGeIZ<N=EqaK@Vu^07+Sy3jyd&2j~P3;Q6u?Dj3V` zd|fQwc>?j*Ud##4tc&3F#@wIjyhr(phW(RWm=?-b$De)EP|u~-;dV5H;nR^<SSMk# zP|UxPmzhY|;KjKSXBBUZIbk0u6M(MIqlG!5#G~tWZyyu8U#wy8SL^0tChoTz{qR!+ zU|^hia!#rclr#y1(c$Nk5U2Vt48?{7A5<!VyqN-ICM-S9;^u5MU994$?K&GH=W}30 zTcSghG6ziX-s_POJe|xApL?vJ2tFH+1KvvRpAlZWw6o;0Beba=uT?J*CZVsN2mH4W zhEK{VQ(OR+?||a@g@<63SX`(DEe;sp?F+)y7!!#G-yp%AtS_YhMYPCX@|$q}BiV1L zKtS;R3kFS`4Bh_^2Qs<Y+KB#l%)<C@e4(WqSICL(Tf6vez=Cw2AY1n;okp(2yVTR1 zye{oqrbY*q73mHl0W@KKKliRhz>#d6V-In2OOk{D8(GbAt_9SyZJE&oGB6EJPibL! zYK+L_h#fMcIg-e<D-3Fw?NZ8AU-wo@PZ(1i8j%<gW^7IJZ>Vn$wz9}qbMu+g4m#mf zd*Jd<Hj`$?+b>P+MJl|p>S0M?R%@wCqLy2jDOwj*VmF_q$#a!RKX&bjlCzyJEQbN8 zhHcFLYD|#Q3}Xp2Vn5n?=k+pxxmE$kHJ~!Y)c$JCRN2I5oQKk6W^ARX#{^TA3kx{i zX-0La0B9C7V20rqSiyexM~#V~Py`)IyQ>P!(1#*Lo95hW81kB{fl9dg;}WbsmaHWv zD0v&9IPa@gdY7QA&}qU_m=KCT7t5|h?HH9dlzAY(G9SWWbbK9*XqRuCWhGR`&t-zh z^GAv@h#bZ2dW;VbDuwn`!ujJw@$GKJpPg|Qpg)shoN?<!k?+3Rs!i-Fqg01=;N%Z@ z*NzJpS<)+OPE}Px`uuE;@U(x_6bB!DbZ1fpg}FM!?58lvh(JIeb@z_J3R9#dyhG-7 z<j;)-RxolA*bW@|F4%gc?(u%Zv<drlf5TrtBYl?m|E(B^6UryL^PqvdryhU`@v1`@ zGZ6{`bu<j`?LIBjMYY$Vdj<94{6;ROSb;oaE7{t(;d*5^_kct?OLPyANCD8~)W+T> ze6?HG?Ya&J*E$BQA0G$oxBdjiaa?i0ts;m78F*XM)@>lNGmkX~z6*jtZ-SoIi*JY9 zBWNY!)Q)bPqS!}dI*)d(#J={~4f%w+#WXE&ob14^%f)cw3&p$t+PFlq#jcu(Y3aq< zs)2bT+r8`wzi~VFx%w4hoa7#1&yW6)1jxQ86R;a@6)7JFN+3UPgYPk0WQN2>3@Ty` zeG(dP@(b=IH3#$MO_OYjDXN>he0Q6;fV27(YJfkAKT`Xwg9x0W*!7ZaA{ta;XoxgJ z<LVWzm2b^nY!4+>&mkffZp|N)>W>i>zTD==pYxIU$6UsFYw@RYhu=RGwnzbI?h-0x zDGdxsXM1@O;tGhSFw6Wn8xd72yYYC)nZ+V~SWV5BWj;l>Q*S$fKp?Go)4j6yHq5oB zqfRnh89%sJf@oehm>|T_uh7+*^Y0;@4|=W-`TP`jzTYwnh7b6?!Jm8%9j0P74n&Bv zI8W&Rs`T>I9jJM?qAeEY7Jy4j@jN-VxCmDA8JliCoyBVSXBJ0X(_=45_2<ak_46rW zw)1xbuf|~=;KzXn@KxA6%z1p0&d^Fc!JjWK6VRW<&dY258{TZp3<C?q9xTT1<33Jj zv%ZG!nx%0a^}g_@T}Hh4?P&*rF9fI`E2}jK@~vx>%S+k#k~c{WVP1~F)SQ5-O{wWu zv#`=9isAF}S$mab6&k`)-GI|RU+C93w3^aV4}PA9LMsjRV7G_|yTk`VSG@+l`-&f{ zzngOT7QQO@{8fK<4H^xJjS^{+pVU}7AQmLt5teyDRjstYvBybY5~}{5T`Kxg9#gI@ zKhQN(ORQj*|2)cqFzu$zgxl9|rC-jz<=FZ)@SgUlxOsAVheoqg=<e#3=26_wVhJu1 zpvJ>sxvx_+-AZqJHoPxAOm@w(A6$2xaO$li=+;H}C1UVxy`I}N<WgU6KRmMB6h1B` z`^s_23!6oVrIuX-b_~=v*HTXQ7Jo8m5wy&{6n<|N_m;f=y&hbDnQ#ec=-><LeVzW_ z#lXJt5bEV%ARw9lxb6RMbIbZ4bDOMdy}^tg_U{MWNWz_CF$=453r!USCp=!Q2rd1L z;xtjkLsq5)#xykx`n_b6mh*G518!^@k{^yUbCvhE5hl9DFkjK8wsKrChAq^l!<TJJ zV&Ir+b~-}C>(UFY_F<AQy5v=A74h_&Pu3mkReyC!ICOymILz92T2?)wF=Pyh*^}ky z6OVS+>IHj`Xh<lVw^?$hMX70ga~dYL@h%%$s1#0K3&St%j+R+UN+9e570~Gm?A($2 zWxYy4xExm(#>@;#&gw@^Q7$)JBx>Q>RN)5>R4aYi<gIzpuG%pb$xaP0=_s<paDQC~ zLlZ_GiwAV&vUlo0j7@fZa=2}8EY$(OK;xG1)r?ASeWax{Uaq7(Ng@>y)VLwTZi^wi z9q%@=RBz3yr!*#$jT>L*eAGJc`;YR>#{Y6{((Ov8<V?M_Dpj~f&Y4S;y9-1oqcj*x z)vnN1rbyGZniiu9XWUN=y!7lowgG~<dp$(eq|d?F*|}BCApa?3tDS){#H%B&{_`iI zFMbd_c4W=Lx;qRWoh~Gp_YfB%A@8-AV@6qq0Mo|cu0((?mIINIvln$w({Rb(h_TtX z(ZEku{9ZFbrSNFRd@_+*-)Z63SHr<k9p0mRd{1_m@enKbKtF=FrFCY(pgGdayw<o$ zN(MZWLGMw<W~w9b1UUgFd&N&}z&MoyxG(wns)WU%`{9I4c{2i}_r3bTJI;M5jkJZZ z5qYE}>mNz{zwD$`DCyZ^2iUf2MaJQQ)b4?gI)<HC>N?tM&=>UiW3j%0?_wAHUxwyg zq7Sra7lV!NuXVnzyFZY|dL!O%iS}}aMgxIZ9=Vi+*MU#u0tn&8;ofgV&;A7NAwF@u z1+g!H+xKV7You9g%Uqz~-QGXM@INYl+I%cyx&IaH*(3k|t#J1Ltnf5N`*mjIfE}t! zL5*#2rXun{5Dplnm2&%yBHMYwk_y;J5<3E3Y1tp2kGvcC2+G>F!Tc^@k*l}aOzy`_ zqmgyv+(ETvk_5Z1G?K|C>~VBd=2A3t?J9Dl&PvDly}k=p6RGTnPUEQ37dLICuQB{h zQq-f)9&*%D$-|I`w`w&Q6)?6cA(j(KDAwFE_id|*dF*PAA>e$>lu+KBM~1e)D{OzE z$dLIgwB58T_M96hNotm0be5t#b=RLiO=FPNDfT5cBesP`X_%tSyl+m(R8C&b&MtF% zB+wobh_xm}O51X4Os5wVKf(vnts`qCJ5H1-P^x%h85R?zNs?F;H07c&8n*l?Me7q` z-;%88yXTL-7;*RB)ui$Ytx7XeiuX{I>QQqk37*(@8V(seRYGBMfa9=MPoj@-$Lt#j zE2bs;qW*dRXe+)(t4#j0we=Yhjm9>@d}^u62yP^OJt4Z72BWbnVWuP+ZnW)J()Ol^ z%!nsxJ3ZOdkWx=vVy7z23aOZ2$yKvd>u=*q@}fNAyBgO7;#!qM)=>0eYYHON+L!^u zAk}Zn#xMk5u=ocRyDGbo&Jirt21Xj>d9$zGE<bL(?&Dig4C0P>1S{ogW4Uwzfw5LJ zXR)SS{@;|sh7)Jm1>4Y0EvTZoLiU7+LRJq0zkX{qKJ}2K>VV&a3Q6Q*wHkb8_8%mw zF!)g0j7exNT+5oKCf!d!Ri*ZVZOhZ^N028LyTZ5o9)N$Eg;)}K+RaLc@P`G5xw!W4 zL58?M;^gWNcNK^mwi@l*jA6B>XI6VPXm>%pI2o>Nv}0O=$APY=wuOFa0(ZCge|YTV zYDDB33!H!m-aZ9B%lCal_aW&5wS#vg$lBzlm<0r`a&|WgzFXUnkcNDJw8^v20e}MV ze026QU0m5(2hD?>Jc3!-R(u~N&w7?IE2;!h9Ta|67S9ILf&#ZA2Z4}&F4VVZ-==$8 zpCv~=HcOp-adLg<lusDof~R0eb#Z~;`}YET!l6ir2fhd0I=i3s_F8W(eexxk&$_lq zatzbf_LeTK+ceig0!<v-f)*Cy{z2~`U#4l_5dTZ#2(85TYl#0Tqtm}U-~U2*HzPv@ zCwnVXV;4pvOS=nQ7z351c7f*^avWhWwCD&3*B{a-l3DyPl*}@e0hFjx|MCi1r6>)_ z({U1DMtnhh=hxqcw9b_)fp<mA*~=aIku8`T>jlFd56v~c!_@=ci^XuBJN{SvZ0?ka z%rG^+SvS-C&w}sV->=u2U(TzKnL#Tp+k%g-DMm+Z1iZ{YV@<73P%j|tqNB&_P6EJ# zgM)+M&kgS}DgymUYoXY;{{|3_ft5`=2RI4|#KT2H-vz+*sThzNK)?n_^`QunFyp-7 zI-oN{{=}sQi-j%$9R(Q$8wI5_;8B2|0+9p3?~^sa6!LK;&1nqe?$gPGpoWwKn+kxe z1LN#N{JG>{dUnsu7OMj;qlxXW*c2KUdUIiRoPWg1mDtZh3b~j6SGK$7Gfi;kHAs|o z(6waGnfspkxlBp+$tFM)^%M<I!bYllyZB~lQmGmLScg#-YS$9heDkJ}=onkVs#zdl z3_(ql4)qa13t9qt3Vnlhn8o)Nl@{<&jugc+Mh-P`9<#V+D47GtowNQe#TXAi`{bVh zO%7i6t^6%}`8@<q4&L;cdjL5e@b^t{Sd2Bu1mw{_frK1f^qb24?E3+~{5$c37)R){ z7rzYT2$TvuSqOufNF#4Vtm_!GONavqZWRz(@ox0!Nt@FwMBsp12T=ws2b}_jAVgLN zuYp(vq6TJOwhI_&gn8j{;H-nSfcg-En}uxx1<fNk0Z{}0#lh!jLNyB+62OcMx(v1@ zgoF&~73kX6E<_~M{pte*0}2C4fkr@7VkpT(SJ`#ANjzZCmIo>qGFx|ydGbWj5IHYh zZCj>}q=G#dVTkTtA!z>-ESwU0A(uEr{!kg6A{XWbz0^?1gLo*(dxZF@zpKTNs2O_3 z*hS^P<{BkKiO5pE+wyyM=^Tr+OH}nx#+z<2WlHTz^^u!0Kz_FM`Kj-!7l!c7vOaai z{Io<Hz2poBjGhxq3Bh(rf*nQL1QmJiC`BYZeD?AsJIV}p7as}N)sStW9p|$~($R+^ zKM8r*%Cwg<mh$Cs2Hom7QT5G)e_a43|JpeBzt6erbV>vCVvbGWQ^0%kz^Idp4keoV zuR7z}3+Jv3WbBh5WheP7<?h}4l}u^vbYko$qa2<yQXE#>N`4A*J!=#C$C3L_`UI}o zv}Bx5GGI_lSuEppaaSqikQ3N*Ze(Ls19}|{&D(}|Xf9bfML(+NmpH!UeC3I2r#A?0 zBoMw*X2{|yuMwMcF)S`08ZU-zFZ@J|-R*luSFm6Pjdz<ygw%@?Ukd}zgxnpn<#Jq_ zOIxXMuVBq`YipP1BBX3_9=>hfV6!5xo-r-NAKri0>BDs;E`AdQZsNCpOL!bx(dd48 zFI?}q?>s01@jQF0UjB>=5n+$KC8JAqP~uc<Af@iS081{ve_rOa*Wdc<&OLRYG{GPI zvZ`AwVoQZS^*p&RG9yLG8it*WAxDy=jn*ZkHl_|;A`&Kx7ANcqNRds3z7rSAQTU@= zg|hu1-3^-&PcRfVD#V(|w^nDGBI`o&00eSQz!@Ps|DIzFN)-KBz!cTEl0Q!HcF@LN zmGJ^`Id?r*dW>MC6B);c-cCjkcPfJx_T%zw$A(=Vu7vv%c|*?UYw~t`Mo-mA+L6;* z3IER{@qpDv-9Uw*Yad30mZbFM?a3#<Q<6-%UuFBI6{8_#w_y;LuvG?W(ZQdP_zcER zl^--BrN~MnWyne;oFqlCt4N{otRbQCM%fkd6M9rcwXo$~1`YAJ`M4damf8_~QO_zo zzR6>=5vpk+{NY=v1}SZwL!+ap8?p=fOI1tJ!cgok{Ya3w?%*R`(?6cG!rXc}-Ej3P z#V6%Yd;gqrU(J?PGNampcrWiu?H|qIl8HmB{6xLI_4A3BwD{33Q*CSe#%@L=@+LBZ z4*e&18FfaH|6jY(Nc%mDJE41~Hw8rzGc7-@LF9WT+gHm^BlyEq3qNt^R@nvj{@Su{ zKu{@b3|cU}z_GBK1rIA_qUZO&q{ej7Tmaz!_qr6+^aNrKTZW+NPrr`O@~IK|60Yx@ zl?BcP=TB>ZU`DU;V{&?QE}@{Ibm6azlb-t-G2P_K9Y@^1_wXa~<^4*D0(^e+u|=%u z@&pqvk*k*rc!JYkD+GeCy+iv{=ih)^Wz0$<!w(Z2vyX%CiMZUJvmN=%m*8bN32Tdw zR6Zwu4|sM#Jf5$H7!p^lKZL$-rNqTZi^VZRFZT3>|8o5FMGCv?Wk>a!9$DBw>leRP zC*)oAf5E*nhqL<4`%CYa!ms2~`F$Jgm&)Q9P?zHT;_3Q|ek)x16Q)=#NVM*iSfDS2 z!1eZMZZ$G^4v^(pO^Cge;UmWx{#tm{AmQ)&eke!y$*O5TWAc3pck%5Fa3{F>hQ9uO zV)FbxKBsDAaW45R{jDMOfp-B_MSE48%z&Mtc~PYk#2+C0l)K&jGWxW-wO0B$^mS(Z z{tjlEUX{B!%W~5{`K-3IrNC3;w3N$$W%iR$5+NnS<>;Jh+R|NE5U@6ksnRetgX<`$ zbs_i{F%8?MxVgN(Pk#WDNL7|&WN<|HaAL{l9@(=@r3Gbt8LO0rL?G9fi40Tu9=6Ko z6Jl_{%;?t~Dy#^lUkq@XW;F>FS7cy^aEJxQ-$=<Prh+HEBex~eih!*w*(zz>#e(C` z;QSY_=T8ci4j;;=g(<1?UrF_{&kIi#v?SpZF-^@MxbvUbd2MbL0FlU6oI*F0ZARWh zY0sRCzt)F`PKmPbsEnV1{mGK<ac4A0J-=zU#9kFo_t_sb^R8w0b0t4%xAJa#A-`IS zZ@-aGwC9J2@)@4kzoZ_-=4agldz{if+QRZvOQqbA_E?Hv^y4~KL2MuxAsNBH5ZdsG zw?Ck6egk2F+I_8obb<cRHDH|J>fk$|eft-kea(H(0p|hVAT`i-Aa)SU;l^B`O>v!& zfa*SmzUn^PzRNy;Fn5sn=L_h*>^`KxPLLDgdhY=5z;}@K-`*grh#x;*SZ;`K(10XB zLm&!J5y)Zyee?4MVipM=Bnhi(x-qxnsr2@PTgdku2Z+_j6sQjh0PBVF1HS;v0ndT< zgNuJ#gzLKoiMprY(}&Un^+Wo=x`779ft(h0L3P!<`J&0ceSE^a*aHE5gg{#$0Z<pn z5A+2x^+EUP84%qF_I-oB2lj#yff0ccf#-twgMT4@fWJ^VfB>ZxPW%@IEsm@u2PKq+ zb}H@qtlwjB9#ZFLu>!Yq!?vN{TF`keTXGQm7I71jlX5G8-xch)HsIY%2;D>vN9+%& zLi$NzJy<TzxIC($<4&vvT_#*)`bJb5nflEWZ5NnsPVy2|>x^xRxyWG$yyMi1cM)$8 ztk;%gCur<Qnn<b973+bJ0G!#!h6Z#coJ9m>NyL)`sl)^@tvrpGJdwyim4a-#%;J|O z;c?LMBB<7&w1QjWp2YRH=|wqQBnhGHxjy%N8kJiwaJ9pR<aJ&9S#x8eZArwsV_|=v zC@VqCj`&8EP8=tPzeZdh#B61*W4-X*&wj~z(XKKHuz7eTzjmB!ty)RxV<3b8ned-; zVo?N9@@rL^^<+6VEn3ZqxxXv;u`sX{FW2I$wGCFtzFD1mV@_Daq=X!j$@wA9YWTz` z?zGllLLECVy14iUWT$8VMu&HTB@bPz4M4*S_|EG4O-eNydGq^2FpE7vBE^y8+{(54 zJ&`yM;rZa~ZS}MJdneeOkUQXZ3Y>__!GA%sS=Mn8%DX@_9!RN1cWqW0KuM-YQL}Zy z{)+Nxr-npp6gD+&DD!!&_Fm{$R5GMGxS&HWvXQjnhvx+|e%yJJGviP3+>7q1;$&?! z>4{61sPxo^Yl_1nD}kS$MD~Y6n1Lc7rA?7QGebjJI@+An6zaFb7}%rb&^=0bnR>o_ z#p1Wa?iQn1M<@!CBS`bM^dpRA12Aw`1qy=!e>By6qJ<n4VPt?Lk!uqpaJ^|94gR$g zFbSBfme`+R<#=dmlP%g?CEIG6IIxR>TB*Nnlw&_Q&d5U`h;qpsm=m&AN!^c4BhBtt z>t(NO3T$?R2Fo}7+X9<{b`}B2q)U#$@YtvzVW7s+mfZg5H`{NS!Q+M;CHpombpAQj zvrG7wPOcdp`c{3tHMUCZP1^d_2<<=e^aSqZ<Hj+j4A1UQ((}5*LzTlN-7%HmB-xgq zi={HU;%v9ZQ$&;ZVduQwox&{L;isEbw%LM<xE7T(@NONktO!ZrBKPpBzS%t8T=?7o zwpM?`gRYUAFyo^`HC3WxHv+DeCi+DpoejJ+_2zm*gLR)=eCJ8<s(y9_?n&1AZ<7^8 zxMOU})a7??4fa`Nm?o&tBROULW2mh#Pk!h3-!WTX3V^;5d;CU4TCdg?CjoTPEj*^$ zN^QZ@^u05Xm%G-g+D1dyI{psXSo*AUxaM@s*50?ID;&?Ip>p*D_C<5uV*zQWt&H5v z+L`EjKO$uap^uf~^&#b-NtG#JG2dVopX3kzRWvN7chc{IXu!WMXc41qFwSgzTLnHZ z3G=NnEWRW-Y`?#yibRYRStCDl>n#Ln-UXW?xh>HTzi~doRNHl-e}8dWT)GS2@36Mq zjgI1_r+kYQIr~qn%8W`hl7h3ZllZ2;bG%tXHVqeNJf2;lbYo(pJp>8I<dSi?eD$kt zLr3VUG|CS~FO<Wo-fxP==!&GqO3yOjcIq4NRMjI{Tr9ANB9GM8UoRY(H--u1)%Vs! zj#iIvXvO$ejxdV<u}hw$(XHMFsv+n)T$ko^QM0i5cBR#Dm))hZPOdz~`!jTgtWgrE zDn?J4nbWiol^<TmTZWg5<b>^XYwineQO#)y{pHKsU{bVo9V>SDV?tNMHQ7WeWOsN1 z&UM@RIiG@`&R;{8ikL13*rZ?R(i?hER70dD&C;YA=*|+W{BaKa)+!e}n|<7`*3+si zYzCw7O_U2~r>;qD7~tFIZL585c<?a|(%BBo-pf}YxPyA(i5IJchEgM^d;J+NMu9uJ zLxWX&=Cqdnc5Bn?l3Ns)YDk5+LT$+pxK9qhvz;l?xJMzYW|B1!Oq4erEdF#(U$*Bl z$@uM`Aa*ZS;WhzY+%8|hhcPa#M;RpNDtQQZ<8{8Qh!b_tE6FZG?w?ffH@vht3S1NA z19NCHNeikX@-smR>h9Yo1+KTL=c(u_PO}~R``qE{;a})n_>pLTM{ySWA;7Q`<e3YE zUI+h9TPaLZijyW}t(2(fF7zODGyNN}u)Jm5mCbt~uBWfkK|n;!VX(*HSm&keagC5< zi<UBI+QZM_anpDZRM-#so6Mm<rVG=0uh=Uqm9}xdD{1MEgHjh0jj!_Sr;8Cg+1GO; z6g1oAbW=*Re1WF_-!AzARZy?_3Q&Ga)xeylGa-4c<UFcsvsvZcLtoVDuvy;fl-Wct zX#V7sZ7P)71ztx`{^Xv5WKonUMcy|2!X?FF>{Y_Z1MHYhb|Uzzve#MI>kpm6JFqZC zgmL`HNmPPs3?A`M;ZcgHDRl-@@otf$(pHoo`rT$Nnq-P_?E0vQsOJ!r+2ba`Jj89} z>FY%=;|hh(3J>xYpFD=tni-9F3i`F9=kc%w)+1UOJ$&O_p=hoySyq-fsxue>pKo0J zX{}iX3(&X_as~s)Kmrh)iRqaMMinW__m$eioI9KTxGoNN`Z4`p%R>>2^b7sV2R%n| zIkB-$4wG%Wv~H)EzflOYEAeT+R2Fvcle)G9JAKI$Or>3Almu<Uz5JqZW8eJGxM_QM zia9GK6(<)WeRGeBd~@Ha8AA$9XZ{W->P&ahKcHR9iFjqTK*M%&>J<zNX`si+JSdTx z_{ByoPmbCk-byCjSLO+_;W>qK5?g5grQHg8o>XC}8OxQggBC8(XVD*ey-eaho_+eD z&E`^Q)8PhVk`~5bGF0WS$1J!^T)!(h=SrVPEr;o^IiT?MIO}Hjdbt4n!b0+Hs$-1i zh@An+W{|af(0?V$V_G7fp88>b3K^`_^BAv68HuRCD9iI`H_DH|5$-UN>?SH*U(W*m zj>hKjJjhq15a*3P(B__x7tY$526GGjS|N+H+2wFLppWg&MBcgBuTF=w2EO9vmk3); z^ZE<Uz_b*m_17!{)t+XU3&U9^QmabR9)t?RmsuatmdN~{_-x6+U~{Cn%=yAjDcJ~$ zpk6t-N3vpKfV{j?jNET!Z5N3E+q8A9;v;*EiC*+pt*oGo$Gc2=*2I>!Z<8^uT?3Mx zfN&?*rDBBQlNtK*rXGu2(aGm@QB?@f7hP3M%N4<C<1y%J|M>wV226yN73IS0$cI`L zc>#-s)@X)L#Z-Ts`BHup*7PjP4A3u{L;m3`MSH0gs+603z@j?X$Im9Ga|N`Yv54Yy zE^u;roPpnJr0N0>kxF*YB-1V|Z_MkZt7WYUV@el=>B#zUPbXjw`Mfb!D9&)#bLG#t z{3EU0%EeKa@iQ}k0MUG%DWF2p1ndDZ1UK?X3_ZCNt^*mI8+802V^8z?(9T%4S^$j? zdJE=IJw*HSt=<kC?s%?f?h5Nfl}B8waLS^=lk&9v^G(4hFR_FbF_*N0h*=+aCiIm+ z8#UF&rh}>-czvY{ogJj2v{H^r0zg`T@fV`?s?qilEc9Z2N{R*fe*cRl{iVE<$UBvx zyaTF<Fp1kYCAhLdZqzY5P(vU{ZVD@ziRRIbUqE+)W<HCgQAPk$2siS*9-W(_XDyP& zwSdE24Y!V$XER$R)Z*8rnm<nOyQ#y?zqidC9PDC?#6^DF_0o~xB;vCXJi^qp8!0AJ z0^JZ`Lk!IfvLB#hz!wjF0FD+A8Gu3!LJ2AX6598-&jd(q0PDccjFI;`!p;{58VeZ< zF5Bmlhn)tBkCc`7_KBO*7=RYI)`xGvr~sA(Mg`jBWer*tfKUgH*awpYdME_RjPw(y z45A|7wGXThvJDgo2w4G2qySD2VS(k(;sEPFs{pqMS_&!?@Gs@ffMyni9ZY5B(4z*# zqBnILFmzA9^ySJ2fKmgU4j{4s{eRl~%AmNmZd)`E9D)W2?(S|ug1fuByA#~qf(3UA z?jAh2ySux?>o51#%{iRh`@UQC{=Mujx~i*LW3L{2?(VtfoNG)~Z)RMeYEWz-6X=3{ z(m;B3EILZ%0@|Z+<a$@B>C?B5uyEtBwE6{hwBuPKx<tY>A*Q=)KRP<tL}XMDBjEUB z>GC570tzS*f`Z-X!KxAo1r7M+?Wu!)WL3jNM1K&R`bu0B8<SN@K+dkqo5$uM52x9l zInT_vFg_g*Did8xpTr%q&ky4kdoa*%SHmt*_9!C5h8Tew(Wfpb!_{5v0m2<TAI?v% zODmyIWm`1T!IumXGbc2r;f>?ZM`&WS$Q@_k93Hm%oJ6;qk4`tWViYa4z9znrTgLAh zOa+G;fPOzQOv7|9Y;$dC)P{Kznp?KubU@v_72Tvo(1yTIX}^%N;I#9Yaol@6jaf8& ze$J3CFrsBM0)2j{q#l{nq-3DwdFI>s(h0l6^|X)#oxpkhoyZBx?a>X6kLP;s{@__( zZFhdT=CN^^VYp3Y^kk_e|2EsraoW7W{30U7v`r#KkaH=&7W=F{V#2B{RO4bX*Q;gF z%R6J`{&D%}u1$4qS-$o3ghgVBUS@PpC4#<(HIySK_LQk<DC|V_Yk4`-C`O8&eq;02 zhnr7U=qd@Gg0%WJaOx$h7SpL>^2HO&5XF*OJ6222pQ8Z^e(TL{0j3m95Vh(dd957o zaMgeiT4!G+EZ11Awb3fakiORV@&|KRtT{7%7ngc)adnW|?M{|(b1Zt6AoyeQ6?Ws2 z)%z4O8Uo=Dj9%7H!(R(+ny!xgmaZi`ExlNDoSGpX2)UHnGf#`_-@vV>R(*JD1CCYZ zC0E}b#G$@~SmyH%Hh{lwgn9kc+*kVh*(cx3h`kLN_X0~k_?>vpxOhJBe4_(;d|M{0 z3wBE)iBy(!LD+B3m0Xb?;*vSu@0<_h38qt9bjGXpg5{;@Ei;YZ9r_JO4&IUCXNLvw zdnt#!_n8)hb;7Gj5tB*WKundla#o~XkyuXj7e3qWS06xS<A$|zkybrGOlD4tBv^76 z=r-Xb;i5s6BeJeZfwaVOtw`6!mq;;Q&E&R4+Gmb_PKz^PBq5_{nf+h|y7N%5lFCV` zJ$f$)9fu9mTxpZJkjx~9Aegle0oN{Lp*jY&(j|h#r>1~yd~@R{G~$93E2K(-j~$ok zVFMDNH;XY;e?pNh_wI7_4Mb?B)~5}=aWNTPZo~=+ZM#)deIjcrZIXUh<6`Aa${#9E zv)DBHrM7U6l7ZX$^fSQG_1A^x{lbs3D^ko_FoN@sykXP9#)$o#Y6<1<d-46}Mg^8} zt7_{mU1|3iiP`t_6nAOo=}u#zNBR@&iZ=lfNi}jteFLfYBvR{L_h!mp-9>X#Ig&({ z=yLf-5+fs|#&sFR27B*^-=5Ejvox9D9iN4+;qx@Xvgl((j2*O4Dt*P3=w$Nz;W7qe zR_O9(+jrSTB8u_LX4=BUq_w8hB=T;hN{T%&g`ry-6+^=0{kQsw(1p8E=zM>E^i45L zk0{)!j%hKEFL=Gm)+}aAt4(Bz))L~XCwkNS$>*k3*m9hK`8^uCbiIJ)VNx|Ch!;nW zO?sMcl!J}^XiX6gy!jJ;t=%A864I*&P=k0LjvC);IKsAyt_uovYH$wZEhNNDO;Pp6 zEo_*lm?+7ZlxC)f>lf=QhV^aO#D?;Ss16{+E9SMuSQuYs3%Q9U7QZdDN6cm{6qa+k zb}>`0ez6@|OZN9S7d;H>nT4-N%{8LXUAU`h9K6r9ZezTf*(fa!OLugAA%vOXjOm`( zq0&i-^bawfUL?uavbj0bjns+w#6EAQ)3r}U2#dz8GP#;ufw!tEg0?BC5ocF=Rh_zm ze6*x++4s04SXDZLs*+BYrFHGg5^CvX>J)u6hrhreI5VNAG-e0;uC1+USMs4s@Y*P3 z_og(f=fLP_#Hl&`#F=A-FS6whr?K^})vCYN;Nc?eWQSvAyRGhI$jh|6{2<2Z!cNo7 zo-8e$Nz&P_WpX7XUVLQgSje|^Vt{!6q_n}uqRDNZq>hVy$)6}qO)ZQud1TRjXL1gH zcXGj(EbjYE6>>_YNB$^s4QfcTkOh>)?DuCFe+8sO+|byU-JBN)(_wYfYuWUHWt1l6 zV5i&s#sJJ4lI3o>G@|cLpF^iEq1LDH$>wMGo75ztG>zEPO<wM3+|-qpPcJbdBOmIu z){8p_2EMJQ5vJ<W2JUR^XU2q*Zl)LpcA7ut@$4Pv*O(k@61bH7=pPRpuUg|Toyl?d z=5B%Qzv8Vkh?_qY{+&{Us0$-jCWqK$-lrCJeC9}+d}_vCLYeD=1oiw~c*gwK&;w)` z8{L?J0Dg@O{zPeW^z{;<;JH}c>TGHzJ$1A&_4rIh)gmUUsf2{0g2imH5XuW@Y?CK; za0O&OY}p|5(D3N*8L|#z9-p&8_3z_tr+6J>qPa}@&Kx@G1%+MlcKHdf&$1xh@DA|y zqk%F@CVc|SOOPuy^(Cs&Q^Q3!8_bw2$W_xE4T8FBO&1rFr^}QW5}bp;<hM9`de^hu zKdpZ8t93UWx9zsx5<?cjFkqpA4a}AZ!PcT(5iKvTova%t8HgIy7CB|ugo<1)UHk5l zOfL`+$EQD&I*wIj7(rPrJQ>+?JAEP5LLc=lVcjf!04pT6uSjbWAqH)V%G9N*=g1-a zsYs_vwdr6ZrY%iDfANeYHM97cAYgv64?3V&g<AkU!R{89FAR~<oS<II1UR0{j>YWW zCW}<m0@ip@a^`aHA%JQ?j7%I~mEh`oO-4TX5TxsVumvAYfu!NbF)Dx3a6lDv>Yypz zJGqS8A7WS0T73#d$b#Vl4Z?HChCY%`u<`|!2Bq^>oEKy)Kv1sx4jw*VU5_d<@--t% zY{!Wh1V3z2Z+1jcL(}9;I&lV3hv1&qPcXe-HMQL+v{+IG|9%IXSxy|iOofqJ)82B) zy54bPK0mpnFb}b}o$chjb`un%b<aa)Tqh9DBAKr-3wbRi;7Ul3iD(49!XZ=4<yLBp zQ9lk<L3<&Pc@-(!ZIn-O?Hn@ii(ioD>W;U3aA!zu0X5+@!HW@LB*9{WUhB^0U)Gc8 zFI7wIrpDwrCOD}Zmj5)CD}1MJJ~-MF$)kz=Q7f$|`yFPFM!}?jLGVa=ij_B++IQ3a zas3)<0YPT`H=+c77AHV2R%oj_K9Xxt7TaL_Enj3ykPVqr%f{}|N=dlQgZ*SR5(uE> zqj;MQ$hxYXK5-ai)$EBZ=F*=zySlE&<A=1s_`hBE63m$0ZKk4;_LloZ_0_yuJ$*H1 zY#!-8l$yzSUvB86dxzrqOs2y4;d=q=lX0le+il_7p|LOnIj1jiFY){$)O4znI^S&O z>AqZh>gqK^5kB+CtvY~dNFI@z&4yFZt%}~QNZuOv7vZT!qms-s3Z&r9GYK@|vNH+< z;<C#{Tov$_#xdDx4tGx}bfd+&4N^vg(2I#UsDzh)x2#AJBR@o0MoO7Fc3fy?F78sd zL$p^a>MCu~@RqI?Rojm>r0&#sUNm+mCHR53SVB-gCuCg}xw?(sF_g5KHz};3!K!kw zq5k|58kPr-e12!D%N->$%YYl7F3lAoaX=UiwWk#oMZ!mA=m9av{|>}>GE7KrvD&ky zDw99L90n0f^xj>Cfd}nHPcWm1mF@<piU!JFsCiAAaTG0AQumDsTZEeKP=ch7ZY~sr ziN$o0uA2XlX^qdzT8CQQx1!mD(}Z}JfkTuj$<T!?aEDA+DTRvOy^55<L5t7XdCiN3 zKxNP^3f;%_bzQPDI;CfEDjrh?%s&2>3PU4YW-{oZMy;*JxRa>6%v5kSfo)&YF;Bb_ zilEK6n9XY5J<X`e=}=cFH!~4!aW|EIH+M{FK|u|jJG9oar{CrxF>|BR-QM7MRKf%n zC;<E<Y&Q@;78h6#GtHGIG~_+1=b9<S$B!@}Qf||8ERvJ2TTYXQVi0aYhzv8al|<AR zmk#KHiriHwBbIkX2hC_2*nG};RK`-QPl)gCVrMK=jmj7bjxs7zu~k$VSkt<zke*I8 z!~|F9PG{Dr!p5jf#H)4+5?w~?=nR~A^f6uiy^wHy3JQ57pkU7Oe=L||_*qV*AOk3v zL-52o<bgcGqY=?_lSnCY%sALxV4r>ChMvosfK(js(9>8i>GdR;Fj1Y2X`l!13r>=7 z!tK#`a?BvoG(dZ=N-`g2^`1gqrik1rRum1ERe3GvvuJhvTDcai@GT0GB83JiqU8hE z67+RJnREcnZYr`L#f55;J2a+#sWi5sx=6Uk{7y~D`MpR$Kh{MU&bf&cokwSd@B(k} zf+`&qgXfaqWkV04C@Tff`@PGS0cRu}!Lj<6tuR(o_Bp)*O2r*z!{+yZ3TKdMKehaT z^<6^qRmBR-(do5m1hZTjrc8RFTxa}ROWVTuZcXw<L#!unaOGmCO9#W|%qhtc0$_O4 z4-__mWBqhSIvXh@=N}QX-=_+}DsGBF82au{h8CS0mL7dCKnZcNkNRQCiLDD+69n(R zDk(-Oh&5!QLEb_w^FtEk+B2(EJ}t<tR*|D+QfNalPF2+4lX-49Svzm6jfK)FHI3yA zta*w4%JvwFQKlMbr;h@1Wchu-L~5y=vUiT*CZ=B?s+;V4%nIU>iEHZ9mb`1xdg$4( zlRFGbjS3TF7oFX+E`850J}Xt}*fAK1qcU@zmDvLw0;W8q`-g-qr0!`1Cwx7!g$vlD zC)v%KYHMy>syT84*Y)=EA>PcKPV0-#X3ZpGRgr_3x-a9H>qyq$FvK{{fK{&)^jylq zZ1}#xdl3`xyyLDCd@Dg3)~$y+(nI;aDyPlxlb)k+h5dYRFK-qH0T<fjakrpl>OOFj zFv480lyzN@dk>9chuyid0l&NQLtNm|l)sb1>9v#x*9*T^Z`H<fPus@9W?T6So>pwi zO)D{f8&T{Fw$@GgfS3nYF3()I6QBg=pMuy^tQ=gYfOj+iwoLFJMLK`!0;6SXVE)%$ zFx3BgT@uqF3Q7l%wD&nd*tV=J#Kj-Fmi5!D;pY7@d>mj&j2ArB+9C($Al5Brm2}t7 zg40E*grkGx_M`g~eTD%jNPqY6g!ZbzgR>ocmb_S68h)lXD5`SJjM5-Us`-|b!DS(2 zo-`fpbn=Hda?9`9p<!=>W4X7rDH2_z`y11+rKys)CAU^!jE%%A_Fy(8oS=2xggNiL z91q=2tg)BWc#DjLNYT3J$bVG3w43q{oa7}OLhU#AtjD~6pe?QP*ufqVXe`&1;b};C zd>$0h8YeT5Lb+g=9?-miP1E1cnLeMIkdYmV$a?y1Rg-iYAHu&y*iZTY|0HbWUxWqz zCt=~Q2z!pgDE4ssFT$d&xC2zmxH}22mHYCIwec_!&^7kmv0$3bajcH{5Ef?V$OIX! zd3_YzJZjuK`)!{DvPTHXiDJB?*?33=0fdR<g;fxHkayTp4Q5n2aByUdOW5TPSiPxc z$eXo)GmzH=Uct4=4j^nFKv>$MT6-UWu*6Qftbu9+US`}aC%M6W0AT^GIX&Fs2f8Z< z5k~k+N|BRNWd<HETh;YBa@5@P2}b&5zX<Ei&Pm1=d;o0T{ilUZ02@R9TZH`-UjGeY zh6~_`Cbfz9mqk~=y7Wz$i;sBk{y|tN&m-Y=++26YzX(GXucn59`sxJF;)Gkjt9^wG zO2&%N)!5Xk$~B%N_VhhljFOz_n@>fYePo{@VYy0=Kqj%r56wi?ZQxkua#vh+<{_%> z6_{N^@roL#^XXb#$3S_KmyrmM6BvsumkZ>9Bt+vZ9wGH4c6fPY;kzYs0|#)<@`~o1 zo@J?&lAnZK7x&)A)>XE?*s7`nEDQ~5i-bJab%@12W=&*i01Yi7hLt1a7$6_@Pr_8E zt3}WLom=v!Lgw!XQ%l4jH>VNYRX70set2lPyA9v3e6)~)Z(Jwee3_7hbGNPcWov&T zTvv*0VHM?vKDbZw+tRBRpP8eHCD5QB^L4H;Uy-p$yNViGRJbP!8r%nb!&Q(&Mxfqo z-v)gd;*~VQz`u*J$737H73UJNaBL(=Rhj~sAOlT)D`&0>ZZfg#JNrR5>5+%SK8>jw zZhI+u2by3M%qh%PXMge?Mbvq|%HTSbvx=UlM0gA4jEGkpwv@p9iDu7AbHdhB;>2OH zB%R7goMupivjmHM!kUW4Aj?5y3@cm6HbA@opNtL8$C1(i)I`G=fq?#1{R?PiWewOe zZf{^KXJBvtbI;JL*K|MY4^^bXFsCqHm@~yPAZ1iDkr7lgi?Sf2oT@a&Nk3N1ENjp_ zGhBs<BrDpIT&uawo`o?nXf~5&KxUFd$SZs%hr_a6mA45R{UH0c1}y@zMh2eCI!gh5 ztMl!Zbj_28l_BSd%EFmXrx2}P_|DVv&5p;)$rE<d$M+~EZ`1Slz_xaKt9c{G-k72l z2oT>9_GL@BLlltlf6YX>N9>+i#QwCylUGcl_`>L<AF;iw<<!)x>b%V6B)M+)DFP%u z#{<^`RsvAKMP0eMgSfaNxaL(q9eQ86&!|;0iIY2c@xpq&wgElD%>8A072mI7<;gWW zrr%YIyKb^;jm|pCvxQSCCox<i;+p&YkC%`HI$uv8w(5`IKeA!Cf*WE4@4!R{Z9Pk~ z?2&Rgsj}j_jWD)QcF)p$nY`M-A_enj;B6E59kh(UZh>`92qGv;rkH|rX_3n;RdBA^ z6>c(1s7ZI2EF`R$@FD~$^27h&z0HC6X8BO;xDpr%gL2_YkWxGiOORqFU%Ax@`#v;S zNRQvFM6o7)D&zpO2fv@zbDqV5Vm24Hcxi36Ir&th2}ik%U&LnU(ziQTT&OTBUOaEA zij>B&Bv0AFVj)fluZ%-Ng+LW7`XGKOA*1r3ekfwc5!}QO9!wbxJRT}Qra)$@%Fh?3 zL8DZa!E`Q*&tR-lT(6FVH)7?TgxR1^S5PayaRtI-cfpcrF7KV^2r6iib%Qy3A63y? zZ%c<<_lcHohH6|x=Ao{Uwc-nCaZV{u)E^5n45Lre>{=2|-HR+OGDX&U*S5wahN5t+ zQ8yQHt(#JbAQ9wCaqDUe8Wz<;Nf_j943aKBR#6sLYF)<cRg0H1nb(hkamC~V8mJw3 zA9Q?%&Blj_#f>8;o(=mt=ebHQdn`c~0bT4d!jBpPRj|mH+Gu>uWgimKQ!9_I+*)na zhIDIP#yUBKotf170@6o-@~!F1dX{e=NM0ud7kLPiZS9`<0XCjbX9k(;k8Av`?+FK) zQw8|d>ADq+OCRB=yPv{YuQ^rQkl_Y0v9~=LtpgMs(jUYld(F^Eya<Z7Yw_LyS2}`~ z&V-KXC&*Qop}kjCq?ml~dd3-GyXiQ!gt(l%leI!#29}+L<$fHC?BY=Y&Q74u%7@jf zV&5>s)M63S5hiaEBy-A{q2R8Q8X+vUUADRvSOhu1?DQ#DXI*7MQg&KhGnY0q;7Om2 z$`Lo8^k7f<AOHmhEFgr&r^KQLn<>g_fct)JLtunz2T*BuwRe7%|3T>7YFX?O=BrW~ z%k+!qbN3>J*PVw-5YMCg=ld5dqDIh4F(?u95eJ09dvqKKPZXyxub^`S@EsXXi4Lkb zuTV^h(ZT3L3BlID^P5%5>l>-n4WNo=Us9fuRjdIkiux{Ftb5&um2c`DUSO>tSPu^{ zx7gr3pU>fO!?yB2uK))<)hNT^>pr5M>a|xqe`=Q^y<KTPA$JB2k#eV~VEMeFaC*?L z<gMp5BQq;^%IQwEVU9Q`1#%+z#dlukE-m;GYC})N&C9?8d<A}G=Chs`)9MY%quJ_& z6VbUDWh8Cg4cp)e+*TMG!r-12$x)h)yOd8l1b6o83h)lf<A>8Xw7p${H$0$;$agzo zbdY+_Qg})d9f_Y|S7$xpZ&UGK{LVarpFsPr(ArD9c#b(wX$XFVJGXmhdHJ6?Aw3+l zgMz=H-3o&9fR|ANSG~aKW3MWq+^My@ttdSLYk?X*iPfKg*PT$@BCbxr2iS6wcAZ1o zHIsP1psq1QS<!ovPHF*J0^{_8SOe#a!Q;d_%D8=1cW)#$alUs>+^A#Ef|@D}WXOk} zjBzmP|45cGtW4I{%QpFybYzgFC&GP2>~eZ_$sgPCEvJ({hk<+DXo3iPD-8(ZP+q<e zYfAT<X0@P3Ur5X85a3qVciCk4p&zc+NFXBNnzy3(Vr#8p<|tK7HLG?b4Vi}U`%T)W zpTp{jv2d83X$`#YoG(^Zw5m*tl7>tfl9z@`VqkjBB1Y&REeP5<RMA@AXY4R-E(qnt zUT}RG3QxlM;mQM1H5XY)H?+OljBNb*jG9m@A$7{(8%kxsjbd1u$x%G5ne!FK)Z{nr zQ*$V|_7`*J*xU3WSNfNUquvXiVMF$p4NAiuuT+Q5HgE2WO=b9=;iEm<_?f+AmNMDJ zk>=xz<8RbT<IXhe9vR`T?=~_bXHNszQp20A{Phk?;x3M-;xSe&-i%OJvuqwSi--}Z zl}&!!KUXduN9E>J1(m%>S{%2>$xf2I*6q_s+7Xd5Z09hQJ071f2iw-*eUr7dbE&RZ zT4OvT;W(R87ZuOdMe)Koyn*V>Q@6FJgL=gJbuO!XU`%JQf?kc$J)(lH(J(j{5>>)j zZPgR?Y5o*UCj*5W!lK{7oIYES?V}5eiZ2xdJJ*ST)+ceSlrJ{w5mi{&;VD6Poz{5> zk{)a1*dy468IUf0e(-zHV@U{Yw(E=Ece;RN&iY9z`ay5>Q$b7c$~gHtyEC;oL2jS; zTAW)AvPKF?U+PbK(kv_$khs2ANaVM(lz&>oI4>rowB(_97#BQ&n)8)lOSbYLG~w~X z%73CA;0N?sR7UA^xTBrP?NUXvEWgVbTmec}POr1RNGW$kZNa{-FIR)9TWF0M2Gxul zcZnLlit~>1?lRJTi+B340QHP5S#CuP`uOlx%WRCqGkvGg-aIVg?Y16U8pWP$G=0>6 zJ!t)bpGW*oE=4iw6=34dA`^XI(r(*^Y-WD1Rx?N%Zte>3M(zlKfDrz}=g_mVurRQ+ zxBDppx~X~)vCxC^jDveZix4&aCLu;x2D#>peYA!NfkR^2Wb7F;n!u)!FmNrZI6jn{ z&+D4&{*6Zp;5iLLHgTzjpjRlSGgL#PZ!Iaw?rzv=_^_({xH4`jL{=<b%VE*7txMoT zq4IW0@$gM~qU5Yj<GRZax%klHVaeTLr4gQ&LzSFL@y&Qcgs6jJNm+$%gYx`r0ey8R z1w7kDm16o0l-7iluhgO)dTEd{1|}m1H@=!Z2Fn#l)K$A!co;NR`<)$!;N_d5d6n}R zD%c5vi8Ra&<agxVA5R7pPc`m%t%<W`O=`*-5hmphYQnzzUNv}@BiF?+PPjaW%_wAv zE_ph8wuo5~8F`7+yG~cJEV#8iWw2ODTC!~5=8P@R$SOV1lpr4z4^JM@QXC7qPkFQ1 zcS4W3w^GoU1rK^bw`lKQ<!4r1-e+z!eLWC=nwwP#cy?#ohhMVZd&Vv<`;JXadiDU( zT=I?nDPZaRvWVe$RV+tVfj+5sU6R*)cf8xS@O%Kyb?9BZ4FmmJF6xj0Ly3`qVsA*j z`*fE%S8`d>ncj#{t8F3lhmKBb5~=~UF>ek{=#C<RZzs2o!xV729fRo2k=<GA^ixYu z<F1M%Ej?g2Q#`DSZ4andHz_y}C%`S2X{}-arE>Vo3+bnm+_0DD1Duz;8y=h@&((<< z#k+T?Qz(x!-V$aVrc}C$WT=Jji3$BT9cqz?<jiVOn~O5(zohW-aZ<>R8PHs_^sOVD zrBvUa4|r-0R3mNqM^A7L<rt#0AtGdAXqn&wsU?y#=)F(cF`a+WmH_4Vd`e!9E20j| zk~-3<$qnyJ!eAsQR>{Sj*(x0xX7NrtRELJ)nFaCHngu!iSY?p#_G^^|NN?Pf6_eN$ zs-_lEy|8K>-?%~fyIUI1x9Gd<kuRzTA?II$5#oC&I~g!_1kZ$=KMqT(*naXkQYTg& z7K&f>WjcBZOMX^8RU?(iDKwn;%y#Gmtd`5!reizgSY)Usd71@j74|`|V%x#A5)?OQ z<p{p%5R<fiQ8X-r1YbF;fI<=NQoS82RQXM~Br&^MZZ^E)-f5!jUXVPyD^9#CPu|E? z^ogn!k78-ssri)8#S<$yhHKZ|V!HG2B=%hd4j3qjdVn6u<vkyQb5#c}x|cju<#v*| zQ^FVLUXjzF5D&&T=e7OILcH6dW@lYQ6RUXN{g?_)g*M6GQ2>Pf<{Zwx2D5@p#&I%L z$0o*y$!5A)wjDI;j4qU5!5?<mXh4+bF2E}%wji&j(Rn9GZV?w|c*9(}Q$_F9;0G1Y zV#*XTe{~MdDJR#zBQ4c~`PfKO@UAf?V;#2&%t5sk(R?TyDz5h7&M4iBc}ZhqoahV~ zY9>p*jvf_O;#(K&y2-WQwJ=gStbcoBkQo;fn`v(?y~JFR4*q9}_qToeG=h{ao^MHW zZvtGkHGDpaaO`==Aa(ad7Z}d1p!t)wl>!UXci#g`tLESr3Ah9Qc+PAX)Bo%?ILMbk z%29ySTp-P+z;$$(0Q=PENSY1Fyfk;D%~D#*`wa?F<SWSS=g+fq_D+?uqO{+TCEBzI z=IT{Ad=&zpXi6SYHEuxh%Cf{Azm%_F!(znD8wVL_cd0|pK6jSswqzd~%r1$_t9cX7 z8IkUNfd<N_!XyGZnxiyN%N(Yjb4VK-z-)Q}J7|`T*FGlNK(H7;62&W@RUpVV)iBBd zIjkY#XP04$p~ZqEKqMujj40R<o*7^d%&OMC7WORS<&y#VCQ+oAV;<Yg^tEom{JUmU zJ7g5GLg^Z)!pA{x4QVPv&Uh<xbw+3g(Ljp~?E;DYd)m$f6bXUC#;mZ~fR(X*q-a!P zM_u#}nV_H(vO`4ly)FI1FI@B(qaDD+oAc8HkkYKwmbEIe`D(2bB_$SE7peIi;WbGX zK#<d?1o|7F>OLN1SuYVI+4zGz2YU|G7Bx3yB-yFsg6btDH9%#af^2~F*e`#8N5ThN z8HRH#!W;3c^LZaES?X-hFAqa8f@0u-tO?0qm3GH+IK@(~CR};JB+99AE!r*`AQz?V zL+^_NJBluq89F`B+u3z3fd~?>gZIYT0L;nwhnMAYN?QC4>P0lBctR$F;{!u)Iy{Kc zx>P^}PpWv&nl5MhFn4Bio!{OGALg2+oa1s}3Q0kDo3m*jF8(-gj)?vnXy7INBh(|9 z94_Dtu3X`qI8<mVF~oEUy`G)XFYVY8Eo=72w_<49SR)C;_H@>1?Km_E%eIHG#5i@a z=f!z+!q{G@_deZ|0w%F5Z$AR}DS`mQ^A|EZW$-Cy#v16BfsC=XEv#pCIAkOndvMJk zRdOy}D}AU|^GT2F=8=3BmkIs=Q`4&&fA*ySZM06pgJ?xHmwA9dDo~Tl5$Z|U9A!{C zR<;L(dq)>^YI$A=cB~m%9oL{f($a_%In)Ilhe?h7^XTA!UyCXe;w+&DxHPQ;l5 zqt@`wg*KKe4|JQC`ay3_IwU3V7&kV;8Isc^Zq6|k1GkI!W&moj#Y4BiH`$KNU2?K$ zyvRzx7IfsYu)rV=OriNp>t}%{aM68cmO?DDAlrg1#k_JuAs;92LtiU?Fk(Ud&wacl z<P_R^a`%-eeZF`$xt~*OiDV`X`1;{p=D_fTajxq5ECl?9of4V)<HmFRSn<gM*m6!5 zx$|ZUMk8A`LEWlP)bJGGlGet_t5?SviNYI$24{=#h%Ial1K-D^U$G2*;up`-jWvXO zmJBV0geg*EFHn1in${V{tks?3&%X_@lGYOxOioI)M9kiDF^1u%!_<7gZ6&}Wfj_B? zPjvAmfjpc=FyY4x8K-Yik?3YJ48P)j<v~jaY*L<}Y}y-)3LiA-`w)3J$tFat`Ry1T zhTWr$46x5)QrKTjcH$e#tPORUEo;pKf>R6V#0;o;!-gm|m_=mVCR6;jumn;18-a-Y z^E-k2l_EZ=9S3nC27-lq`1#_`>FNkkLE`e4e6~VCc;SWbI2551vIwY57J)UfLSJ}i zgBIM$yiWHn!TXE7^;le<6l_<MR)K0MB4|M8=LlgE!^I*Yj(KtCRoHFxf;iQ!`_aXw zPf(={9lc@9;y)DC3^?2r&ppB0!B}xRul+zy)D$55Vdx0QSb6Lxf+Odd@eO*?=r(Pi zHRQHZhZGj8B0RP)f;m$CzA`A^>d_|jw<k$<jV48Q?P+skG=^rr>|Lr8n_)?D$|`cH z$}kZXSm@115Zs|(*`<J>m7?z}Wncs6qxv<)wUez=wq}oJ-}SAUso>kS^DvQfNy|=G z7~x^dact*z=m}GuktJnna#Ay;d==@57DDKACP>Xtotms*KIZFZFPHBPvk>KVJ03R& z8EnWz-<Z$bc0^NO(~p%Gs3Bo4Cq8vPYz*ghbDT(s)B0fMTdT4&^zjG38=-k`3Tfun z)7T~IVShIVn&~%S#yx@ZS`aRobIo(H>MOhoJxoe$S-G7XPK~NW`b3qq7CUBq{ULFC zPU|uI6!D`$m5;i~4Ga{ePo}Og5-%wT$9zci%t@(D3dT?FSgdSb?h9uptjf-V;2>V; z5kl?W6$fEEUF{fQ_te8aOa67vn~!Ip5neBVXiu$;?<!+jjGAWh@krjL;J6_)dwj+! zP+4|5U*F31cEW}M=Sqtej?Vs9tdf?Yph0gQ===+Rh7VV05i2(yMu*{Jyy<>AyDpvx zt`kxaKNbe6j=+$f36TWy!$zCe3qgm0=moUJ!afwFuJUqhXnzlk&D#i6Fv*<T)*qE{ zlfKn`{?OrUYktGbA-UKwp56crwlQk=Xie|bd1reVme}=NZ~~h+cy~Baweh?H>luS5 z^)(`wTDFxV)!?*RnAKH`{^Q5b@hjeCM{O>%FV!Y(dpm#*_y12~gY<u5gMVU!e`14w zVuOETgMVU!e`14wVuMli|I64QU7oS-B|zi43{XMvDmKXQvr&Y+L@1yls0HmF$L}}- zt5?Ka&rO7xj`C7^c1g@i)deebDlmF;(sIwMyje1#650Ts3pFgs<^3||2qr}vMh>qL zZGxG?>K6WYNaP77rF;+`Zh1bRXb1Tzd1G+!4;|}(=G^(ICS-hUp`6~8;Ff6z^dBkB z_+*<#u}QtgjR^`b3fT}BFfwW^fK!nPp}JtB<0v>8^kY2y1Y|E3oxGJ!3O1zFtyx;F zbqG})on+@nN!Bp2f5bb}o&*sY38^jB`;*`|ny~shJyGfDuKvjHU!~7)F{`PV$!Nyu zPLjJU&1EOK&~brvpEBUos&^=VOFxS2T=(6$tgRT(;KM>~zJut4__Q7Ti~5n3l8y81 z8xbrkA7`g9Uq&v%4rbN-1%sM-g{CK(l?2}oX!IBNt8x<A?P~~jIhj@D2S1Zc*!7|S zU`~t5R9-ogt?ev`c~cRRZ-7!+8=>sCxHMtMBk)7E>)m}hj);%Fx)bF3wb?%~ck8-+ z*g&rr99MuhhJxN_l-;n}QaKn;a^07?$q;m;1G9y{A94+?R&8b$-3E>h#j$Fh_38`N ztK=%lwtk0|?Sb-L!4+nFfgHBB*BP5=jhy?i=tZ;n3giv({d;M>>Pt`WHYF^U&jne! z)r-&F*~MxXD|SQXcK|uH{}&l+#y?)^AFuR}SNg{*{o|GX@k;-ArGLEA|Fu_|Y>U4@ z2$0myh5LtB`m<JbQ%S;VffczKZPyFeevhc~(PA!0?v`H@jMNdj+HJOF49Ng7Q`kEP zsj2(<23s~}qk6A1umv0vIT#fO$7p*f+mh#l7~h8^X7zpp+TAPAPW>W!YWTs`k@6;5 zxS$GoLwVBvgpS~()s3YGr@BING2&7@NZcWM(vnhC8U1me_{mH8ZbHxgc<RFLY;J+P z@i4<Y?H>vF@2?*w3UbunNcagRX}-(Jxz8;7SfVL2c%yFZc3uG@n-}`vE1!G}j>NT3 zeJb;UT)$bUD6~}@$QwOi_q8(p95aky=tBcPaeY&64OJZT($G?3tn+5vqlUBtCKJ!H z1B-L9|IRy>>qlP3!Nd+zgsAPO2XhP{E4k5-nww881#AIr@DIh$){N-^P5svB5&8YZ z>ZRM+-*~>3eOl1)s!A(;=2K0;uK3)uVx)dOwp4m21jx+4xd4JML=t^R&ETjHi3z`$ zu}o%?lGr-_F@6$uyLxkdG7N5scava#Fr3<wfD*M~Y99P$r#lrHH7|<{(oYp;B9ij+ zJU6rutS>!w9U)tucvK8KT(Hb{)JY`%L8!JRD_D_2Lbq++o@Gj@!uFoev3Z3C82vTm z{UTGz#`gmBrY0xN<pmuKrh5<GjaLej;XMmF^caVoQ_hS~z5VV8!!!MyixGhbGV~vY z?tPCIGLKnG_lvxnGIxFGp&i3Wf|LzciKy%_O7d==0oSXIVi0wco1}W{#upC)=hRcI z9?y+$s|$85fzt*LYwz$&2>XF@t^d35t@ox;@~JQ^Gqx=aAqOs;DVKF3JgHI;#R%U= zI`2%815?h}=Sq4&)u*ij&6M1+5owA~+6gBn6>vt<&>^Ajqcd&!@Hva9Whx>DKQtsB z70@^*qiPsDv_lwX`iJeXP2eURvcVj$iehc6zFkz_0X8~|>c~>7u!io0AI+dgiVpdR zl0@Ux=^FTf%uJ`DDvD-VWKomUeq*U~4|2cTT>;Z^!$})TAvA3DuA8~U)^UMwJ?{=< zIKkCI8zP#DHXqjIt!m=3Tt|*s6Hi>=VDHo$R6#FdA2A^&#lA1R3RCSNAjA8(qJ3Nh zw=Xc?Pl*r`a{b+dQ%tN+bGPm6v_GdKEvs>q3`44D(r5S*VZ}=!3urI2D^2+F`^K{z zldD<Jhd=vJv^%`jj(qh0LRBx{zKrSDgd4%zKKDacpmSyfN|(%mhT?5o+b#sn*|WlX zOAaI1>m_=)b9s?221q!XW)QVVcksk=kb76UYF@l#d<`AmsTi&rm76{sOK9qpoenz# zoltJM)I8h4%bI}NmVr@4kr?y5SUy$^LSs#vX6TuXBiz)Z^()gB)o|ummdR%BmSfqp zT<1&|EgQi(<2Svf;)sRPNZ^U^Lj}dC%8I5bLc7CQBe4DbpBmu;ldnn8tNHNSZN52@ zgIOffLyaizn_;$$O=;^{rup89&qH~fQxn~S5*BI4sHEK^SxgvkmtxwkT7TL;texx0 z?eCN-=fovaGBIbEHzzqWM(>$Gm$36|x}*IXcxG8EY53&8xooj4IkCt`*E+tGAG5lm zeJ{E$z7hNs)v<-SpIRm3BWLra0&6&(k9K7=-MT=~Z}`b0r*V-%LSE8Z5_S8{#HpZk zwiE5G;NH&hS`-R?ptY9a#PX0?_>W12sk^Aixi76cog#BoHQXLHBCEv*Y)*_Y?&O&h z99l-33{%zL8+WCTCr)^BCes>IBJ+;>7oUDxcdUqJGN21s=lcORPapz80_$7p0pv}s z?P&C@Yz_Xk(8tXJ118M?0{qMW--mmQmQ^P%a`#H&AqQh~9j}9sh+3@0nMf>Aqxo5Q zheE%q@#njpg-FFC)7nGC;Rc@J1`B2`MT2UHgSauathaoK#yUCU3s@~>hG*!4!**0I zCXAcAU~m$l4^85e;^@#c)05gQAZXw*j^dMS^DbHpVQI5`-5N+(>_M1%6vtg-mNrA( zBr712Z!y-OI1?A`7J)Vz39_JNL*P`r<|@I3y)*rws-3<NBfmIv%?>=VjgoU2<vtd| zrBT~Fn3E0b?r^hDVZaCS+oD&~ee?`jLAx+eV^VCKiwBa|m!;M5cGC;>uTW&m-gzYG zYI!gWBy#dZoNkf9#gXbhmp?{{oaV5ke^j$1raV-*s0p#3Rc*9$uC#FWMAJv6ud?eZ zK4*4z)f|$JMq-a33+^`-#F^FZgH&M|v)oRB+X>xAI`vv7%rg-WD2}72*V(K8)C>~B ze0Jm|=l<KISLhR|NGZVT5&`cb{M9gZtgZjp<*$1&H5TCV(82_qfL`OkEMYb*A;5TN z=?bG=e$(ZPbu^O9)61+S6b?%~oQbR*n^f_b4__B+_Po_`a5LX8>HbpUBxb`kC>St& z+u8Ue%Xu`dq$Q95sSaG@RF6y|zOLx?n87>-5vEUPV7@jkD=(mfXF!j_9A!w!p(u)j ziR;<Za%o@!yX<0TJR05Hf6?%rAd})-J(LkY9V@l!z1Bo$N(2qpiKMr8&Hek27O)6T zg*YX_+;=gRR-b%pV)SI_mQaWARN42``_JSoU?=N!^)x=FN%z<JHGbL+s0GH1RivU) zW!%(DIK=s|5Z>p>*9m7trW%BbMW?95%z0lA7o7FHw;@)9Kt|RN&R+@NTeC6iIepRk z*>|vWSQ|)Z4avr+u;(&+kvHe<`rVg$BXnno)KACPz%z>!;dmqUk*5!Eki<;SxQxQ? zcCqIJEZ#;j%_-9}koG~WRF;`<7pIy!u9{X6A15xI4m8K*W?cfhXL2=(r)qtIz*jn4 z_ICUlcRYF;8B`Qjua%ppo6#EIxqo<O@Zd_{O1Et-lacJN;}LxU_^khopusAF^nn6Q z+zasc?yn}U=U``VW%0+#|6=8%Qc}LO=x<J&f>$#4uvKyUaQy|fo190#H*lB?PWwHJ zs7dQS9y`j^HTANLq-cG;AzFZE^Y&LGsV>q%(u1teN`C85w6>swV}Wb!DQ*6xdYwtV z!FOT%IigwFqW;Q<0|9c&)cJZ@9+J&<+t%79NI_W0Vr#X_SV}lDJEM8VK0)0jPWSQi z?qPF6bc5|_Z!zkNR*Z?O$xcMzlwTw4UAGdgD2*~K&qMRU1J;gABSjNwhX`i8ii&`< zug!ekQmgx-2cI_eAQ_lx)Y@hK%uah{(!C?<I}TgS-EO5=7pG<Bx5JgZ-D`n1nyWq; zYC@x0{s9>aJIPi67XE2Fu^AwdHh|yk5rAdAiqF&k%WyX37fdnf5Io>Fao%v%jMhv? z9I^H|o)W3S<RZtkUP??4y9(N?5YAPo-Jhgg&6@z!hK|&9RC2QC7_!-njO<$&-q-H> z+B*K2FO36M-Lo~A(YwB)7JCEp5t$inUInX)>y0B98xN+45zIXLXGnAOD*7JNK9?w> z_lT-xAH_)gFN6f<@0sA^KSc*3(Y+tyFqz*GA{~4mkQr@-N7uPxRpeKxgkKatPJRrp z4ikn&Ju~f6F)qJ5oF3%#j8Y*p-n1Gri;sUC@=?THL91zy%Ti$pgBc@Of-((0zkAcM zyiT;ePPQJ=3XUWnTo8NVLd5^J9yqf5i^Kv<@~r{4;pB;Q2vp=j8jdA$2aoGsGfN8G zI3Yj+I)`JZQs0rGN%+muV5rOVo_^@_BS}1p&8N4VTw1#|>#;JoaB7*ohaqD^>_?oW zDI0LXj-zFK>1>!O85}a~1g75*g!Vjq9Zf&upl?is!hK4OIb&a-TjP;bRMSBx#}nbT z72|Ca6tZlhKqN#tN(-oQN`X@619E9EUOj9&Z3V}HzqJK6e22d(*(ZF1chis5*La1N zqs<G6BC0_8qmTv+ZxKq)1cneyrRvK0)AD5kPOdxN#-Rkgq*2;n%H0-YgXqIVj@!fm zk#C24G)$$75nY=~&QkG<c#kkY?sV)%UDSGT9F;O^&Bx16F&}Xv_`J&Z@8<3>=*p|V z-QcX62r28h^O~2OF6uL01DUWMcSK7)R(yfTW%2-HYba8<@N`5Blx9uyrS#VzNd1^I zYMU)cSRM@Y#-D5UR1o@hqrd9X%X+;(cF39SsB;tTx90i3-_FMHCpqvZIq)Yr@FzL& zCpqvZIq)Yr@FzL&ZzTtSL8t(Y*53dDcFqIg0cEU+du{_3Q$L{owV--^sPS6?_S0M1 zYstG<8`x?5%S)CP`#Zop%~k^t0Kf{svj&9A{&L^|>R|puO5~?n!Rvtxt?v^T0g(zB zfPMs^8R%aQS%A6@U^D-}6iIZgtjvCQ)}db79ZtYmuK~XN`@akY0#XA^0x<YLPW)Gh z-%bM7@X=xfFwcCzf&l$*I5V*SH3>Uo0|Wd26h!&C%kT9^kp}U|eOCpb(+>O<$|BNl zQGTvBUZbe^2c$v*K6)o$B7dI?QO|EtepVH|MnSHh_w@$Q+6AEX?<jghzeV|(1$vG0 z<<K@(7Qm#Y#A}rIW4}fDDdX`P#pLXyHUMB}h(xbZvKD`f@>AR5HOfBmNcte)Ix}Ou zHc`#3pD6am1{MZ1Kc9blc`EISTj&79QqTYa0j}y_4%LI70Cx5+<_31Zv$L&8HIM}W zlTUziBmZR}5K!&;e+)I&u{F?_vj+%I82uLEbxYe{AYG3?5&pXdx36h<-K_B!5Ix}j z1DO18fPS`ad=2!vC*Ch0WbpqE=x2w#*FdjZef<Iofcy#c_r_nZhrbrn`ZZjS?7y4G zPid{!K(Ez|egRQX{sj8JX&b!;e69BL3s9Zve*pZY{qq{|wHnSZz;NpS0q~bL&TGKe zYB9e6TWS6Wz+c)iuK{1LK7Rq8(*6W|C8hFOYUQW2%0EC~?Uw!ps?YWl^ly8o|1tWN zzRR!CFoOTilz#5$cn$ixNb?uyXNkW}>9=35|2<1tQvc4BewOpSp3>_K{jVvZD*tUt zzx}5Edp7FVSRf!{4Pc;u<fdMq{O>t~zYYi1{>$+H%pklT{r7O}Uq^2OGGxC_{$Ih_ z*8s0$a=##O4gU7ofBTF1-(QY4qkrde{0yK-i-EuS*J<$p4<x`>J%csCuK@Z#Sa7g( literal 0 HcmV?d00001 diff --git a/V1.7.1/run_analysis.sh b/V1.7.1/run_analysis.sh new file mode 100755 index 0000000..bb285c4 --- /dev/null +++ b/V1.7.1/run_analysis.sh @@ -0,0 +1,17 @@ +#!/bin/bash +if [[ $# < 2 ]] +then + echo "Usage:" + echo "sh this_file.sh input_path/input_file_name.xlsx output_path/output_file_name.xlsx optional_BioHPC_queue_name" + exit +fi +module load R/3.5.1-gccmkl +DIR=`dirname $0` +INPUT_FILE=$1 +OUTPUT_FILE=$2 +QUEUE=$3 +if [ "$QUEUE" == "" ] +then + QUEUE=256GB +fi +srun -p $QUEUE Rscript $DIR/ODA_1.7.1.R $INPUT_FILE $OUTPUT_FILE -- GitLab