diff --git a/runs.ipynb b/runs.ipynb index 51afd6a06..8483ab88b 100644 --- a/runs.ipynb +++ b/runs.ipynb @@ -9,7 +9,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -17,6 +17,7 @@ "import pandas as pd\n", "import anndata as ad \n", "import numpy as np\n", + "import scanpy as sc \n", "from src.exp_analysis.helper import plot_cumulative_density\n", "def extract_data(data, reg='reg1', dataset_id='scgen_pearson'):\n", " i = 0\n", @@ -59,53 +60,45 @@ ] }, { - "cell_type": "code", - "execution_count": 10, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "from sklearn.preprocessing import StandardScaler \n", - "adata = ad.read_h5ad('resources/grn-benchmark/multiomics_rna.h5ad')" + "# Add to process data pipeline" ] }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 85, "metadata": {}, "outputs": [], "source": [ - "X_tran = StandardScaler().fit_transform(adata.X.todense().A)" + "# process multiomics\n", + "# !pip install --user scikit-misc\n", + "from sklearn.preprocessing import StandardScaler \n", + "adata = ad.read_h5ad('resources/grn-benchmark/multiomics_rna.h5ad')\n", + "def antoine_filtering(adata):\n", + " threshold = 0.1\n", + " mask = adata.X!=0\n", + " mask_obs = (np.sum(mask, axis=1).A.flatten()/mask.shape[1])>threshold\n", + " mask_var = (np.sum(mask, axis=0).A.flatten()/mask.shape[0])>threshold\n", + " adata.obs['high_coverage'] = mask_obs\n", + " adata.var['high_coverage'] = mask_var\n", + "antoine_filtering(adata)\n", + "\n", + "# raw hvgs\n", + "var = sc.pp.highly_variable_genes(adata, flavor='seurat_v3', n_top_genes=7000, inplace=False)\n", + "adata.var['hvg_counts'] = var.highly_variable\n" ] }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 5, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "(
,\n", - " )" - ] - }, - "execution_count": 27, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFzCAYAAADSc9khAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCCklEQVR4nO3deXxTVd4/8E+SJmlLV+hKKbRlK/tOLSqLtiw6jDCjIvoAg4rjUn9KwcEqUNFnLOOCuKCMKCIzo6jzqONIRWppkaWy75ssLWXpSkv3Nmlyfn+kjcaWkpSkN7n5vF8vXpCbc2++h7YfLueee65CCCFARERuQSl1AURE1HEY+kREboShT0TkRhj6RERuhKFPRORGGPpERG6EoU9E5EYY+kREbsRD6gI6mtFoxOXLl+Hr6wuFQiF1OUREN0wIgaqqKnTt2hVKZdvn8m4X+pcvX0ZkZKTUZRAR2d2FCxfQrVu3Ntu4Xej7+voCMP3l+Pn5SVxN++j1emzevBkTJ06EWq2WupwbJrf+APLrk9z6A8irT5WVlYiMjDTnW1vcLvSbh3T8/PxcOvS9vb3h5+fn8t+sgPz6A8ivT3LrDyDPPlkzZM0LuUREboShT0TkRhj6RERuhKFPRORGGPpERG6EoU9E5EYY+kREbkTS0P/xxx8xdepUdO3aFQqFAl9//fV198nOzsbw4cOh1WrRq1cvrFu3zuF1EhHJhaShX1NTgyFDhmDVqlVWtc/NzcWdd96JCRMm4ODBg3j66afx8MMP4/vvv3dwpURE8iDpHblTpkzBlClTrG6/evVqREdH4/XXXwcA9OvXD9u3b8cbb7yBSZMmOapMInIRQggIAQgARiFgbH4tTK+btwsj0KDXoVoPXKlugFA2wmAUMBgFGpt/NwjoDUbU6gxQqxQQ5s9o+pymPwOAgOlDLdrgl1rMtZkao+lP5tpE0/vbz5RiVFRn3NIrCIGdNA75O3KpZRhycnKQkJBgsW3SpEl4+umnr7lPQ0MDGhoazK8rKysBmG7B1uv1DqnT0ZrrdtX6f0tu/QFcq09Go0B9owH1eiNqdI2o1xlRrWtERZ0eSoUCeoMRDbpG7C9V4OKPZxHgrYXeYITOIHChvA7eahWUSphDs/mXzmDEz0XViOriDYMRMIhf3jOKX9qcLKxCr2AfGIUpcI1G0dQWMBiNOF1cg0BvNbzUKotjNAd0rc4ApQIWIWw9Dzy/d6sD/lbbb33OeXz+yGgMiwyweh9bvs9cKvQLCwsRGhpqsS00NBSVlZWoq6uDl5dXi33S0tKwbNmyFts3b94Mb29vh9XaETIyMqQuwa7k1h/Afn0yGIGaRqBKD5TrTAnXKIC6RqC2EVAogLpGBUrrTa/VSqDBCDQYFGg0mtoW15naeakAowAMAtAZbVleXAWcPmtz7QcuVFy3ze6a8jbfL6/VoxzXDjajzWHfOq1SQKkAlApA1fR7pc70e4AWaP7bMv+usHxt/rPC9LvF9mu0/e32YE+Bi4d2ouCI9XXX1tZa3dalQr89UlJSkJycbH7dvBrdxIkTXXrBtYyMDCQmJspioSi59Qewrk+1ukYUVjTg4tU65JfV4mRhFYqrGqBSKFBc3YDDFysR6K1Gea39/rcghOkfj2tRqxTQGwS6+nuiuqERSoUC3Tt7QaUAqiorENylM85dqcWQbv7QeiihVilRUtWAMH9P+Gg9oFIq4KFUQKVUQNWUZuW1OnQL9IKy6T2lwvS+6XfT59bpDAjy0Zq2KxVQKQCVUgmVEk3/2xDw9fSAR/Mxfv05TccyHdu06JhCYdrP/Bq/ea0ADI2NyPzhB0yc6Prfd80jGNZwqdAPCwtDUVGRxbaioiL4+fm1epYPAFqtFlqttsV2tVrt8l9oOfTh1+TUHyEEyhuA3fmVOHypCo0GI0qqG3ClWoeiqgZcKq9FabXuusf5beB37qSBp4cS5bV6DIrwh8ZDCV2jEVAAfUN9EeCthlKhgFqlQNcAL3hrPOClUUGjUkLjYQphX081PJQKqFWm0FYpFfDWqOCpVkGlbP3MX6/XIz09HXfcMUo2XyO9XgmFQh7fd7bU71KhHx8fj/T0dIttGRkZiI+Pl6gicmdCCBRW1mN3bhkOXahAbmk1Omk9UFRZjzPF1Siv9QD277vucTQqJYb3CEDnThoYjUCvEB/Ehvsi1M8Tnh4qBHir0UnrAX8v9TVDmchakoZ+dXU1zpw5Y36dm5uLgwcPonPnzujevTtSUlJw6dIlrF+/HgDw6KOP4p133sFf/vIXPPjgg9iyZQs+//xzbNy4UaoukBtoNBiRW1qDzJPFKK/R4djlSmw/U3rd/RQQiA7qhB5dOqFWZ0B8zy7o3EmDEF9PdAv0QmRnb/h5evCxndShJA39vXv3YsKECebXzWPvc+bMwbp161BQUID8/Hzz+9HR0di4cSPmz5+PN998E926dcMHH3zA6ZpkN1eqG3D4YgWOXKpA5sliXCirRVW9HnrDta8U9grxgYdSge6dvdEt0BuDu/kjMkCLcwd2YNrUW1x+6IDkRdLQHz9+PEQbc6xau9t2/PjxOHDggAOrInchhMD+/KvYm1eG3bllOHjhKq7UtD7OrlYp0DfMFxV1esyJj0KvEB+E+nmiZ7APNB4t73HU6/W4eNjRPSCynUuN6RPdqMKKemSeLMLOM1ew93wZiiobWrSJ6uKNQd0C0CvYB0G+GsTHdEFUl05QcjydZIChT7JmNArsySvDN4cuY3duGc6WVFvM6fZSq3BL7yAM6x6AUVGd0S/cDz5a/liQfPG7m2Snedhm2+kSfLH3Ii5drbN4f1j3AIzvE2IOei+NSqJKiToeQ59k49jlCnxz6DIyjhXhXGmNeXsnjQoTYkMQF9MFwyIDMDDCX8IqiaTF0CeXVlRZj4935mHHmVIcuvjL7f4aDyWiu3TCvLEx+N3gcHiqeTZPBDD0yQUZjAI7zpRi4+ECfHXwkumO1CaxYb54bHxP3BYbAl9PTpUk+i2GPrmMilo91ufk4Z+7zlvMuhkU4Y/E/qG4e0Q3dA1ofTkOIjJh6JPTK6lqQOo3R/HDiWLzWb2vpwfuGBiOPwyPwOjozryrlchKDH1yWuU1Ovz9x3NYn5OHWp0BABDiq8WzU2Jx5+BwaD04Tk9kK4Y+OZ1aXSP++dN5vLPlDCrrTesAD4zww6LJsbilVxDP6oluAEOfnEZJVQPW7zqLf+TkoabpzD42zBcLJvZFQr8Qhj2RHTD0SXJV9Xr8+5wSC3b9iMam22UjO3vhsXG9cO/IbvBQtVzbhojah6FPkhFCYMOeC3hl00mU1yoBCAzvHoBHxvbExP6hXOuGyAEY+iSJPXlleO7LIzhdXA0A6KwVWDFzJMbHhklcGZG8MfSpQ5VUNeDV70/i870XAQAqpQLzb++F8MoTuLlnF4mrI5I/hj51CF2jEWt35OLdrF9m5PxxeDcsmtIXgZ4qpKefkLhCIvfA0CeHK6yox5Of7seevHIApgd4L7trAG6KMZ3Z6/X6tnYnIjti6JPDCCHw6e4LePHbY6jXG+GpVmLR5Fj8z009oOaMHCJJMPTJIU4VVuGFb44h59wVAEDvEB+smT0SUUGdJK6MyL0x9Mnusk4W4/F/7Ued3gCNSokFE/tg3q0xnIJJ5AQY+mQ3Qgh8uD0Xf00/ASGAMT27YPkfBqN7F2+pSyOiJgx9souGRgMW/fswvj54GQBwz4hueGnaQD68hMjJMPTphv1cVIUFnx/CkUumJ1ctmhyLR8fFcK0cIifE0Kcb8s2hy3jmi0NoaDTCz9MDb80chvF9Q6Qui4iugaFP7aI3GPG3707ig+25AIBRUYF45/7hCPXzlLgyImoLQ59sVlWvx//79ACyTpUAAB66JRqLJsdC48G590TOjqFPNjl/pQZz1u5G3pVaqJQKvHr3YPxheDepyyIiKzH0yWpHL1XgTx/tRmm1DkE+Wqy6fxjiYrhIGpErYeiTVXaeKcWfPtoDncGIvqG++PjB0Qjz5/g9kath6NN1bT9digc/NgX+qKhAfDB7FPy91VKXRUTtwNCnNm39uQTz1u+FrtGI22JD8O4Dw3nDFZELY+jTNW0/XYpHmgJ/XJ9gvPc/w6H1YOATuTLOsaNW7Ttfjj//Yy8aGo24tXcQ3p89goFPJAMMfWrhVGEV5n60GzU6A0b0CMT7s0Yy8IlkgqFPFipq9Uj6ZD8q6xsxvHsAPpo7Cl4aBj6RXDD0yayh0YBH/7kPp4urEeSjwTv3D4efJ2fpEMkJQ58AAEajwDNfHEbOuSvwVCux/sE4dA3wkrosIrIzhj4BAN7JOoNvDl2GUgG898AI9O/qJ3VJROQADH3Cjz+XYEXGzwCAF+8aiAmxXBqZSK4Y+m7uaq0Oi/7vMADgvlGReCCuu8QVEZEjMfTd3PLvTqKgoh4RAV5YOrU/n3ZFJHMMfTd2/HIlvth3EQDw2j1D4K3hDdpEcsfQd1O6RiMW/d9hGIwCkwaEIr4nl0gmcgcMfTe1PicPRy5VwM/TAy/8foDU5RBRB2Hou6GrtTqsyjoDAHh2Sj+E+3M+PpG7kDz0V61ahaioKHh6eiIuLg67d+9us/3KlSvRt29feHl5ITIyEvPnz0d9fX0HVSsPq7eeQ3mtHn1CfXDvSD7qkMidSBr6n332GZKTk5Gamor9+/djyJAhmDRpEoqLi1tt/8knn+DZZ59FamoqTpw4gQ8//BCfffYZnnvuuQ6u3HVdKKvF+pw8AMDCiX3hoZL8330i6kCS/sSvWLEC8+bNw9y5c9G/f3+sXr0a3t7eWLt2bavtd+7ciZtvvhn3338/oqKiMHHiRMycOfO6/zugX7ycfgK1OgNGR3VGQr9Qqcshog4m2Rw9nU6Hffv2ISUlxbxNqVQiISEBOTk5re4zZswY/POf/8Tu3bsxevRonDt3Dunp6Zg1a9Y1P6ehoQENDQ3m15WVlQAAvV4PvV5vp950rOa6ba3/wIWr+O5oIZQKYOmdfWEwNMJgcESFtmlvf5yZ3Pokt/4A8uqTLX2QLPRLS0thMBgQGmp5thkaGoqTJ0+2us/999+P0tJS3HLLLRBCoLGxEY8++mibwztpaWlYtmxZi+2bN2+Gt7f3jXVCYhkZGVa3FQJ4+5gKgAIjg4w4u38bzjqutHaxpT+uQm59klt/AHn0qba21uq2LnU3TnZ2Nl5++WW8++67iIuLw5kzZ/DUU0/hpZdewpIlS1rdJyUlBcnJyebXlZWViIyMxMSJE+Hn55qLiun1emRkZCAxMRFqtXVLH+87X46zP+2BxkOJV2ff6lQraLanP85Obn2SW38AefWpeQTDGpKFflBQEFQqFYqKiiy2FxUVISwsrNV9lixZglmzZuHhhx8GAAwaNAg1NTV45JFH8Pzzz0OpbHmJQqvVQqvVttiuVqtd/gttSx/W77oAAJg2tCt6BDvnP3Zy+Jr8ltz6JLf+APLoky31S3YhV6PRYMSIEcjMzDRvMxqNyMzMRHx8fKv71NbWtgh2lcr0VCchhOOKdXFnS6qx6WghAODBW6IlroaIpCTp8E5ycjLmzJmDkSNHYvTo0Vi5ciVqamowd+5cAMDs2bMRERGBtLQ0AMDUqVOxYsUKDBs2zDy8s2TJEkydOtUc/tTSqqwzMArg9tgQxIY551k+EXUMSUN/xowZKCkpwdKlS1FYWIihQ4di06ZN5ou7+fn5Fmf2ixcvhkKhwOLFi3Hp0iUEBwdj6tSp+Otf/ypVF5zehbJa/OfgZQBA0m29JK6GiKQm+YXcpKQkJCUltfpedna2xWsPDw+kpqYiNTW1AyqTh7U7cmEwCozp2QXDugdKXQ4RSYy3Y8pYRa0en+8xXcB98GaO5RMRQ1/Wvj1yGTU6A3oGd8JtfAQiEYGhL1tCCKzfeR4AcPeISCiVfCIWETH0ZevAhas4VVQFT7US9/O5t0TUhKEvU5/sygcA3DEwHP5ern3jCRHZD0Nfhirr9fj2sGmaJs/yiejXGPoy9N9Dl1GvN6J3iA9G9OA0TSL6BUNfhr4+cAkA8Ifh3aBQ8AIuEf2CoS8z56/UYE9eOZQK4K6hXaUuh4icDENfZr5qOssf0zPIqZZPJiLnYHPoR0VF4cUXX0R+fr4j6qEbIITAfw+ZLuBOGxYhcTVE5IxsDv2nn34aX375JWJiYpCYmIgNGzZYPI6QpHO2pBpnS2qgVikwcQCff0tELbUr9A8ePIjdu3ejX79+ePLJJxEeHo6kpCTs37/fETWSlf57qAAAcHOvIPh5cm4+EbXU7jH94cOH46233sLly5eRmpqKDz74AKNGjcLQoUOxdu1aPtREAt80De1M59AOEV1Du5dW1uv1+Oqrr/DRRx8hIyMDN910Ex566CFcvHgRzz33HH744Qd88skn9qyV2nCupBq5paahHS6uRkTXYnPo79+/Hx999BE+/fRTKJVKzJ49G2+88QZiY2PNbaZPn45Ro0bZtVBqW9apEgBAXHQX+HJoh4iuwebQHzVqFBITE/Hee+9h2rRprT6QNzo6Gvfdd59dCiTrZJ8qBgDc0jtI4kqIyJnZHPrnzp1Djx492mzTqVMnfPTRR+0uimxTWt2AnWevAAAm9uesHSK6Npsv5E6YMAFXrlxpsf3q1auIiYmxS1Fkm01HC2EwCgzu5o+YYB+pyyEiJ2Zz6Ofl5cFgMLTY3tDQgEuXLtmlKLJNxvEiAMCkAWESV0JEzs7q4Z1vvvnG/Ofvv/8e/v7+5tcGgwGZmZmIioqya3F0ffV6A3KahnYm8YYsIroOq0N/2rRpAACFQoE5c+ZYvKdWqxEVFYXXX3/drsXR9e3JK4POYESQjxY9ObRDRNdhdegbjUYAppk5e/bsQVAQZ4k4g61NUzXH9gniMspEdF02z97Jzc11RB3UDkIIbG4az0/ox6EdIro+q0L/rbfewiOPPAJPT0+89dZbbbb9f//v/9mlMLq+k4XVyC+rhdZDiXF9gqUuh4hcgFWh/8Ybb+CBBx6Ap6cn3njjjWu2UygUDP0OlHmy6YasXkHopG33ihpE5EasSopfD+lweMd5ZJwwhT6XUSYia93wk7MMBgMOHjyI8vJye9RDViqpA44XVMFDqcDtHM8nIiu1az39Dz/8EIAp8MeOHYvhw4cjMjIS2dnZ9q6PruFkhWmmzogegQjy0UpcDRG5CptD/9///jeGDBkCAPjvf/+LvLw8nDx5EvPnz8fzzz9v9wKpdaebQv/mXpw6S0TWszn0S0tLERZmut0/PT0d99xzD/r06YMHH3wQR44csXuB1FKjwYifm0Kfq2oSkS1sDv3Q0FAcP34cBoMBmzZtQmJiIgCgtrYWKpXK7gVSS/vyr6LOoECgtxpDugVIXQ4RuRCb5/nNnTsX9957L8LDw6FQKJCQkAAA2LVrl8WDVMhxdp4tAwCM6dkFKiXvwiUi69kc+i+88AIGDhyICxcu4J577oFWa7qIqFKp8Oyzz9q9QGqp+SlZt/bqInElRORq2nVHz913391i228XYSPHKKlqwInCKgDA+L68C5eIbNOu0M/MzERmZiaKi4vNC7E1W7t2rV0Ko9Y1PxaxWyeBLp00EldDRK7G5tBftmwZXnzxRYwcOdI8rk8dZ+vPpqGd/gFC4kqIyBXZHPqrV6/GunXrMGvWLEfUQ20wGoX5gSl9A4zXaU1E1JLNUzZ1Oh3GjBnjiFroOo5cqsCVGh06aVWI4vNSiKgdbA79hx9+GJ988okjaqHraB7aGRPTBR43vGoSEbkjm4d36uvr8f777+OHH37A4MGDoVarLd5fsWKF3YojS9tPlwIAxvYOAkouS1wNEbkim0P/8OHDGDp0KADg6NGjFu/xoq7j1OoaceCCaSXT+JjOOFYicUFE5JJsDv2srCxH1EHXsTu3DHqDQLi/J7p39sIxqQsiIpfU7pHhM2fO4Pvvv0ddXR0A0/NayXF2nGke2gnm/6iIqN1sDv0rV67g9ttvR58+fXDHHXegoKAAAPDQQw9hwYIFdi+QTLY1jeffzFU1iegG2Bz68+fPh1qtRn5+Pry9vc3bZ8yYgU2bNtlcwKpVqxAVFQVPT0/ExcVh9+7dbba/evUqnnjiCYSHh0Or1aJPnz5IT0+3+XNdSWl1A04VmZZeiI/hejtE1H42j+lv3rwZ33//Pbp162axvXfv3jh//rxNx/rss8+QnJyM1atXIy4uDitXrsSkSZNw6tQphISEtGiv0+mQmJiIkJAQ/Pvf/0ZERATOnz+PgIAAW7vhUnaevQIhgH7hfgj21UKv10tdEhG5KJtDv6amxuIMv1lZWZl5xU1rrVixAvPmzcPcuXMBmO723bhxI9auXdvqip1r165FWVkZdu7caZ4qGhUVZWsXXM5P50x34d4U01niSojI1dkc+rfeeivWr1+Pl156CYBpmqbRaMQrr7yCCRMmWH0cnU6Hffv2ISUlxbxNqVQiISEBOTk5re7zzTffID4+Hk888QT+85//IDg4GPfffz8WLVp0zQe4NDQ0oKGhwfy6srISAKDX613ijFkIgR+bFlmLiwqwqNsV6reG3PoDyK9PcusPIK8+2dIHm0P/lVdewe233469e/dCp9PhL3/5C44dO4aysjLs2LHD6uOUlpbCYDAgNDTUYntoaChOnjzZ6j7nzp3Dli1b8MADDyA9PR1nzpzB448/Dr1ej9TU1Fb3SUtLw7Jly1ps37x5c6v/Y3E2RXXAxase8FAIVJ7ei/Rzv7yXkZEhXWEOILf+APLrk9z6A8ijT7W1tVa3tTn0Bw4ciJ9//hnvvPMOfH19UV1djT/84Q/mi6uOZDQaERISgvfffx8qlQojRozApUuX8Oqrr14z9FNSUpCcnGx+XVlZicjISEycOBF+fn4Ordce/rkrHzh4EsN7dMb0qaMAmP5Vz8jIQGJiYos7ol2R3PoDyK9PcusPIK8+NY9gWKNd6+n7+/vj+eefb8+uZkFBQVCpVCgqKrLYXlRUZH7w+m+Fh4dDrVZbDOX069cPhYWF0Ol00Ghari+v1WpbvdagVqtd4gudedI0VXNCbGiLel2lD9aSW38A+fVJbv0B5NEnW+q3OfRPnz6N//znP8jLy4NCoUBMTAymTZuG6Ohom46j0WgwYsQIZGZmYtq0aQBMZ/KZmZlISkpqdZ+bb74Zn3zyCYxGI5RK02zTn3/+GeHh4a0Gvqur0xmwO8/0PNzE/i1nMxER2cqmefppaWno378/Fi1ahP/7v//DF198gYULF6Jv37547bXXbP7w5ORkrFmzBh9//DFOnDiBxx57DDU1NebZPLNnz7a40PvYY4+hrKwMTz31FH7++Wds3LgRL7/8Mp544gmbP9sV7DxbCl2jEREBXugZzLWUiejGWX2mn5WVhcWLF2PJkiV46qmnEBgYCMA0VXPlypV49tlnMXr0aIwdO9bqD58xYwZKSkqwdOlSFBYWYujQodi0aZP54m5+fr75jB4AIiMj8f3332P+/PkYPHgwIiIi8NRTT2HRokVWf6Yr2XjYdLdzYv9QLr1ARHZhdeivXr0aDz/8MF544QWL7Z07d8aLL76IwsJCvPfeezaFPgAkJSVdczgnOzu7xbb4+Hj89NNPNn2GKxJCmOfnJ/QLvU5rIiLrWD28s3v37jYfkThr1iy3COOOkl9Wi8sV9VCrFBjeI0DqcohIJqwO/aKiojbvfo2OjkZhYaE9aiIA+86b1s7v39Uf3pp2TbIiImrB6tCvr69vc4aMWq2GTqezS1H0y1OyuMAaEdmTTaeQH3zwAXx8Wp9FUlVVZZeCyDSev71p/fxbuZQyEdmR1aHfvXt3rFmz5rpt6MadLalBcVUDNB5KjOgRKHU5RCQjVod+Xl6eA8ugX8tpmrUzonsgPNWtLyRHRNQe7X5cIjnO/qaLuKOjuZQyEdkXQ98J7Wo+0+fQDhHZGUPfyVwsN83P91AqMDKKoU9E9sXQdzI5Z01n+QMiOD+fiOyPoe9kdjRN1RzLqZpE5ADtCv2zZ89i8eLFmDlzJoqLTY/y++6773Ds2DG7FueODly4CgAYzvF8InIAm0N/69atGDRoEHbt2oUvv/wS1dXVAIBDhw5d8+lVZJ0LZbU4f6XWNJ7P0CciB7A59J999ln87//+LzIyMiyWZbjtttu44NoN2nve9MCUgRH+8PV07Sf5EJFzsjn0jxw5gunTp7fYHhISgtLSUrsU5a62nzZdxI3j/HwichCbQz8gIAAFBQUtth84cAARERF2KcodCSHMF3Fv7sWLuETkGDaH/n333YdFixahsLAQCoUCRqMRO3bswMKFCzF79mxH1OgW8stqUVhpWj9/VBTP9InIMWwO/ZdffhmxsbGIjIxEdXU1+vfvj7Fjx2LMmDFYvHixI2p0Cz82LaU8rHsgvDRcb4eIHMPmu380Gg3WrFmDJUuW4OjRo6iursawYcPQu3dvR9TnNvbmmS7icv18InIkm0N/+/btuOWWW9C9e3cupWxHhzg/n4g6gM3DO7fddhuio6Px3HPP4fjx446oye2U1+iQd6UWADA4wl/iaohIzmwO/cuXL2PBggXYunUrBg4ciKFDh+LVV1/FxYsXHVGfW9jbtJRyz+BOCOx07UdSEhHdKJtDPygoCElJSdixYwfOnj2Le+65Bx9//DGioqJw2223OaJG2Wu+KWt4dw7tEJFj3dCCa9HR0Xj22WexfPlyDBo0CFu3brVXXW7lwPmrAIA4XsQlIgdrd+jv2LEDjz/+OMLDw3H//fdj4MCB2Lhxoz1rcwsGo8CJgkoAwICufhJXQ0RyZ/PsnZSUFGzYsAGXL19GYmIi3nzzTdx1113w9vZ2RH2yd7KwElUNjfDReqB3iI/U5RCRzNkc+j/++COeeeYZ3HvvvQgK4nIBN2pf00XcoZEB8FDx8QZE5Fg2h/6OHTscUYfb+qnpebg3xXDpBSJyPKtC/5tvvsGUKVOgVqvxzTfftNn297//vV0KcxcH868CAEb0YOgTkeNZFfrTpk1DYWEhQkJCMG3atGu2UygUMBgM9qpN9sprdLhcUQ8AGBDBi7hE5HhWhb7RaGz1z3RjDjYtvRAT3Al+fGgKEXUAm68crl+/Hg0NDS2263Q6rF+/3i5FuYuThVUAgAFdufQCEXUMm0N/7ty5qKioaLG9qqoKc+fOtUtR7uJ0kSn0+3CqJhF1EJtDXwgBhULRYvvFixfh788zVls0n+n3CfOVuBIichdWT9kcNmwYFAoFFAoFbr/9dnh4/LKrwWBAbm4uJk+e7JAi5UjXaMSZ4moAQL8wXsQloo5hdeg3z9o5ePAgJk2aBB+fX4YkNBoNoqKi8Mc//tHuBcpVflkNdAYjOmlUiOzsJXU5ROQmrA791NRUAEBUVBRmzJgBT09PhxXlDn4uMp3l9wr1bXW4jIjIEWy+I3fOnDmOqMPtnG0a2ukZ3EniSojIndgc+gaDAW+88QY+//xz5OfnQ6fTWbxfVlZmt+Lk7GxJ05k+Z+4QUQeyefbOsmXLsGLFCsyYMQMVFRVITk7GH/7wByiVSrzwwgsOKFGecpsejxjdhWf6RNRxbA79f/3rX1izZg0WLFgADw8PzJw5Ex988AGWLl2Kn376yRE1yo4QAmea5ujzTJ+IOpLNoV9YWIhBgwYBAHx8fMw3av3ud7/jQ1SsVFTZgBqdASqlAlFBPNMnoo5jc+h369YNBQUFAICePXti8+bNAIA9e/ZAq9XatzqZah7P797ZG2quoU9EHcjmxJk+fToyMzMBAE8++SSWLFmC3r17Y/bs2XjwwQftXqAcNS+/0DOYQztE1LFsnr2zfPly859nzJiB7t27IycnB71798bUqVPtWpxcNS+/EMvlF4iog93w2EJ8fDySk5NvKPBXrVqFqKgoeHp6Ii4uDrt377Zqvw0bNkChULS5xr8zOldSAwDoHcozfSLqWFY/Octatj4567PPPkNycjJWr16NuLg4rFy5EpMmTcKpU6cQEhJyzf3y8vKwcOFC3HrrrTZ9njM4XczhHSKShtVPzrJGe56ctWLFCsybN8+8LPPq1auxceNGrF27Fs8++2yr+xgMBjzwwANYtmwZtm3bhqtXr9r0mVKqrNejvFYPAIjmzB0i6mA2PznLnnQ6Hfbt24eUlBTzNqVSiYSEBOTk5FxzvxdffBEhISF46KGHsG3btjY/o6GhweKhL5WVlQAAvV4PvV5/gz2w3dlC0+d36aSBRinaVUPzPlLU7why6w8gvz7JrT+AvPpkSx9svpBrT6WlpTAYDAgNDbXYHhoaipMnT7a6z/bt2/Hhhx/i4MGDVn1GWloali1b1mL75s2b4e3tbXPNN+rAFQUAFXwVDUhPT7+hY2VkZNinKCcht/4A8uuT3PoDyKNPtbW1Vre1OfRffPHFNt9funSprYe0WlVVFWbNmoU1a9YgKCjIqn1SUlKQnJxsfl1ZWYnIyEhMnDgRfn4dv479+a3ngJ/PYGjPrrjjjkHtOoZer0dGRgYSExOhVrv+s3Xl1h9Afn2SW38AefWpeQTDGjaH/ldffWXxWq/XIzc3Fx4eHujZs6dNoR8UFASVSoWioiKL7UVFRQgLC2vR/uzZs8jLy7OYKdQ89OTh4YFTp06hZ8+eFvtotdpWbxpTq9WSfKHPl9UDAHqG+N7w50vVB0eRW38A+fVJbv0B5NEnW+q3OfQPHDjQYltlZSX+9Kc/Yfr06TYdS6PRYMSIEcjMzDRfLDYajcjMzERSUlKL9rGxsThy5IjFtsWLF6OqqgpvvvkmIiMjbfp8KeRdMU3X5PILRCQFu4zp+/n5YdmyZZg6dSpmzZpl077JycmYM2cORo4cidGjR2PlypWoqakxz+aZPXs2IiIikJaWBk9PTwwcONBi/4CAAABosd1Z5ZeZxt6iuLomEUnAbhdyKyoqzIuv2WLGjBkoKSnB0qVLUVhYiKFDh2LTpk3mi7v5+flQKuWxPk293oCSKtNMoohAPiKRiDqezaH/1ltvWbwWQqCgoAD/+Mc/MGXKlHYVkZSU1OpwDgBkZ2e3ue+6deva9ZlSuHS1DgDgrVEh0Nu1xxCJyDXZHPpvvPGGxWulUong4GDMmTPHYr49tZTf9OCU7p29+VxcIpKEzaGfm5vriDrcwrnSpou4HM8nIonIY7DcReQ3zdyJ5sPQiUgiNp/p19fX4+2330ZWVhaKi4tbLNGwf/9+uxUnNwUVpjn6Xf09Ja6EiNyVzaH/0EMPYfPmzbj77rsxevRojk3boLDSFPqhfgx9IpKGzaH/7bffIj09HTfffLMj6pG1y02zd7oGcLomEUnD5jH9iIgI+PryiU+2qmloRGm1DgAQ2bnjF3ojIgLaEfqvv/46Fi1ahPPnzzuiHtkqqDCd5ftqPeDvxTn6RCQNm4d3Ro4cifr6esTExMDb27vFQj9lZWV2K05OLl1tuojLoR0ikpDNoT9z5kxcunQJL7/8MkJDQ3kh10oXmtbciezM0Cci6dgc+jt37kROTg6GDBniiHpk60K5KfS7BXI8n4ikY/OYfmxsLOrq6hxRi6xdbhre6caF1ohIQjaH/vLly7FgwQJkZ2fjypUrqKystPhFrSuq4Bx9IpKezcM7kydPBgDcfvvtFtuFEFAoFDAYDPapTGZ4YxYROQObQz8rK8sRdciaEALFVc2h3/LRjUREHcXm0B83bpwj6pC1yvpG1OtNaxSF+PJMn4ikY3Po//jjj22+P3bs2HYXI1fFTUM7fp4e8NKoJK6GiNyZzaE/fvz4Ftt+PVefY/otNa+uGcbVNYlIYjbP3ikvL7f4VVxcjE2bNmHUqFHYvHmzI2p0eZyjT0TOwuYzfX9//xbbEhMTodFokJycjH379tmlMDkpMC/BwDN9IpKW3Z6cFRoailOnTtnrcLLSPHMnjNM1iUhiNp/pHz582OK1EAIFBQVYvnw5hg4daq+6ZOVK05LKXXw4XZOIpGVz6A8dOhQKhQJCCIvtN910E9auXWu3wuSktLoBANClk0biSojI3dkc+rm5uRavlUolgoOD4enJoYtrKakyhX6wL8/0iUhaNod+jx49HFGHbBmNAsVNoc8pm0QkNasv5G7ZsgX9+/dvdVG1iooKDBgwANu2bbNrcXJwpUaHRqOAQgEEcUyfiCRmdeivXLkS8+bNg5+fX4v3/P398ec//xkrVqywa3FyUNh0Y1aQjxZqld0mSxERtYvVKXTo0CHzCputmThxIufot6L52bhdObRDRE7A6tAvKipq8TzcX/Pw8EBJSYldipKToqbxfC6pTETOwOrQj4iIwNGjR6/5/uHDhxEeHm6XouSkeeZOEGfuEJETsDr077jjDixZsgT19fUt3qurq0Nqaip+97vf2bU4OSi4ahreCeeZPhE5AaunbC5evBhffvkl+vTpg6SkJPTt2xcAcPLkSaxatQoGgwHPP/+8wwp1VSVNN2aF8OEpROQErA790NBQ7Ny5E4899hhSUlLMd+QqFApMmjQJq1atQmhoqMMKdVVFlc2hzzN9IpKeTTdn9ejRA+np6SgvL8eZM2cghEDv3r0RGBjoqPpcXvMDVEL5xCwicgI235ELAIGBgRg1apS9a5GdRoMRZbWmxdY4vENEzoB3CznQlRodhACUCiDQm4utEZH0GPoO1Dxds4uPFiql4jqtiYgcj6HvQFdqmtbR55LKROQkGPoOxCWVicjZMPQdqKhp5k4IZ+4QkZNg6DtQ85k+Z+4QkbNg6DvQL2f6DH0icg4MfQdqDn2usElEzoKh70DNSzAw9InIWTD0HcT0bFzTmT6fjUtEzsIpQn/VqlWIioqCp6cn4uLisHv37mu2XbNmDW699VYEBgYiMDAQCQkJbbaXypUaHfQG07NxOaZPRM5C8tD/7LPPkJycjNTUVOzfvx9DhgzBpEmTUFxc3Gr77OxszJw5E1lZWcjJyUFkZCQmTpyIS5cudXDlbWt+TGIwn41LRE5E8jRasWIF5s2bh7lz56J///5YvXo1vL29sXbt2lbb/+tf/8Ljjz+OoUOHIjY2Fh988AGMRiMyMzM7uPK2NT8QPZxDO0TkRNq1yqa96HQ67Nu3DykpKeZtSqUSCQkJyMnJseoYtbW10Ov16Ny5c6vvNzQ0oKGhwfy6srISAKDX66HX62+g+rZdLq8BYBrasffnNB/PkfV3JLn1B5Bfn+TWH0BefbKlD5KGfmlpKQwGQ4uHr4SGhuLkyZNWHWPRokXo2rUrEhISWn0/LS0Ny5Yta7F98+bN8Pb2tr1oK+3MVwJQor6sEOnp6Q75jIyMDIccVypy6w8gvz7JrT+APPpUW1trdVtJQ/9GLV++HBs2bEB2djY8PVsfRklJSUFycrL5dWVlpfk6gJ+fn8Nq2/bVMeDSJYwY2Ad3jI+x67H1ej0yMjKQmJgItVpt12NLQW79AeTXJ7n1B5BXn5pHMKwhaegHBQVBpVKhqKjIYntRURHCwsLa3Pe1117D8uXL8cMPP2Dw4MHXbKfVaqHVtpw9o1arHfqFLmpagqFroLfDPsfRfehocusPIL8+ya0/gDz6ZEv9kl7I1Wg0GDFihMVF2OaLsvHx8dfc75VXXsFLL72ETZs2YeTIkR1Rqs2a78YN441ZROREJB/eSU5Oxpw5czBy5EiMHj0aK1euRE1NDebOnQsAmD17NiIiIpCWlgYA+Nvf/oalS5fik08+QVRUFAoLCwEAPj4+8PHxkawfv9U8e4c3ZhGRM5E89GfMmIGSkhIsXboUhYWFGDp0KDZt2mS+uJufnw+l8pf/kLz33nvQ6XS4++67LY6TmpqKF154oSNLv6aqej0q6xsBcMomETkXyUMfAJKSkpCUlNTqe9nZ2Rav8/LyHF/QDWpec8fX0wO+nq49VkhE8iL5zVlyZH5ilg+XXyAi58LQd4Dmhdb4mEQicjYMfQe4UGa6UaJboONu/iIiag+GvgNcLDctthbZ2UviSoiILDH0HeBy03TNrgEMfSJyLgx9B7h81XSmH8HQJyInw9B3AN6YRUTOiqFvZ5X1elQ38MYsInJODH07u9R0ETfAWw1vjVPc+0ZEZMbQt7Pm0O8WyPF8InI+DH07u1DeNEc/gHP0icj5MPTt7CLP9InIiTH07aygomm6JkOfiJwQQ9/OmqdrcuYOETkjhr6dXb7aPEefZ/pE5HwY+nbU0GhAUdMKmxzTJyJnxNC3o8KKeggBaD2U6NJJI3U5REQtMPTt6NfLLygUComrISJqiaFvR5eaFlrryvF8InJSDH074t24ROTsGPp2dLlpjn44l1QmIifF0Lej3NIaAEBUFy7BQETOiaFvRxfKTGf63Tsz9InIOTH07aRebzAP70QFdZK4GiKi1jH07SS3tAZCAL6eHpyjT0ROi6FvJ2eKqwEAPYN9OEefiJwWQ99OThdVAQBiw3wlroSI6NoY+nZyuulMv1eIj8SVEBFdG0PfTo4XVAIA+oX7SVwJEdG1MfTtQNdoNN+NGxPMmTtE5LwY+nZwtqQajUYBP08PhPnx4SlE5LwY+nbwc9NF3F4hnLlDRM6NoW8HJwpMod+/K8fzici5MfTt4NjlCgBAbBhDn4icG0P/BgkhcOyyaebOoAh/iashImobQ/8GXSirQ1mNDhqVErHhvDGLiJwbQ/8GHS8wDe30DvWB1kMlcTVERG1j6N+ggxdMoc+hHSJyBQz9G3Tk0lUAwJDIAEnrICKyBkP/BugNRhzIvwoAGNyNZ/pE5PwY+jfg+OVK1OoM8PdSox+naxKRC2Do34Btp0sAACN6BEKp5J24ROT8GPo34PtjRQCASQNCJa6EiMg6DP12ulBWiyOXTDN3JvQNkbgaIiLrMPTbafXWswCAW3oFIYQraxKRi3CK0F+1ahWioqLg6emJuLg47N69u832X3zxBWJjY+Hp6YlBgwYhPT29gyo1qajV4+sDlwAA88bGdOhnExHdCMlD/7PPPkNycjJSU1Oxf/9+DBkyBJMmTUJxcXGr7Xfu3ImZM2fioYcewoEDBzBt2jRMmzYNR48e7ZB6hRD4343HUaMzoE+oD8b2DuqQzyUisgfJQ3/FihWYN28e5s6di/79+2P16tXw9vbG2rVrW23/5ptvYvLkyXjmmWfQr18/vPTSSxg+fDjeeecdh9Z5oqASX+6/iFkf7sYX+y4CAFKnDuD6+UTkUjyk/HCdTod9+/YhJSXFvE2pVCIhIQE5OTmt7pOTk4Pk5GSLbZMmTcLXX3/davuGhgY0NDSYX1dWmlbE1Ov10Ov1Vtf68sbj2HbmCgBAoQCem9IXo3v423QMe2n+TCk+2xHk1h9Afn2SW38AefXJlj5IGvqlpaUwGAwIDbWc8hgaGoqTJ0+2uk9hYWGr7QsLC1ttn5aWhmXLlrXYvnnzZnh7e1tdq3e9Ej19FejaSWBMiBEh5ceQnn7M6v0dISMjQ9LPtze59QeQX5/k1h9AHn2qra21uq2kod8RUlJSLP5nUFlZicjISEycOBF+ftbfRXuHI4prJ71ej4yMDCQmJkKtVktdzg2TW38A+fVJbv0B5NWn5hEMa0ga+kFBQVCpVCgqKrLYXlRUhLCwsFb3CQsLs6m9VquFVqttsV2tVrv8F1oOffg1ufUHkF+f5NYfQB59sqV+SS/kajQajBgxApmZmeZtRqMRmZmZiI+Pb3Wf+Ph4i/aA6b9n12pPRES/kHx4Jzk5GXPmzMHIkSMxevRorFy5EjU1NZg7dy4AYPbs2YiIiEBaWhoA4KmnnsK4cePw+uuv484778SGDRuwd+9evP/++1J2g4jIJUge+jNmzEBJSQmWLl2KwsJCDB06FJs2bTJfrM3Pz4dS+ct/SMaMGYNPPvkEixcvxnPPPYfevXvj66+/xsCBA6XqAhGRy5A89AEgKSkJSUlJrb6XnZ3dYts999yDe+65x8FVERHJj+Q3ZxERUcdh6BMRuRGGPhGRG2HoExG5EYY+EZEbYegTEbkRp5iy2ZGEEABsW6vC2ej1etTW1qKystLlbx8H5NcfQH59klt/AHn1qTnPmvOtLW4X+lVVVQCAyMhIiSshIrKvqqoq+Pv7t9lGIaz5p0FGjEYjLl++DF9fX5d9AErzSqEXLlywaaVQZyW3/gDy65Pc+gPIq09CCFRVVaFr164WKxi0xu3O9JVKJbp16yZ1GXbh5+fn8t+svya3/gDy65Pc+gPIp0/XO8Nvxgu5RERuhKFPRORGGPouSKvVIjU1tdWHw7giufUHkF+f5NYfQJ59sobbXcglInJnPNMnInIjDH0iIjfC0CciciMMfSIiN8LQdzGrVq1CVFQUPD09ERcXh927d0tdEgDghRdegEKhsPgVGxtrfr++vh5PPPEEunTpAh8fH/zxj39EUVGRxTHy8/Nx5513wtvbGyEhIXjmmWfQ2Nho0SY7OxvDhw+HVqtFr169sG7dOrvU/+OPP2Lq1Kno2rUrFAoFvv76a4v3hRBYunQpwsPD4eXlhYSEBJw+fdqiTVlZGR544AH4+fkhICAADz30EKqrqy3aHD58GLfeeis8PT0RGRmJV155pUUtX3zxBWJjY+Hp6YlBgwYhPT3dIX3605/+1OJrNnnyZKftU1paGkaNGgVfX1+EhIRg2rRpOHXqlEWbjvw+c9afxesS5DI2bNggNBqNWLt2rTh27JiYN2+eCAgIEEVFRVKXJlJTU8WAAQNEQUGB+VdJSYn5/UcffVRERkaKzMxMsXfvXnHTTTeJMWPGmN9vbGwUAwcOFAkJCeLAgQMiPT1dBAUFiZSUFHObc+fOCW9vb5GcnCyOHz8u3n77baFSqcSmTZtuuP709HTx/PPPiy+//FIAEF999ZXF+8uXLxf+/v7i66+/FocOHRK///3vRXR0tKirqzO3mTx5shgyZIj46aefxLZt20SvXr3EzJkzze9XVFSI0NBQ8cADD4ijR4+KTz/9VHh5eYm///3v5jY7duwQKpVKvPLKK+L48eNi8eLFQq1WiyNHjti9T3PmzBGTJ0+2+JqVlZVZtHGmPk2aNEl89NFH4ujRo+LgwYPijjvuEN27dxfV1dXmNh31febMP4vXw9B3IaNHjxZPPPGE+bXBYBBdu3YVaWlpElZlkpqaKoYMGdLqe1evXhVqtVp88cUX5m0nTpwQAEROTo4QwhRQSqVSFBYWmtu89957ws/PTzQ0NAghhPjLX/4iBgwYYHHsGTNmiEmTJtm1L78NSKPRKMLCwsSrr75q0SetVis+/fRTIYQQx48fFwDEnj17zG2+++47oVAoxKVLl4QQQrz77rsiMDDQ3B8hhFi0aJHo27ev+fW9994r7rzzTot64uLixJ///Ge79kkIU+jfdddd19zH2ftUXFwsAIitW7cKITr2+8yZfxavh8M7LkKn02Hfvn1ISEgwb1MqlUhISEBOTo6Elf3i9OnT6Nq1K2JiYvDAAw8gPz8fALBv3z7o9XqL2mNjY9G9e3dz7Tk5ORg0aBBCQ0PNbSZNmoTKykocO3bM3ObXx2hu4+j+5+bmorCw0OKz/f39ERcXZ1F/QEAARo4caW6TkJAApVKJXbt2mduMHTsWGo3Gov5Tp06hvLzc3KYj+5idnY2QkBD07dsXjz32GK5cuWJ+z9n7VFFRAQDo3LkzgI77PnOFn8W2MPRdRGlpKQwGg8U3KwCEhoaisLBQoqp+ERcXh3Xr1mHTpk147733kJubi1tvvRVVVVUoLCyERqNBQECAxT6/rr2wsLDVvjW/11abyspK1NXVOahnv3x+W3/3hYWFCAkJsXjfw8MDnTt3tksfHfE1njx5MtavX4/MzEz87W9/w9atWzFlyhQYDAan75PRaMTTTz+Nm2++GQMHDjR/Tkd8nzn7z+L1uN0qm+QYU6ZMMf958ODBiIuLQ48ePfD555/Dy8tLwsroWu677z7znwcNGoTBgwejZ8+eyM7Oxu233y5hZdf3xBNP4OjRo9i+fbvUpbgcnum7iKCgIKhUqhYzEYqKihAWFiZRVdcWEBCAPn364MyZMwgLC4NOp8PVq1ct2vy69rCwsFb71vxeW238/Pwc+g9L8+e39XcfFhaG4uJii/cbGxtRVlZmlz52xNc4JiYGQUFBOHPmjLkWZ+xTUlISvv32W2RlZVksk95R32eu9rP4Wwx9F6HRaDBixAhkZmaatxmNRmRmZiI+Pl7CylpXXV2Ns2fPIjw8HCNGjIBarbao/dSpU8jPzzfXHh8fjyNHjliETEZGBvz8/NC/f39zm18fo7mNo/sfHR2NsLAwi8+urKzErl27LOq/evUq9u3bZ26zZcsWGI1GxMXFmdv8+OOP0Ov1FvX37dsXgYGB5jZS9BEALl68iCtXriA8PNxcizP1SQiBpKQkfPXVV9iyZQuio6Mt3u+o7zNX+1lsQeoryWS9DRs2CK1WK9atWyeOHz8uHnnkEREQEGAxE0EqCxYsENnZ2SI3N1fs2LFDJCQkiKCgIFFcXCyEME2l6969u9iyZYvYu3eviI+PF/Hx8eb9m6fSTZw4URw8eFBs2rRJBAcHtzqV7plnnhEnTpwQq1atstuUzaqqKnHgwAFx4MABAUCsWLFCHDhwQJw/f14IYZqyGRAQIP7zn/+Iw4cPi7vuuqvVKZvDhg0Tu3btEtu3bxe9e/e2mN549epVERoaKmbNmiWOHj0qNmzYILy9vVtMb/Tw8BCvvfaaOHHihEhNTW33lM22+lRVVSUWLlwocnJyRG5urvjhhx/E8OHDRe/evUV9fb1T9umxxx4T/v7+Ijs722KaaW1trblNR32fOfPP4vUw9F3M22+/Lbp37y40Go0YPXq0+Omnn6QuSQhhmtIWHh4uNBqNiIiIEDNmzBBnzpwxv19XVycef/xxERgYKLy9vcX06dNFQUGBxTHy8vLElClThJeXlwgKChILFiwQer3eok1WVpYYOnSo0Gg0IiYmRnz00Ud2qT8rK0sAaPFrzpw5QgjTtM0lS5aI0NBQodVqxe233y5OnTplcYwrV66ImTNnCh8fH+Hn5yfmzp0rqqqqLNocOnRI3HLLLUKr1YqIiAixfPnyFrV8/vnnok+fPkKj0YgBAwaIjRs32r1PtbW1YuLEiSI4OFio1WrRo0cPMW/evBah5Ux9aq0vACy+Bzry+8xZfxavh0srExG5EY7pExG5EYY+EZEbYegTEbkRhj4RkRth6BMRuRGGPhGRG2HoExG5EYY+kZ1lZ2dDoVC0WAOmLS+88AKGDh3qsJqImjH0ya2tXr0avr6+Fo/Lq66uhlqtxvjx4y3aNof52bNn2zzmmDFjUFBQAH9/f7vWOn78eDz99NN2PSa5H4Y+ubUJEyaguroae/fuNW/btm0bwsLCsGvXLtTX15u3Z2VloXv37ujZs2ebx9RoNAgLC4NCoXBY3UTtxdAnt9a3b1+Eh4cjOzvbvC07Oxt33XUXoqOj8dNPP1lsnzBhAoxGI9LS0hAdHQ0vLy8MGTIE//73vy3a/XZ4Z82aNYiMjIS3tzemT5+OFStWtHjYBwD84x//QFRUFPz9/XHfffehqqoKgOkh5lu3bsWbb75pfoh5Xl6evf86yA0w9MntTZgwAVlZWebXWVlZGD9+PMaNG2feXldXh127dmHChAlIS0vD+vXrsXr1ahw7dgzz58/H//zP/2Dr1q2tHn/Hjh149NFH8dRTT+HgwYNITEzEX//61xbtzp49i6+//hrffvstvv32W2zduhXLly8HALz55puIj4/HvHnzUFBQgIKCAkRGRjrgb4Pkjk/OIrc3YcIEPP3002hsbERdXR0OHDiAcePGQa/XY/Xq1QBMz01taGjA+PHj0b9/f/zwww/mtdNjYmKwfft2/P3vf8e4ceNaHP/tt9/GlClTsHDhQgBAnz59sHPnTnz77bcW7YxGI9atWwdfX18AwKxZs5CZmYm//vWv8Pf3h0ajgbe3t0s8qIOcF0Of3N748eNRU1ODPXv2oLy8HH369EFwcDDGjRuHuXPnor6+HtnZ2YiJiUF1dTVqa2uRmJhocQydTodhw4a1evxTp05h+vTpFttGjx7dIvSjoqLMgQ8A4eHhLZ5cRXSjGPrk9nr16oVu3bohKysL5eXl5rP1rl27IjIyEjt37kRWVhZuu+02VFdXAwA2btyIiIgIi+NotdobqkOtVlu8VigUMBqNN3RMot9i6BPBNMSTnZ2N8vJyPPPMM+btY8eOxXfffYfdu3fjscceQ//+/aHVapGfn9/qUE5r+vbtiz179lhs++1ra2g0GhgMBpv3I/o1hj4RTKH/xBNPQK/XW4T5uHHjkJSUBJ1OhwkTJsDX1xcLFy7E/PnzYTQaccstt6CiogI7duyAn58f5syZ0+LYTz75JMaOHYsVK1Zg6tSp2LJlC7777jubp3RGRUVh165dyMvLg4+PDzp37gylknMxyDb8jiGCKfTr6urQq1cvhIaGmrePGzcOVVVV5qmdAPDSSy9hyZIlSEtLQ79+/TB58mRs3LixxYO6m918881YvXo1VqxYgSFDhmDTpk2YP38+PD09bapx4cKFUKlU6N+/P4KDg5Gfn9/+DpPb4uMSiSQwb948nDx5Etu2bZO6FHIzHN4h6gCvvfYaEhMT0alTJ3z33Xf4+OOP8e6770pdFrkhnukTdYB7770X2dnZqKqqQkxMDJ588kk8+uijUpdFboihT0TkRnghl4jIjTD0iYjcCEOfiMiNMPSJiNwIQ5+IyI0w9ImI3AhDn4jIjTD0iYjcCEOfiMiN/H9Cri4JPFdI7wAAAABJRU5ErkJggg==", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ - "plot_cumulative_density(np.sum(X_tran, axis=1))" + "# subset to donor 0\n", + "adata = ad.read_h5ad('resources/grn-benchmark/multiomics_rna.h5ad')\n", + "adata[adata.obs.donor_id=='donor_0', :].write('resources/grn-benchmark/multiomics_rna_0.h5ad')" ] }, { @@ -120,9 +113,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "RUN_ID=\"benchmark_donor_0_default\" -> grnboost2, check how this compares to when there was 5000 genes \n" - ] + "source": [] }, { "cell_type": "markdown", @@ -350,35 +341,20 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "upload: resources/grn-benchmark/multiomics_rna_qc.h5ad to s3://openproblems-data/resources/grn/grn-benchmark/multiomics_rna_qc.h5ad\n", - "upload: resources/grn-benchmark/multiomics_rna_0.h5ad to s3://openproblems-data/resources/grn/grn-benchmark/multiomics_rna_0.h5ad\n", - "upload: resources/grn-benchmark/multiomics_atac_0.h5ad to s3://openproblems-data/resources/grn/grn-benchmark/multiomics_atac_0.h5ad\n" + "delete: s3://openproblems-data/resources/grn/grn-benchmark/multiomics_rna_qc.h5ad\n", + "delete: s3://openproblems-data/resources/grn/grn-benchmark/multiomics_rna_qc_bc.h5ad\n" ] } ], "source": [ - "!aws s3 sync resources/grn-benchmark s3://openproblems-data/resources/grn/grn-benchmark " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import yaml\n", - "import pandas as pd \n", - "import matplotlib.pyplot as plt\n", - "import anndata as ad\n", - "controls = ['negative_control','positive_control']\n", - "grn_models = ['collectri','granie', 'figr', 'celloracle', 'scglue', 'scenicplus']" + "!aws s3 sync resources/grn-benchmark s3://openproblems-data/resources/grn/grn-benchmark --delete" ] }, { diff --git a/src/api/comp_method.yaml b/src/api/comp_method.yaml index f156ccce6..fbe58e220 100644 --- a/src/api/comp_method.yaml +++ b/src/api/comp_method.yaml @@ -51,11 +51,18 @@ functionality: type: boolean direction: input default: false + description: normalize rna seq data before inference. currently, it's only applicable to baseline models + - name: --only_hvgs + type: boolean + direction: input + default: false + description: subset rna seq data to only 7000 hvgs to reduce dimensionality - + resources: + - path: ../utils/util.py test_resources: - type: python_script diff --git a/src/control_methods/pearson/script.py b/src/control_methods/pearson/script.py index 045198a95..f2396c957 100644 --- a/src/control_methods/pearson/script.py +++ b/src/control_methods/pearson/script.py @@ -6,6 +6,7 @@ from tqdm import tqdm from scipy.stats import spearmanr from sklearn.preprocessing import StandardScaler +from utils.util import process_data, create_corr_net ## VIASH START par = { @@ -15,62 +16,26 @@ 'cell_type_specific': True, 'max_n_links': 50000, 'prediction': 'resources/grn_models/donor_0_default/pearson.csv', - "seed": 32 + "seed": 32, + 'only_hvgs': True, + 'normalize': False } ## VIASH END print(par) -def process_links(net, par): - net = net[net.source!=net.target] - net_sorted = net.reindex(net['weight'].abs().sort_values(ascending=False).index) - net = net_sorted.head(par['max_n_links']).reset_index(drop=True) - return net - -def corr_net(X, gene_names, par): - X = StandardScaler().fit_transform(X) - net = np.dot(X.T, X) / X.shape[0] - net = pd.DataFrame(net, index=gene_names, columns=gene_names) - net = net.sample(len(tf_all), axis=1, random_state=par['seed']) - net = net.reset_index() - index_name = net.columns[0] - net = net.melt(id_vars=index_name, var_name='source', value_name='weight') - - net.rename(columns={index_name: 'target'}, inplace=True) - net = process_links(net, par) - - return net - -def create_corr_net(X, gene_names, groups, par): - if par['cell_type_specific']: - i = 0 - for group in tqdm(np.unique(groups), desc="Processing groups"): - X_sub = X[groups == group, :] - net = corr_net(X_sub, gene_names, par) - net['cell_type'] = group - if i==0: - grn = net - else: - grn = pd.concat([grn, net], axis=0).reset_index(drop=True) - i += 1 - else: - grn = corr_net(X, gene_names, par) - return grn print('Read data') multiomics_rna = ad.read_h5ad(par["multiomics_rna"]) - - +process_data(multiomics_rna) + gene_names = multiomics_rna.var_names.to_numpy() tf_all = np.loadtxt(par['tf_all'], dtype=str) groups = multiomics_rna.obs.cell_type tf_all = np.intersect1d(tf_all, gene_names) -if par['normalize']: - print('Noramlize data') - sc.pp.normalize_total(multiomics_rna) - sc.pp.log1p(multiomics_rna) - sc.pp.scale(multiomics_rna) +n_tfs = len(tf_all) + print('Create corr net') -net = create_corr_net(multiomics_rna.X, multiomics_rna.var_names, groups, par) +net = create_corr_net(multiomics_rna.X, multiomics_rna.var_names, groups, par, tf_all=None) print('Output GRN') net.to_csv(par['prediction']) diff --git a/src/methods/single_omics/genie3/script.py b/src/methods/single_omics/genie3/script.py index c4a6f6af1..7057bbcee 100644 --- a/src/methods/single_omics/genie3/script.py +++ b/src/methods/single_omics/genie3/script.py @@ -3,6 +3,7 @@ import anndata import numpy as np import pandas as pd +import scanpy as sc from arboreto.algo import genie3 from distributed import Client @@ -12,7 +13,8 @@ 'multiomics_rna': 'resources/grn-benchmark/multiomics_rna_0.h5ad', "tf_all": 'resources/prior/tf_all.csv', 'prediction': 'output//genie3_donor0.csv', - 'max_n_links': 50000 + 'max_n_links': 50000, + 'only_hvgs': True } ## VIASH END @@ -27,6 +29,7 @@ tfs = set(list(df['gene_name'])) tf_names = [gene_name for gene_name in gene_names if (gene_name in tfs)] + print('GRN inference') client = Client(processes=False) network = genie3(X, client_or_address=client, gene_names=gene_names, tf_names=tf_names) diff --git a/src/utils/util.py b/src/utils/util.py new file mode 100644 index 000000000..bac97c72e --- /dev/null +++ b/src/utils/util.py @@ -0,0 +1,55 @@ +import pandas as pd +import anndata as ad +import numpy as np +from tqdm import tqdm +import scanpy as sc + +def process_links(net, par): + net = net[net.source!=net.target] + net_sorted = net.reindex(net['weight'].abs().sort_values(ascending=False).index) + net = net_sorted.head(par['max_n_links']).reset_index(drop=True) + return net + +def corr_net(X, gene_names, par, tf_all=None): + X = StandardScaler().fit_transform(X) + net = np.dot(X.T, X) / X.shape[0] + net = pd.DataFrame(net, index=gene_names, columns=gene_names) + if tf_all is None: + net = net.sample(n_tfs, axis=1, random_state=par['seed']) + else: + net = net[tf_all] + net = net.reset_index() + index_name = net.columns[0] + net = net.melt(id_vars=index_name, var_name='source', value_name='weight') + + net.rename(columns={index_name: 'target'}, inplace=True) + net = process_links(net, par) + + return net + +def create_corr_net(X, gene_names, groups, par, tf_all): + if par['cell_type_specific']: + i = 0 + for group in tqdm(np.unique(groups), desc="Processing groups"): + X_sub = X[groups == group, :] + net = corr_net(X_sub, gene_names, par) + net['cell_type'] = group + if i==0: + grn = net + else: + grn = pd.concat([grn, net], axis=0).reset_index(drop=True) + i += 1 + else: + grn = corr_net(X, gene_names, par) + return grn +def process_data(adata, par): + if par['normalize']: + print('Noramlize data') + sc.pp.normalize_total(adata) + sc.pp.log1p(adata) + sc.pp.scale(adata) + adata.X = adata.X.todense() + if par['only_hvgs']: + print('Subsetting data to hvgs') + adata = adata[:, adata.var.hvg_counts] + print('New dimension of data: ', adata.shape) \ No newline at end of file