diff --git a/openproblems/tasks/regulatory_effect_prediction/methods/__init__.py b/openproblems/tasks/regulatory_effect_prediction/methods/__init__.py index 37894113cd..683566546a 100644 --- a/openproblems/tasks/regulatory_effect_prediction/methods/__init__.py +++ b/openproblems/tasks/regulatory_effect_prediction/methods/__init__.py @@ -1,3 +1,5 @@ from .beta import archr_model21 from .beta import beta from .beta import marge +from .maestro import rp_enhanced +from .maestro import rp_simple diff --git a/openproblems/tasks/regulatory_effect_prediction/methods/beta.py b/openproblems/tasks/regulatory_effect_prediction/methods/beta.py index cc348e6002..e8ede814f8 100644 --- a/openproblems/tasks/regulatory_effect_prediction/methods/beta.py +++ b/openproblems/tasks/regulatory_effect_prediction/methods/beta.py @@ -1,5 +1,7 @@ from ....patch import patch_datacache from ....tools.decorators import method +from .maestro import _rp_enhanced +from .maestro import _rp_simple import numpy as np import pandas as pd @@ -28,6 +30,7 @@ def _chrom_limit(x, tss_size=2e5): return [gene_end - tss_size // 2, gene_end + tss_size // 2] +# included 2-3 lines to make it run with exons. def _get_annotation(adata, retries=3): """Insert meta data into adata.obs.""" from pyensembl import EnsemblRelease @@ -55,12 +58,7 @@ def _get_annotation(adata, retries=3): ) gene = data.gene_by_id(i) genes.append( - [ - "chr%s" % gene.contig, - gene.start, - gene.end, - gene.strand, - ] + ["chr%s" % gene.contig, gene.start, gene.end, gene.strand, gene.exons] ) except ValueError: try: @@ -76,6 +74,7 @@ def _get_annotation(adata, retries=3): gene.start, gene.end, gene.strand, + gene.exons, ] ) except (IndexError, ValueError) as e: @@ -86,7 +85,7 @@ def _get_annotation(adata, retries=3): [adata.var, pd.DataFrame(genes, index=adata.var_names)], axis=1 ) adata.var.columns = np.hstack( - [old_col, np.array(["chr", "start", "end", "strand"])] + [old_col, np.array(["chr", "start", "end", "strand", "exons"])] ) @@ -135,7 +134,7 @@ def _filter_has_chr(adata): return adata -def _atac_genes_score(adata, top_genes=2000, threshold=1, method="beta"): +def _atac_genes_score(adata, top_genes=2000, threshold=1, method="beta", **kwargs): """Calculate gene scores and insert into .obsm.""" import pybedtools @@ -191,6 +190,8 @@ def _atac_genes_score(adata, top_genes=2000, threshold=1, method="beta"): axis=1, ) + extend_tss["gene_short_name"] = adata.var["gene_short_name"] + # peak summits peaks = pd.DataFrame( { @@ -211,8 +212,10 @@ def _atac_genes_score(adata, top_genes=2000, threshold=1, method="beta"): ) # overlap TSS bins with peaks + x = pybedtools.BedTool.from_dataframe(summits) y = pybedtools.BedTool.from_dataframe(extend_tss) + tss_to_peaks = x.intersect(y, wb=True, wa=True, loj=True).to_dataframe() # remove non-overlapped TSS and peaks @@ -226,6 +229,18 @@ def _atac_genes_score(adata, top_genes=2000, threshold=1, method="beta"): _archr_model21(tss_to_peaks, adata) elif method == "marge": _marge(tss_to_peaks, adata) + elif method == "rp_simple": + _rp_simple(tss_to_peaks, adata) + elif method == "rp_enhanced": + _rp_enhanced(tss_to_peaks, adata) + + # the genes and gene_scores have to have the same dimensions + same_shape = adata.shape == adata.obsm["gene_score"].shape + if not same_shape: + print(adata.shape, adata.obsm["gene_score"].shape) + print("dimensions are not the same. Check calculation/filters") + assert same_shape + return adata diff --git a/openproblems/tasks/regulatory_effect_prediction/methods/maestro.py b/openproblems/tasks/regulatory_effect_prediction/methods/maestro.py new file mode 100644 index 0000000000..9389d7f34d --- /dev/null +++ b/openproblems/tasks/regulatory_effect_prediction/methods/maestro.py @@ -0,0 +1,200 @@ +from ....tools.decorators import method + +import pandas as pd + + +def _rp_simple(tss_to_peaks, adata, log_each=500): + import numpy as np + import scipy + + # the coordinates of the current peaks + peaks = adata.uns["mode2_var"] + + decay = 10000 + + def Sg(x): + return 2 ** (-x) + + # gene_distance = 15 * decay + + weights = [] + + tss_to_peaks = tss_to_peaks.drop_duplicates("itemRgb") + + # print(tss_to_peaks.shape) + # print(tss_to_peaks.head()) + + for ri, r in tss_to_peaks.iterrows(): + wi = 0 + summit_chr, tss_summit_start, tss_summit_end = r[:3] + tss_extend_chr, tss_extend_start, tss_extend_end = r[4:7] + + # print(summit_chr, summit_start, summit_end, + # extend_chr, extend_start, extend_end) + sel_chr = [pi for pi in peaks if pi[0] == tss_extend_chr] + sel_peaks = [ + pi + for pi in sel_chr + if int(pi[1]) >= tss_extend_start and int(pi[2]) <= tss_extend_end + ] + + # print('# peaks in chromosome', len(sel_chr), + # '# of peaks around tss', len(sel_peaks)) + # if len(sel_peaks) > 0: + # print(sel_peaks) + + # if peaks then this is take them into account, one by one + for pi in sel_peaks: + summit_peak = int((int(pi[2]) + int(pi[1])) / 2) + distance = np.abs(tss_summit_start - summit_peak) + # print(pi, distance, Sg(distance / decay)) + wi += Sg(distance / decay) + + if log_each is not None and len(weights) % log_each == 0: + if len(weights) > 0: + print( + "# weights calculated so far", + len(weights), + "out of", + tss_to_peaks.shape[0], + ) + weights.append(wi) + + tss_to_peaks["weight"] = weights + + gene_peak_weight = scipy.sparse.csr_matrix( + ( + tss_to_peaks.weight.values, + (tss_to_peaks.thickEnd.astype("int32").values, tss_to_peaks.name.values), + ), + shape=(adata.shape[1], adata.uns["mode2_var"].shape[0]), + ) + + adata.obsm["gene_score"] = adata.obsm["mode2"] @ gene_peak_weight.T + + +def _rp_enhanced(tss_to_peaks, adata, log_each=500): + import numpy as np + import scipy + + # prepare the exonic ranges + exon_ranges = [] + for exons in adata.var["exons"]: + exon_ranges.append([e.start, e.end] for e in exons) + adata.var["exon.ranges"] = exon_ranges + exon_coordinates_by_gene = adata.var["exon.ranges"].to_dict() + + tss_to_peaks["exon.ranges"] = tss_to_peaks["itemRgb"].map(exon_coordinates_by_gene) + tss_to_peaks = tss_to_peaks.drop_duplicates("itemRgb").reset_index(drop=True) + + # the coordinates of the current peaks + peaks = adata.uns["mode2_var"] + + decay = 10000 + + def Sg(x): + return 2 ** (-x) + + print("calculating weights per gene...") + weights = [] + for ri, r in tss_to_peaks.iterrows(): + wi = 0 + summit_chr, tss_summit_start, tss_summit_end = r[:3] + tss_extend_chr, tss_extend_start, tss_extend_end = r[4:7] + + # gene_name = r[-2] + exon_ranges = r[-1] + + # print(summit_chr, tss_summit_start, tss_summit_end, + # tss_extend_chr, tss_extend_start, tss_extend_end, + # gene_name, exon_ranges) + sel_chr = [pi for pi in peaks if pi[0] == tss_extend_chr] + sel_peaks = [ + pi + for pi in sel_chr + if int(pi[1]) >= tss_extend_start and int(pi[2]) <= tss_extend_end + ] + # check whether the peak overlaps with a given exon + if not pd.isnull(exon_ranges): + sel_peak_summits = [(int(pi[1]) + int(pi[2])) / 2.0 for pi in sel_peaks] + peak_in_exons = [ + np.sum([ps >= ex[0] and ps <= ex[1] for ex in exon_ranges]) >= 1 + for ps in sel_peak_summits + ] + else: + peak_in_exons = [False for pi in sel_peaks] + + # if sum(peak_in_exons) > 0: + # print ('exon / peak overlap found!') + # print(ri, peak_in_exons) + # if peaks then this is take them into account, one by one + for pi, peak_in_exon in zip(sel_peaks, peak_in_exons): + # the current peak is part of an exon + # if peak_in_exon: + # print(pi) + summit_peak = int((int(pi[2]) + int(pi[1])) / 2) + distance = np.abs(tss_summit_start - summit_peak) + # print(pi, distance, Sg(distance / decay)) + wi += Sg(distance / decay) if not peak_in_exon else 1.0 + + if log_each is not None and len(weights) % log_each == 0: + if len(weights) > 0: + print( + "# weights calculated so far", + len(weights), + "out of", + tss_to_peaks.shape[0], + ) + + weights.append(wi) + + out = tss_to_peaks.copy() + out["weight"] = weights + + gene_peak_weight = scipy.sparse.csr_matrix( + ( + out.weight.values, + (out.thickEnd.astype("int32").values, out.name.values), + ), + shape=(adata.shape[1], adata.uns["mode2_var"].shape[0]), + ) + + adata.obsm["gene_score"] = adata.obsm["mode2"] @ gene_peak_weight.T + + +@method( + method_name="RP_simple", + paper_name="""Integrative analyses of single-cell transcriptome\ +and regulome using MAESTRO.""", + paper_url="https://pubmed.ncbi.nlm.nih.gov/32767996", + paper_year=2020, + code_version="1.0", + code_url="https://github.com/liulab-dfci/MAESTRO", + image="openproblems-python-extras", +) +def rp_simple(adata, n_top_genes=2000): + from .beta import _atac_genes_score + + adata = _atac_genes_score( + adata, + top_genes=n_top_genes, + method="rp_simple", + ) + return adata + + +@method( + method_name="RP_enhanced", + paper_name="""Integrative analyses of single-cell transcriptome\ +and regulome using MAESTRO.""", + paper_url="https://pubmed.ncbi.nlm.nih.gov/32767996", + paper_year=2020, + code_version="1.0", + code_url="https://github.com/liulab-dfci/MAESTRO", + image="openproblems-python-extras", +) +def rp_enhanced(adata, n_top_genes=2000): + from .beta import _atac_genes_score + + adata = _atac_genes_score(adata, top_genes=n_top_genes, method="rp_enhanced") + return adata diff --git a/openproblems/tasks/regulatory_effect_prediction/tests/snare_chrompotential_maestro.ipynb b/openproblems/tasks/regulatory_effect_prediction/tests/snare_chrompotential_maestro.ipynb new file mode 100644 index 0000000000..fb37fadad5 --- /dev/null +++ b/openproblems/tasks/regulatory_effect_prediction/tests/snare_chrompotential_maestro.ipynb @@ -0,0 +1,568 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "from openproblems.tasks.regulatory_effect_prediction import datasets, methods" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "tags": [] + }, + "source": [ + "## I. Use the utility function to prepare a subset of genes that are HVG and also mappable to chromosomes/TSS" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [] + }, + "outputs": [], + "source": [ + "from openproblems.patch import patch_datacache\n", + "import numpy as np\n", + "import pandas as pd\n", + "import scanpy as sc\n", + "import warnings\n", + "\n", + "def _chrom_limit(x, tss_size=2e5):\n", + " \"\"\"Extend TSS to upstream and downstream intervals.\n", + "\n", + " Parameters\n", + " ----------\n", + " x : pd.Series\n", + " a pd.Series containing [start, end, direction]\n", + " where start and end are ints and direction is {'+', '-'}.\n", + " tss_size: int\n", + " a int that defines the upstream and downstream regions around TSS\n", + " \"\"\"\n", + " y = x.values\n", + " gene_direction = y[-1]\n", + " gene_start = y[-3]\n", + " gene_end = y[-2]\n", + " if gene_direction == \"+\":\n", + " return [gene_start - tss_size // 2, gene_start + tss_size // 2]\n", + " else:\n", + " return [gene_end - tss_size // 2, gene_end + tss_size // 2]\n", + "\n", + "\n", + "def _get_annotation(adata, retries=3):\n", + " \"\"\"Insert meta data into adata.obs.\"\"\"\n", + " from pyensembl import EnsemblRelease\n", + "\n", + " data = EnsemblRelease(\n", + " adata.uns[\"release\"],\n", + " adata.uns[\"species\"],\n", + " )\n", + " for _ in range(retries):\n", + " try:\n", + " with patch_datacache():\n", + " data.download(overwrite=False)\n", + " data.index(overwrite=False)\n", + " break\n", + " except TimeoutError:\n", + " pass\n", + "\n", + " # get ensemble gene coordinate\n", + " genes = []\n", + " for i in adata.var.index.map(lambda x: x.split(\".\")[0]):\n", + " try:\n", + " with warnings.catch_warnings():\n", + " warnings.filterwarnings(\n", + " action=\"ignore\", message=\"No results found for query\"\n", + " )\n", + " gene = data.gene_by_id(i)\n", + " genes.append(\n", + " [\n", + " \"chr%s\" % gene.contig,\n", + " gene.start,\n", + " gene.end,\n", + " gene.strand,\n", + " gene.exons\n", + " ]\n", + " )\n", + " except ValueError:\n", + " try:\n", + " with warnings.catch_warnings():\n", + " warnings.filterwarnings(\n", + " action=\"ignore\", message=\"No results found for query\"\n", + " )\n", + " i = data.gene_ids_of_gene_name(i)[0]\n", + " gene = data.gene_by_id(i)\n", + " genes.append(\n", + " [\n", + " \"chr%s\" % gene.contig,\n", + " gene.start,\n", + " gene.end,\n", + " gene.strand,\n", + " gene.exons \n", + " ]\n", + " )\n", + " except (IndexError, ValueError) as e:\n", + " # print(e)\n", + " genes.append([np.nan, np.nan, np.nan, np.nan])\n", + " old_col = adata.var.columns.values\n", + " adata.var = pd.concat(\n", + " [adata.var, pd.DataFrame(genes, index=adata.var_names)], axis=1\n", + " )\n", + " adata.var.columns = np.hstack(\n", + " [old_col, np.array([\"chr\", \"start\", \"end\", \"strand\", 'exons'])]\n", + " )\n", + "\n", + "\n", + "def _filter_mitochondrial(adata):\n", + " if adata.uns[\"species\"] in [\"mus_musculus\", \"homo_sapiens\"]:\n", + " adata.var[\"mt\"] = adata.var.gene_short_name.str.lower().str.startswith(\n", + " \"mt-\"\n", + " ) # annotate the group of mitochondrial genes as 'mt'\n", + " sc.pp.calculate_qc_metrics(\n", + " adata, qc_vars=[\"mt\"], percent_top=None, log1p=False, inplace=True\n", + " )\n", + "\n", + " adata_filter = adata[adata.obs.pct_counts_mt <= 10]\n", + " if adata_filter.shape[0] > 100:\n", + " adata = adata_filter.copy()\n", + " return adata\n", + "\n", + "\n", + "def _filter_n_genes_max(adata):\n", + " adata_filter = adata[adata.obs.n_genes_by_counts <= 2000]\n", + " if adata_filter.shape[0] > 100:\n", + " adata = adata_filter.copy()\n", + " return adata\n", + "\n", + "\n", + "def _filter_n_genes_min(adata):\n", + " adata_filter = adata.copy()\n", + " sc.pp.filter_cells(adata_filter, min_genes=200)\n", + " if adata_filter.shape[0] > 100:\n", + " adata = adata_filter\n", + " return adata\n", + "\n", + "\n", + "def _filter_n_cells(adata):\n", + " adata_filter = adata.copy()\n", + " sc.pp.filter_genes(adata_filter, min_cells=5)\n", + " if adata_filter.shape[1] > 100:\n", + " adata = adata_filter\n", + " return adata\n", + "\n", + "\n", + "def _filter_has_chr(adata):\n", + " adata_filter = adata[:, ~pd.isnull(adata.var.loc[:, \"chr\"])].copy()\n", + " if adata_filter.shape[1] > 100:\n", + " adata = adata_filter\n", + " return adata" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**This step takes like five minutes for the whole object. Please wait.**\n", + "## In the context of testing, the pre-analysis here check for genes that are\n", + "i. Highly variable.\n", + "ii. Chromosome+exons mappable." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "obtaining annotation...\n", + "done...\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "... storing 'chr' as categorical\n", + "... storing 'strand' as categorical\n" + ] + } + ], + "source": [ + "# test = False does not work as tss_to_peaks need to be sub-sampled, and it seems that everything is blended.\n", + "adata = datasets.snare_p0_braincortex(test=False)\n", + "\n", + "top_genes = 2000\n", + "sc.pp.normalize_total(adata, target_sum=1e4)\n", + "sc.pp.log1p(adata)\n", + "\n", + "if top_genes <= adata.shape[1]:\n", + " sc.pp.highly_variable_genes(adata, n_top_genes=top_genes)\n", + " adata = adata[:, adata.var.highly_variable].copy()\n", + " \n", + "# get annotation for TSS\n", + "print('obtaining annotation...')\n", + "_get_annotation(adata)\n", + "print('done...')\n", + "\n", + "# basic quality control\n", + "adata = _filter_has_chr(adata)\n", + "adata = _filter_mitochondrial(adata)\n", + "adata = _filter_n_genes_max(adata)\n", + "adata = _filter_n_genes_min(adata)\n", + "adata = _filter_n_cells(adata)\n", + "\n", + "# regress out and scale\n", + "sc.pp.regress_out(adata, [\"total_counts\", \"pct_counts_mt\"])\n", + "sc.pp.scale(adata, max_value=10)\n", + "\n", + "sel_genes = set(adata.var.index)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## II. Once an annotation has been prepared, we test the main methods by just sampling the selected genes from before" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "before filtering (5081, 19322)\n", + "after filtering (5081, 1234)\n" + ] + } + ], + "source": [ + "# test = False does not work as tss_to_peaks need to be sub-sampled, and it seems that everything is blended.\n", + "adata = datasets.snare_p0_braincortex(test=False)\n", + "print('before filtering', adata.shape)\n", + "adata = adata[:, adata.var.index.isin(sel_genes)]\n", + "print('after filtering', adata.shape)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## RP-basic" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "... storing 'chr' as categorical\n", + "... storing 'strand' as categorical\n", + "/home/icb/ignacio.ibarra/miniconda3/envs/openproblems/lib/python3.7/site-packages/pandas/core/arrays/categorical.py:2487: FutureWarning: The `inplace` parameter in pandas.Categorical.remove_unused_categories is deprecated and will be removed in a future version.\n", + " res = method(*args, **kwargs)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "# weights calculated so far 500 out of 1232\n", + "# weights calculated so far 1000 out of 1232\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/mnt/znas/icb_zstore01/groups/ml01/workspace/ignacio.ibarra/SingleCellOpenProblems/openproblems/tasks/regulatory_effect_prediction/methods/maestro.py:63: SettingWithCopyWarning: \n", + "A value is trying to be set on a copy of a slice from a DataFrame.\n", + "Try using .loc[row_indexer,col_indexer] = value instead\n", + "\n", + "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n", + " tss_to_peaks[\"weight\"] = weights\n" + ] + } + ], + "source": [ + "adata = methods.rp_simple(adata, n_top_genes=2000) # log_each=10)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(0.010924862357515821, -0.010644082225420766)\n" + ] + } + ], + "source": [ + "%autoreload 2\n", + "import seaborn as sns\n", + "from openproblems.tasks.regulatory_effect_prediction import metrics\n", + "cors = metrics.spearman_correlation(adata), metrics.pearson_correlation(adata)\n", + "print(cors)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING: You’re trying to run this on 1233 dimensions of `.X`, if you really want this, set `use_rep='X'`.\n", + " Falling back to preprocessing with `sc.pp.pca` and default params.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(5081, 1233) (5081, 1233)\n" + ] + }, + { + "data": { + "text/plain": [ + "Text(0.5, 0, 'correlations')" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEICAYAAABYoZ8gAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8+yak3AAAACXBIWXMAAAsTAAALEwEAmpwYAAAsJUlEQVR4nO3deXxU9b3/8ddnZrKQkJB9gQTCjuxKUBEFV9zX2kqr1lp/te3tYpfbXtve7u29Xe5tb3t7e61trXhbl1rUUrVuoGgVlYCGsENCgJCdQDay5/P7YyYYIZDJMnNm+Twfjzwyc2bmnM+ZTN45+Z7v+X5FVTHGGBM9XE4XYIwxJrgs+I0xJspY8BtjTJSx4DfGmChjwW+MMVHGgt8YY6JMwIJfRB4QkVoR2TrAY18WERWRjEBt3xhjzMA8AVz3g8CvgIf6LxSRfGAFcMDfFWVkZGhBQcFo1maMMRFv06ZN9aqaeeLygAW/qr4qIgUDPPRz4KvAX/1dV0FBAUVFRaNVmjHGRAUR2T/Q8qC28YvI9cAhVS3247l3i0iRiBTV1dUFoTpjjIkOQQt+EUkAvg58y5/nq+r9qlqoqoWZmSf9p2KMMWaYgnnEPxWYDBSLSDmQB2wWkZwg1mCMMVEvkCd330dVS4Csvvu+8C9U1fpg1WCMMSaw3TkfATYAM0WkQkTuCtS2jDHG+C+QvXo+PMjjBYHatjHGmFOzK3eNMSbKWPAbv/T02oQ9xkSKoJ3cNeGnobWT/3xhF8+UVHH0WBeLC1L54mUzOG+qjbRhTDizI34zoK2HGrn0Z+t5dONBLp6ZxSeXT+FgQxu3/u4t/lZc6XR5xpgRsCN+c5JtlY3c+ru3GBvn4eFPnMOsnGQAvnjpDD76+7f50p/fJTMpjnOnpDtcqTFmOOyI37zPkdZO7n5oEwmxbh69+9zjoQ8QH+Pmtx8tJC81gX9ZvYWO7h4HKzXGDJcFvzlOVfny48XUNXdw322LyE9LOOk54xJi+O51c9h/+BgP/KM8+EUaY0bMgt8c99jGg6zbWcvXr5rFgvyUUz5v2YxMLj0jm1+t20NDa2fwCjTGjAoLfgNAbVM7//bsDs6dksYd5xUM+vx/uWImrZ09PPzWgKO+GmNCmAW/AeAHz+ygvbuXf79pPiIy6POnZyexfEYmqzbst7Z+Y8KMBb9h0/4jrCmu5FPLpjA5I9Hv1911/mTqmjt4urgqgNUZY0abBX+UU1V+8Mx2spLi+OTyqUN67QXTM5iWNZY/WXOPMWHFgj/KvbC9hncOHOWfV8wkMW5ol3WICB9clMfmA0cpq2sJUIXGmNFmwR/FenuVn7+4mymZidx01oRhrePGMyfgEli9uWKUqzPGBIoFfxR7bls1O6ubueeS6Xjcw/soZCXHs2xGJk9uPkSvDeRmTFiw4I9Sqsov1+5hWtZYrpk/fkTruumsPCob23m7vGGUqjPGBJIFf5Rav7uOndXNfGr5VNyuwbtvns4ls7KI9bh4flv1KFVnjAkkC/4o9Zv1ZeQkx3PdgpEd7QMkxnlYNj2T57dWo2rNPcaEOgv+KLT1UCMbyg5z1/mTifWMzkfgirk5VDa2U3KocVTWZ4wJHAv+KLTqjXISYt3ccnb+qK3z0jOycLuE57Zac48xoc6CP8ocPdbJmuJKbjhzAsnxMaO23pSEWJZMSec5a+4xJuRZ8EeZx4sq6Oju5fZzJ436ui+fm0NZfSt7a+1iLmNCWcCCX0QeEJFaEdnab9lPRWSniGwRkSdFJCVQ2zcnU1UeefsAhZNSOSM3efAXDNHls7MRgb9bc48xIS2QR/wPAlecsOxFYK6qzgd2A18L4PbNCTYfOEpZfSsfWjx6bfv9ZSXHs2hiqrXzGxPiAhb8qvoq0HDCshdUtdt3900gL1DbNydbvbmCMTFurpqXG7BtXDE3h+1VTRw4fCxg2zDGjIyTbfwfB/5+qgdF5G4RKRKRorq6uiCWFZnau3r4W3ElV87NYewQB2MbistmZwPw8q7agG3DGDMyjgS/iHwD6Ab+dKrnqOr9qlqoqoWZmZnBKy5CrdtZS3N7NzedFdh/sialJzIlI5F1Oy34jQlVQQ9+EfkYcA1wq1q/v6B5pqSKjLGxLJmaHvBtXTgziw1lh2nrtJm5jAlFQQ1+EbkC+CpwnapaI3CQtHX2sG5HLZfPyRnxuDz+uGhWJp3dvWwoqw/4towxQxfI7pyPABuAmSJSISJ3Ab8CkoAXReRdEbkvUNs371m/u5a2rp6AntTt7+zJaSTEunl5p52bMSYUBewsn6p+eIDFvw/U9sypPVNSTVpiLOdMTgvK9uI8bpZOy2Ddzlq+p+rX5O3GmOCxK3cjXFdPL6/srOXSM7KGPdnKcFw0M4tDR9vsKl5jQpAFf4TbtP8IzR3dXDwrO6jbvXCmtyeWdes0JvRY8Ee4l3fVEuMWlk4LfG+e/sanjGFWTpK18xsTgiz4I9wrO+tYXJBG0iiOxOmvi2ZlsbG8gab2rqBv2xhzahb8EezQ0TZ21TRz0cwsR7Z/0cwsunuV1/dYt05jQokFfwR7bbe3maWvvT3YzpyYQlKch1f3WHOPMaHEgj+CbSg7TGZSHNOyxjqy/Ri3i6XTMli/q84mZzEmhFjwRyhVZUPpYZZMSXe0H/3ymZlUNrZbt05jQogFf4Qqq2+ltrmDc6cEtzfPiZbN8DYzrd9tzT3GhAoL/gi1ofQwQFAGZTudCSljmJ411oLfmBBiwR+hNpQdJic5noL0BKdLYfmMTN7a12CjdRoTIiz4I5Cq8va+Bs6dkhYS4+Qsm+EdrfPNfYedLsUYgwV/RKo40kZdcweLCoIzKNtgzp6cRnyMi/W7rLnHmFBgwR+BNh84AsBZE1OcLcQnPsbNuVPSedXa+Y0JCRb8EWjz/iMkxLqZmZ3kdCnHLZ+RSVl9q03CbkwIsOCPQJsPHGVBXkpQh2EezPK+bp12Fa8xjgudZDCjoq2zhx1VTZw1KcXpUt5nckYi+WljrJ3fmBBgwR9htlQcpbtXOWtiqtOlvI+IsGx6Jm+U1tPZ3et0OcZENQv+CPPuwaMALMxPcbSOgSyfkcmxzh6K9jc4XYoxUc2CP8KUHGpkQsoY0sfGOV3KSc6bloHHJby624ZpNsZJFvwRZuuhRuZNGOd0GQMaG+ehsCDVhm8wxmEW/BGkqb2L8sPHmJcXmsEPsHxGFjuqmqhpane6FGOiVsCCX0QeEJFaEdnab1maiLwoInt830PrDGSY23qoEYA545MdruTU+rp12sVcxjgnkEf8DwJXnLDsXmCtqk4H1vrum1HSF/yh2tQDcEZuEplJcdbcY4yDAhb8qvoqcGL3jeuBVb7bq4AbArX9aFRyqInx4+JD8sRun75una/tqaen12blMsYJwW7jz1bVKt/taiD7VE8UkbtFpEhEiurq7OjQH1sPNTI3hI/2+yyfmUljWxfFFUedLsWYqOTYyV31TsJ6ykM+Vb1fVQtVtTAz05nJwsPJsc5uyg+3Mmd86Af/BdMycAms21HrdCnGRKVgB3+NiOQC+L7bb/4o2V3TgirMzAmdgdlOJTUxlsUFaby4vcbpUoyJSsEO/jXAHb7bdwB/DfL2I9bu6mYAZoVB8AOsmJPDrppmyutbnS7FmKgTyO6cjwAbgJkiUiEidwE/Ai4TkT3Apb77ZhTsrG4mPsZFfprzUy36Y8Vs7+kdO+o3Jvg8gVqxqn74FA9dEqhtRrNdNU3MyE7C7XJ+qkV/5KclMCsniRe2V/OJZVOcLseYqGJX7kaIXdUtITXxij8un5ND0f4jdhWvMUFmwR8BDrd0UN/SERYndvu7Zn4uqvBsSdXgTzbGjBoL/giwy3diN9yCf3p2ErNyknh6iwW/McFkwR8BdtWEZ/CD96h/0/4jVB5tc7oUY6KGBX8E2FXdTGpCDJkhPFTDqVwzfzwAa4orHa7EmOhhwR8BdlY3MzMnCZHw6NHTX0FGIosmpfLnooN4L+Y2xgSaBX+Y6+1V9tQ0MysndIdiHswti/Mpq2ulaP8Rp0sxJipY8Ie5Q0fbaO3sYUaYdeXs7+p5uSTGunls40GnSzEmKljwh7mdYdqjp7/EOA/XLRzPM1uqaG7vcrocYyKeBX+Y2+3r0TMje6zDlYzMLYsn0tbVw9+KrWunMYFmwR/mdlY3MyFlDEnxMU6XMiIL8sYxMzuJxzYecLoUYyKeBX+Y21XdFDYjcp6OiHDL4nyKKxrZUdXkdDnGRDQL/jDW1dNLWV0rMyIg+AFuPHMCsW4Xj7xtR/3GBJIFfxjbf/gY3b3K9Kzwbt/vk5oYyzXzc3li8yFaOrqdLseYiGXBH8ZK61oAmJoZGcEPcPuSSbR0dPPk5gqnSzEmYlnwh7G+4J+SmehwJaNnYX4K8/PGsWrDfruS15gAseAPY3trW8hJjg/7Hj39iQi3nzuJvbUtbCg77HQ5xkQkC/4wVlrXytSsyDna73PtgvGkJsTw0Bv7nS7FmIhkwR+mVJWy2paIat/vEx/j5kOL83lxR40N12xMAFjwh6na5g6aO7qZFiE9ek502zmT6FXl4besa6cxo82CP0yV1kZej57+8tMSuGRWNo9uPEBHd4/T5RgTUSz4w1QkduU80UeXTKK+pZO/l1Q7XYoxEcWR4BeRL4rINhHZKiKPiEi8E3WEs721LYyN85CdHH6zbvnr/GkZTMlI5KEN5U6XYkxE8Sv4ReQJEblaREb8h0JEJgCfBwpVdS7gBlaOdL3RprSulamZiWE565a/XC7htnMnsfnAUbZX2vg9xowWf4P818BHgD0i8iMRmTnC7XqAMSLiARIAm3B1iErrIrNHz4n6xu/5c5FN0mLMaPEr+FX1JVW9FTgLKAdeEpE3ROROERnS1UOqegj4D+AAUAU0quoLJz5PRO4WkSIRKaqrqxvKJiJeS0c3VY3tTI3QHj39pSbGsmJONk+9e8hO8hozSvxuuhGRdOBjwP8D3gF+gfcPwYtD2aCIpALXA5OB8UCiiNx24vNU9X5VLVTVwszMzKFsIuKVRcGJ3f4+VJjP0WNdvLS91ulSjIkI/rbxPwm8hrdZ5lpVvU5VH1PVzwFDTZ9LgX2qWqeqXcATwHlDXEdU6+vRMy0Cr9odyNJpGYwfF89j1txjzKjw94j/t6o6W1X/XVWrAEQkDkBVC4e4zQPAuSKSIN4zk5cAO4a4jqhWWtuK2yVMTIuO4He7hJsX5fHanjq7kteYUeBv8P9ggGUbhrNBVX0L+AuwGSjx1XD/cNYVrfbWtjApPYFYT/RchnHzonxUYfUmG67ZmJHynO5BEckBJuDtgXMm0Nd3MBlvs8+wqOq3gW8P9/XRLlp69PQ3MT2B86am8/imCj5z0TRcrsjtxmpMoA12yHg53h44ecDPgP/0fX0J+HpgSzMD6e7ppfxwa9QFP8AHC/M40HCMjeUNTpdiTFg77RG/qq4CVonIB1R1dZBqMqdxoOEYXT0asYOznc6K2TnEx5Twty2VnDMl3elyjAlbgzX13KaqfwQKRORLJz6uqj8LWGVmQKV1rQBMjaBZt/yVGOfh0jOyebakmm9fO4cYd/Sc4zBmNA32m9OXLmOBpAG+TJAdH5wtCo/4Aa5bMJ6G1k5e31vvdCnGhK3Bmnp+4/v+3eCUYwZTWttCVlIcyRE03eJQLJ+ZSXK8hzXFlVw4M8vpcowJS/5ewPUTEUkWkRgRWSsidQNdbWsCb28U9ujpL87j5oq5ObywrYb2LhvCwZjh8LeRdIWqNgHX4B2rZxrwlUAVZQamqpTWtkTkPLtDcd2CCbR0dPPyThvCwZjh8Df4+5qErgYeV9XGANVjTqO+pZOm9m6mRfERP8CSqelkjI1jTbEN6mrMcPgb/E+LyE5gEbBWRDKB9sCVZQaytza6T+z2cbuEa+bnsnZnLc3tXU6XY0zY8XdY5nvxDqRW6BtYrRXvCJsmiKJhukV/Xbsgl87uXtbusOYeY4bqtL16TjALb3/+/q95aJTrMadRWtdCQqyb3HE2U+WZ+ankJMfzbEkVN5w5welyjAkrfgW/iPwfMBV4F+jrSqFY8AfV3lpvj55Inm7RXy6XcMXcHB5++wAtHd2MjRvKMYwx0c3f35ZCYLaqaiCLMadXVtfK4oJUp8sIGVfPz+XBN8pZt7OW6xaMd7ocY8KGvyd3twI5gSzEnF5rRzeHjrZF5Rg9p7JoYipZSXE8u6XK6VKMCSv+HvFnANtF5G2go2+hql4XkKrMSfbV943RY8Hfx+USrpybw6MbD9La0U2iNfcY4xd/f1O+E8gizODem27Rgr+/K+flsmrDfl7eVcs18625xxh/+Nudcz3eK3ZjfLc34p1BywTJ3toW73SL6cOe/yYiLS5II2NsHH8vqXa6FGPChr9j9XwC73SJv/EtmgA8FaCazABK61qYmJZAnMftdCkhxe0SrpibzbqdtbR12tg9xvjD35O7nwGWAk0AqroHsKERg6ivK6c52VVzc2nr6uGVXXYxlzH+8Df4O1S1s++O7yIu69oZJN09vZTXH4v6wdlO5ezJaaQnxvLsVmvuMcYf/gb/ehH5Ot5J1y8DHgf+FriyTH8Hj7TR2dMb9YOznYrH7WLFnBzW7rChmo3xh7/Bfy9QB5QAnwSeBf51uBsVkRQR+YuI7BSRHSKyZLjrigalNjjboK6el8uxzh7W765zuhRjQp5f3TlVtVdEngKeUtXR+M36BfCcqt4sIrGAdVU5jb02ONugzp2SRmpCDM+WVHH5HLvW0JjTOe0Rv3h9R0TqgV3ALt/sW98a7gZFZBywDPg9gKp2qurR4a4vGpTWtpCZFMe4MdE53aI/PG4Xl8/JYe2OWmvuMWYQgzX1fBFvb57FqpqmqmnAOcBSEfniMLc5GW+z0R9E5B0R+Z2InHTWUkTuFpEiESmqq4vuf99L61qYmmkndgdz5bxcWjq6eW2PTcRuzOkMFvy3Ax9W1X19C1S1DLgN+Ogwt+kBzgL+V1XPxDu2/70nPklV71fVQlUtzMzMHOamwp+qsre2xa7Y9cN5U9MZNyaGv5fY2D3GnM5gwR+jqicdPvna+Yfb7lABVKjqW777f8H7h8AMoG+6RWvfH1yM28WK2dm8uKOGjm5r7jHmVAYL/s5hPnZKqloNHBSRmb5FlwDbh7OuaNA33aId8fvnqnm5NLd38/pea+4x5lQG69WzQESaBlguwEimgfoc8Cdfj54y4M4RrCui2XSLQ7N0WgZJ8R6eLanm4lnZTpdjTEg6bfCrakAGhlHVd/FO7mIGsbfWplsciliPi8tmZ/PCtmo6b5xHrMffS1WMiR72WxHivD16bLrFobh6Xi5N7d28UWrNPcYMxII/xJXVtVpXziE6f3oGSXEenrXePcYMyII/hNl0i8MT53Fz6exsXtheQ1dPr9PlGBNyLPhDmE23OHxXzs3h6LEuNpQedroUY0KOBX8Is66cw7dsRiaJsW7+vtWae4w5kQV/CCuts+kWhys+xs0lZ2Tz/LYauq25x5j3seAPYXtrbbrFkbhqXg4NrZ28WdbgdCnGhBQL/hDW15XTDM+FM7NIivewenOF06UYE1Is+EOUTbc4cvExbm5YOIFnS6poPNbldDnGhAwL/hDVN92iHfGPzMqz8+no7uXJd+yo35g+FvwhandNMwAzspMcriS8zRk/jnkTxvHoxoOoqtPlGBMSLPhD1B5f8FtXzpFbeXY+O6ubKa5odLoUY0KCBX+I2l3TwoSUMYyN82taZHMa1y0Yz5gYN49tPOB0KcaEBAv+ELW7ppnp2Xa0PxqS4mO4Zn4ua96tpLWj2+lyjHGcBX8I6u7ppayu1dr3R9HKs/Np7ezhr+9WOl2KMY6z4A9B+xuO0dnTy3Rr3x81Z01MZXZuMg+8vo/eXjvJa6KbBX8I2mM9ekadiHD3sinsrW3hld21TpdjjKMs+EPQ7hobnC0Qrp6fS+64eO5/tczpUoxxlAV/CNpd00xe6hgSrUfPqIpxu/j40sm8WdZAiXXtNFHMgj8E7a1tsWaeAFl5dj5JcR5++5od9ZvoZcEfYvp69FhXzsBIio9h5dn5PFNSRcWRY06XY4wjLPhDTPlhb4+eGVl2xB8ody6djAC/tbZ+E6UcC34RcYvIOyLytFM1hCLr0RN441PG8MHCPB55+yCHjrY5XY4xQefkEf89wA4Htx+Sdte0IGI9egLtsxdPB+BX6/Y4XIkxwedI8ItIHnA18Dsnth/Kdtc2k5+awJhYm3UrkCakjOEj50zkz0UVlPsmtTcmWjh1xP9fwFeBU06GKiJ3i0iRiBTV1dUFrTCn7alpZoad2A2Kf7pwKjFu4Zdr7ajfRJegB7+IXAPUquqm0z1PVe9X1UJVLczMzAxSdc7q6ullX30r0619PyiykuO5Y0kBT7576Pi5FWOigRNH/EuB60SkHHgUuFhE/uhAHSGnvL6Vrh61I/4g+uTyqSTGevjxc7ucLsWYoAl68Kvq11Q1T1ULgJXAOlW9Ldh1hKId1d6jzpnZyQ5XEj3SEmP5zEXTeGlHDa/tiZ4mRRPdrB9/CNle2USMW6xHT5B9/PwCJqUn8N2/baer55SnnYyJGI4Gv6q+oqrXOFlDKNle1cT0rCRiPfb3OJjiPG6+efVs9ta22ABuJipYwoSQHVVNnJFrzTxOuHR2NlfOzeEXa/ewz7p3mghnwR8iapvbqWvuYPZ4C36nfPe6OcR5XHz1L8X02GQtJoJZ8IeIHVXeE7uz7YjfMVnJ8Xzv+jlsLD/CfetLnS7HmICx4A8R2yubAAt+p92wcALXzM/l5y/uZmN5g9PlGBMQFvwhYkdVExNSxjAuIcbpUqKaiPDDG+eRn5bAp/+4mapGG8TNRB4L/hCx3U7shoxxY2K4//ZFtHV286k/bqa9q8fpkowZVRb8IaCts4eyuhY7sRtCpmcn8Z8fWkjxwaN886mtqNrJXhM5LPhDwK6aZnoVZufaGD2h5Iq5OXz+4mk8vqmCP7xe7nQ5xowaC/4QsKOq78TuOIcrMSf6wqUzuGx2Nj94Zjuv7rYhHUxksOAPAdsrm0iK85CXOsbpUswJXC7hv25ZyIzsJD7z8GZK61qcLsmYEbPgDwF9J3ZdLnG6FDOAxDgPv/1oITFuF59YVURjW5fTJRkzIhb8DuvtVd9QDda+H8ry0xK477ZFHGg4xtefLLGTvSasWfA77EDDMY519liPnjBw9uQ0vnjZDJ7ZUsXjRRVOl2PMsFnwO2xrZSNgJ3bDxaeWT2XJlHS+vWabtfebsGXB77Dig0eJ9biYmWNNPeHA7RJ+fstC4mNcfP6Rd+jotou7TPix4HdY8cFG5oxPtjH4w0jOuHh+/IH5bKts4n9etsHcTPixtHFQd08vJYcaWZCX4nQpZohWzMnh+oXjue+VUmvyMWHHgt9Be2pbaOvqYWF+itOlmGH416tnEx/j4hvWy8eEGQt+BxUfPArAAgv+sJSZFMe9V57Bm2UNrN58yOlyjPGbBb+DiiuOkhzvoSA9welSzDCtXJzPokmp/PCZ7TS0djpdjjF+seB30DsHjrIgPwURu2I3XLlcwg9vnEtzezc/fX6n0+UY4xcLfoc0tnWxq6aZwklpTpdiRmhWTjJ3Li3g0Y0HeefAEafLMWZQQQ9+EckXkZdFZLuIbBORe4JdQyjYvP8IqrC4INXpUswouOfSGWQlxfGtv26zidpNyHPiiL8b+LKqzgbOBT4jIrMdqMNRb5c34HEJZ0604I8EY+M8fOPq2ZQcauSRtw84XY4xpxX04FfVKlXd7LvdDOwAJgS7DqcVlTcwd8I4xsS6nS7FjJJr5+eyZEo6P31+F4dbOpwux5hTcrSNX0QKgDOBtwZ47G4RKRKRorq6yJoAo72rh+KDjZw92dr3I4mI8L3r59Da0c1PntvldDnGnJJjwS8iY4HVwBdUtenEx1X1flUtVNXCzMzM4BcYQFsqGuns6WVxgQV/pJmencRd50/msaKDbLYTvSZEORL8IhKDN/T/pKpPOFGDk94orUfETuxGqs9dMp2c5Hi++dRWunt6nS7HmJM40atHgN8DO1T1Z8Hefih4bU898/NSSEmIdboUEwBj4zx8+9rZbKts4r71NoibCT1OHPEvBW4HLhaRd31fVzlQhyMa27p49+BRlk3PcLoUE0BXzsvlmvm5/GLtHrZXntSSaYyjnOjV8w9VFVWdr6oLfV/PBrsOp2woPUxPr3LB9Mg6b2FO9v3r55KSEMtnH9lMa0e30+UYc5xduRtkr+2pIzHWzZkTU5wuxQRYamIsv1x5JuX1rXztCRvB04QOC/4gUlXW765jydR0Ytz21keDJVPT+fKKmawpruS/1+11uhxjAAv+oNpe1UTFkTYum53tdCkmiP7pwqncdNYEfvbibh61q3pNCPA4XUA0eW5rNS6BS8+w4I8mIsKPbprP4ZZO7n2ihO5e5bZzJzldlolidsQfRM9treacyemkj41zuhQTZLEeF7+5fREXz8riX5/ayk+e20mvDeZmHGLBHyR7a1vYU9vCFXNznC7FOCQ+xs19ty3iw2fn8+tXSrnjD29T3djudFkmClnwB8ma4kpEYMUca+aJZrEeF/924zx+eONcisqPsOLn63nynQrr8WOCyoI/CHp6lb8UHeSC6ZnkjhvjdDnGYSLCredM4tl7LmB6dhJffKyYu1YVsa++1enSTJSw4A+C1/bUUdnYzsrF+U6XYkLI5IxE/vzJJXzjqjN4q+wwK36+nn97dgdN7V1Ol2YinAV/EDy28SBpibHWm8ecxO0SPrFsCi9/5UJuWDiB+18t4+L/eIXHNh6wmbxMwFjwB9jBhmO8sL2GmxflEeuxt9sMLCspnp9+cAF//cxSJqYl8C+rS7j8v17l6S2V1vvHjDpLogC7b30pbhE+vnSy06WYMLAgP4XVnz6P//nIWQB89uF3uOqXr/HE5grau3ocrs5ECgv+AKpubOfxogpuLswjZ1y80+WYMCEiXD0/l+e/sIxfrFxIV08vX/pzMef9aB0/eW4n+w/bSWAzMnblbgD9Yu1uelT51LKpTpdiwpDbJVy/cALXLRjPG6WHWfVGOfetL+XXr5SyuCCVm87K4+r5uSTHxzhdqgkzFvwB8s6BIzy68SB3LZ3MxPQEp8sxYUxEWDotg6XTMqhqbOPJdw6xelMFX3uihO+s2caKOTncdNYElk3PxO0Sp8s1YUDC4cKRwsJCLSoqcroMv3V293Ljr1+nvqWDtV++kLFx9vfVjC5VZUtFI6s3V7CmuJKjx7rISx3Dx84r4JbF+STZfwEGEJFNqlp40nIL/tH3nTXbePCNcu67bZEN0WACrrO7l5d21PDg6+W8Xd7A2DgPHyrM586lBeSn2X+b0exUwW+HoqPsic0VPPhGOXedP9lC3wRFrMfFVfNyuWpeLlsqjvL7f+zjoQ3lPPjGPi6fk8Nd509m0aRUvNNdG2NH/KPqb8WV3PPoO5w7JZ1VHz/bJlsxjqlqbOOhDft5+K0DNLZ1MW/COK6al8vSaenMykm2a0qihDX1BJCq8r/rS/mP53dRWJDGg3cuJiHW/pkyzjvW2c3qTRU8uvEg23yTvrtdQsbYWDLGxpE+No44j4sYt+BxuYhx+2777sd6XHhcgsftIjHWTV5qAvlpY5iUnsi4MXYeIdRZU0+A7Ktv5Vt/3cpre+q5dsF4fvyBeRb6JmQkxHq4fUkBty8poLqxnU37j7Cjqona5nbqmjtoaO2ko7uX7l6lu6eXrh6lq8d7v6unl+4epbvXu/xEE1LGMHt8MrNzk5kzPpkzcpPJSx1jTUphwBJqmHbXNPPAP/axenMF8R43379hLredM9E+9CZk5YyL5+r5uVw9P3fIr1VVWjq6OdjQxsEjx9hX38q2yia2Vzby0o4a+hoOkuI8zMxJOv6VkxxPdnI8mUlxpCTEMCbGbb8jIcCR4BeRK4BfAG7gd6r6IyfqGIrGY11sq2rkzbIG1u6oYVtlE7EeFysXT+RzF08jK9muzDWRS0RIio9h9vgYZo9Pft9jxzq72VHVzM7qJnZWNbOrupk1xZU0v9V90npi3EJyfAzjxsSQ7PsaNyaG9MRYpmQmUpCeyOSMRManjLFrEgIo6MEvIm7gf4DLgApgo4isUdXtgd62qtKr3vHxe1Xp6VXaunpo7eimtaOH1s5umtq6qGnqoLqpnZrGdqqb2imta6HiSBsALoGF+Sl865rZXL9wvE2jaKJeQqyHRZNSWTQp9fgyVaW2uYOapnZqmjqoa+6gsa2LpvYu7/c27/fGti4ONhyjtqmd1s73xiKK9bgoSE9gUnoiealjyEtNYPy4+ON/MJLiPSTHxxAX48Ll+w9CBFwiCN4/VL3q/T1XxXfb+117Of5Yr76XC72qx9fh/fKeDxER3C7vfZcIIhxf5/HvcHy9ynuP9a1b0ePrdbsEtwhut/e7y4X3vm9bweDEEf/ZwF5VLQMQkUeB64FRD/7vP72dP721n95e6O7tZSiDHIpAemIcOePiWJCfwkfOmcjc8eOYnzeOlITY0S7VmIgiImT7mnn8oarUNXdQVt9KeX0r++pbj99+fW89xzqjY4C6vj82fX8kROA3ty/igumZo7odJ4J/AnCw3/0K4JwTnyQidwN3++62iMiuIW4nA6gfVoU+5SN58egZ8X6ECNuP0GL7ETpOuw/Lvj+idU8aaGHIntxV1fuB+4f7ehEpGqgbU7ix/Qgtth+hJRL2w4l9cOIqjkNA/zkI83zLjDHGBIETwb8RmC4ik0UkFlgJrHGgDmOMiUpBb+pR1W4R+SzwPN7unA+o6rYAbGrYzUQhxvYjtNh+hJZI2I+g70NYDNlgjDFm9NhITcYYE2Us+I0xJsqEdfCLSJqIvCgie3zfU0/xvOdE5KiIPH3C8ski8paI7BWRx3wnm4NuCPtxh+85e0Tkjn7LXxGRXSLyru8rK4i1X+Hb9l4RuXeAx+N87+1e33td0O+xr/mW7xKRy4NV80CGux8iUiAibf3e+/uCXvz76xxsP5aJyGYR6RaRm094bMDPlxNGuB89/X4ejnYc8WM/viQi20Vki4isFZFJ/R4L3M9DVcP2C/gJcK/v9r3Aj0/xvEuAa4GnT1j+Z2Cl7/Z9wKdDdT+ANKDM9z3VdzvV99grQKEDdbuBUmAKEAsUA7NPeM4/Aff5bq8EHvPdnu17fhww2bcet0Pv/0j2owDY6kTdw9yPAmA+8BBwsz+fr3DaD99jLU7/LIawHxcBCb7bn+73uQrozyOsj/jxDvWwynd7FXDDQE9S1bVAc/9l4h0U42LgL4O9Pgj82Y/LgRdVtUFVjwAvAlcEp7xTOj78hqp2An3Db/TXf9/+Alzie++vBx5V1Q5V3Qfs9a3PCSPZj1Ay6H6oarmqbgF6T3htKH2+RrIfocSf/XhZVY/57r6J97omCPDPI9yDP1tVq3y3q4HsIbw2HTiqqn1DCFbgHU7CCf7sx0BDXfSv9w++f22/GcRAGqym9z3H91434n3v/XltsIxkPwAmi8g7IrJeRC4IdLGnMZL3NNx+HqcTLyJFIvKmiNwwqpUNzVD34y7g78N87ZCE7JANfUTkJWCgyWu/0f+OqqqIhGzf1ADvx62qekhEkoDVwO14/wU2gVcFTFTVwyKyCHhKROaoapPThUWxSb7fhynAOhEpUdVSp4s6HRG5DSgElgdjeyEf/Kp66akeE5EaEclV1SoRyQVqh7Dqw0CKiHh8R3ABHTpiFPbjEHBhv/t5eNv2UdVDvu/NIvIw3n8xgxH8/gy/0fecChHxAOPwvvehNHTHsPdDvQ2yHQCquklESoEZgBNzhY7kPT3l58sBI/ps9Pt9KBORV4Az8ba1B5tf+yEil+I9AFyuqh39XnvhCa99ZbQKC/emnjVA39nuO4C/+vtC3y/sy0Bfj4AhvX6U+bMfzwMrRCTV1+tnBfC8iHhEJANARGKAa4CtQagZ/Bt+o/++3Qys8733a4CVvt4yk4HpwNtBqvtEw94PEckU7xwT+I4wp+M9EeeEkQyHMuDnK0B1DmbY++GrP853OwNYSgCGfPfToPshImcCvwGuU9X+B3yB/Xk4feZ7JF9421jXAnuAl4A03/JCvDN79T3vNaAOaMPbVna5b/kUvGGzF3gciAvx/fi4r9a9wJ2+ZYnAJmALsA3fzGZBrP0qYDfeI6pv+JZ9D+8HGSDe997u9b3XU/q99hu+1+0CrnT4szSs/QA+4Hvf3wU2A9eG+H4s9v0OtOL9z2vb6T5f4bYfwHlACd4eNCXAXSG+Hy8BNb7Pz7vAmmD8PGzIBmOMiTLh3tRjjDFmiCz4jTEmyljwG2NMlLHgN8aYKGPBb4wxUcaC35ghEJHviMg/D/KcG0Rkdr/73/NdpGNMSLDgN1HPdyXuKe8Pww14Rx8FQFW/paovjXCdxowaC34TUUTko76xzYtF5P/EO17+un7jnU/0Pe9BEblPRN4CfjLA/anincdhk4i8JiKzBtjWJ0Rko29bq0UkQUTOA64DfuobNG+qb903+15ziW9AtxIReaDfVablIvJd8Y4xX9K3PRFZLu+NLf+ObzwmY0bEgt9EDBGZA/wrcLGqLgDuAf4bWKWq84E/Ab/s95I84DxV/dIA9+8HPqeqi4B/Bn49wCafUNXFvm3twHuV6Bt4L8v/iqou1H6Dg4lIPPAgcIuqzsM7Vtan+62vXlXPAv7Xt0183z+jqguBC/BefW7MiFjwm0hyMfC4qtYDqGoDsAR42Pf4/wHn93v+46rac+J9ERmL99L/x0XkXbxjqeQOsL25vv8GSoBbgTmD1DcT2Kequ333VwHL+j3+hO/7JrwTjQC8DvxMRD4PpOh7w4gbM2whPzqnMQHUeor7LrxzNSwc5PUPAjeoarGIfIz3j6Y4HH0jM/bg+91U1R+JyDN4x3x5XUQuV9WdI9yOiXJ2xG8iyTrggyKSDt65jIE38I6KCN6j8tcGW4l6x9LfJyIf9K1HRGTBAE9NAqp8o6Le2m95s++xE+0CCkRkmu/+7cD609UiIlNVtURVf4x3tMeTzjUYM1QW/CZiqOo24IfAehEpBn4GfA64U0S24A3ae/xc3a3AXb71bOPkqRgBvgm8hbc5pv9R+KPAV3wnY6f2q68duBNvE1IJ3mkDB5uc/QsistVXfxfvzdBkzLDZ6JzGGBNl7IjfGGOijAW/McZEGQt+Y4yJMhb8xhgTZSz4jTEmyljwG2NMlLHgN8aYKPP/AVIbjEo7EjaoAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "sc.pp.neighbors(adata, n_neighbors=300)\n", + "adata.layers['X_knn'] = adata.obsp['connectivities'].dot(adata.X)\n", + "adata.layers['gene_score_knn'] = adata.obsp['connectivities'].dot(adata.obsm[\"gene_score\"])\n", + "print(adata.layers['X_knn'].shape, adata.layers['gene_score_knn'].shape)\n", + "from scipy.stats import pearsonr\n", + "cors = []\n", + "for i in range(adata.layers['X_knn'].shape[0]):\n", + " x = adata.layers['X_knn'][i,:]\n", + " y = adata.layers['gene_score_knn'][i,:].toarray().flatten()\n", + " cors.append(pearsonr(x, y))\n", + "cors_cor = list(map(lambda x: x[0], cors))\n", + "sns.kdeplot(np.array(list(cors_cor)))\n", + "plt.xlabel('correlations')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## RP-enhanced" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "before filtering (5081, 19322)\n", + "after filtering (5081, 1234)\n" + ] + } + ], + "source": [ + "# test = False does not work as tss_to_peaks need to be sub-sampled, and it seems that everything is blended.\n", + "adata = datasets.snare_p0_braincortex(test=False)\n", + "print('before filtering', adata.shape)\n", + "adata = adata[:, adata.var.index.isin(sel_genes)]\n", + "print('after filtering', adata.shape)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "... storing 'chr' as categorical\n", + "... storing 'strand' as categorical\n", + "/home/icb/ignacio.ibarra/miniconda3/envs/openproblems/lib/python3.7/site-packages/pandas/core/arrays/categorical.py:2487: FutureWarning: The `inplace` parameter in pandas.Categorical.remove_unused_categories is deprecated and will be removed in a future version.\n", + " res = method(*args, **kwargs)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "calculating weights per gene...\n", + "# weights calculated so far 500 out of 1232\n", + "# weights calculated so far 1000 out of 1232\n" + ] + } + ], + "source": [ + "adata = methods.rp_enhanced(adata, n_top_genes=2000) # log_each=10)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(0.010924862357515821, -0.010644082225420766)\n" + ] + } + ], + "source": [ + "%autoreload 2\n", + "import seaborn as sns\n", + "from openproblems.tasks.regulatory_effect_prediction import metrics\n", + "cors = metrics.spearman_correlation(adata), metrics.pearson_correlation(adata)\n", + "print(cors)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING: You’re trying to run this on 1233 dimensions of `.X`, if you really want this, set `use_rep='X'`.\n", + " Falling back to preprocessing with `sc.pp.pca` and default params.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(5081, 1233) (5081, 1233)\n" + ] + }, + { + "data": { + "text/plain": [ + "Text(0.5, 0, 'correlations')" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEICAYAAABYoZ8gAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8+yak3AAAACXBIWXMAAAsTAAALEwEAmpwYAAAsJUlEQVR4nO3deXxU9b3/8ddnZrKQkJB9gQTCjuxKUBEFV9zX2kqr1lp/te3tYpfbXtve7u29Xe5tb3t7e61trXhbl1rUUrVuoGgVlYCGsENCgJCdQDay5/P7YyYYIZDJMnNm+Twfjzwyc2bmnM+ZTN45+Z7v+X5FVTHGGBM9XE4XYIwxJrgs+I0xJspY8BtjTJSx4DfGmChjwW+MMVHGgt8YY6JMwIJfRB4QkVoR2TrAY18WERWRjEBt3xhjzMA8AVz3g8CvgIf6LxSRfGAFcMDfFWVkZGhBQcFo1maMMRFv06ZN9aqaeeLygAW/qr4qIgUDPPRz4KvAX/1dV0FBAUVFRaNVmjHGRAUR2T/Q8qC28YvI9cAhVS3247l3i0iRiBTV1dUFoTpjjIkOQQt+EUkAvg58y5/nq+r9qlqoqoWZmSf9p2KMMWaYgnnEPxWYDBSLSDmQB2wWkZwg1mCMMVEvkCd330dVS4Csvvu+8C9U1fpg1WCMMSaw3TkfATYAM0WkQkTuCtS2jDHG+C+QvXo+PMjjBYHatjHGmFOzK3eNMSbKWPAbv/T02oQ9xkSKoJ3cNeGnobWT/3xhF8+UVHH0WBeLC1L54mUzOG+qjbRhTDizI34zoK2HGrn0Z+t5dONBLp6ZxSeXT+FgQxu3/u4t/lZc6XR5xpgRsCN+c5JtlY3c+ru3GBvn4eFPnMOsnGQAvnjpDD76+7f50p/fJTMpjnOnpDtcqTFmOOyI37zPkdZO7n5oEwmxbh69+9zjoQ8QH+Pmtx8tJC81gX9ZvYWO7h4HKzXGDJcFvzlOVfny48XUNXdw322LyE9LOOk54xJi+O51c9h/+BgP/KM8+EUaY0bMgt8c99jGg6zbWcvXr5rFgvyUUz5v2YxMLj0jm1+t20NDa2fwCjTGjAoLfgNAbVM7//bsDs6dksYd5xUM+vx/uWImrZ09PPzWgKO+GmNCmAW/AeAHz+ygvbuXf79pPiIy6POnZyexfEYmqzbst7Z+Y8KMBb9h0/4jrCmu5FPLpjA5I9Hv1911/mTqmjt4urgqgNUZY0abBX+UU1V+8Mx2spLi+OTyqUN67QXTM5iWNZY/WXOPMWHFgj/KvbC9hncOHOWfV8wkMW5ol3WICB9clMfmA0cpq2sJUIXGmNFmwR/FenuVn7+4mymZidx01oRhrePGMyfgEli9uWKUqzPGBIoFfxR7bls1O6ubueeS6Xjcw/soZCXHs2xGJk9uPkSvDeRmTFiw4I9Sqsov1+5hWtZYrpk/fkTruumsPCob23m7vGGUqjPGBJIFf5Rav7uOndXNfGr5VNyuwbtvns4ls7KI9bh4flv1KFVnjAkkC/4o9Zv1ZeQkx3PdgpEd7QMkxnlYNj2T57dWo2rNPcaEOgv+KLT1UCMbyg5z1/mTifWMzkfgirk5VDa2U3KocVTWZ4wJHAv+KLTqjXISYt3ccnb+qK3z0jOycLuE57Zac48xoc6CP8ocPdbJmuJKbjhzAsnxMaO23pSEWJZMSec5a+4xJuRZ8EeZx4sq6Oju5fZzJ436ui+fm0NZfSt7a+1iLmNCWcCCX0QeEJFaEdnab9lPRWSniGwRkSdFJCVQ2zcnU1UeefsAhZNSOSM3efAXDNHls7MRgb9bc48xIS2QR/wPAlecsOxFYK6qzgd2A18L4PbNCTYfOEpZfSsfWjx6bfv9ZSXHs2hiqrXzGxPiAhb8qvoq0HDCshdUtdt3900gL1DbNydbvbmCMTFurpqXG7BtXDE3h+1VTRw4fCxg2zDGjIyTbfwfB/5+qgdF5G4RKRKRorq6uiCWFZnau3r4W3ElV87NYewQB2MbistmZwPw8q7agG3DGDMyjgS/iHwD6Ab+dKrnqOr9qlqoqoWZmZnBKy5CrdtZS3N7NzedFdh/sialJzIlI5F1Oy34jQlVQQ9+EfkYcA1wq1q/v6B5pqSKjLGxLJmaHvBtXTgziw1lh2nrtJm5jAlFQQ1+EbkC+CpwnapaI3CQtHX2sG5HLZfPyRnxuDz+uGhWJp3dvWwoqw/4towxQxfI7pyPABuAmSJSISJ3Ab8CkoAXReRdEbkvUNs371m/u5a2rp6AntTt7+zJaSTEunl5p52bMSYUBewsn6p+eIDFvw/U9sypPVNSTVpiLOdMTgvK9uI8bpZOy2Ddzlq+p+rX5O3GmOCxK3cjXFdPL6/srOXSM7KGPdnKcFw0M4tDR9vsKl5jQpAFf4TbtP8IzR3dXDwrO6jbvXCmtyeWdes0JvRY8Ee4l3fVEuMWlk4LfG+e/sanjGFWTpK18xsTgiz4I9wrO+tYXJBG0iiOxOmvi2ZlsbG8gab2rqBv2xhzahb8EezQ0TZ21TRz0cwsR7Z/0cwsunuV1/dYt05jQokFfwR7bbe3maWvvT3YzpyYQlKch1f3WHOPMaHEgj+CbSg7TGZSHNOyxjqy/Ri3i6XTMli/q84mZzEmhFjwRyhVZUPpYZZMSXe0H/3ymZlUNrZbt05jQogFf4Qqq2+ltrmDc6cEtzfPiZbN8DYzrd9tzT3GhAoL/gi1ofQwQFAGZTudCSljmJ411oLfmBBiwR+hNpQdJic5noL0BKdLYfmMTN7a12CjdRoTIiz4I5Cq8va+Bs6dkhYS4+Qsm+EdrfPNfYedLsUYgwV/RKo40kZdcweLCoIzKNtgzp6cRnyMi/W7rLnHmFBgwR+BNh84AsBZE1OcLcQnPsbNuVPSedXa+Y0JCRb8EWjz/iMkxLqZmZ3kdCnHLZ+RSVl9q03CbkwIsOCPQJsPHGVBXkpQh2EezPK+bp12Fa8xjgudZDCjoq2zhx1VTZw1KcXpUt5nckYi+WljrJ3fmBBgwR9htlQcpbtXOWtiqtOlvI+IsGx6Jm+U1tPZ3et0OcZENQv+CPPuwaMALMxPcbSOgSyfkcmxzh6K9jc4XYoxUc2CP8KUHGpkQsoY0sfGOV3KSc6bloHHJby624ZpNsZJFvwRZuuhRuZNGOd0GQMaG+ehsCDVhm8wxmEW/BGkqb2L8sPHmJcXmsEPsHxGFjuqmqhpane6FGOiVsCCX0QeEJFaEdnab1maiLwoInt830PrDGSY23qoEYA545MdruTU+rp12sVcxjgnkEf8DwJXnLDsXmCtqk4H1vrum1HSF/yh2tQDcEZuEplJcdbcY4yDAhb8qvoqcGL3jeuBVb7bq4AbArX9aFRyqInx4+JD8sRun75una/tqaen12blMsYJwW7jz1bVKt/taiD7VE8UkbtFpEhEiurq7OjQH1sPNTI3hI/2+yyfmUljWxfFFUedLsWYqOTYyV31TsJ6ykM+Vb1fVQtVtTAz05nJwsPJsc5uyg+3Mmd86Af/BdMycAms21HrdCnGRKVgB3+NiOQC+L7bb/4o2V3TgirMzAmdgdlOJTUxlsUFaby4vcbpUoyJSsEO/jXAHb7bdwB/DfL2I9bu6mYAZoVB8AOsmJPDrppmyutbnS7FmKgTyO6cjwAbgJkiUiEidwE/Ai4TkT3Apb77ZhTsrG4mPsZFfprzUy36Y8Vs7+kdO+o3Jvg8gVqxqn74FA9dEqhtRrNdNU3MyE7C7XJ+qkV/5KclMCsniRe2V/OJZVOcLseYqGJX7kaIXdUtITXxij8un5ND0f4jdhWvMUFmwR8BDrd0UN/SERYndvu7Zn4uqvBsSdXgTzbGjBoL/giwy3diN9yCf3p2ErNyknh6iwW/McFkwR8BdtWEZ/CD96h/0/4jVB5tc7oUY6KGBX8E2FXdTGpCDJkhPFTDqVwzfzwAa4orHa7EmOhhwR8BdlY3MzMnCZHw6NHTX0FGIosmpfLnooN4L+Y2xgSaBX+Y6+1V9tQ0MysndIdiHswti/Mpq2ulaP8Rp0sxJipY8Ie5Q0fbaO3sYUaYdeXs7+p5uSTGunls40GnSzEmKljwh7mdYdqjp7/EOA/XLRzPM1uqaG7vcrocYyKeBX+Y2+3r0TMje6zDlYzMLYsn0tbVw9+KrWunMYFmwR/mdlY3MyFlDEnxMU6XMiIL8sYxMzuJxzYecLoUYyKeBX+Y21XdFDYjcp6OiHDL4nyKKxrZUdXkdDnGRDQL/jDW1dNLWV0rMyIg+AFuPHMCsW4Xj7xtR/3GBJIFfxjbf/gY3b3K9Kzwbt/vk5oYyzXzc3li8yFaOrqdLseYiGXBH8ZK61oAmJoZGcEPcPuSSbR0dPPk5gqnSzEmYlnwh7G+4J+SmehwJaNnYX4K8/PGsWrDfruS15gAseAPY3trW8hJjg/7Hj39iQi3nzuJvbUtbCg77HQ5xkQkC/4wVlrXytSsyDna73PtgvGkJsTw0Bv7nS7FmIhkwR+mVJWy2paIat/vEx/j5kOL83lxR40N12xMAFjwh6na5g6aO7qZFiE9ek502zmT6FXl4besa6cxo82CP0yV1kZej57+8tMSuGRWNo9uPEBHd4/T5RgTUSz4w1QkduU80UeXTKK+pZO/l1Q7XYoxEcWR4BeRL4rINhHZKiKPiEi8E3WEs721LYyN85CdHH6zbvnr/GkZTMlI5KEN5U6XYkxE8Sv4ReQJEblaREb8h0JEJgCfBwpVdS7gBlaOdL3RprSulamZiWE565a/XC7htnMnsfnAUbZX2vg9xowWf4P818BHgD0i8iMRmTnC7XqAMSLiARIAm3B1iErrIrNHz4n6xu/5c5FN0mLMaPEr+FX1JVW9FTgLKAdeEpE3ROROERnS1UOqegj4D+AAUAU0quoLJz5PRO4WkSIRKaqrqxvKJiJeS0c3VY3tTI3QHj39pSbGsmJONk+9e8hO8hozSvxuuhGRdOBjwP8D3gF+gfcPwYtD2aCIpALXA5OB8UCiiNx24vNU9X5VLVTVwszMzKFsIuKVRcGJ3f4+VJjP0WNdvLS91ulSjIkI/rbxPwm8hrdZ5lpVvU5VH1PVzwFDTZ9LgX2qWqeqXcATwHlDXEdU6+vRMy0Cr9odyNJpGYwfF89j1txjzKjw94j/t6o6W1X/XVWrAEQkDkBVC4e4zQPAuSKSIN4zk5cAO4a4jqhWWtuK2yVMTIuO4He7hJsX5fHanjq7kteYUeBv8P9ggGUbhrNBVX0L+AuwGSjx1XD/cNYVrfbWtjApPYFYT/RchnHzonxUYfUmG67ZmJHynO5BEckBJuDtgXMm0Nd3MBlvs8+wqOq3gW8P9/XRLlp69PQ3MT2B86am8/imCj5z0TRcrsjtxmpMoA12yHg53h44ecDPgP/0fX0J+HpgSzMD6e7ppfxwa9QFP8AHC/M40HCMjeUNTpdiTFg77RG/qq4CVonIB1R1dZBqMqdxoOEYXT0asYOznc6K2TnEx5Twty2VnDMl3elyjAlbgzX13KaqfwQKRORLJz6uqj8LWGVmQKV1rQBMjaBZt/yVGOfh0jOyebakmm9fO4cYd/Sc4zBmNA32m9OXLmOBpAG+TJAdH5wtCo/4Aa5bMJ6G1k5e31vvdCnGhK3Bmnp+4/v+3eCUYwZTWttCVlIcyRE03eJQLJ+ZSXK8hzXFlVw4M8vpcowJS/5ewPUTEUkWkRgRWSsidQNdbWsCb28U9ujpL87j5oq5ObywrYb2LhvCwZjh8LeRdIWqNgHX4B2rZxrwlUAVZQamqpTWtkTkPLtDcd2CCbR0dPPyThvCwZjh8Df4+5qErgYeV9XGANVjTqO+pZOm9m6mRfERP8CSqelkjI1jTbEN6mrMcPgb/E+LyE5gEbBWRDKB9sCVZQaytza6T+z2cbuEa+bnsnZnLc3tXU6XY0zY8XdY5nvxDqRW6BtYrRXvCJsmiKJhukV/Xbsgl87uXtbusOYeY4bqtL16TjALb3/+/q95aJTrMadRWtdCQqyb3HE2U+WZ+ankJMfzbEkVN5w5welyjAkrfgW/iPwfMBV4F+jrSqFY8AfV3lpvj55Inm7RXy6XcMXcHB5++wAtHd2MjRvKMYwx0c3f35ZCYLaqaiCLMadXVtfK4oJUp8sIGVfPz+XBN8pZt7OW6xaMd7ocY8KGvyd3twI5gSzEnF5rRzeHjrZF5Rg9p7JoYipZSXE8u6XK6VKMCSv+HvFnANtF5G2go2+hql4XkKrMSfbV943RY8Hfx+USrpybw6MbD9La0U2iNfcY4xd/f1O+E8gizODem27Rgr+/K+flsmrDfl7eVcs18625xxh/+Nudcz3eK3ZjfLc34p1BywTJ3toW73SL6cOe/yYiLS5II2NsHH8vqXa6FGPChr9j9XwC73SJv/EtmgA8FaCazABK61qYmJZAnMftdCkhxe0SrpibzbqdtbR12tg9xvjD35O7nwGWAk0AqroHsKERg6ivK6c52VVzc2nr6uGVXXYxlzH+8Df4O1S1s++O7yIu69oZJN09vZTXH4v6wdlO5ezJaaQnxvLsVmvuMcYf/gb/ehH5Ot5J1y8DHgf+FriyTH8Hj7TR2dMb9YOznYrH7WLFnBzW7rChmo3xh7/Bfy9QB5QAnwSeBf51uBsVkRQR+YuI7BSRHSKyZLjrigalNjjboK6el8uxzh7W765zuhRjQp5f3TlVtVdEngKeUtXR+M36BfCcqt4sIrGAdVU5jb02ONugzp2SRmpCDM+WVHH5HLvW0JjTOe0Rv3h9R0TqgV3ALt/sW98a7gZFZBywDPg9gKp2qurR4a4vGpTWtpCZFMe4MdE53aI/PG4Xl8/JYe2OWmvuMWYQgzX1fBFvb57FqpqmqmnAOcBSEfniMLc5GW+z0R9E5B0R+Z2InHTWUkTuFpEiESmqq4vuf99L61qYmmkndgdz5bxcWjq6eW2PTcRuzOkMFvy3Ax9W1X19C1S1DLgN+Ogwt+kBzgL+V1XPxDu2/70nPklV71fVQlUtzMzMHOamwp+qsre2xa7Y9cN5U9MZNyaGv5fY2D3GnM5gwR+jqicdPvna+Yfb7lABVKjqW777f8H7h8AMoG+6RWvfH1yM28WK2dm8uKOGjm5r7jHmVAYL/s5hPnZKqloNHBSRmb5FlwDbh7OuaNA33aId8fvnqnm5NLd38/pea+4x5lQG69WzQESaBlguwEimgfoc8Cdfj54y4M4RrCui2XSLQ7N0WgZJ8R6eLanm4lnZTpdjTEg6bfCrakAGhlHVd/FO7mIGsbfWplsciliPi8tmZ/PCtmo6b5xHrMffS1WMiR72WxHivD16bLrFobh6Xi5N7d28UWrNPcYMxII/xJXVtVpXziE6f3oGSXEenrXePcYMyII/hNl0i8MT53Fz6exsXtheQ1dPr9PlGBNyLPhDmE23OHxXzs3h6LEuNpQedroUY0KOBX8Is66cw7dsRiaJsW7+vtWae4w5kQV/CCuts+kWhys+xs0lZ2Tz/LYauq25x5j3seAPYXtrbbrFkbhqXg4NrZ28WdbgdCnGhBQL/hDW15XTDM+FM7NIivewenOF06UYE1Is+EOUTbc4cvExbm5YOIFnS6poPNbldDnGhAwL/hDVN92iHfGPzMqz8+no7uXJd+yo35g+FvwhandNMwAzspMcriS8zRk/jnkTxvHoxoOoqtPlGBMSLPhD1B5f8FtXzpFbeXY+O6ubKa5odLoUY0KCBX+I2l3TwoSUMYyN82taZHMa1y0Yz5gYN49tPOB0KcaEBAv+ELW7ppnp2Xa0PxqS4mO4Zn4ua96tpLWj2+lyjHGcBX8I6u7ppayu1dr3R9HKs/Np7ezhr+9WOl2KMY6z4A9B+xuO0dnTy3Rr3x81Z01MZXZuMg+8vo/eXjvJa6KbBX8I2mM9ekadiHD3sinsrW3hld21TpdjjKMs+EPQ7hobnC0Qrp6fS+64eO5/tczpUoxxlAV/CNpd00xe6hgSrUfPqIpxu/j40sm8WdZAiXXtNFHMgj8E7a1tsWaeAFl5dj5JcR5++5od9ZvoZcEfYvp69FhXzsBIio9h5dn5PFNSRcWRY06XY4wjLPhDTPlhb4+eGVl2xB8ody6djAC/tbZ+E6UcC34RcYvIOyLytFM1hCLr0RN441PG8MHCPB55+yCHjrY5XY4xQefkEf89wA4Htx+Sdte0IGI9egLtsxdPB+BX6/Y4XIkxwedI8ItIHnA18Dsnth/Kdtc2k5+awJhYm3UrkCakjOEj50zkz0UVlPsmtTcmWjh1xP9fwFeBU06GKiJ3i0iRiBTV1dUFrTCn7alpZoad2A2Kf7pwKjFu4Zdr7ajfRJegB7+IXAPUquqm0z1PVe9X1UJVLczMzAxSdc7q6ullX30r0619PyiykuO5Y0kBT7576Pi5FWOigRNH/EuB60SkHHgUuFhE/uhAHSGnvL6Vrh61I/4g+uTyqSTGevjxc7ucLsWYoAl68Kvq11Q1T1ULgJXAOlW9Ldh1hKId1d6jzpnZyQ5XEj3SEmP5zEXTeGlHDa/tiZ4mRRPdrB9/CNle2USMW6xHT5B9/PwCJqUn8N2/baer55SnnYyJGI4Gv6q+oqrXOFlDKNle1cT0rCRiPfb3OJjiPG6+efVs9ta22ABuJipYwoSQHVVNnJFrzTxOuHR2NlfOzeEXa/ewz7p3mghnwR8iapvbqWvuYPZ4C36nfPe6OcR5XHz1L8X02GQtJoJZ8IeIHVXeE7uz7YjfMVnJ8Xzv+jlsLD/CfetLnS7HmICx4A8R2yubAAt+p92wcALXzM/l5y/uZmN5g9PlGBMQFvwhYkdVExNSxjAuIcbpUqKaiPDDG+eRn5bAp/+4mapGG8TNRB4L/hCx3U7shoxxY2K4//ZFtHV286k/bqa9q8fpkowZVRb8IaCts4eyuhY7sRtCpmcn8Z8fWkjxwaN886mtqNrJXhM5LPhDwK6aZnoVZufaGD2h5Iq5OXz+4mk8vqmCP7xe7nQ5xowaC/4QsKOq78TuOIcrMSf6wqUzuGx2Nj94Zjuv7rYhHUxksOAPAdsrm0iK85CXOsbpUswJXC7hv25ZyIzsJD7z8GZK61qcLsmYEbPgDwF9J3ZdLnG6FDOAxDgPv/1oITFuF59YVURjW5fTJRkzIhb8DuvtVd9QDda+H8ry0xK477ZFHGg4xtefLLGTvSasWfA77EDDMY519liPnjBw9uQ0vnjZDJ7ZUsXjRRVOl2PMsFnwO2xrZSNgJ3bDxaeWT2XJlHS+vWabtfebsGXB77Dig0eJ9biYmWNNPeHA7RJ+fstC4mNcfP6Rd+jotou7TPix4HdY8cFG5oxPtjH4w0jOuHh+/IH5bKts4n9etsHcTPixtHFQd08vJYcaWZCX4nQpZohWzMnh+oXjue+VUmvyMWHHgt9Be2pbaOvqYWF+itOlmGH416tnEx/j4hvWy8eEGQt+BxUfPArAAgv+sJSZFMe9V57Bm2UNrN58yOlyjPGbBb+DiiuOkhzvoSA9welSzDCtXJzPokmp/PCZ7TS0djpdjjF+seB30DsHjrIgPwURu2I3XLlcwg9vnEtzezc/fX6n0+UY4xcLfoc0tnWxq6aZwklpTpdiRmhWTjJ3Li3g0Y0HeefAEafLMWZQQQ9+EckXkZdFZLuIbBORe4JdQyjYvP8IqrC4INXpUswouOfSGWQlxfGtv26zidpNyHPiiL8b+LKqzgbOBT4jIrMdqMNRb5c34HEJZ0604I8EY+M8fOPq2ZQcauSRtw84XY4xpxX04FfVKlXd7LvdDOwAJgS7DqcVlTcwd8I4xsS6nS7FjJJr5+eyZEo6P31+F4dbOpwux5hTcrSNX0QKgDOBtwZ47G4RKRKRorq6yJoAo72rh+KDjZw92dr3I4mI8L3r59Da0c1PntvldDnGnJJjwS8iY4HVwBdUtenEx1X1flUtVNXCzMzM4BcYQFsqGuns6WVxgQV/pJmencRd50/msaKDbLYTvSZEORL8IhKDN/T/pKpPOFGDk94orUfETuxGqs9dMp2c5Hi++dRWunt6nS7HmJM40atHgN8DO1T1Z8Hefih4bU898/NSSEmIdboUEwBj4zx8+9rZbKts4r71NoibCT1OHPEvBW4HLhaRd31fVzlQhyMa27p49+BRlk3PcLoUE0BXzsvlmvm5/GLtHrZXntSSaYyjnOjV8w9VFVWdr6oLfV/PBrsOp2woPUxPr3LB9Mg6b2FO9v3r55KSEMtnH9lMa0e30+UYc5xduRtkr+2pIzHWzZkTU5wuxQRYamIsv1x5JuX1rXztCRvB04QOC/4gUlXW765jydR0Ytz21keDJVPT+fKKmawpruS/1+11uhxjAAv+oNpe1UTFkTYum53tdCkmiP7pwqncdNYEfvbibh61q3pNCPA4XUA0eW5rNS6BS8+w4I8mIsKPbprP4ZZO7n2ihO5e5bZzJzldlolidsQfRM9treacyemkj41zuhQTZLEeF7+5fREXz8riX5/ayk+e20mvDeZmHGLBHyR7a1vYU9vCFXNznC7FOCQ+xs19ty3iw2fn8+tXSrnjD29T3djudFkmClnwB8ma4kpEYMUca+aJZrEeF/924zx+eONcisqPsOLn63nynQrr8WOCyoI/CHp6lb8UHeSC6ZnkjhvjdDnGYSLCredM4tl7LmB6dhJffKyYu1YVsa++1enSTJSw4A+C1/bUUdnYzsrF+U6XYkLI5IxE/vzJJXzjqjN4q+wwK36+nn97dgdN7V1Ol2YinAV/EDy28SBpibHWm8ecxO0SPrFsCi9/5UJuWDiB+18t4+L/eIXHNh6wmbxMwFjwB9jBhmO8sL2GmxflEeuxt9sMLCspnp9+cAF//cxSJqYl8C+rS7j8v17l6S2V1vvHjDpLogC7b30pbhE+vnSy06WYMLAgP4XVnz6P//nIWQB89uF3uOqXr/HE5grau3ocrs5ECgv+AKpubOfxogpuLswjZ1y80+WYMCEiXD0/l+e/sIxfrFxIV08vX/pzMef9aB0/eW4n+w/bSWAzMnblbgD9Yu1uelT51LKpTpdiwpDbJVy/cALXLRjPG6WHWfVGOfetL+XXr5SyuCCVm87K4+r5uSTHxzhdqgkzFvwB8s6BIzy68SB3LZ3MxPQEp8sxYUxEWDotg6XTMqhqbOPJdw6xelMFX3uihO+s2caKOTncdNYElk3PxO0Sp8s1YUDC4cKRwsJCLSoqcroMv3V293Ljr1+nvqWDtV++kLFx9vfVjC5VZUtFI6s3V7CmuJKjx7rISx3Dx84r4JbF+STZfwEGEJFNqlp40nIL/tH3nTXbePCNcu67bZEN0WACrrO7l5d21PDg6+W8Xd7A2DgPHyrM586lBeSn2X+b0exUwW+HoqPsic0VPPhGOXedP9lC3wRFrMfFVfNyuWpeLlsqjvL7f+zjoQ3lPPjGPi6fk8Nd509m0aRUvNNdG2NH/KPqb8WV3PPoO5w7JZ1VHz/bJlsxjqlqbOOhDft5+K0DNLZ1MW/COK6al8vSaenMykm2a0qihDX1BJCq8r/rS/mP53dRWJDGg3cuJiHW/pkyzjvW2c3qTRU8uvEg23yTvrtdQsbYWDLGxpE+No44j4sYt+BxuYhx+2777sd6XHhcgsftIjHWTV5qAvlpY5iUnsi4MXYeIdRZU0+A7Ktv5Vt/3cpre+q5dsF4fvyBeRb6JmQkxHq4fUkBty8poLqxnU37j7Cjqona5nbqmjtoaO2ko7uX7l6lu6eXrh6lq8d7v6unl+4epbvXu/xEE1LGMHt8MrNzk5kzPpkzcpPJSx1jTUphwBJqmHbXNPPAP/axenMF8R43379hLredM9E+9CZk5YyL5+r5uVw9P3fIr1VVWjq6OdjQxsEjx9hX38q2yia2Vzby0o4a+hoOkuI8zMxJOv6VkxxPdnI8mUlxpCTEMCbGbb8jIcCR4BeRK4BfAG7gd6r6IyfqGIrGY11sq2rkzbIG1u6oYVtlE7EeFysXT+RzF08jK9muzDWRS0RIio9h9vgYZo9Pft9jxzq72VHVzM7qJnZWNbOrupk1xZU0v9V90npi3EJyfAzjxsSQ7PsaNyaG9MRYpmQmUpCeyOSMRManjLFrEgIo6MEvIm7gf4DLgApgo4isUdXtgd62qtKr3vHxe1Xp6VXaunpo7eimtaOH1s5umtq6qGnqoLqpnZrGdqqb2imta6HiSBsALoGF+Sl865rZXL9wvE2jaKJeQqyHRZNSWTQp9fgyVaW2uYOapnZqmjqoa+6gsa2LpvYu7/c27/fGti4ONhyjtqmd1s73xiKK9bgoSE9gUnoiealjyEtNYPy4+ON/MJLiPSTHxxAX48Ll+w9CBFwiCN4/VL3q/T1XxXfb+117Of5Yr76XC72qx9fh/fKeDxER3C7vfZcIIhxf5/HvcHy9ynuP9a1b0ePrdbsEtwhut/e7y4X3vm9bweDEEf/ZwF5VLQMQkUeB64FRD/7vP72dP721n95e6O7tZSiDHIpAemIcOePiWJCfwkfOmcjc8eOYnzeOlITY0S7VmIgiImT7mnn8oarUNXdQVt9KeX0r++pbj99+fW89xzqjY4C6vj82fX8kROA3ty/igumZo7odJ4J/AnCw3/0K4JwTnyQidwN3++62iMiuIW4nA6gfVoU+5SN58egZ8X6ECNuP0GL7ETpOuw/Lvj+idU8aaGHIntxV1fuB+4f7ehEpGqgbU7ix/Qgtth+hJRL2w4l9cOIqjkNA/zkI83zLjDHGBIETwb8RmC4ik0UkFlgJrHGgDmOMiUpBb+pR1W4R+SzwPN7unA+o6rYAbGrYzUQhxvYjtNh+hJZI2I+g70NYDNlgjDFm9NhITcYYE2Us+I0xJsqEdfCLSJqIvCgie3zfU0/xvOdE5KiIPH3C8ski8paI7BWRx3wnm4NuCPtxh+85e0Tkjn7LXxGRXSLyru8rK4i1X+Hb9l4RuXeAx+N87+1e33td0O+xr/mW7xKRy4NV80CGux8iUiAibf3e+/uCXvz76xxsP5aJyGYR6RaRm094bMDPlxNGuB89/X4ejnYc8WM/viQi20Vki4isFZFJ/R4L3M9DVcP2C/gJcK/v9r3Aj0/xvEuAa4GnT1j+Z2Cl7/Z9wKdDdT+ANKDM9z3VdzvV99grQKEDdbuBUmAKEAsUA7NPeM4/Aff5bq8EHvPdnu17fhww2bcet0Pv/0j2owDY6kTdw9yPAmA+8BBwsz+fr3DaD99jLU7/LIawHxcBCb7bn+73uQrozyOsj/jxDvWwynd7FXDDQE9S1bVAc/9l4h0U42LgL4O9Pgj82Y/LgRdVtUFVjwAvAlcEp7xTOj78hqp2An3Db/TXf9/+Alzie++vBx5V1Q5V3Qfs9a3PCSPZj1Ay6H6oarmqbgF6T3htKH2+RrIfocSf/XhZVY/57r6J97omCPDPI9yDP1tVq3y3q4HsIbw2HTiqqn1DCFbgHU7CCf7sx0BDXfSv9w++f22/GcRAGqym9z3H91434n3v/XltsIxkPwAmi8g7IrJeRC4IdLGnMZL3NNx+HqcTLyJFIvKmiNwwqpUNzVD34y7g78N87ZCE7JANfUTkJWCgyWu/0f+OqqqIhGzf1ADvx62qekhEkoDVwO14/wU2gVcFTFTVwyKyCHhKROaoapPThUWxSb7fhynAOhEpUdVSp4s6HRG5DSgElgdjeyEf/Kp66akeE5EaEclV1SoRyQVqh7Dqw0CKiHh8R3ABHTpiFPbjEHBhv/t5eNv2UdVDvu/NIvIw3n8xgxH8/gy/0fecChHxAOPwvvehNHTHsPdDvQ2yHQCquklESoEZgBNzhY7kPT3l58sBI/ps9Pt9KBORV4Az8ba1B5tf+yEil+I9AFyuqh39XnvhCa99ZbQKC/emnjVA39nuO4C/+vtC3y/sy0Bfj4AhvX6U+bMfzwMrRCTV1+tnBfC8iHhEJANARGKAa4CtQagZ/Bt+o/++3Qys8733a4CVvt4yk4HpwNtBqvtEw94PEckU7xwT+I4wp+M9EeeEkQyHMuDnK0B1DmbY++GrP853OwNYSgCGfPfToPshImcCvwGuU9X+B3yB/Xk4feZ7JF9421jXAnuAl4A03/JCvDN79T3vNaAOaMPbVna5b/kUvGGzF3gciAvx/fi4r9a9wJ2+ZYnAJmALsA3fzGZBrP0qYDfeI6pv+JZ9D+8HGSDe997u9b3XU/q99hu+1+0CrnT4szSs/QA+4Hvf3wU2A9eG+H4s9v0OtOL9z2vb6T5f4bYfwHlACd4eNCXAXSG+Hy8BNb7Pz7vAmmD8PGzIBmOMiTLh3tRjjDFmiCz4jTEmyljwG2NMlLHgN8aYKGPBb4wxUcaC35ghEJHviMg/D/KcG0Rkdr/73/NdpGNMSLDgN1HPdyXuKe8Pww14Rx8FQFW/paovjXCdxowaC34TUUTko76xzYtF5P/EO17+un7jnU/0Pe9BEblPRN4CfjLA/anincdhk4i8JiKzBtjWJ0Rko29bq0UkQUTOA64DfuobNG+qb903+15ziW9AtxIReaDfVablIvJd8Y4xX9K3PRFZLu+NLf+ObzwmY0bEgt9EDBGZA/wrcLGqLgDuAf4bWKWq84E/Ab/s95I84DxV/dIA9+8HPqeqi4B/Bn49wCafUNXFvm3twHuV6Bt4L8v/iqou1H6Dg4lIPPAgcIuqzsM7Vtan+62vXlXPAv7Xt0183z+jqguBC/BefW7MiFjwm0hyMfC4qtYDqGoDsAR42Pf4/wHn93v+46rac+J9ERmL99L/x0XkXbxjqeQOsL25vv8GSoBbgTmD1DcT2Kequ333VwHL+j3+hO/7JrwTjQC8DvxMRD4PpOh7w4gbM2whPzqnMQHUeor7LrxzNSwc5PUPAjeoarGIfIz3j6Y4HH0jM/bg+91U1R+JyDN4x3x5XUQuV9WdI9yOiXJ2xG8iyTrggyKSDt65jIE38I6KCN6j8tcGW4l6x9LfJyIf9K1HRGTBAE9NAqp8o6Le2m95s++xE+0CCkRkmu/+7cD609UiIlNVtURVf4x3tMeTzjUYM1QW/CZiqOo24IfAehEpBn4GfA64U0S24A3ae/xc3a3AXb71bOPkqRgBvgm8hbc5pv9R+KPAV3wnY6f2q68duBNvE1IJ3mkDB5uc/QsistVXfxfvzdBkzLDZ6JzGGBNl7IjfGGOijAW/McZEGQt+Y4yJMhb8xhgTZSz4jTEmyljwG2NMlLHgN8aYKPP/AVIbjEo7EjaoAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "sc.pp.neighbors(adata, n_neighbors=300)\n", + "adata.layers['X_knn'] = adata.obsp['connectivities'].dot(adata.X)\n", + "adata.layers['gene_score_knn'] = adata.obsp['connectivities'].dot(adata.obsm[\"gene_score\"])\n", + "print(adata.layers['X_knn'].shape, adata.layers['gene_score_knn'].shape)\n", + "from scipy.stats import pearsonr\n", + "cors = []\n", + "for i in range(adata.layers['X_knn'].shape[0]):\n", + " x = adata.layers['X_knn'][i,:]\n", + " y = adata.layers['gene_score_knn'][i,:].toarray().flatten()\n", + " cors.append(pearsonr(x, y))\n", + "cors_cor = list(map(lambda x: x[0], cors))\n", + "sns.kdeplot(np.array(list(cors_cor)))\n", + "plt.xlabel('correlations')" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "openproblems", + "language": "python", + "name": "openproblems" + }, + "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.7.10" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}