{ "cells": [ { "cell_type": "markdown", "id": "ede3416b", "metadata": { "user_expressions": [] }, "source": [ "# Homework: Introductory GWAS Analysis\n", "# Due Date: Feb 16, 2026\n", "\n", "**Goal:** Run a reproducible GWAS analysis on simulated genotype data, evaluate basic QC, perform a per-SNP association test with covariates, visualize results, and write a short interpretation.\n", "\n", "**Learning objectives**\n", "- Understand genotype simulation and basic quality control (minor allele frequency, MAF, missingness).\n", "- Run PCA to capture population structure and include PCs as covariates.\n", "- Perform per-SNP linear regression GWAS with covariates.\n", "\n", "**Instructions**\n", "- Run each cell in order in JupyterLab. The notebook will create a `gwas_outputs/` folder with CSVs and PNGs.\n", "- Modify parameters (e.g., `n_samples`, `n_snps`, or number of causal SNPs) if you want to explore effects on power, but stick to the approximate time budget.\n", "\n", "**Grading rubric (100 pts)**\n", "- Notebook runs without errors and outputs saved: 40 pts\n", "- Correct GWAS implementation and inclusion of covariates: 25 pts\n", "- QC, PCA, and reasonable plotting: 15 pts\n", "- Short written interpretation, discussing results & limitations: 20 pts\n", "\n", "Good luck — ask if you get stuck!" ] }, { "cell_type": "code", "execution_count": 1, "id": "19842d85", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Environment ready. Outputs will be written to gwas_outputs/\n" ] } ], "source": [ "# Setup and imports\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from scipy import stats\n", "from sklearn.decomposition import PCA\n", "import os\n", "\n", "np.random.seed(2025)\n", "os.makedirs('gwas_outputs', exist_ok=True)\n", "print('Environment ready. Outputs will be written to gwas_outputs/')" ] }, { "cell_type": "markdown", "id": "f65dfff0-d7c8-4620-9d68-b9e8c6a4b86f", "metadata": { "user_expressions": [] }, "source": [ "## Part 1:\n", "Generate Synthetic SNP Genotype Data" ] }, { "cell_type": "code", "execution_count": 2, "id": "7eb60e71-4dfe-4b44-a479-d663bb3a2a05", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Environment ready. Outputs will be written to gwas_outputs/\n" ] } ], "source": [ "# Setup and imports\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from scipy import stats\n", "from sklearn.decomposition import PCA\n", "import os\n", "\n", "np.random.seed(2025)\n", "os.makedirs('gwas_outputs', exist_ok=True)\n", "print('Environment ready. Outputs will be written to gwas_outputs/')" ] }, { "cell_type": "code", "execution_count": 6, "id": "331504cb-7663-4f0e-88a9-95eb924f2b5a", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "array([[0., 2., 0., 0., 0.],\n", " [1., 0., 0., 0., 0.],\n", " [2., 0., 0., 0., 1.],\n", " [2., 1., 1., 1., 0.],\n", " [2., 0., 0., 0., 1.]])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Generate synthetic genotype data\n", "n_samples = 50\n", "n_snps = 500\n", "samples = [f\"S{i+1:02d}\" for i in range(n_samples)]\n", "pop_labels = np.array(['A'] * (n_samples//2) + ['B'] * (n_samples - n_samples//2))\n", "\n", "baseline_af = np.random.uniform(0.05, 0.5, size=n_snps)\n", "popB_shift = np.random.normal(loc=0.0, scale=0.08, size=n_snps)\n", "af_A = np.clip(baseline_af, 0.001, 0.999)\n", "af_B = np.clip(baseline_af + popB_shift, 0.001, 0.999)\n", "\n", "genotypes = np.full((n_samples, n_snps), np.nan)\n", "for v in range(n_snps):\n", " pA = af_A[v]\n", " pB = af_B[v]\n", " probsA = [ (1-pA)**2, 2*pA*(1-pA), pA**2 ]\n", " probsB = [ (1-pB)**2, 2*pB*(1-pB), pB**2 ]\n", " genotypes[:n_samples//2, v] = np.random.choice([0,1,2], size=n_samples//2, p=probsA)\n", " genotypes[n_samples//2:, v] = np.random.choice([0,1,2], size=n_samples - n_samples//2, p=probsB)\n", "genotypes[0:5, 0:5]" ] }, { "cell_type": "code", "execution_count": 7, "id": "e2e89805-af2a-4e1d-8580-9cc5d7c2345f", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Simulate covariates: age and sex (binary)\n", "age = np.random.normal(50, 10, size=n_samples).astype(int)\n", "sex = np.random.binomial(1, 0.5, size=n_samples)\n", "\n", "# Create phenotype: a polygenic trait with 10 causal SNPs\n", "causal_idx = np.random.choice(n_snps, size=10, replace=False)\n", "beta_snps = np.zeros(n_snps)\n", "beta_snps[causal_idx] = np.random.normal(0.3, 0.1, size=10)\n", "\n", "# Genetic effect\n", "# mean-impute genotypes for calculation\n", "G_imp = np.where(np.isnan(genotypes), np.nanmean(genotypes, axis=0), genotypes)\n", "polygenic = G_imp.dot(beta_snps)" ] }, { "cell_type": "code", "execution_count": 8, "id": "70cdb6c9-ecf1-4c5f-a2b5-bdd36256981a", "metadata": { "tags": [] }, "outputs": [], "source": [ "# Covariate effects\n", "beta_age = 0.02\n", "beta_sex = -0.1\n", "beta_pop = 0.4\n", "\n", "# Environmental noise\n", "env = np.random.normal(0, 1.0, size=n_samples)\n", "## convert population label to numeric \n", "pop_numeric = np.array([0 if x == 'A' else 1 for x in pop_labels])\n", "\n", "y = 0.5 + polygenic + beta_age * age + beta_sex * sex + beta_pop * pop_numeric + env" ] }, { "cell_type": "code", "execution_count": 10, "id": "dc535f3a", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Phenotypes saved: gwas_outputs/phenotypes.csv\n", "Generated genotype matrix shape: (50, 500)\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
rs100000rs100001rs100002rs100003rs100004rs100005rs100006rs100007rs100008rs100009...rs100490rs100491rs100492rs100493rs100494rs100495rs100496rs100497rs100498rs100499
S010.02.00.00.00.00.02.01.00.00.0...0.02.01.00.0NaN1.00.01.0NaN1.0
S021.00.00.00.00.01.00.00.00.00.0...0.01.01.00.00.00.01.00.00.00.0
S03NaN0.00.00.01.01.01.00.02.01.0...0.00.01.00.02.01.00.01.01.0NaN
S042.01.01.01.00.00.02.0NaN1.00.0...1.01.01.01.00.00.01.00.00.00.0
S052.00.00.00.01.01.00.01.00.00.0...0.01.01.00.01.01.01.01.01.00.0
\n", "

5 rows × 500 columns

\n", "
" ], "text/plain": [ " rs100000 rs100001 rs100002 rs100003 rs100004 rs100005 rs100006 \\\n", "S01 0.0 2.0 0.0 0.0 0.0 0.0 2.0 \n", "S02 1.0 0.0 0.0 0.0 0.0 1.0 0.0 \n", "S03 NaN 0.0 0.0 0.0 1.0 1.0 1.0 \n", "S04 2.0 1.0 1.0 1.0 0.0 0.0 2.0 \n", "S05 2.0 0.0 0.0 0.0 1.0 1.0 0.0 \n", "\n", " rs100007 rs100008 rs100009 ... rs100490 rs100491 rs100492 \\\n", "S01 1.0 0.0 0.0 ... 0.0 2.0 1.0 \n", "S02 0.0 0.0 0.0 ... 0.0 1.0 1.0 \n", "S03 0.0 2.0 1.0 ... 0.0 0.0 1.0 \n", "S04 NaN 1.0 0.0 ... 1.0 1.0 1.0 \n", "S05 1.0 0.0 0.0 ... 0.0 1.0 1.0 \n", "\n", " rs100493 rs100494 rs100495 rs100496 rs100497 rs100498 rs100499 \n", "S01 0.0 NaN 1.0 0.0 1.0 NaN 1.0 \n", "S02 0.0 0.0 0.0 1.0 0.0 0.0 0.0 \n", "S03 0.0 2.0 1.0 0.0 1.0 1.0 NaN \n", "S04 1.0 0.0 0.0 1.0 0.0 0.0 0.0 \n", "S05 0.0 1.0 1.0 1.0 1.0 1.0 0.0 \n", "\n", "[5 rows x 500 columns]" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_meta = pd.DataFrame({'sample': samples, 'age': age, 'sex': sex, 'pop': pop_labels, 'y': y}).set_index('sample')\n", "df_meta.to_csv('gwas_outputs/phenotypes.csv')\n", "print('Phenotypes saved: gwas_outputs/phenotypes.csv')\n", "\n", "# ~2% missing\n", "mask_missing = np.random.rand(n_samples, n_snps) < 0.02\n", "genotypes[mask_missing] = np.nan\n", "\n", "snps_ids = [f'rs{100000+i}' for i in range(n_snps)]\n", "df_geno = pd.DataFrame(genotypes, index=samples, columns=snps_ids)\n", "\n", "print('Generated genotype matrix shape:', df_geno.shape)\n", "df_geno.head()" ] }, { "cell_type": "markdown", "id": "98f904dd-e631-4cbe-a675-48180c1264da", "metadata": { "user_expressions": [] }, "source": [ "## Questions\n", "1. What does the variable `pop_labels' represents in the simulation code above? How does it affect the genotype data being generated?\n", "2. How is the phenotype (y) generated in the simulation above? What are the important factors contribute to the value of y in the simulation?\n", "3. What is the structure of the genotype matrix (df_genotype)? what does each row represent? How about columns? What do the values in the cells mean? \n", "4. How are the genotypes generated in the code above? Explain what is Hardy–Weinberg equilibrium and how does it related to how the genotypes being generated in the simulation? " ] }, { "cell_type": "code", "execution_count": 11, "id": "e030b933-4568-4857-a95d-74a63ec069ef", "metadata": { "tags": [] }, "outputs": [], "source": [ "import numpy as np\n", "\n", "def compute_maf(genotypes):\n", " \"\"\"\n", " Compute minor allele frequency (MAF) for a single SNP.\n", "\n", " Parameters\n", " ----------\n", " genotypes : array-like\n", " Genotypes coded as 0, 1, 2 (alt allele count).\n", " Missing values can be np.nan.\n", "\n", " Returns\n", " -------\n", " maf : float\n", " Minor allele frequency, or np.nan if no valid genotypes.\n", " \"\"\"\n", " g = np.asarray(genotypes, dtype=float)\n", " g = g[~np.isnan(g)] # remove missing\n", "\n", " if g.size == 0:\n", " return np.nan\n", "\n", " # total number of alternate alleles\n", " alt_allele_count = g.sum()\n", " # total number of alleles\n", " total_alleles = 2 * g.size\n", "\n", " af = alt_allele_count / total_alleles\n", " maf = min(af, 1 - af)\n", "\n", " return maf\n" ] }, { "cell_type": "code", "execution_count": 95, "id": "ddc40307-46aa-4760-b0bd-dd025be9967b", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "rs100000 0.061224\n", "rs100001 0.377551\n", "rs100002 0.340000\n", "rs100003 0.180000\n", "rs100004 0.153061\n", " ... \n", "rs100495 0.360000\n", "rs100496 0.071429\n", "rs100497 0.438776\n", "rs100498 0.010417\n", "rs100499 0.132653\n", "Length: 500, dtype: float64" ] }, "execution_count": 95, "metadata": {}, "output_type": "execute_result" } ], "source": [ "maf = df_geno.apply(compute_maf, axis=0)\n", "maf" ] }, { "cell_type": "markdown", "id": "ca600060-abfd-4eaa-8306-ec0f9a2c3251", "metadata": { "user_expressions": [] }, "source": [ "## Questions:\n", "5. What is an allele? Explain how allele frequencies are calculated. What is minor allele frequency (MAF)?\n", "6. Generate a plot for the MAF distribution. " ] }, { "cell_type": "code", "execution_count": 14, "id": "cfbc8f54-caa0-4f24-a2d7-efaac578f8ed", "metadata": { "tags": [] }, "outputs": [], "source": [ "# TASK: Compute QC metrics and run PCA\n", "# 1) Compute missingness for each SNP and save as gwas_outputs/variant_qc.csv\n", "# 2) Mean-impute genotypes and run PCA on a subset of SNPs to get PCs (save to gwas_outputs/pcs.csv)\n", "# 3) Plot PC1 vs PC2 colored by population and save as gwas_outputs/pca.png\n", "#\n", "# --- TODO: Computing missingness for each snp ---\n", "### YOUR CODE HERE ###\n", "\n", "# Save QC (once implemented)\n", "# qc = pd.DataFrame({'MAF': maf, 'missingness': missingness})\n", "# qc[1:10, 1:10]" ] }, { "cell_type": "markdown", "id": "4631ee8b-837e-4be0-b378-c6562c1f06cb", "metadata": { "user_expressions": [] }, "source": [ "### Questions\n", "7. Report the first5 rows and first 5 columns of the object `qc'. " ] }, { "cell_type": "code", "execution_count": 15, "id": "fb299bc1-558d-47dd-8f48-b7401103cf36", "metadata": { "tags": [] }, "outputs": [], "source": [ "### Impute Missing Values and Center genotype for PCA analysis\n", "G_imp = df_geno.fillna(df_geno.mean(axis=0))\n", "G_centered = G_imp - G_imp.mean(axis=0)" ] }, { "cell_type": "code", "execution_count": 16, "id": "b375944d-20ad-4c6a-b6e9-776b36b60379", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
PC1PC2PC3PC4
S01-1.888524-3.0093120.7076430.659718
S02-2.3751504.459011-3.7931631.617335
S03-5.630537-2.389041-3.736545-2.093097
S04-1.9665140.735002-4.7221710.654252
S05-2.1903210.6364153.232525-0.872838
\n", "
" ], "text/plain": [ " PC1 PC2 PC3 PC4\n", "S01 -1.888524 -3.009312 0.707643 0.659718\n", "S02 -2.375150 4.459011 -3.793163 1.617335\n", "S03 -5.630537 -2.389041 -3.736545 -2.093097\n", "S04 -1.966514 0.735002 -4.722171 0.654252\n", "S05 -2.190321 0.636415 3.232525 -0.872838" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# PCA on mean-imputed genotypes (for population structure)\n", "# perform PCA on a random subset of SNPs for speed\n", "pca = PCA(n_components=4)\n", "pcs = pca.fit_transform(G_centered.values)\n", "\n", "pc_df = pd.DataFrame(pcs, index=samples, columns=['PC1','PC2','PC3','PC4'])\n", "pc_df.head()" ] }, { "cell_type": "code", "execution_count": 17, "id": "5eb4be4b-7ba7-4aee-bf6d-c6d5d15c77a6", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
agesexpopy
sample
S01450A1.592461
S02451A3.078932
S03481A0.203718
S04451A2.332965
S05540A2.585327
\n", "
" ], "text/plain": [ " age sex pop y\n", "sample \n", "S01 45 0 A 1.592461\n", "S02 45 1 A 3.078932\n", "S03 48 1 A 0.203718\n", "S04 45 1 A 2.332965\n", "S05 54 0 A 2.585327" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_meta.head()" ] }, { "cell_type": "code", "execution_count": 18, "id": "95df38e0-9ad4-45c1-8f32-59be847cd17b", "metadata": { "tags": [] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAHqCAYAAADyPMGQAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAS0VJREFUeJzt3Xt01OW97/HPzGRuQDIh3EKAcBNv3KJREBWFg1pYRVEQpaIFb1srsLVslbrpObpdsDlWObtbLLQqAtZrOSCxsXqk1gsWKIgNUKoWCBQWSQgYyASSuWTyO3/8SCQk4CSZzPX9WitrzO+ZzHwzePn4PM/v+1gMwzAEAACA72WNdQEAAACJguAEAAAQJoITAABAmAhOAAAAYSI4AQAAhIngBAAAECaCEwAAQJgITgAAAGEiOAEAAISJ4ASkkJUrV8pisTR8paWlqXfv3rr77rt16NChJs8vLi7W7Nmzdf7558vtdqtDhw4aPHiwfv7znzf7fEmaPHmyLBaLZs+e3eZ6n3rqqUb1OhwO9e/fXw8//LCOHz/e5Pk7duzQ3Xffrf79+8vlcqlTp0669NJL9Ytf/EIVFRUNz/v888913333KT8/X06nUxaLRfv3729zvZFy+u9ssVjk8Xg0ZswYvffee02e6/f79cILL+jqq69W586d5XA41KtXL91222369NNPG563bds2zZo1S0OHDlV6erp69Oih6667Tn/605+i+asBCY/gBKSgFStWaNOmTVq/fr3uv/9+vfnmmxo9erROnjzZ8JzCwkINGzZMhYWF+pd/+RcVFhY2/PXvf/97TZw4scnrlpeXq7CwUJL0+uuvy+fzRaTeDz74QJs2bdJ7772nm2++WUuWLNGECRN0+olRL730kvLz87V161Y99thj+uCDD/TOO+9o6tSp+vWvf61777234bkfffSR/vjHPyo3N1dXXnllRGqMtFtvvVWbNm3Sn//8Z/3qV79SWVmZbrzxxkbh6ejRo7rqqqs0d+5cDRkyRCtXrtRHH32kxYsXy2azady4cdq+fbsk6c0339SWLVt0zz33qKCgQC+//LKcTqfGjRunV199NVa/JpB4DAApY8WKFYYkY+vWrY2u/8//+T8NScZrr71mGIZhFBcXGx07djQuueQS4/jx401ep66uzlizZk2T688++6whyfjhD39oSDJef/31NtX75JNPGpKMI0eONLp+1113GZKMzz//3DAMw9i4caNhs9mM8ePHGz6fr8nr+P1+o6CgoOH7UCjUpOZ9+/a1qdZIkmTMmjWr0bU9e/YYkozrrruu4dqECROMtLQ046OPPmr2dbZs2WL885//NAzDMA4fPtxkvLa21hg2bJgxcODACFYPJDdmnADoiiuukCT985//lCT9n//zf3Ty5EktXbpUHo+nyfMtFosmT57c5Porr7yiHj16aNWqVXK73XrllVeiUu9//ud/ymKx6MUXX5TT6WzyfIfDoZtuuqnhe6u19f/qu+SSSzR69Ogm10OhkHr16tXoc1m2bJmGDx+uTp06KT09XRdeeKH+/d//vVXvO3DgQHXr1q3hd962bZvef/993Xvvvfof/+N/NPszl19+uXJzcyVJ3bt3bzJus9mUn5+vgwcPtqomIBURnABoz549kqRu3bpJkj788EP16NGjIaCEY+PGjfrqq6/04x//WF26dNGUKVP0pz/9Sfv27Wv0vP3798tisWjmzJkRqTcUCulPf/qT8vPz1adPn1a/Zrjuvvtuff7559q9e3ej6x9++KFKSkp09913S5LeeustPfTQQ7r22mv1zjvvaN26dfrpT3/aaDm0JY4dO6Zvv/220Z+RJN18882t/l1qa2u1YcMGDR48uNWvAaQaghOQgkKhkGpra3XixAm99957WrBggdLT0xtmZQ4cOKD+/fu36DWXL18uSbrnnnskSffee68Mw9CKFSsaPc9ischms8lms7W43uPHj+v111/Xr3/9a/Xp00ejR4/W0aNHVV1d3eJ6W2v69OlyOBxauXJlo+srV65Ujx49NGHCBEnSn//8Z2VmZur555/X9ddfr3HjxumBBx7Qf//3f4f1PoZhqLa2VsFgUF9//bWmT5+uuro6TZ8+XZL5ZySpTb/3U089pT179ujJJ59s9WsAqYbgBKSgK664Qna7Xenp6Zo4caKys7P1/vvvq0ePHq16vRMnTuh3v/udrrzySl144YWSpGuvvVYDBw7UypUrVVdX1/Dcvn37qra2tiFohSM7O1t2u12dO3fWnXfeqUsvvVQffPCBXC5Xq+ptiy5duujGG2/UqlWrGn6vY8eOqaCgQD/+8Y+VlpYmSRoxYoSOHz+uH/3oRyooKNDRo0db9D5Lly6V3W6Xw+HQRRddpI0bN+rpp5/WQw89FJHf4+WXX9bChQv1b//2b5o0aVJEXhNIBWmxLgBA9L366qu66KKLlJaWph49eqhnz56NxnNzc5sssZ3L22+/rRMnTui2225r1Cbgtttu06JFi7R+/Xr94Ac/aHW9f/zjH+XxeGS329W7d2916dKlYaxr167q0KFDi+ptq3vuuUdr1qxp+L3efPNN+f3+RsuPd911l2pra/XSSy9pypQpqqur0+WXX64FCxbo+uuv/973uO222/TYY4/JYrEoPT1dAwcObDRLV793ad++fbrgggtaVP+KFSv0wAMP6F/+5V/07LPPtuhngZQX693pAKLnbHfVnWnOnDmGJGPTpk1hve6oUaMMSWf9mjp1aqvqPdtddWe68cYbjbS0NOPgwYMtfo/W3FVXW1tr5OTkGLfffrthGIZx2WWXGSNHjjzr80+cOGH84Q9/MC6//HLD4XAY+/fvP+frq5m76s60bds2Q5LxwAMPhF23YRjGK6+8YlitVuPuu+826urqWvSzALirDkAzfvrTn6pjx4566KGHVFlZ2WTcMAy98847kqSvvvpKmzZt0pQpU/Txxx83+Ro3bpwKCgr07bfftlu9TzzxhAzD0P33369AINBkPBgM6ve//33E3s9ms+muu+7SunXrtGHDBn3xxRcNe7ua07FjR02YMEHz589XIBDQrl272lzDpZdeqgkTJmj58uVnbWL5xRdfNOyFksx9WPfdd5/uvPNOvfzyy7JYLG2uA0g1LNUBaKJ///566623dPvttysvL0+zZ8/WJZdcIkn6+9//rldeeUWGYeiWW25p2Kv0+OOPa8SIEU1eq6qqSh999JFee+01Pfzww/rnP/+pgQMHasaMGS3a53Quo0aN0rJly/TQQw8pPz9fP/nJTzR48GAFg0H99a9/1YsvvqghQ4boxhtvlCQdOXKkoav2zp07JUnvv/++unXrpm7duunaa6/93ve855579Mwzz+iOO+6Q2+3W7bff3mj8/vvvl9vt1lVXXaWePXuqrKxMixYtksfj0eWXXx6R3/vVV1/V+PHjNWHCBN1zzz2aMGGCOnfurNLSUv3+97/Xm2++qW3btik3N1erV6/Wvffeq7y8PD3wwAPasmVLo9e65JJLmm3lAOAMMZ7xAhBF4S7V1du7d6/x0EMPGeedd57hdDoNt9ttXHzxxcbcuXONffv2GYFAwOjevbuRl5d31teora01evfubQwdOtQwDMPYt2+fIcmYMWPG975/uEt19YqKiowZM2YYubm5hsPhaGji+b/+1/8yysvLG5738ccfn3VZ8dprrw3rvQzDMK688kpDkjF9+vQmY6tWrTLGjh1r9OjRw3A4HEZOTo5x2223GTt27Pje11UYS3X1ampqjOeff94YNWqUkZGRYaSlpRk5OTnG5MmTjffee6/heTNmzDjncmo8NQAF4pnFME47swAAAABnxR4nAACAMBGcAAAAwkRwAgAACBPBCQAAIEwEJwAAgDARnAAAAMKUUg0w6+rqVFJSovT0dDrmAgAASeZpCFVVVcrJyZHVeu45pZQKTiUlJerTp0+sywAAAHHo4MGD6t279zmfk1LBKT09XZL5wWRkZMS4GgAAEA+8Xq/69OnTkBPOJaWCU/3yXEZGBsEJAAA0Es42HjaHAwAAhIngBAAAECaCEwAAQJgITgAAAGEiOAEAAISJ4AQAABAmghMAAECYCE4AAABhIjgBAACEieAEAAAQJoITAABAmAhOAAAAYUqpQ34BIGEFqqWdq6VdayVviZSRIw2eLA2dKjk6xLo6IGUQnAAg3gWqpXfnSPs3SLJIaS6p/GupfKG07zPppiWEJyBKWKoDgHi3c7UZmtydpfRsyZ1pPro6m9d3ro51hUDKIDgBQLzbtVYNM02ns7ski/XUOIBoIDgBQLzzljQNTfVsTslbGt16gBRGcAKAeJeRI9X6mh8L+aWMntGtB0hhBCcAiHeDJ0sypOAZ4Snok4y6U+MAooG76gAg3g2dat49t3+D5Leay3Mhvxma+o02xwFEBcEJAOKdo4PZcqChj1OplNWPPk5ADBCcACARODpI+TPMLwAxwx4nAACAMBGcAAAAwkRwAgAACBPBCQAAIEwEJwAAgDARnAAAAMJEcAIAAAgTwQkAACBMBCcAAIAwEZwAAADCRHACAAAIE8EJAAAgTAQnAACAMBGcAAAAwkRwAgAACFNarAsAACSmmkBIBUWHVLijVGWVNcr2uDVxWE9Nyuslt8MW6/KAdkFwAgC0WE0gpHlrtmtTcYWskpx2m3YfrtLi9VXauPeonpkynPCEpMRSHQCgxQqKDmlTcYUy3XZ1z3DJc+ox023X5uIKFRQdinWJQLsgOAEAWqxwR6msklz2xrNKLrtNFos5DiQjghMAoMXKKmvktDe/FOdMs6nMWxPlioDoIDgBAFos2+OWPxhqdsxfG1J2hjvKFQHRQXACALTYxGE9VSfJd0Z48gVDMgxzHEhG3FUHAGixSXm9tHHvUW0urpDXF5QzzSZ/rRmarhiQpUl5vWJdItAuCE4AgBZzO2x6Zsrw7/o4eWuUm5VOHyckPYthGEasi4gWr9crj8ejyspKZWRkxLocAAAQB1qSD9jjBAAAECaCEwAAQJgITgAAAGEiOAEAAISJu+qAOMOJ8wAQvwhOQBzhxHkAiG8s1QFxhBPnASC+EZyAOMKJ8wAQ3whOQBzhxHkAiG8EJyCOcOI8AMQ3ghMQRzhxHgDiG3fVAXGEE+cBIL4RnIA4wonzABDfLIZhGLEuIlpacvoxAABIDS3JB8w4AQCAyAtUSztXS7vWSt4SKSNHGjxZGjpVcnSIdXWtRnACAACRFaiW3p0j7d8gySKluaTyr6XyhdK+z6SbliRseOKuOgAAEFk7V5uhyd1ZSs+W3Jnmo6uzeX3n6lhX2GoEJwAAEFm71qphpul0dpdksZ4aT0wEJwAAEFnekqahqZ7NKXkT9/goghMAAIisjByp1tf8WMgvZSRuM1+CEwAAiKzBkyUZUvCM8BT0SUbdqfHElLDBadGiRbJYLHrkkUdiXQoAADjd0KlSv9GS75h04rBUc9x89B0zrw+dGusKWy0h2xFs3bpVL774ooYNGxbrUgAAwJkcHcyWAw19nEqlrH70cYqFEydOaPr06XrppZe0YMGCWJcDAACa4+gg5c8wv5JIwi3VzZo1Sz/84Q913XXXxboUAACQYhJqxumtt97Sl19+qa1bt4b1fL/fL7/f3/C91+ttr9IAAEAKSJgZp4MHD+rhhx/Wa6+9JpfrLL0hzrBo0SJ5PJ6Grz59+rRzlQAAIJlZDMMwYl1EONatW6dbbrlFNput4VooFJLFYpHVapXf7280JjU/49SnT5+wTj8GAACpwev1yuPxhJUPEmapbty4cdq5c2eja3fffbcuvPBCzZs3r0lokiSn0ymn0xmtEgEAQJJLmOCUnp6uIUOGNLrWsWNHdenSpcl1AACA9pAwe5wAAABiLWFmnJrzySefxLoEAACQQphxAgAACBPBCQAAIEwEJwAAgDAl9B4nJKeaQEgFRYdUuKNUZZU1yva4NXFYT03K6yW3o2nbCQAAooXghLhSEwhp3prt2lRcIaskp92m3YertHh9lTbuPapnpgwnPAEAYoalOsSVgqJD2lRcoUy3Xd0zXPKcesx027W5uEIFRYdiXSIAIIURnBBXCneUyirJZW88q+Sy22SxmOMAAMQKwQlxpayyRk5780txzjSbyrw1Ua4IAIDvEJwQV7I9bvmDoWbH/LUhZWe4o1wRAADfITghrkwc1lN1knxnhCdfMCTDMMcBAIgV7qpDXJmU10sb9x7V5uIKeX1BOdNs8teaoemKAVmalNcr1iUCAFIYwQlxxe2w6Zkpw7/r4+StUW5WOn2cgPYUqJZ2rpZ2rZW8JVJGjjR4sjR0quToEOvqgLhiMQzDiHUR0eL1euXxeFRZWamMjIxYlwMAsReolt6dI+3fIMkipbmkWp8kQ+o3WrppCeEJSa8l+YA9TgCQynauNkOTu7OUni25M81HV2fz+s7Vsa4QiCsEJwBIZbvWqmGm6XR2l2SxnhoHUI/gBACpzFvSNDTVszklL01ngdMRnAAglWXknNrT1IyQX8qgBQhwOoITAKSywZMlGVLwjPAU9ElG3alxAPVoRwAAqWzoVGnfZ+ZGcL/VXJ4L+c3Q1G+0OQ6gAcEJAFKZo4PZcqChj1OplNWPPk7AWRCcACDVOTpI+TPMLwDnRHBCSqgJhL7rRl5Zo2yPm27kAIAWIzgh6dUEQpq3Zrs2FVfIKslpt2n34SotXl+ljXuP6pkpwwlPAICwcFcdkl5B0SFtKq5Qptuu7hkueU49Zrrt2lxcoYKiQ7EuEQCQIAhOSHqFO0plleSyN55VctltsljMcQAAwkFwQtIrq6yR0978UpwzzaYyb02UKwIAJCqCE5JetsctfzDU7Ji/NqTsDHeUKwIAJCqCE5LexGE9VSfJd0Z48gVDMgxzHACAcHBXHZLepLxe2rj3qDYXV8jrC8qZZpO/1gxNVwzI0qS8XrEuEQCQIAhOSHpuh03PTBn+XR8nb41ys9Lp4wQAaDGLYRhGrIuIFq/XK4/Ho8rKSmVkZMS6HAAAEAdakg/Y4wQAABAmghMAAECYCE4AAABhIjgBAACEieAEAAAQJoITAABAmAhOAAAAYSI4AQAAhInO4QASQk0g9F3398oaZXvcdH9PZYFqaedqaddayVsiZeRIgydLQ6dKjg6xrg5JjM7hAOJeTSCkeWu2a1NxhaySnHab/MGQ6iSNGpClZ6YMJzylkkC19O4caf8GSRYpzSXV+iQZUr/R0k1LCE9oETqHA0gqBUWHtKm4Qpluu7pnuOQ59ZjptmtzcYUKig7FukRE087VZmhyd5bSsyV3pvno6mxe37k61hUiiRGcAMS9wh2lskpy2RvPKrnsNlks5jhSyK61aphpOp3dJVmsp8aB9kFwAhD3yipr5LQ3vxTnTLOpzFsT5YoQU96SpqGpns0peQnSaD8EJwBxL9vjlj8YanbMXxtSdoY7yhUhpjJyTu1pakbIL2X0jG49SCkEJwBxb+KwnqqT5DsjPPmCIRmGOY4UMniyJEMKnhGegj7JqDs1DrQP2hEAiHuT8npp496j2lxcIa8vKGeaTf5aMzRdMSBLk/J6xbpERNPQqdK+z8yN4H6ruTwX8puhqd9ocxxoJ7QjAJAQGvVx8tYoO4M+TimtUR+nUnN5jj5OaKWW5AOCEwAASGn0cQIAAGgHBCcAAIAwEZwAAADCRHACAAAIE+0IgCTX6G60yhple7gbDQBai+AEJLGaQEjz1mzXpuIKWSU57TbtPlylxeurtHHvUT0zZTjhCQBagKU6IIkVFB3SpuIKZbrt6p7hkufUY6bbrs3FFSooOhTrEgEgoTDjBCSxwh2lskpynXFArstuk9cXVOGOUk0bkRub4oC2atQEs8Q8w44mmGhnBCcgiZVV1shpb34pzplmU5m3JsoVARESqJbenWMeuyKLlOaSyr+Wyheax7HctITwhHbBUh2QxLI9bvnPOBi3nr82pOwMd5QrAiJk52ozNLk7S+nZkjvTfHR1Nq/vXB3rCpGkCE5AEps4rKfqJPnOCE++oHlA7sRhPWNTGNBWu9aqYabpdHaXZLGeGgcij6U6IIlNyuuljXuPanNxhby+oJxpNvlrzdB0xYAsTcrrFesSgdbxljQNTfVsTvPgX6AdEJyAJOZ22PTMlOHf9XHy1ig3K50+Tkh8GTnmnqbmhPxSVr+oloPUQXACkpzbYdO0EbncPYfkMniyuRE86DOX5+oFfZJRZ44D7YDgBABIPEOnmnfP7d8g+a3m8lzIb4amfqPNcaAdEJwAAInH0cFsOdDQx6nUXJ6jjxPaGcEJAJCYHB2k/BnmFxAlCdOOYNGiRbr88suVnp6u7t276+abb9Y333wT67IAAEAKSZjg9Omnn2rWrFnavHmz1q9fr9raWt1www06efJkrEsDAAApwmIYhhHrIlrjyJEj6t69uz799FNdc801Yf2M1+uVx+NRZWWlMjIy2rlCAACQCFqSDxJmxulMlZWVkqSsrKwYVwIAAFJFQm4ONwxDc+fO1dVXX60hQ4ac9Xl+v19+v7/he6/XG43yAABAkkrIGafZs2drx44devPNN8/5vEWLFsnj8TR89enTJ0oVAgCAZJRwe5zmzJmjdevW6bPPPlP//v3P+dzmZpz69OnDHicAANCgJXucEmapzjAMzZkzR++8844++eST7w1NkuR0OuV0OqNQHQAASAUJE5xmzZqlN954QwUFBUpPT1dZWZkkyePxyO12x7g6AACQChJmqc5isTR7fcWKFZo5c2ZYr0E7AgAAcKakXaoDAACIpYS8qw4AACAWCE4AAABhIjgBAACEKWH2OAEAEkygWtq5Wtq1VvKWSBk50uDJ0tCpkqNDrKsDWoXgBACIvEC19O4caf8GSRYpzSWVfy2VL5T2fSbdtITwhITEUh0AIPJ2rjZDk7uzlJ4tuTPNR1dn8/rO1bGuEGgVghMAIPJ2rVXDTNPp7C7JYj01DiQeghMAIPK8JU1DUz2bU/KWRrceIEIITgCAyMvIkWp9zY+F/FJGz+jWA0QIwQkAEHmDJ0sypOAZ4Snok4y6U+NA4uGuOiS0mkBIBUWHVLijVGWVNcr2uDVxWE9Nyuslt8MW6/KA1DV0qnn33P4Nkt9qLs+F/GZo6jfaHAcSUMIc8hsJHPKbXGoCIc1bs12biitkleS02+QPhlQnadSALD0zZTjhCYilRn2cSs3lOfo4IQ4l5SG/wJkKig5pU3GFMt12ueynApLbLl8wpM3FFSooOqRpI3JjWySQyhwdpPwZ5heQJNjjhIRVuKNUVum70HSKy26TxWKOAwAQSQQnJKyyyho57c0vxTnTbCrz1kS5IgBAsiM4IWFle9zyB0PNjvlrQ8rOcEe5IgBAsiM4IWFNHNZTdZJ8Z4QnXzAkwzDHAQCIJDaHI2FNyuuljXuPanNxhby+oJxpNvlrzdB0xYAsTcrrFesSAaSyRncVlphNQbmrMOHRjgAJrVEfJ2+NsjNi18eJnlIAGgSqpXfnmH2s6s/sq/VJMsw+VjctITzFkZbkA4ITEAH0lALQyLZV0scLJXfnxmf2BX2S75g0dj5tGuJIS/IBe5yACDi9p1T3DJc8px4z3faGnlIAUsiutWqYaTqd3SVZrKfGkYjY4wREwLl6Snl9QRXuKKUZJxAPorXvyFvSNDTVsznNTupISAQnIALoKQUkgOb2HZV/LZUvNM/Vi+S+o4wc87WbE/JLWf0i8z6IOpbqgAigpxSQAHauNkOTu7OUni25M81HV2fz+s7VkXuvwZMlGeaeptMFfeZBx4MnR+69EFUEJyAC6CmFdhOoNjcavzpJeuFy83HbKvM6Wiaa+46GTjXvnvMdk04clmqOm4++Y+b1oVMj916IKpbqgAigpxTaRTSXllJBNPcdOTqYfz4N+6lKzeU5+jglPIITEAFuh03PTBneqKdUblY6fZzQNqcvLZ15S3v90hK3tIcv2vuOHB3MPx/+jJIKwQmIELfDpmkjcrl7DpFzrqUl/6mlJf6jHL7Bk83ZuqDP/Azrse8ILUBwAoB4xS3tkTV0qrnEuX+DGTxtTnOmyahj3xHCRnACgHjFLe2Rxb4jRADBCQDiFUtLkce+I7QRwQkA4hVLS0DcITgBQLxiaQmIOwQnAIhnLC0BcYXO4QAAAGFqcXDavn27FixYoKVLl+ro0aONxrxer+65556IFQcAABBPLIZhGOE++cMPP9SNN96oQYMGqaqqStXV1frd736nsWPHSpIOHz6snJwchULNH3Yaa16vVx6PR5WVlcrIyIh1OQAAIA60JB+0aMbpqaee0qOPPqq//e1v2r9/vx5//HHddNNN+uCDD9pUMAAAQCJo0ebwXbt26be//a0kyWKx6LHHHlPv3r1166236s0339SIESPapUgAAIB40KLg5HQ6dfz48UbXfvSjH8lqtWratGlavHhxJGsD0EY1gdB3Bw9X1ijb4+bgYQBogxYFp7y8PH388cfKz89vdP32229XXV2dZszgdlkgXtQEQpq3Zrs2FVfIKslpt2n34SotXl+ljXuP6pkpwwlPANBCLdrj9JOf/ESHDh1qduxHP/qRVq1apWuuuSYihQFom4KiQ9pUXKFMt13dM1zynHrMdNu1ubhCBUXN/7MMADi7Ft1Vl+i4qw7JIpwluDtf/ot2H65S9wxXk58vr/JpUPd0vXbfyGiXDgBxp93uqjt27JiWLFkir9fbZKyysvKsYwAip34JbvH6f2j34SoFQsapJbh/aN6a7aoJmO1Ayipr5LQ3vxTnTLOpzFsTzbIBICm0KDi98MIL+uyzz5pNYx6PRxs2bNCSJUsiVhyApsJdgsv2uOUPNt9TzV8bUnaGO5plA4AUqJa2rZJenSS9cLn5uG2VeT1BtCg4rVmzRg8++OBZxx944AH93//7f9tcFICzK9xRKqsk1xmzSS67TRaLOS5JE4f1VJ0k3xnhyRcMyTDMcQCImkC19O4c6eOFUvnXUm3AfPx4oXk9QcJTi+6q27t3rwYNGnTW8UGDBmnv3r1tLgrA2YW7BDcpr5c27j2qzcUV8vqCcqbZ5K81Q9MVA7I0Ka9XNMsGkEgC1dLO1dKutZK3RMrIkQZPloZONQ+ebo2dq6X9GyR3ZynttL2XQZ95fefqhDjMukUzTjabTSUlJWcdLykpkdXKucFAewp3Cc7tsOmZKcM19/rzNah7uhxpFg3qnq65159PKwIAZ9deM0O71kqyNA5NkmR3SRbrqfH416IZp0suuUTr1q3TFVdc0ez4O++8o0suuSQihQFo3sRhPbV4fZV8wVCj5brmluDcDpumjcjVtBG5sSgVQCJqr5khb0nT0FTP5pS8pa2rN8paND00e/ZsLV68WC+88EKjg3xDoZCWLFmi//qv/9KsWbMiXiSA70zK66VRA7JUWRNUeZWv0SNLcADarL1mhjJypFpf82Mhv5SRGPsuWzTjNGXKFD3++OP613/9V82fP18DBgyQxWLR3r17deLECT322GO69dZb26tWAPpuCa6hj5O3RrlZ6RylAiAy2mtmaPBkqXyhOXNlP2Mmy6gzxxNAqxpgbt26Va+//rp2794twzB0/vnn64477oj7Q35pgAkAwPd4dZK5pyk9u+nYicNStwukHxe0/HXr907t32DOXNmc5kyTUSf1Gy3dtKT1G8/bqCX5oEUzTtXV1Xrssce0bt06BYNBjRs3TkuWLFHXrl3bVDCA5MMBw0CCaq+ZIUcHMxw13K1XKmX1a/vdelHWohmnxx57TEuXLtX06dPldrv1xhtvaMyYMVq9enV71hgxzDgB0dHcAcP+YEh1kkYNyOKuPiCexfHMUHtptxmntWvXavny5Zo2bZokafr06brqqqsUCoVks/EvQQCm07ubN9z557bLFww1dDfnTj8gTiXJzFB7aVFwOnjwoEaPHt3w/YgRI5SWlqaSkhL16dMn4sUB0cTSUuScq7u51xdU4Y5SghMQzxwdzJYDCdCQMtpaFJxCoZAcDkfjF0hLU21tbUSLAqKtuaUl8+DcKm3ce5SlpRbigGEAyapFwckwDM2cOVNOp7Phms/n04MPPqiOHTs2XFu7NjG6fwL1WFqKrGyPW7sPV0lue5Mxf21IuVnpMagKANquRcFpxoymU3Z33nlnxIoBYoWlpchqSXdzAEgkLQpOK1asaK86gJhiaSmyOGAYQLJqUXACkhVLS5FFd3MAyYrgBIilpfbAAcMAkhHBCRBLSwCA8BCcALG0BAAIT6sO+U1UHLkCAADO1JJ8YI1STQAAAAkv4YLT0qVL1b9/f7lcLuXn52vDhg2xLgkAAKSIhNrj9Pbbb+uRRx7R0qVLddVVV+k3v/mNJkyYoL///e/KzeXOnWTFGXIAgHiRUHucRo4cqUsvvVTLli1ruHbRRRfp5ptv1qJFi77359njlHiaO0POHwypTtKoAVmcIQcAaLOk3OMUCAS0bds23XDDDY2u33DDDdq4cWOMqkJ7O/0Mue4ZLnlOPWa67Q1nyAEAEC0JE5yOHj2qUCikHj16NLreo0cPlZWVNfszfr9fXq+30RcSy7nOkLNYzHEAAKIlofY4SZLFYmn0vWEYTa7VW7Rokf7jP/4jGmWhnXCGHAAkkEC1tHO1tGut5C2RMnKkwZOloVMlR4dYVxcRCTPj1LVrV9lstiazS+Xl5U1moeo98cQTqqysbPg6ePBgNEpFBGV73PIHQ82O+WtDys5wR7kiAECzAtXSu3OkjxdK5V9LtQHz8eOF5vVAdawrjIiECU4Oh0P5+flav359o+vr16/XlVde2ezPOJ1OZWRkNPpCYpk4rKfqZJ4ZdzrOkAOAOLNztbR/g+TuLKVnS+5M89HV2by+c3WsK4yIhAlOkjR37ly9/PLLeuWVV/TVV1/ppz/9qQ4cOKAHH3ww1qWhnUzK66VRA7JUWRNUeZWv0SNnyAFAHNm1VpJFSnM1vm53SRbrqfHEl1B7nG6//XZ9++23evrpp1VaWqohQ4boD3/4g/r27Rvr0tBOOEMOABKEt6RpaKpnc0re5LiZJ6H6OLUVfZwAAGgnr04y9zSlZzcdO3FY6naB9OOC6NcVhqTs4wQAAOLY4MmSDCnoa3w96JOMulPjiS+hluoAAECcGjpV2veZuRHcbzWX50J+MzT1G22OJwGCEwBEGecvIik5Okg3LTmtj1OplNUv6fo4sccJAKKI8xeB+NOSfMCMEwBE0ennLzYcJeS2yxcMNZy/OG1EbmyLRPtJgc7ayY7N4QAQRZy/mMJSpLN2siM4AUAUcf5iCkuRztrJjuAEAFHE+YspLEU6ayc7ghMARBHnL6awFOmsnewITgAQRZy/mMIycqRaX/NjIb+UQWhOBAQnAIii+vMX515/vgZ1T5cjzaJB3dM19/rzaUWQ7FKks3ayox0BAESZ22HTtBG5tB1INSnSWTvZEZwAAOGhB1HbpEhn7WRH53AAwPer70G0f4Ma7gyr9UkyzNmSm5bwH34krJbkA/Y4AQC+Hz2IAEkEJwBAOOhBBEgiOAEAwkEPIkASwQkAEA56EAGSCE4AgHDQgwiQRDsCAEA46EEESCI4AQDCQQ8iQBLBCVBNIKSCokMq3FGqssoaZXvcmjispybl9eL4C+B0jg5S/gzzC0hRBCektJpASPPWbNem4gpZJTntNu0+XKXF66u0ce9Rzg4DADTC5nCktIKiQ9pUXKFMt13dM1zynHrMdNu1ubhCBUWHYl0iACCOEJyQ0gp3lMoqyWVvPKvksttksZjjAADUY6kOKa2sskb2NKsqTgZ0vDqg2jpDaVaLMjs4ZLdZVeatiXWJAIA4wowTUlq3dJeOVPl12OuTv7ZOhiH5a+t02OvTkSq/unU6S6dkAEBKYsYJKa2nx6VgqE42q0VpVvP/I2yyqLauTsFQnXp6CE4AgO8QnJDSSit9stusCtUZCobqZLVYVGcYkiSHzarSyrMcMQEASEkEJ6S0I1U+dUt3yjCkypqggqE6udJs8rjtslikIycITkgAgerTGlOWmOfK0ZgSaBcEJ6S0bI9buw9XqXuGS1kdHY3Gyqt86telY4wqA8IUqJbenWMehSKLlOaSyr+WyheaR6TctITwBEQQm8OR0iYO66k6Sb5gqNF1XzAkwzDHgbi2c7UZmtydpfRsyZ1pPro6m9d3ro51hUBSITghpU3K66VRA7JUWRNUeZWv0eMVA7I0Ka9XrEsEzm3XWjXMNJ3O7pIs1lPjACKFpTqkNLfDpmemDP/urDpvjXKz0jmrDonDW9I0NNWzOc3DeAFEDMEJKc/tsGnaiFxNG5Eb61LajAOLU1BGjrmnqTkhv5TVL6rlAMmOpTogSdQfWLx4/T+0+3CVAiHj1IHF/9C8NdtVEwh9/4sg8QyeLMmQgmfcARr0SUbdqXEAkUJwApIEBxanqKFTpX6jJd8x6cRhqea4+eg7Zl4fOjXWFQJJheAEJAkOLE5Rjg5my4Gx86VuF0hpTvNx7HxaEQDtgD1OQJIoq6yR0978PiZnmo0Di5OZo4OUP8P8AtCuCE5AHGvJZu/6Zp5y25u8jr82pNys9GiVDQBJi6U6IE61dLM3zTwBoP0x44SEley33p++2bth35LbLl8w1LDZ+/QWCpPyemnj3qPaXFwhry8oZ5pN/lozNNHMEwAig+CEhFQ/G7OpuEJWSU677dRsTJU27j2qZ6YMT/jwdK7N3l5fUIU7ShsFJ5p5AkD7IzghIbV0NiYRtWazdzI18wSAeMQeJySkVLj1Ptvjlj/YfNNKf21I2RnuKFcEACA4ISGlwq33bPYGgPjDUh0SUirces9mbyABBaqlnaulXWvNA5gzcsxjb4ZOpRlpkiA4ISFNHNZTi9dXyRcMNVquS6bZGDZ7AwkmUC29O0fav0GSRUpzmQcwly+U9n1GJ/ckQXBCQkqV2Rg2ewMJZOdqMzS5O5uhqV7QZ17fuZru7kmAPU5ISPWzMXOvP1+DuqfLkWbRoO7pmnv9+UnRigBAAtq1Vg0zTaezuySL9dQ4Eh0zTkhYzMYgKbFHJnF5S5qGpno2p+RN/Lt9QXACgPjBHpnElpFj/nk1J+SXsvpFtRy0D5bqACBenL5HJj1bcmeaj67O3+2RQfwaPFmSYe5pOl3QJxl1p8aR6AhOABAv2COT2IZOlfqNlnzHpBOHpZrj5qPvmHl96NRYV4gIYKkOAOIFe2QSm6ODuZzasEet1FyeY49aUiE4AUC8YI9M4nN0MFsO0HYgabFUBwDxgj0yQNxjxgkA4sXQqebdc/s3SH6ruTwX8puhiT0yQFwgOAFtVBMIfXcsSmWNsj1ujkVB67BHBoh7FsMwjFgXES1er1cej0eVlZXKyMiIdTlIAjWBkOat2a5NxRWySnLabfIHQ6qTNGpAFl3MASABtCQfMOMEtEFB0SFtKq5Qptv+3WHDbrt8wZA2F1eooOgQnc0BoLXisJM+m8OBNijcUSqr9F1oOsVlt8liMccBAK1Q30n/44Xm3aa1AfPx44Xm9UB1TMoiOAFtUFZZI6e9+aU4Z5pNZd6aKFcEAEkiTjvpE5yANsj2uOUPhpod89eGlJ3hjnJFAOJSoFratkp6dZL0wuXm47ZVMZs1SQhx2kmf4AS0wcRhPVUnyXdGePIFQzIMcxxAiovTJae4F6ed9AlOQBtMyuulUQOyVFkTVHmVr9HjFQOyNCmvV6xLBBBrcbrkFPcycqRaX/NjIb+UEZv/MU2I4LR//37de++96t+/v9xutwYOHKgnn3xSgUAg1qUhxbkdNj0zZbjmXn++BnVPlyPNokHd0zX3+vNpRQDAFKdLTnEvTjvpJ0Q7gq+//lp1dXX6zW9+o/POO09/+9vfdP/99+vkyZN67rnnYl0eUpzbYdO0Ebm0HQDQvDhdcop7cdpJP2EbYD777LNatmyZiouLw/4ZGmACAKLu1Unmnqb07KZjJw5L3S6QflwQ/boSQaM+TqXm8lw79HFKiQaYlZWVysrKinUZAACc2+DJUvlCc4nJftrMU7SWnOKwiWTYHB2k/BnmV5xIyOC0d+9eLVmyRIsXLz7n8/x+v/x+f8P3Xq+3vUsDAKCxWC451d/Rt3+DGvZZlX9tBrl9n5lnI8Z7eIozMd0c/tRTT8lisZzz64svvmj0MyUlJRo/frymTp2q++6775yvv2jRInk8noavPn36tOevAwBAU/WHN4+dby7LpTnNx7Hz2z+4cEdfxMV0j9PRo0d19OjRcz6nX79+crnMqc2SkhKNHTtWI0eO1MqVK2W1njv3NTfj1KdPH/Y4AQBSA/urwpIwe5y6du2qrl27hvXcQ4cOaezYscrPz9eKFSu+NzRJktPplNPpbGuZAAAkJu7oi7iE2ONUUlKiMWPGKDc3V88995yOHDnSMJad3UyKBgAA5kbw8q+bHwv5pax+US0nGSREcPrwww+1Z88e7dmzR7179240lqDdFAAAaH+xvqMvCSVsH6fWoI8TgJSSyLehIzJOv6vO0swdfdxVJ6ll+YDgBADJqLnb0Gt9kgz+g5lqotREMpElzOZwAEA7Of029LQzlmjqb0OPo6aCaEdx2EQykSXEIb8AgBbiYFmgXRCcACAZcRs60C5YqgOAZJRst6Gz0R1xguDURjWBkAqKDqlwR6nKKmuU7XFr4rCempTXS26HLdblAUhVyXQbOuetIY4QnNqgJhDSvDXbtam4QlZJTrtNuw9XafH6Km3ce1TPTBlOeAIQG7E8WDbS2OiOOMIepzYoKDqkTcUVynTb1T3DJc+px0y3XZuLK1RQdCjWJQJIVbE8WDbS2OiOOMKMUxsU7iiVVZLL3nhWyWW3yesLqnBHqaaNyI1NcQCQLLehs9EdcYQZpzYoq6yR0978UpwzzaYyb02UKwKAJJSRc6p5ZzNCfrOhIxAlBKc2yPa45Q+Gmh3z14aUneGOckUAkIQGT5ZkmHuaTpeIG92R8AhObTBxWE/VSfKdEZ58wZAMwxwHALTR0KnmhnbfMenEYanmuPnoO5Z4G92R8Njj1AaT8npp496j2lxcIa8vKGeaTf5aMzRdMSBLk/J6xbpEAEh89RvdTz9vLasffZwQExzy20aN+jh5a5SdQR8nAAASSUvyAcEJAACktJbkA/Y4AQAAhIngBAAAECaCEwAAQJgITgAAAGGiHQEAAHGsrq5OgUAg1mUkNLvdLpstMne6E5wAAIhTgUBA+/btU11dXaxLSXiZmZnKzs6WxWJp0+sQnAAAiEOGYai0tFQ2m019+vSR1crumtYwDEPV1dUqLy+XJPXs2bZTPQhOAADEodraWlVXVysnJ0cdOtAdvS3cbvPs2PLycnXv3r1Ny3bEVwAA4lAoZJ6D6nA4YlxJcqgPn8FgsE2vQ3ACACCOtXVPDkyR+hwJTgAAoFVmzpypm2++OW5eJxoITgAAJKCZM2fKYrHIYrHIbrdrwIABevTRR3Xy5MlYl3ZW+/fvl8ViUVFRUaPr//3f/62VK1fGpKaWYnM4AAAJavz48VqxYoWCwaA2bNig++67TydPntSyZctiXVqLeDyeWJcQNmacAABIUE6nU9nZ2erTp4/uuOMOTZ8+XevWrZPf79e//uu/qnv37nK5XLr66qu1devWhp/75JNPZLFY9N5772n48OFyuVwaOXKkdu7c2fCcp556Snl5eY3e75e//KX69et31no++OADXX311crMzFSXLl00ceJE7d27t2G8f//+kqRLLrlEFotFY8aMkdR0qS7c+j/66CNddtll6tChg6688kp98803rfgUW4bgBKSImkBIb205oDtf/ouuW/yJ7nz5L3prywHVBEKxLg1AhLjdbgWDQT3++ONas2aNVq1apS+//FLnnXeefvCDH6iioqLR8x977DE999xz2rp1q7p3766bbrqpTXednTx5UnPnztXWrVv10UcfyWq16pZbbmlo4LllyxZJ0h//+EeVlpZq7dq1zb5OuPXPnz9fixcv1hdffKG0tDTdc889ra49XAQnIAXUBEKat2a7Fq//h3YfrlIgZGj34SotXv8PzVuznfAEJIEtW7bojTfe0NixY7Vs2TI9++yzmjBhgi6++GK99NJLcrvdWr58eaOfefLJJ3X99ddr6NChWrVqlQ4fPqx33nmn1TVMmTJFkydP1qBBg5SXl6fly5dr586d+vvf/y5J6tatmySpS5cuys7OVlZWVpPXqF9qDKf+hQsX6tprr9XFF1+sn/3sZ9q4caN8Pl+r6w8HwQlIAQVFh7SpuEKZbru6Z7jkOfWY6bZrc3GFCooOxbpEAK1QWFioTp06yeVyadSoUbrmmms0Z84cBYNBXXXVVQ3Ps9vtGjFihL766qtGPz9q1KiGv87KytIFF1zQ5DktsXfvXt1xxx0aMGCAMjIyGpbmDhw40KLXCLf+YcOGNfx1fUfw+g7h7YXN4Ug5NYGQCooOqXBHqcoqa5TtcWvisJ6alNdLbkdkDoGMN4U7SmWV5LI3/v1cdpu8vqAKd5Rq2ojc2BQHoNXqZ5fsdrtycnJkt9u1fft2SU37FhmGEVYvo/rnWK1WGYbRaOz7lvFuvPFG9enTRy+99JJycnJUV1enIUOGtOiQ4vr3DKd+u93epO72PtePGSeklFRdsiqrrJHT3nwodKbZVOatiXJFACKhY8eOOu+889S3b9+GEHHeeefJ4XDo888/b3heMBjUF198oYsuuqjRz2/evLnhr48dO6Z//OMfuvDCCyWZy2plZWWNwtOZbQRO9+233+qrr77Sz3/+c40bN04XXXSRjh071ug59V3Q67uiN6cl9ccCM05IKacvWTXMvrjt8gVDDUtWyTjzku1xa/fhKsltbzLmrw0pNys9BlUBaA8dO3bUT37yEz322GPKyspSbm6ufvGLX6i6ulr33ntvo+c+/fTT6tKli3r06KH58+era9euDXe3jRkzRkeOHNEvfvEL3Xrrrfrggw/0/vvvKyMjo9n37dy5s7p06aIXX3xRPXv21IEDB/Szn/2s0XO6d+8ut9utDz74QL1795bL5WrSiqAl9ccCM05IKedasrJYzPFkNHFYT9VJ8gUb/1+eLxiSYZjjAJLH//7f/1tTpkzRXXfdpUsvvVR79uzR//t//0+dO3du8ryHH35Y+fn5Ki0t1bvvvtswK3TRRRdp6dKl+tWvfqXhw4dry5YtevTRR8/6nlarVW+99Za2bdumIUOG6Kc//ameffbZRs9JS0vT888/r9/85jfKycnRpEmT2lR/LFiMMxcwk5jX65XH41FlZeVZEzOS23WLP1EgZMjTzMxLZU1QjjSL/jh3TPQLa2f1S5SbiytksZjLc/5aMzRdMSBLz0wZnrT7u4BE5fP5tG/fPvXv318ulyuir/3JJ59o7NixOnbsmDIzMyP62vHqXJ9nS/IBS3VIKam6ZOV22PTMlOHfbYr31ig3Kz3pN8UDQKQRnJBSJg7rqcXrq+QLhhot16XCkpXbYdO0EblJuYcLAKKF4ISUMimvlzbuParNxRXy+oJNlqwm5fWKdYkA0O7GjBnTpNUAwkNwQkphyQoA0BYEJ6QclqwAAK1FOwIAAIAwEZwAAADCRHACAAAIE8EJAAAgTAQnAACAMBGcAACIhUC1tG2V9Ook6YXLzcdtq8zrKWrjxo2y2WwaP358rEs5K4ITAADRFqiW3p0jfbxQKv9aqg2Yjx8vNK+naHh65ZVXNGfOHH3++ec6cOBArMtpFsEJAIBo27la2r9BcneW0rMld6b56OpsXt+5OtYVRt3Jkyf1u9/9Tj/5yU80ceJErVy5MtYlNYvgBABAtO1aK8kipbkaX7e7JIv11Hhs1ARCemvLAd358l903eJPdOfLf9FbWw6oJhBq1/d9++23dcEFF+iCCy7QnXfeqRUrVsTlsTAEJwAAos1b0jQ01bM5JW9pdOs5pSYQ0rw127V4/T+0+3CVAiFDuw9XafH6f2jemu3tGp6WL1+uO++8U5I0fvx4nThxQh999FG7vV9rEZwAAIi2jByp1tf8WMgvZfSMbj2nFBQd0qbiCmW67eqe4ZLn1GOm267NxRUqKDrULu/7zTffaMuWLZo2bZokKS0tTbfffrteeeWVdnm/tuCsOgAAom3wZKl8oRT0mctz9YI+yagzxyXzr6uPSSdOSnUByeow90O5O0vWyB9KXrijVFZJLnvj13bZbfL6gircUdou53wuX75ctbW16tWrV8M1wzBkt9t17Ngxde7cOeLv2VrMOAEAEG1Dp0r9Rku+Y9KJw1LNcfPRd8y8PnSqFKiRao5JJ4+Ys1OGYT5WlUrHD0p1kV82K6uskdPefCBzptlU5q2J+HvW1tbq1Vdf1eLFi1VUVNTwtX37dvXt21evv/56xN+zLZhxAgAg2hwdpJuWmHfP7Vpr7mnK6mfONA2dao7vLJBqu0nWzpLttDBjGFKgygxVHbtGtKxsj1u7D1dJbnuTMX9tSLlZ6RF9P0kqLCzUsWPHdO+998rj8TQau/XWW7V8+XLNnj074u/bWsw4AQAQC44OUv4M6ccF0uwt5mP+DPO6JO1ebz5aLI1/zmKRZDFnqSJs4rCeqpPkCzaezfIFQzIMczzSli9fruuuu65JaJKkKVOmqKioSF9++WXE37e1mHECgFQXqD5t5qPE3Lh8+swHYuNEuSRL82MWi1QXjPhbTsrrpY17j2pzcYW8vqCcaTb5a83QdMWALE3K6/X9L9JCv//97886dumll8ZdSwKCEwCksvoO1vs3qKGvUPnX5sblfZ+Zy0mEp9jo1F3SWUKDYUi2pstpbeV22PTMlOEqKDqkwh2lKvPWKDcrXROH9dSkvF5yOyK/IT3REJwAIJWd3sE67Yy7u+o7WOfPiF19qWzQ9VKtzJB0OsOQZJh317UDt8OmaSNy2+XuuWTAHicASGVx3ME65Q0ab/651IWkUO1pj7WSI90Mu4g6ZpwAIJXFaQdrSHK4zXDUMUuqO2nuabLZ27WPE74fwQkAUllGjrmnqTkhv3mLPGLHYpU6dJZcsekkjqZYqgOAVDZ4siTD3NN0ujM7WAOQxIwTAKS2oVPNu+f2b5D8VnN5LuQ3Q1N9B2sADQhOAJDKwulgDaABwQkAUl19B2vaDgDfK+H2OPn9fuXl5clisaioqCjW5QAAgBSScMHp8ccfV05OTqzLAAAAETJz5kxZLJaGry5dumj8+PHasWNHrEtrIqGC0/vvv68PP/xQzz33XKxLAaKqJhDSW1sO6M6X/6LrFn+iO1/+i97ackA1gdD3/zAAJIDx48ertLRUpaWl+uijj5SWlqaJEyfGuqwmEmaP0+HDh3X//fdr3bp16tAhvM2Kfr9ffr+/4Xuv19te5QHtpiYQ0rw127WpuEJWSU67TbsPV2nx+ipt3HtUz0wZzvlRABKe0+lUdna2JCk7O1vz5s3TNddcoyNHjqhbt24xru47CRGcDMPQzJkz9eCDD+qyyy7T/v37w/q5RYsW6T/+4z/atzignRUUHdKm4gpluu1y2U8FJLddvmBIm4srVFB0KOXPlKoJhL47lLSyRtkeN4eSAq0VqD7tLssSs0lqlO+yPHHihF5//XWdd9556tKlS1TeM1wxXap76qmnGq1pNvf1xRdfaMmSJfJ6vXriiSda9PpPPPGEKisrG74OHjzYTr8J0H4Kd5TKKn0Xmk5x2W2yWMzxVFY/I7d4/T+0+3CVAiHj1IzcPzRvzXaWM4GWCFRL786RPl5odpSvDZiPHy80rweq2+2tCwsL1alTJ3Xq1Enp6el699139fbbb8tqja9dRTGdcZo9e7amTZt2zuf069dPCxYs0ObNm+V0OhuNXXbZZZo+fbpWrVrV7M86nc4mPwMkmrLKGjntzc+aONNsKvPWRLmi+MKMHBBBO1ebzVDdnRufYRj0mdd3rm63thVjx47VsmXLJEkVFRVaunSpJkyYoC1btqhv377t8p6tEdPg1LVrV3Xt2vV7n/f8889rwYIFDd+XlJToBz/4gd5++22NHDmyPUsEYi7b49buw1WS295kzF8bUm5Wegyqih/nmpHz+oIq3FFKcALCtWutJEvTg5/tLrOz/K617RacOnbsqPPOO6/h+/z8fHk8Hr300kuNMkCsJcQep9zcxv/S69SpkyRp4MCB6t27dyxKAqJm4rCeWry+Sr5gqFE48AVDMgxzPJUxIwdEkLekaWiqZ3OaneWjxGKxyGq1qqYmvv4ZTojgBKSySXm9tHHvUW0urpDXF5QzzSZ/rRmarhiQpUl5vWJdYkwxIwdEUEaOuaepOSG/eRxPO/H7/SorK5MkHTt2TC+88IJOnDihG2+8sd3eszUSMjj169dPhmHEugwgKtwOm56ZMvy7u8a8NcrNSueusVOYkQMiaPBkqXyhuafJfsYeJ6POHG8nH3zwgXr2NP95TU9P14UXXqjVq1drzJgx7faerZGQwQlINW6HTdNG5LJXpxnMyAERNHSqtO8zcyO432ouz4X8ZmjqN9ocbwcrV67UypUr2+W1I43gBCChMSMHRJCjg3TTktP6OJWay3NR7uMUzwhOABIeM3JABDk6mHfOtdPdc4kuvrpKAQAAxDGCEwAAQJgITgAAAGEiOAEAEMdovxMZkfocCU4AAMQhm828IzQQCMS4kuRQXW0eUGy3N22W2xLcVQcAQBxKS0tThw4ddOTIEdntdlmtzHW0hmEYqq6uVnl5uTIzMxsCaWsRnAAAiEMWi0U9e/bUvn379M9//jPW5SS8zMxMZWdnt/l1CE4AAMQph8OhQYMGsVzXRna7vc0zTfUITgAAxDGr1SqXy/X9T0RUsGAKAAAQJoITAABAmAhOAAAAYUqpPU71za+8Xm+MKwEAAPGiPheE0yQzpYJTVVWVJKlPnz4xrgQAAMSbqqoqeTyecz7HYqRQL/e6ujqVlJQoPT1dFoul0ZjX61WfPn108OBBZWRkxKjC1MBnHT181tHF5x09fNbRleyft2EYqqqqUk5Ozvc2Gk2pGSer1arevXuf8zkZGRlJ+TdFPOKzjh4+6+ji844ePuvoSubP+/tmmuqxORwAACBMBCcAAIAwEZxOcTqdevLJJ+V0OmNdStLjs44ePuvo4vOOHj7r6OLz/k5KbQ4HAABoC2acAAAAwkRwAgAACBPBCQAAIEwEp7N47733NHLkSLndbnXt2lWTJ0+OdUlJze/3Ky8vTxaLRUVFRbEuJynt379f9957r/r37y+3262BAwfqySefVCAQiHVpSWHp0qXq37+/XC6X8vPztWHDhliXlJQWLVqkyy+/XOnp6erevbtuvvlmffPNN7EuKyUsWrRIFotFjzzySKxLiSmCUzPWrFmju+66S3fffbe2b9+uP//5z7rjjjtiXVZSe/zxx5WTkxPrMpLa119/rbq6Ov3mN7/Rrl279F//9V/69a9/rX//93+PdWkJ7+2339Yjjzyi+fPn669//atGjx6tCRMm6MCBA7EuLel8+umnmjVrljZv3qz169ertrZWN9xwg06ePBnr0pLa1q1b9eKLL2rYsGGxLiX2DDQSDAaNXr16GS+//HKsS0kZf/jDH4wLL7zQ2LVrlyHJ+Otf/xrrklLGL37xC6N///6xLiPhjRgxwnjwwQcbXbvwwguNn/3sZzGqKHWUl5cbkoxPP/001qUkraqqKmPQoEHG+vXrjWuvvdZ4+OGHY11STDHjdIYvv/xShw4dktVq1SWXXKKePXtqwoQJ2rVrV6xLS0qHDx/W/fffr9/+9rfq0KFDrMtJOZWVlcrKyop1GQktEAho27ZtuuGGGxpdv+GGG7Rx48YYVZU6KisrJYm/j9vRrFmz9MMf/lDXXXddrEuJCwSnMxQXF0uSnnrqKf385z9XYWGhOnfurGuvvVYVFRUxri65GIahmTNn6sEHH9Rll10W63JSzt69e7VkyRI9+OCDsS4loR09elShUEg9evRodL1Hjx4qKyuLUVWpwTAMzZ07V1dffbWGDBkS63KS0ltvvaUvv/xSixYtinUpcSNlgtNTTz0li8Vyzq8vvvhCdXV1kqT58+drypQpys/P14oVK2SxWLR69eoY/xaJIdzPesmSJfJ6vXriiSdiXXJCC/fzPl1JSYnGjx+vqVOn6r777otR5cnFYrE0+t4wjCbXEFmzZ8/Wjh079Oabb8a6lKR08OBBPfzww3rttdfkcrliXU7cSIt1AdEye/ZsTZs27ZzP6devn6qqqiRJF198ccN1p9OpAQMGsNEzTOF+1gsWLNDmzZubtPC/7LLLNH36dK1atao9y0wa4X7e9UpKSjR27FiNGjVKL774YjtXl/y6du0qm83WZHapvLy8ySwUImfOnDl699139dlnn6l3796xLicpbdu2TeXl5crPz2+4FgqF9Nlnn+mFF16Q3++XzWaLYYWxkTLBqWvXruratev3Pi8/P19Op1PffPONrr76aklSMBjU/v371bdv3/YuMymE+1k///zzWrBgQcP3JSUl+sEPfqC3335bI0eObM8Sk0q4n7ckHTp0SGPHjm2YSbVaU2bSud04HA7l5+dr/fr1uuWWWxqur1+/XpMmTYphZcnJMAzNmTNH77zzjj755BP1798/1iUlrXHjxmnnzp2Nrt1999268MILNW/evJQMTVIKBadwZWRk6MEHH9STTz6pPn36qG/fvnr22WclSVOnTo1xdcklNze30fedOnWSJA0cOJD/g2wHJSUlGjNmjHJzc/Xcc8/pyJEjDWPZ2dkxrCzxzZ07V3fddZcuu+yyhpm8AwcOsH+sHcyaNUtvvPGGCgoKlJ6e3jDT5/F45Ha7Y1xdcklPT2+yd6xjx47q0qVLSu8pIzg149lnn1VaWpruuusu1dTUaOTIkfrTn/6kzp07x7o0oNU+/PBD7dmzR3v27GkSTA3O+m6T22+/Xd9++62efvpplZaWasiQIfrDH/7ALHU7WLZsmSRpzJgxja6vWLFCM2fOjH5BSDkWg39jAgAAhIUNDgAAAGEiOAEAAISJ4AQAABAmghMAAECYCE4AAABhIjgBAACEieAEAAAQJoITAABAmAhOAAAAYSI4AUgaM2fOlMVikcVikd1u14ABA/Too4/q5MmTDc9Zs2aNxowZI4/Ho06dOmnYsGF6+umnVVFRIUkqLS3VHXfcoQsuuEBWq1WPPPJIjH4bAPGI4AQgqYwfP16lpaUqLi7WggULtHTpUj366KOSpPnz5+v222/X5Zdfrvfff19/+9vftHjxYm3fvl2//e1vJUl+v1/dunXT/PnzNXz48Fj+KgDiEGfVAUgaM2fO1PHjx7Vu3bqGa/fff78KCwtVUFCgkSNH6pe//KUefvjhJj97/PhxZWZmNro2ZswY5eXl6Ze//GX7Fg4gYTDjBCCpud1uBYNBvf766+rUqZMeeuihZp93ZmgCgOYQnAAkrS1btuiNN97QuHHjtHv3bg0YMEB2uz3WZQFIYAQnAEmlsLBQnTp1ksvl0qhRo3TNNddoyZIlMgxDFosl1uUBSHBpsS4AACJp7NixWrZsmex2u3JychpmmM4//3x9/vnnCgaDzDoBaDVmnAAklY4dO+q8885T3759GwWkO+64QydOnNDSpUub/bnjx49HqUIAiYwZJwApYeTIkXr88cf1b//2bzp06JBuueUW5eTkaM+ePfr1r3+tq6++uuFuu6KiIknSiRMndOTIERUVFcnhcOjiiy+O4W8AIB7QjgBA0miuHcGZfve73+lXv/qV/vrXv6qurk4DBw7Urbfeqjlz5jTcWdfcXqi+fftq//797VM4gIRBcAIAAAgTe5wAAADCRHACAAAIE8EJAAAgTAQnAACAMBGcAAAAwkRwAgAACBPBCQAAIEwEJwAAgDARnAAAAMJEcAIAAAgTwQkAACBMBCcAAIAw/X9JETBD5wYHIQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(6, 5))\n", "for label in df_meta['pop'].unique():\n", " idx = df_meta['pop'] == label\n", " plt.scatter(\n", " pc_df.loc[idx, 'PC1'],\n", " pc_df.loc[idx, 'PC2'],\n", " label=label,\n", " s=30,\n", " alpha=0.8\n", " )\n", "\n", "plt.xlabel('PC1')\n", "plt.ylabel('PC2')\n", "plt.title('PCA: PC1 vs PC2')\n", "plt.legend(title='Population', loc='lower right')\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "185eb587-ac46-4226-bb99-e51ec90724f0", "metadata": { "user_expressions": [] }, "source": [ "## Questions\n", "8. Explain what is PCA (principal component analysis) in simple terms and explain what does the PCA plot above demonstrate?" ] }, { "cell_type": "markdown", "id": "5feda329", "metadata": { "user_expressions": [] }, "source": [ "## Notes & Next steps\n", "\n", "- This assignment uses **simulated data** \n", "- For real GWAS, use PLINK, Hail, or SAIGE for large-scale analyses." ] }, { "cell_type": "code", "execution_count": null, "id": "0c7c5717-9594-4491-a2da-e2bf81bc452f", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python (base)", "language": "python", "name": "base" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.7" } }, "nbformat": 4, "nbformat_minor": 5 }