diff --git a/runs.ipynb b/runs.ipynb index ee42d7f03..4151fa8a3 100644 --- a/runs.ipynb +++ b/runs.ipynb @@ -87,7 +87,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -110,7 +110,7 @@ "\n", "sys.path.append('../')\n", "from grn_benchmark.src.commons import surragate_names\n", - "from src.utils.helper import *\n", + "from src.helper import *\n", "par = {\n", " # 'methods': [ 'collectri', 'negative_control', 'positive_control', 'pearson_corr', 'portia', 'ppcor', 'genie3', 'grnboost2', 'scenic', 'scglue', 'celloracle'],\n", " 'methods': [ 'collectri', 'negative_control', 'positive_control', 'pearson_corr', 'portia', 'ppcor', 'grnboost2', 'scenic', 'scglue', 'celloracle'],\n", @@ -185,14 +185,14 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Submitted batch job 7757950\n" + "Submitted batch job 7758859\n" ] } ], @@ -201,12 +201,525 @@ " run_consensus(par)\n", "\n", "if True: # run metrics/script_all.py\n", - " # !bash scripts/sbatch/calculate_scores.sh #includes both reg1 and 2. #inside the script, set the reg_type, read and write dirs, and methods\n", - " !sbatch scripts/sbatch/calculate_scores.sh #includes both reg1 and 2. #inside the script, set the reg_type, read and write dirs, and methods\n", - "\n", - "if False: # check the scores\n", - " df_scores = check_scores()\n", - " df_scores.style.background_gradient()" + " calculate_scores()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "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", + " \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", + "
S1S2static-theta-0.0static-theta-0.5rank
scenicplus0.2450330.4034940.7605830.5392094
collectri-0.100238-0.2111820.4855060.45725911
negative_control-0.039305-0.0410040.2746590.44038312
positive_control0.1971290.5788220.8720030.5954892
pearson_corr0.2693790.5092970.7351560.5170563
portia0.1489410.2272480.4736070.4676079
ppcor0.0228460.0941070.4307760.44914410
grnboost20.3810320.4598600.7481750.6157901
scenic0.1446960.2065710.6850340.5564857
scglue0.0783090.2388590.5305310.4834238
celloracle0.2168970.3114510.7115490.5641606
scenicplus0.2450330.4034940.7605830.5392094
\n", + "
" + ], + "text/plain": [ + " S1 S2 static-theta-0.0 static-theta-0.5 rank\n", + "scenicplus 0.245033 0.403494 0.760583 0.539209 4\n", + "collectri -0.100238 -0.211182 0.485506 0.457259 11\n", + "negative_control -0.039305 -0.041004 0.274659 0.440383 12\n", + "positive_control 0.197129 0.578822 0.872003 0.595489 2\n", + "pearson_corr 0.269379 0.509297 0.735156 0.517056 3\n", + "portia 0.148941 0.227248 0.473607 0.467607 9\n", + "ppcor 0.022846 0.094107 0.430776 0.449144 10\n", + "grnboost2 0.381032 0.459860 0.748175 0.615790 1\n", + "scenic 0.144696 0.206571 0.685034 0.556485 7\n", + "scglue 0.078309 0.238859 0.530531 0.483423 8\n", + "celloracle 0.216897 0.311451 0.711549 0.564160 6\n", + "scenicplus 0.245033 0.403494 0.760583 0.539209 4" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "df_scores = pd.read_csv(f\"resources/scores/hvg/skeleton_False/scgen_pearson-ridge.csv\", index_col=0)\n", + "df_all_n = (df_scores-df_scores.min(axis=0))/(df_scores.max(axis=0)-df_scores.min(axis=0))\n", + "df_scores['rank'] = df_all_n.mean(axis=1).rank(ascending=False).astype(int)\n", + "df_scores" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "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", + " \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", + "
 S1S2static-theta-0.0static-theta-0.5rank
collectri-0.100238-0.2111820.4855060.45725911
negative_control-0.044574-0.0451580.3598050.43845110
positive_control0.1971290.5788220.8720030.5954892
pearson_corr0.2734430.5163430.7829780.5382523
portia0.2633100.3570060.5663650.5075706
ppcor0.0179540.1597540.4680490.4549959
grnboost20.4219360.4893220.7889310.6294711
scenic0.1680060.2189160.7569650.5654345
granie0.0832980.1060120.1941640.36342512
scglue0.0808570.2936300.6603570.4807347
celloracle0.2091510.2914780.6900990.5763434
figr0.1136450.1931310.4280320.4652688
\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "df_scores = pd.read_csv(f\"resources/scores/full/skeleton_True/scgen_pearson-ridge.csv\", index_col=0)\n", + "df_all_n = (df_scores-df_scores.min(axis=0))/(df_scores.max(axis=0)-df_scores.min(axis=0))\n", + "df_scores['rank'] = df_all_n.mean(axis=1).rank(ascending=False).astype(int)\n", + "df_scores.style.background_gradient()" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(
,\n", + " )" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAFzCAYAAADSc9khAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA6ZklEQVR4nO3de1xUdf4/8NfMMAwXB7yhiKJG4iVNNMVLiiAKWLRruGtR8Q0v3+2yrVlmLi75S8q+Ylmbrd3rG+qXStuyMjcRCc2yFKHybnnLK4gIDnIZ5nJ+fwwzOAJ6DgJnzszr+Xj4yDnzOWfefpxeHj7ncz5HJQiCACIi8ghquQsgIqL2w9AnIvIgDH0iIg/C0Cci8iAMfSIiD8LQJyLyIAx9IiIPwtAnIvIgXnIX0N6sVivOnj0LvV4PlUoldzlERDdMEARUVlYiJCQEavW1z+U9LvTPnj2L0NBQucsgImp1p06dQq9eva7ZxuNCX6/XA7B1TkBAgOj9TCYTNm/ejPj4eGi12rYqTzHYHw3YFw3YF87aqz8MBgNCQ0Md+XYtHhf69iGdgIAAyaHv5+eHgIAAfpnB/rgS+6IB+8JZe/eHmCFrXsglIvIgDH0iIg/C0Cci8iAMfSIiD8LQJyLyIAx9IiIPInvo19XVIS0tDV5eXjhx4sR123/33XcYM2YMoqOjMWbMGGzfvr3tiyQichOyztM/ceIE7rvvPvTv3x8Wi+W67X///XckJibiq6++QlRUFLZt24a77roLe/bsQZ8+fdqhYiIiZZP1TP/y5ctYs2YNZs6cKar9ihUrcMsttyAqKgoAEB0djQEDBuC1115ryzKJiNyGrGf6Q4YMAQCcPn1aVPu8vDxMmDDBaVtkZCS2bNnS6rURNUUQBFgFwCoIEARAgACTyYI6C1BTZ4FJUEGwvw9AEADUt7O1tx2jzmKF2dJwDNS3Fa74HMf+ti2O31+53XHcq1431NuwvaLaBG8vced5Vx7jmu3g3NBituC3Syr8eOwiNF6aqxuLPKZ4La2zdY55fRazGQcrVOjw2wVoNJrrts8/dB5jwrpg4sBu8NFev31LKGoZhmPHjmH69OlO24KDg3H8+PFm9zEajTAajY7XBoMBgO32aJPJJPqz7W2l7OPO2rM/BEGAodaMc5dqUVppRHWdBRarAJNVgMVqxYkL1ejk7w2LVXD8MpqtOHahCsEBOlisAsxWARcu1+FiVR06+3ujus6Cn09V4OZu/rBabSFttQqw1If6ZaMZJQYjunbwduxvsQqormtuGNILT+/Ka/O+UAYNVh7YLXcRLkSDtw4WiW696offsWNBNIL0OtH7SPn/UFGhX11dDZ3OuSN0Oh2qq6ub3Wfp0qXIyMhotH3z5s3w8/OTXENubq7kfdyZmP6wCkCtBag0ASYrUGFUocYCnLysgo8GsAjAqSpArwWOV6rQWQccMdjWENGoBFiEtlsCe89pwzXfv3C5rs0+GwC81QJUAFD/R1Q1/Nbx+sr3mtsuZp86K1BtVqG7r7hT2bboddHHlPDhYpu2yZ+nDQ56UwcBu77Nk3Tsa2Xg1RQV+n5+fk5n7YDtTP5a4b1w4ULMmzfP8dq+Gl18fLzkBddyc3MRFxfHhaTQ0B+TJ0/G8YtG5B8uRUWNCWcralFVZ8a3v5Uh0NcLNSYr6sxWSce+eMVf8dWB36ujD8xWAb07+8FLo4KXWgWVSoXfy6oxrFcg1GrbNrVaBUGwDWn069YBmvrtJosVGrUawQE6aNQqmK0CunbwhkZlO45aBWjUKqhUth/zfbUa+Hpr4KVWQVP/y8dLDS+N2ha2KsBsNiM/Px+xsRPhrdVCrVLVB6/K0cb23yteu+mzHPj/ibP26g/7CIYYigr9sLAwlJSUOG0rLi5GWFhYs/vodLpGPx0AgFarbdFfQkv3UzKTxYrN+0vwa0kljl+owv6zl2C1Cjhe5gX8sLXZ/S7VmJvc3k2vw2WjGQOC9SivqkPcLd2h1ahx4bIRQ3oGwmiy4qau/vD2UqNXJ18E+Grh7aWGXuflkmFpMpngowE6+vt63HejOZ74/8m1tHV/SDm2okJ/0qRJ+PHHH5227d69G5MnT5apIvciCAL2nzXgaOllfPfbBWz9tRTBAT7Ye+aSqP0jegViQLAeg3oEoHuAD0I6+qKLvze8vdTo6KeF7uoLe0TU7lw69O+//35oNBqsWbMGADB37ly8//77+P777zFu3Dhs374dhw4dwrp162SuVJmOnK/Em1uP4WjpZfx8qqLJNqWVzsNpiUN74LbenQDBijO/7ceUmLG4tVdn+Hoz0ImUQNbQr6urQ3x8PCoqKgAAycnJCA0NxSeffAIAqK2tdXreY58+ffDVV1/hqaeegre3N4xGI7766ivemCVR4e/lyPz6IApOlDf5foCPF8b16wqrICAsqAMmD+qOW3sGOk33M5lM+M/FfRge2hHaNppaRkStT9bQ9/b2xtatW5t9/7PPPmu0LSoqqtEQD11fTZ0FK/J+w1vbjjZ6r2dHX8xP6I9BPQIQ1rWD6LncRKQ8Lj28Qzfuu98uIOX9nY7ZKHaBvlo8MTkcM8fdJF9xRNTuGPpuquhkOV7cdAg/HrsIoCHw77w1GH8e0QuxA7vLWB0RyYWh72Yqquvw90/3IGe/89TWGbf3xTOJg+Cl4dANkSdj6LuRA2cNuPO1hqWme3XyxYrkYRjRp7OMVRGRK2Hou4kj5yvxl9UN652kjOmN56cOccmbmYhIPgx9N5B7oMQR+IG+Wnz66Fj066aXuSoickUMfQUrvlSLv6ze7bhjtmdHX7z9XyMY+ETULIa+Qn2y+xSe/vcex+v+3Tvgy7+Nb7M1uInIPTD0FWjVjhN49sv9jtf/N3s0xod3lbEiIlIKhr7C/O93x/HcVwcAALf0CMDah8dA78PVDIlIHIa+gqwtOOkI/DtvDcZrycM5756IJGHoK4AgCFj42V58XHAKAHDX0B5YkTwcGjWnYxKRNAx9F2e1Ckj813c4eM72ZJzbb+6C5dMjGPhE1CIMfRdWa7LgT2/ucAT+zHF98ewfBstcFREpGUPfhaW8txP7z9oCf+m0W3HfqN4yV0RESsfQd1HrfzqN3b/bHnLyVsptmDKkh8wVEZE74NQPF1RlNCPt070AgHtHhjLwiajVMPRd0IubDsFotsLPW4NFf7hF7nKIyI0w9F1M/uHzWPXD7wCAF5KGoIOOI3BE1HoY+i7kbEUNHv/oJwDAlMHBSBreS+aKiMjdMPRdhCAIeHhNISprzQjv1gHL74mQuyQickMMfRfxxc9nHUskv/jnoRzWIaI2wdB3AecNtUj7zLZM8oNj+2B4704yV0RE7oqh7wLe2HoUtSYrgvQ6/OPOQXKXQ0RujKEvs8tGM9btti2kNi+uPx+CQkRtiqEvs1U7TqC6zoJenXxx78hQucshIjfH0JdR2WUjVn5zBADw3+NvgporZxJRG2Poy+id7cdQY7Kge4AOD47tK3c5ROQBGPoyOVNRg/e2HwcAzI8fwLN8ImoXDH2ZrC04BYtVwLDQjpjOsXwiaicMfRlUGc14f/sxAEDq7X1kroaIPAlDXwbfHbmAqjoLenb0xR8jespdDhF5EIa+DNbVP+A8Krwrn3VLRO2Kod/ODhUbkHfoPAAgZQyHdoiofTH029mnhacBALEDu2FIz0CZqyEiT8PQb0cWq4BP6kP/T7dxrXwian8M/XaUf+g8KqpN8PfWYNKgbnKXQ0QeiKHfjnIPlAAA/jisJxdWIyJZMPTbiSAI2HLQFvoTBwTJXA0ReSqGfjvZe+YSyqrq4K1RIyqcoU9E8mDot5MNv5wFAEQPCIKvN4d2iEgeDP12YLEK+LI+9O8a2kPmaojIkzH028EvpytQYjDCR6vG5EHd5S6HiDwYQ78d5OwrBgBMHtQd/jovmashIk/mEqG/fv16REZGIioqCtHR0di/f3+zbY1GI5588klEREQgOjoao0ePxvr169uxWul+OlUBAJjQnxdwiUhesp927tq1C6mpqSgsLER4eDhWr16NhIQEHDx4EHq9vlH7JUuW4PPPP8fPP/+MwMBA/PTTTxgzZgx27dqFiIgIGf4E11ZaacTPJysAALf17iRvMUTk8WQ/08/MzERiYiLCw8MBACkpKTCbzcjKymqy/c8//4zIyEgEBtrWrRk+fDgCAwPxzTfftFfJkmzadw51Fitu7RmIm4P85S6HiDyc7KGfl5eHkSNHOl6r1WqMGDECW7ZsabL9n/70J2zfvh0nT54EAOTk5KC0tBTdu7vmBVL70E5UeFeoVFxGmYjkJevwTllZGQwGQ6PADg4ORkFBQZP7zJgxA9XV1Rg6dCh69OiBX3/9FX/+859xzz33NNneaDTCaDQ6XhsMBgCAyWSCyWQSXau9rZR9AGDHkQsAgBG9AyXv68pa2h/uiH3RgH3hrL36Q8rxZQ396upqAIBOp3PartPpHO9d7b333kNmZiYKCwtx880345dffsGWLVugVjf9Q8vSpUuRkZHRaPvmzZvh5+cnuebc3FzRbS8agWKDF1QQcOHgLvznN8kf5/Kk9Ie7Y180YF84a+v+aC4vmyJr6NtD98ozcfvrpgJZEAQsWLAATz31FG6++WYAQEREBObNm4eamho888wzjfZZuHAh5s2b53htMBgQGhqK+Ph4BAQEiK7VZDIhNzcXcXFx0Gq1ovbZefwiULQbfbv4Y9ofx4v+LCVoSX+4K/ZFA/aFs/bqD/sIhhiyhn6XLl0QGBiIkpISp+3FxcUICwtr1L60tBTl5eXo27ev0/abbroJn376aZOhr9PpGv0kAQBarbZFfwlS9hNUtp8+dFqN2/4P0NJ+dEfsiwbsC2dt3R9Sji37hdzY2FgUFhY6XguCgKKiIkyePLlR265du0Kn0+HcuXNO28+dO9eioZq2ZrYIAACtRvZuJiIC4AKhn5aWho0bN+LIkSMAgOzsbGg0GqSmpgIAxo8fj/T0dAC2mT2pqal47733UF5eDgAoKipCbm5usxdy5WSyWAEAWg1n7RCRa5D95qxRo0YhKysLycnJ8PX1hVqtRk5OjuPGrOrqaqcx/3/+859YvHgxJk2aBD8/P1RWViIzMxOPP/64XH+EZpmttjN9L57pE5GLkD30ASApKQlJSUlNvldUVOT02s/PDy+++GJ7lHXDLlbVAQCfkkVELoOnoG1o94mLAIBbeoifJURE1JYY+m3owDnbNKrIvlxzh4hcA0O/jVTXmXHk/GUAwC0hPNMnItfA0G8je09fglUAuul1CA7wkbscIiIADP02c7ikEgAwqEcAF1ojIpfB0G8jh4ttoT+wR+NnAhARyYWh30aOX6gCAIR15Rr6ROQ6GPptwGIVsOf0JQDA4JBAmashImogOfSzs7Pbog63cqKsCpeNZui81BgYzOEdInIdku/IXbBgAUwmE+655x6XXOTMFfxaP57fr1sHLsFARC5FciL16NEDFRUViIuLw6xZs7B9+/a2qEvRDtaH/hAO7RCRi5Ec+mvWrMETTzyB77//HnPmzMEnn3yCMWPGYOnSpThz5kxb1Kg4J8tsF3H7dOVPQkTkWiSH/qBBgxy/Hz58OO6//34MGDAA6enpGDx4MO688058+umnEAShVQtVkqOlttDv24Uzd4jItUgO/cTERJw/fx7Lly/H4MGDMX78eJw+fRqrVq3CuXPn8O6772Lv3r24//7726Jel1dntjrm6HN4h4hcjeQLuVu3bkVoaChCQkKQmpqKGTNmOD2+sGfPnli8eDGGDx/emnUqxpHzl1FnsULv44XQzr5yl0NE5ERy6Hfr1g3vv/8+YmNjm22zdOlSj53Zc7TUtshav24duPwCEbkcycM76enpjQK/rq4Ob775JsrKygAACxcuxPfff986FSqMfTllzs8nIlckOfQ//PDDRttUKhUqKysxffr0VilKyQ45Qp/LKROR62mVO4e0Wi0WLFiAqqqq1jicYgmCgF/ql1+4tRcv4hKR6xE1pr9ixQqsWLECAFBcXIywsLBGbS5duoSRI0e2bnUKU2k0O56Ly+EdInJFokI/JiYGHTt2hCAIWLZsGdLS0pzeV6vVCAoKuubFXU9QfKkWABDg4wU/b5d45jwRkRNRyRQREYGIiAgAgE6nw3333demRSnV6fJqAEDPTp45c4mIXJ/kMf1rBf5DDz10Q8Uo3YkLttDv05mhT0SuSdSZ/ueff47OnTtjwoQJmDVrVrPtNm3a1GqFKdH5SiMAoEdHPhOXiFyTqDP9559/Hm+99RYA4Ouvv4YgCE3+8nTl9RdxO/t5y1wJEVHTRJ3pFxYWOn4fHx+PDz74oMl2qamprVOVQpVetp3pd+mgk7kSIqKmSR7TX7VqVYve8wT2C7khHN4hIhclOfQLCgrw3HPP4dixYwCAN954AxEREZg+fTpKSkpavUClEAQBZytsUzZ7cfYOEbkoyaGfkZEBq9WKTp064aeffsKcOXOQkJCAXr164W9/+1tb1KgINSYLLhvNAIDgQJ7pE5FrknwHUVVVFRYvXgwAePbZZzFx4kS8+OKLAIDx48e3anFKct5gG8/31Wrg762RuRoioqZJPtM3mUwAAKPRiE8++QSzZ892vOft7bmzVooNtqGd4EAfLqlMRC5L8pl+jx49MGPGDBQXFwMApk2bBkEQkJOTA6PR2OoFKsW5SzUAgOAADu0QkeuSfKb/9ttvw9/fHz4+Pli/fj10Oh2++OILLFu2DI899lhb1KgIZ8ptod+zE5+WRUSuS/KZfufOnfH66687bbv77rtx9913Ox6i4olK6sf0eaZPRK6sVdbTt/Pkh6iU1I/pdwvgjVlE5Lokh/7PP/+MmJgYdOrUCRqNxunXtm3b2qJGRThTYRveCQnk8A4RuS7JwzupqamYPHkynnrqKej1esdMFUEQ8OSTT7Z6gUpxrn4t/ZCODH0icl2SQ1+v1+Pll19u8j3707U8jdlidTwxi8M7ROTKJA/vDB06FBcuXGjyvaKiohsuSIkuVtsCX60COnGFTSJyYS060x89ejRiY2MREhICjabh7tOsrCw88cQTrVmfIpRdtoV+Rz9vaNS8MYuIXJfk0H/nnXcwbNgwHDlyBEeOHHF6r6KiorXqUhT7s3E5XZOIXJ3k0B8/fjw2bNjQ5Hue+uxc+8ydHlxojYhcnOQx/eYCHwA++uijGypGqYo5c4eIFKJFN2ft3LkTqampuPfeewEAb731Fufog0sqE5Hrkxz6n3/+OSZPnozy8nIcPHgQADBw4EAsXLgQH3/8cYuKWL9+PSIjIxEVFYXo6Gjs37//mu2PHTuGP/3pT5g4cSIGDx6MMWPGYPfu3S367NZwtj70e3HdHSJycZJD/+WXX8Yvv/yCL7/8El26dAEAxMTEIDc3F2+88YbkAnbt2oXU1FR8+OGH2L59O2bPno2EhARUVlY22b60tBSTJk3C3LlzkZ+fj19++QV+fn6NLiq3J8eyyryQS0QuTnLoazQahIWFAYDTuvH+/v6wWq2SC8jMzERiYiLCw8MBACkpKTCbzcjKymqy/bJlyzB27FhMmDABAODl5YV33nnH8VoOFypti60F6XljFhG5NsmhX1lZiXPnzjXavnfv3mbPzq8lLy8PI0eObChIrcaIESOwZcuWJtt/9tlnjQK+X79+CAkJkfzZraGy1oSqOgsAoBvP9InIxUmesjl37lxEREQgOTkZp06dQkZGBg4fPowvv/wS77zzjqRjlZWVwWAwoHv37k7bg4ODUVBQ0Kh9VVUVjh8/DovFggceeAAnTpxAhw4d8MQTT+COO+5o8jOMRqPTw10MBgMA2xPA7E8BE8Pe9up9zpVXAQD8dRro1IKkYypZc/3hidgXDdgXztqrP6QcX3LoP/jgg+jevTuWLl2K8vJy/Otf/8KQIUOwfv16xMXFSTpWdXU1AECncx4W0el0jveuZL/5a9GiRcjPz0dERATy8vKQkJCAr7/+usnPX7p0KTIyMhpt37x5M/z8/CTVCwC5ublOr48aAMALvjDjP//5j+TjKd3V/eHJ2BcN2BfO2ro/msrL5kgOfQBISEhAQkJCS3Z1Yg/dqx+zaDQamwxk+5IPf/jDHxAREQEAmDRpEmJjY7FixYomQ3/hwoWYN2+e47XBYEBoaCji4+MREBAgulaTyYTc3FzExcVBq9U6tm8+UALs/wW9unXEnXeOFn08pWuuPzwR+6IB+8JZe/WHfQRDDMmhX1hYiIKCAlRUVKBz584YNWoUhg0bJvUwAIAuXbogMDAQJSUlTtuLi4sdF4uvFBQUBJ1Oh549ezpt79OnD3bs2NHkZ+h0ukY/SQCAVqtt0V/C1ftVGm0Xr7v46zzyS97SfnRH7IsG7Atnbd0fUo4tOvSPHj2KlJQU7Nq1C4IgOLarVCqMHTsW2dnZ6NOnj7RKAcTGxqKwsNDxWhAEFBUVIT09vVFbjUaDcePGNbqQXFJSgt69e0v+7NZQVr+kcid/rq5JRK5P1OydsrIyTJw4EZ07d8amTZtQVlYGk8mECxcuYOPGjdDr9YiOjkZ5ebnkAtLS0rBx40bHPPvs7GxoNBqkpqYCsK31c+U/AH//+9/xxRdf4OTJkwCAAwcOYPPmzbI9lP3CZdvQVNcOnK5JRK5P1Jn+yy+/jMTERLz55ptO2zt37owpU6ZgypQpeOSRR7B8+XK88MILkgoYNWoUsrKykJycDF9fX6jVauTk5ECv1wOwXaC4csw/Pj4er732GqZOnYoOHTrAbDZj1apVuOuuuyR9bmuxL6vctQPP9InI9YkK/dzcXOTn51+zzYsvvohJkyZJDn0ASEpKQlJSUpPvNfVglpSUFKSkpEj+nLZQVmX7B6kLQ5+IFEDU8I6vry86dOhwzTYBAQHw9fW8tWfsZ/pd/Dm8Q0SuT1Toe3mJu94rtp07Ka9/VCIfk0hESiAqpQ8ePIhZs2Zdt92hQ4duuCAlEQQBFdW2O+E6+nF6GhG5PlGhX1tbi+PHj4tq50mq6ywwmm3z9DtzyiYRKYCo0B82bNh1L+QCwMSJE2+4ICWxD+14a9Tw89ZcpzURkfxEjelv3rxZ1MHEtnMX5VW2oZ1O/lqnZaaJiFyVqNAXe4uvp912fZEXcYlIYVr0jFyyqay1nekH+HrWP3ZEpFwM/RtQbl93hzN3iEghGPo34GL9mD5n7hCRUrQ49E0mk2PRs5Y8G9cdOJZg4N24RKQQkkPfaDTikUcegb+/v2OK5qxZszB79mzU1NS0eoGuzLEEA9fdISKFkBz6aWlpOHPmDD7++GN069YNAPDee+9h0KBBTk+o8gSllVxWmYiURXLo7969G1988QWmTZvmWGDNy8sL8+fP97hlGOxr6QfpGfpEpAySQ99isUCttu125RO0AODixYutU5VC2Ofpd+GFXCJSCMmhHxgYiHfffRcAHHehVlVV4Zlnnmn07Fp3JggCLtXYZu8Ecp4+ESmE5LWQV6xYgSlTpuDpp5+GxWLBTTfdhHPnzqFXr17IyclpixpdUnWdBfYfdDr4eN6S0kSkTJLTqn///jh06BCys7Oxf/9+AMCQIUNw//33w9vbc4Y57Gf5Wo0KvloutkZEyiA59N9++208/PDDmDlzZlvUoxj2dfQDfbnYGhEph+TQT09PR0VFBVJSUjxqDP9qBq67Q0QKJPlCbnh4OEJDQ/HQQw9hypQpWL16Naqrq9uiNpdmH94J8GHoE5FySD7TX7duHUJDQ3H//fejuLgY2dnZiIuLw80334zU1FRMmjSpLep0OfbQ52MSiUhJJJ/p9+rVy/H74OBgjB07Frfeeiv+/e9/4+67727N2lyagWf6RKRAkkN/0qRJOHr0KBYvXox+/fohOjoaR48exZtvvolz5861RY0u6bLRDIDTNYlIWSQn1vfff4/+/ftj0KBB+Mtf/uKxF3Qv19aHvo6hT0TKITmx+vTpg7Vr12L48OFtUY9iVNVZAAD+3gx9IlIOycM7q1evbjbwP/vssxsuSClq6mxn+v463phFRMoh6jS1rq4OWq3tJqS6ujp8++23TbZbsmQJpk2b1qoFuioDh3eISIFEJVa/fv0wcOBAbN68GTExMc2286Q7U8vrV9js6Oc5S08QkfKJCv3PPvsMer0eABAdHY38/Pwm29mfpOUJKuvP9LnCJhEpiajQHzlypOP3b7zxRpNtqqqqmn3PHVXWL8PA4R0iUhLJF3Jfe+21RtuqqqowevRoZGdnt0pRSlBltM3e4Tx9IlISyYl1+PDhRtv8/f2xb98+REVFtUpRrk4QBFRx9g4RKZCo0N+2bRu2bdsGADhx4gSee+65Rm3Ky8tRVlbWutW5qEqj2fEAFb2OY/pEpByiQv/EiROOi7fl5eWNLuSq1WoEBQU5HqPo7uwXcbUaFXy0kkfIiIhkIyr0U1NTkZqa6vj9qlWr2rQoV2dfgkHvwweoEJGySD5NvVbg5+Xl3VAxSmGfuaPnRVwiUpgWpZbVasXRo0dRXFwMwT64DeDpp59GUVFRqxXnqiqNvBuXiJRJcmodPHgQSUlJ+PXXX6FSqZxC31OGOhqGdxj6RKQskod3nnjiCSxatAg1NTWYMGECrFYramtrkZ2djWeffbYtanQ5VfVn+lxhk4iURnLoG41GPPDAA9DpdI5t3t7euO+++/DTTz+1anGuqsZkuzHL15tz9IlIWSSHvslkcvzeYrE45ubX1NRg3759rVeZC6uuX0vfV8vQJyJlkRz6PXv2RHJyMioqKjBx4kSMHj0af/nLXxAZGYkBAwa0RY0up6Y+9P14pk9ECiN5UPqll17Cvn37oNVqsXDhQly4cAHbt2/HkCFD8Morr7RFjS6nYXiHY/pEpCySz/T79OmDxMRE+Pv7w8fHB6+//jr27NmDjz/+GCEhIZILWL9+PSIjIxEVFYXo6Gjs379f1H4rV66ESqXC1q1bJX/mjeLwDhEpVauuIZCUlCSp/a5du5CamooPP/wQ27dvx+zZs5GQkIDKyspr7nf27Fm89NJLN1LqDbHP3uEKm0SkNKJSKzY2VtTBfv75Z0kfnpmZicTERISHhwMAUlJSsGDBAmRlZWHOnDnN7jdnzhz84x//wCOPPCLp81pLw5RNnukTkbKICv3jx49jxowZ12134sQJSR+el5eH//f//p/jtVqtxogRI7Bly5ZmQ3/Dhg3QarVISEiQ9FmtyTG8w9AnIoURFfr33XefqBuvjEaj6A8uKyuDwWBA9+7dnbYHBwejoKCgyX2qqqqQnp6OnJwc0Z9lNBqd2hoMBgC2qadXTj+9Hntbk8mEqjrb73VqSDqGO7myPzwd+6IB+8JZe/WHlOOLCv3/+Z//EXWwu+66S/QHV1dXA4DTTV721/b3rrZo0SI88sgj6NGjh+ifKpYuXYqMjIxG2zdv3gw/Pz/R9drl5uaitEwDQIU9PxfCeFy47j7uLDc3V+4SXAb7ogH7wllb90dzmdkUyVciT5482ex78+fPx44dO0Qdxx64V5+xG43GJsO4qKgIO3fuxPLlyyVUCyxcuBDz5s1zvDYYDAgNDUV8fDwCAgJEH8dkMiE3NxdxcXFYfuhHoKYGMePGYnjvjpLqcRdX9odW69kPkmFfNGBfOGuv/rCPYIghOfT79u3bKgurdenSBYGBgSgpKXHaXlxcjLCwsEbtN27ciJqaGsdF5draWgC2tYA6duyI9957D/369Wu0n06na/TTBABotdoW/SVotVrUmKwAAL2fzuO/2C3tR3fEvmjAvnDW1v0h5diSQ3/06NH4+OOPHa8tFgtOnz6NtWvXYsKECZKOFRsbi8LCQsdrQRBQVFSE9PT0Rm0XLVqERYsWOV6fOHECN910E1599VXExMRI/WPckFoT5+kTkTK16CEqffr0cfwKCwvDhAkT8Prrrzv9YyBGWloaNm7ciCNHjgAAsrOzodFoHE/pGj9+fJP/AMjNaLaFvo6PSiQihZF8pt+/f/8mt5tMJvz222+SjjVq1ChkZWUhOTkZvr6+UKvVyMnJgV6vB2C7ONHULJ0nnngCP/74o+P3AwcOlPwPTktZrQJMFtvFW28NQ5+IlEVy6M+aNavRtsrKShQVFWHUqFGSC0hKSmr2Tt7mnsL16quvSv6c1lJnsTp+7+3F0CciZZGcWl9//TUEQXD8AoCQkBA888wzyMrKau36XE6duSH0dV4c0yciZZF8pp+cnIx//vOfbVGLItTWh75aBWg1nvF4SCJyH5LP9K8V+KtXr76hYpTAPnPHR6vxmGcCE5H7aNEykb///jt++eUXXLp0yenB6JmZmXjwwQdbrThXdGXoExEpjeTQX7ZsGdLT09G5c2f4+/s7vXf1jVbuqLb+xizO0SciJZIc+u+//z7279/f5KMR5Vz5sr0Y68f0dZy5Q0QKJDm5Bg8e3OyzcNeuXXvDBbk6e+hzuiYRKZHk5Hr88cfx1ltv4ezZs07j+QAwbdq0VivMVTnuxmXoE5ECSU4uvV6PN954A6GhofDy8oJGo3H82rZtW1vU6FLs8/R1HNMnIgWSPKY/c+ZMTJ06FcuWLXNaAlkQBDz55JOtWpwrsl/I5ewdIlIiyaHfqVMnLFmypMn3XnnllRsuyNXxQi4RKZnk5Lr99ttx/PjxJt/Lycm54YJcXQ2XVSYiBZN8pn/u3DmMGjUKw4cPR48ePaDRNITfpk2bkJmZ2aoFuhqupU9ESiY59Ddv3uz0LNyrZ/C4u5q6+tD3ZugTkfJIDv277roL7777bpPvecKFXPvwjh9Dn4gUSPKYfnOBDwAvvPDCDRWjBBzTJyIla9UpKFcO+7irKiOHd4hIuSQP74SFhTX7XnFx8Q0VowTV9WP6/roWLVBKRCQrycml0+mQlpbmeG2xWHDmzBls2LABjz76aKsW54o4pk9ESiY59DMyMnDPPfc02v7kk0/ikUceaZWiXJl9yiYflUhESiR5TL+pwAeADh064MiRIzdckKurqV+GgWf6RKREks/0m3okYmVlJXbs2AG12v2XJuCTs4hIySSH/sMPP4zg4GDHa5VKBb1ej2HDhiE7O7tVi3NFHNMnIiWTHPpjxoxBfn5+W9SiCEbHKpvu/1MNEbkfycnlyYEPNJzp80IuESmRqNAvLS3Fc889h+eeew4HDhxo9P6CBQtQWlra6sW5GqvQsLQy5+kTkRKJCv21a9fihRdewKVLl9CxY8dG7x88eBBjx47FmTNnWrs+l1I/sgOAY/pEpEyiQv+LL77AunXr8PLLLyMkJKTR+xs2bMDcuXORkZHR6gW6kvoVGADwISpEpEyikqu6uhpTp069Zps5c+Zg//79rVKUqzLXryKt81JDpVLJWwwRUQuICn0fHx9RB9PpdDdUjKuzD+/wLJ+IlEpUeplMJlit1mu2sVgsqKura5WiXFX9NVx4c+YOESmUqNCPi4vD3//+92u2SU9PR0JCQqsU5aquHN4hIlIiUfMO58+fj4kTJ2LEiBG47777MHDgQHTo0AFVVVU4cOAA1q1bBz8/P+Tm5rZ1vbIyc3iHiBROVOj7+voiPz8fixYtwpIlS2AwGKBSqSAIAgIDA/Hoo49i8eLF8Pb2but6ZWW22i7eajUMfSJSJtF3GPn6+mL58uVYtmwZDh065JizP3DgQI9YaA1oGN7x5pk+ESmU5NtKNRoNBg8e3Ba1uLyGC7kMfSJSJqaXBPYpm3woOhEpFUNfAs7eISKlY3pJ4Lg5i8sqE5FCMb0ksIe+D2/OIiKFYuhLwAu5RKR0TC8JTPXz9DmmT0RKxfSSgPP0iUjpmF4S2Id3eEcuESmVS6TX+vXrERkZiaioKERHR19zXf5169YhPj4ekyZNQmRkJKZPn44TJ060S50WnukTkcLJnl67du1CamoqPvzwQ2zfvh2zZ89GQkICKisrm2yfkpKCp556Cnl5edi5cyd8fX0xZcoUGI3GNq+VF3KJSOlkT6/MzEwkJiYiPDwcgC3UzWYzsrKymmw/depUxxLOarUajz/+OA4fPoyioqI2r9Uxps/hHSJSKNnTKy8vDyNHjnS8VqvVGDFiBLZs2dJk+08++cTptf2pXu1xpm8f3uGYPhEpleQF11pTWVkZDAYDunfv7rQ9ODgYBQUFoo7xww8/ICQkBOPGjWvyfaPR6PQPgsFgAGB7GpjJZBJdq8lkgqV+eEejEiTt647sf35P7weAfXEl9oWz9uoPKceXNfSrq6sBNH62rk6nc7x3LUajES+99BJWrlwJrVbbZJulS5ciIyOj0fbNmzfDz89PUr1mwXaGf3DfXvzn/B5J+7ord39wjhTsiwbsC2dt3R9i8tJO1tC3h+7VQzNGo1FUID/88MO49957kZSU1GybhQsXYt68eY7XBoMBoaGhiI+PR0BAgOhaTSYTXj+QBwAYcdsw3Dm0h+h93ZHJZEJubi7i4uKa/QfXU7AvGrAvnLVXf9hHMMSQNfS7dOmCwMBAlJSUOG0vLi5GWFjYNfdNS0uDn58fnn/++Wu20+l0jX6SAACtViv5L8FSf0eur7f0fd1VS/rRXbEvGrAvnLV1f0g5tuxXJGNjY1FYWOh4LQgCioqKMHny5Gb3yczMxKlTp7By5UoAQGFhodMx2oqJN2cRkcLJnl5paWnYuHEjjhw5AgDIzs6GRqNBamoqAGD8+PFIT093tH/rrbfwf//3f5gzZw6Kioqwe/dubNiwAXv37m3zWh3r6XNpZSJSKFmHdwBg1KhRyMrKQnJyMnx9faFWq5GTkwO9Xg/AdoHCPuZfWVmJxx57DFarFWPHjnU6zgcffNDmtVo4T5+IFE720AeApKSkZi/GXnnTlV6vh8Viaa+yGnGsvcM7colIoZheEvCOXCJSOqaXBLwjl4iUjuklgf2OXC+NSt5CiIhaiKEvgeNMX81uIyJlYnpJ4Ah9L57pE5EyMfRFEgQBFsEW9l480ycihWJ6iWS2Co7fc/YOESkV00skk/0qLji8Q0TKxdAXyWxpONPn8A4RKRXTSySnM31O2SQihWLoi2SqH9P3UqugUjH0iUiZGPoi2c/0eWMWESkZQ18k+5g+l2AgIiVjgolkD30vNc/0iUi5GPoimay24R2e6RORkjHBRDLxTJ+I3ABDXyQzL+QSkRtg6ItkP9PnEgxEpGRMMJHMVg7vEJHyMfRFsljtwzvsMiJSLiaYSPYpmxqe6RORgjH0ReLwDhG5A4a+SBZ76HP2DhEpGENfJPuCaxoutkZECsbQF6nhQi5Dn4iUi6EvUsPaO+wyIlIuJphIJo7pE5EbYOiLZF+GQcszfSJSMCaYSPbZO5ynT0RKxtAXybHKJod3iEjBGPoi2c/0+VB0IlIyhr5I5vopmxzeISIlY+iLZHaM6bPLiEi5mGAiWbj2DhG5AYa+SAx9InIHDH2RuMomEbkDhr5I9jN9NUOfiBSMoS8Sz/SJyB0w9EXimD4RuQOGvkiOKZu8OYuIFIyhL5J9wTUurUxESsYEE4nDO0TkDhj6Ipm5yiYRuQGGvkgNT85i6BORcrlE6K9fvx6RkZGIiopCdHQ09u/f36rtW4PJ/hAVjUt0GRFRi3jJXcCuXbuQmpqKwsJChIeHY/Xq1UhISMDBgweh1+tvuH1rMfNxiUTkBmQ/bc3MzERiYiLCw8MBACkpKTCbzcjKymqV9q3F5Ji9w9AnIuWSPfTz8vIwcuRIx2u1Wo0RI0Zgy5YtrdK+tZgdD1GRvcuIiFpM1uGdsrIyGAwGdO/e3Wl7cHAwCgoKbrg9ABiNRhiNRsdrg8EAADCZTDCZTKJrNZltZ/oQrJL2c1f2PmBfsC+uxL5w1l79IeX4soZ+dXU1AECn0zlt1+l0jvdupD0ALF26FBkZGY22b968GX5+fqJr7WBW4ya9Cr/uLULdCdG7ub3c3Fy5S3AZ7IsG7Atnbd0fzeVfU2QNfXvoXnkmbn/dVCBLbQ8ACxcuxLx58xyvDQYDQkNDER8fj4CAANG1xplMyM3NRVxcHLRarej93JWJ/eHAvmjAvnDWXv1hH8EQQ9bQ79KlCwIDA1FSUuK0vbi4GGFhYTfcHrD9FHD1TwYAoNVqW/SX0NL93BX7owH7ogH7wllb94eUY8t+VTI2NhaFhYWO14IgoKioCJMnT26V9kRE1ED20E9LS8PGjRtx5MgRAEB2djY0Gg1SU1MBAOPHj0d6erro9kRE1DzZb84aNWoUsrKykJycDF9fX6jVauTk5DhutKqurnYaw79eeyIiap7soQ8ASUlJSEpKavK9oqIiSe2JiKh5sg/vEBFR+2HoExF5EIY+EZEHYegTEXkQhj4RkQdh6BMReRCXmLLZngTBtkSylLUqANsaGtXV1TAYDLy9HOyPK7EvGrAvnLVXf9jzzJ5v1+JxoV9ZWQkACA0NlbkSIqLWVVlZicDAwGu2UQli/mlwI1arFWfPnoVer4dKJf4pWPbVOU+dOiVpdU53xf5owL5owL5w1l79IQgCKisrERISArX62qP2Hnemr1ar0atXrxbvHxAQwC/zFdgfDdgXDdgXztqjP653hm/HC7lERB6EoU9E5EEY+iLpdDo8++yzTT6QxROxPxqwLxqwL5y5Yn943IVcIiJPxjN9IiIPwtAnIvIgDH0iIg/C0Bdp/fr1iIyMRFRUFKKjo7F//365S2p3ixcvxrBhwxATE+P4NW3aNLnLald1dXVIS0uDl5cXTpw40ej9t99+GyNGjMC4ceOQmJiIM2fOtH+R7eRafTFjxgyMGTPG6bvy17/+VZ5C28G6desQHx+PSZMmITIyEtOnT3fqE0EQ8Nxzz+G2227DqFGjkJKSgkuXLslTrEDXtXPnTkGv1wu//vqrIAiCsGrVKqFnz56CwWCQubL29eyzzwr5+flylyGb48ePC2PGjBEefPBBAYBw/Phxp/c//fRToUePHkJpaakgCIKQkZEhDBs2TLBYLDJU27au1xepqamNtrkzrVYrbNq0SRAEQbBYLMJ//dd/CQMGDBBqa2sFQRCEl19+WRg6dKhQXV0tCIIgzJw5U/jDH/4gS6080xchMzMTiYmJCA8PBwCkpKTAbDYjKytL3sKoXV2+fBlr1qzBzJkzm3x/yZIlSE1NRdeuXQEAc+fOxb59+7Bx48b2LLNdXK8vPM3UqVORkJAAwHbX/+OPP47Dhw+jqKgIFosFmZmZ+Otf/wpfX18AwPz587Fhwwbs3bu33Wtl6IuQl5eHkSNHOl6r1WqMGDECW7ZskbEqam9DhgxBv379mnzv4sWL+Omnn5y+J4GBgejfv79bfk+u1Ree6JNPPnF67ePjAwAwGo3Ys2cPSktLnb4bgwYNgr+/vyzfDYb+dZSVlcFgMKB79+5O24ODg3H8+HGZqpLP//7v/yImJgbjxo1Damoqjh49KndJLsH+XeD3pMHSpUsRExOD8ePH47HHHkNJSYncJbWbH374ASEhIRg3bhyOHTsGwPm7oVKp0L17d1m+Gwz966iurgaARnfU6XQ6x3ueonfv3hg+fDi2bNmC7du346abbsKIESPc+mKlWPyeOOvfvz8mTJiAb775Bvn5+TAajRgzZgwuX74sd2ltzmg04qWXXsLKlSuh1Wpd7rvB0L8OPz8/ALa/yCsZjUbHe55i1qxZePLJJ+Hl5QW1Wo1FixbBx8cHb7zxhtylyY7fE2f/+Mc/8MADD0CtVkOr1eKVV17ByZMn8dFHH8ldWpt7+OGHce+99yIpKQmA6303GPrX0aVLFwQGBjb60bS4uBhhYWEyVeUaNBoN+vbtyyEewPFd4PekaQEBAQgKCnL770paWhr8/Pzw/PPPO7Y1990oKSmR5bvB0BchNjYWhYWFjteCIKCoqAiTJ0+Wsar2N3fu3Ebbzp49i969e8tQjWvp1KkThg8f7vQ9MRgM+PXXXz3uewI0/q4YjUaUlZW59XclMzMTp06dwsqVKwEAhYWFKCwsxNChQxEUFOT03Th48CCqqqrk+W7IMlFUYXbu3CkEBAQIv/32myAIgrBmzRqPnKfft29f4YsvvnC8fvfddwUfHx/h4MGDMlbV/vLz85udpx8SEiJcuHBBEARBeP755912nr5dc33h7e0tFBQUOF4/88wzQlBQkHD+/Pl2rrB9vPnmm8LgwYOFH374QSgoKBAKCgqEZ599Vvjggw8EQbDN04+IiHDM0589e7Zs8/Q97slZLTFq1ChkZWUhOTkZvr6+UKvVyMnJgV6vl7u0dvXCCy/g1VdfxSuvvIK6ujrodDps2bIFAwcOlLu0dlFXV4f4+HhUVFQAAJKTkxEaGuqYrjdt2jScP38ecXFx8PHxQadOnbBhw4brPr5Oia7XF8uXL3dc/6murkZQUBDy8/MRFBQkY9Vto7KyEo899hisVivGjh3r9N4HH3wAAHjyySdx+fJljBs3Dl5eXggPD8fq1avlKJdLKxMReRL3OwUhIqJmMfSJiDwIQ5+IyIMw9ImIPAhDn4jIgzD0iYg8CEOfiMiDMPTJo911113Q6XTo3bs35syZ49j+ww8/QKVS4bfffnNse+aZZ9CrVy9ERkbiwIEDzR7zzJkz6N69u6TVR1euXImBAweib9++12z3+eef4/PPPxd9XKKrMfTJo3311VeYMGEChg8fjn/961+O7Xl5eQCAb775xrFtyZIlGDZsGLZu3Ypbbrml2WP6+PhgwIABjqckifG3v/0NaWlp123H0KcbxdAnjxcbG4tvv/0WFovFse27777D7bff7gh/ADCZTDCZTPD397/m8bp06YJvv/0WnTt3brOaiVqKoU8eLzY2FhUVFSgqKgIA1NbWwmw2449//CPy8/NhX6lk586dGD16NADgpZdewrBhwxAdHY3o6Ghs374dgO2xiTExMfDx8XF6hnJxcTHuvPNO9O/fH3FxccjOzoZKpcKwYcPw73//26merKws3HHHHejXrx8yMzMd2xcsWIBNmzZh06ZNiImJwdSpU9uyW8hNccE18ngjR45EQEAA8vLyEBkZiR07dmDs2LGIjY1FWloa9uzZg4iICHzzzTeIjY3FO++8g/fffx8//vgjOnbsiB07dmDSpEk4dOgQ+vTpg61btzYam58xYwZ8fHxw6NAhqNVqx9LDr776KmJiYhztSkpKoFKp8PXXX2Pfvn0YOnQopk+fjptvvhkvvvgizp8/DwBO/6AQScEzffJ4Go3G8Wg/wDaOP2nSJNx2220IDAx0DPH8+OOPGDt2LF544QX893//Nzp27AgAuP3229GvXz+89957TR7/8OHDyMnJwdy5cx0rbj7++ONNthUEAQ888AAA28PHO3bsiD179rTmH5c8HEOfCLYhnu+//x51dXWOcNdoNIiOjkZeXh5qa2uhVqtRV1eHkydP4oMPPkBMTIzjl8lkQmVlZZPHPnToEAA4PSWpuYeJBAUFwcur4QfwgIAAGAyGVvyTkqfj8A4RbKFfXV2N3NxcaLVax0OsY2NjsWjRImzbtg233367o/38+fMxc+bMFn+eSqVqcrtGo2m0jaufU2vimT4RgKFDh6Jr167IyMjAhAkTHNtjY2NRWVmJZcuWITY2Fnq9Hr1798bhw4ed9l+7di0+/fTTJo9tf8jMsWPHHNtOnjzZojqvfCBLdXW104wjIjEY+kSwnXnHxMSgoKAAsbGxju1DhgxBt27dsHv3bowcORIAkJ6ejlWrVjmCu7S0FBkZGRgyZEiTxx4wYAASEhKwYsUKWK1WAMA777zTojqDgoJQXl4OAPjzn//sGDoiEouhT1QvNjYWAQEBjnAHGv4xmDBhgmOs/aGHHsLTTz+NKVOmICoqCtOnT8err76KAQMGOKZsFhcXIzMz0/GQ7KysLBiNRgwcOBBTpkxxPFZPq9U63s/MzERxcTHi4+MBAHfccYfjOGvWrAEAzJw5E8eOHUNUVBS6du2KwYMHt1v/kHvg4xKJ2kFpaanT82HPnj2Lnj174vTp0+jZs6eMlZGn4Zk+UTt49NFHsW3bNsfr119/HTExMQx8anecvUPUDqZOnYr58+ejQ4cOMBqN6NOnDz766CO5yyIPxOEdIiIPwuEdIiIPwtAnIvIgDH0iIg/C0Cci8iAMfSIiD8LQJyLyIAx9IiIPwtAnIvIgDH0iIg/y/wH1IXyEWB+f1AAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plot_cumulative_density(pd.read_csv('resources/grn_models/scglue.csv').weight)" ] }, { diff --git a/scripts/sbatch/figr.sh b/scripts/sbatch/r_script.sh similarity index 51% rename from scripts/sbatch/figr.sh rename to scripts/sbatch/r_script.sh index 5f59662ad..ebabda172 100644 --- a/scripts/sbatch/figr.sh +++ b/scripts/sbatch/r_script.sh @@ -1,5 +1,5 @@ #!/bin/bash -#SBATCH --job-name=figr +#SBATCH --job-name=r_model #SBATCH --time=48:00:00 #SBATCH --output=logs/%j.out #SBATCH --error=logs/%j.err @@ -9,4 +9,5 @@ #SBATCH --mem=250G -singularity run ../../images/figr Rscript src/methods/multi_omics/figr/script.R +# singularity run ../../images/figr Rscript src/methods/multi_omics/figr/script.R +singularity run ../../images/ppcor Rscript src/methods/single_omics/ppcor/script.R diff --git a/src/control_methods/pearson/script.py b/src/control_methods/pearson/script.py index ce895fb28..3fb6af02b 100644 --- a/src/control_methods/pearson/script.py +++ b/src/control_methods/pearson/script.py @@ -54,6 +54,7 @@ from util import create_corr_net par['causal'] = True +par['normalize'] = True net = create_corr_net(par) print('Output GRN') net.to_csv(par['prediction']) diff --git a/src/control_methods/positive_control/script.py b/src/control_methods/positive_control/script.py index de07531d0..e142fd3bc 100644 --- a/src/control_methods/positive_control/script.py +++ b/src/control_methods/positive_control/script.py @@ -24,6 +24,8 @@ print('Create causal corr net') par['causal'] = True +par['normalize'] = True + par['multiomics_rna'] = par['perturbation_data'] net = create_corr_net(par) diff --git a/src/exp_analysis/helper.py b/src/exp_analysis/helper.py index b27461d69..d455853e7 100644 --- a/src/exp_analysis/helper.py +++ b/src/exp_analysis/helper.py @@ -59,7 +59,7 @@ def calculate_feature_distance(sample: pd.DataFrame, control: pd.DataFrame, top_ distance_df = pd.DataFrame({'distance_raw':distance_raw, 'distance_rank':distance_rank}, index=entries_common) return distance_df -def cosine_similarity(nets_dict, col_name='source', weight_col='weight', figsize=(4, 4), title='Cosine Similarity', ax=None): +def cosine_similarity_net(nets_dict, col_name='source', weight_col='weight', figsize=(4, 4), title='Cosine Similarity', ax=None): from itertools import combinations from sklearn.metrics.pairwise import cosine_similarity import seaborn as sns @@ -99,14 +99,14 @@ def cosine_similarity(nets_dict, col_name='source', weight_col='weight', figsize sns.heatmap(cosine_sim_matrix, annot=True, cmap="coolwarm", xticklabels=nets_names, yticklabels=nets_names, ax=ax) ax.grid(False) - ax.set_title(title) + ax.set_title(title, pad=20) # Rotate x labels for readability plt.xticks(rotation=45, ha='right') return cosine_sim_matrix, fig -def jaccard_similarity(nets_dict, col_name='link', figsize=(4, 4), title='jaccard Similarity', ax=None): +def jaccard_similarity_net(nets_dict, col_name='link', figsize=(4, 4), title='jaccard Similarity', ax=None): from itertools import combinations import seaborn as sns import matplotlib.pyplot as plt diff --git a/src/utils/helper.py b/src/helper.py similarity index 93% rename from src/utils/helper.py rename to src/helper.py index f78edab56..c0ed50505 100644 --- a/src/utils/helper.py +++ b/src/helper.py @@ -4,7 +4,6 @@ import anndata as ad import numpy as np import scanpy as sc -from src.exp_analysis.helper import plot_cumulative_density import sys import subprocess import io @@ -41,6 +40,15 @@ def check_scores(par): df_all_n = (df_scores-df_scores.min(axis=0))/(df_scores.max(axis=0)-df_scores.min(axis=0)) df_scores['rank'] = df_all_n.mean(axis=1).rank(ascending=False).astype(int) return df_scores + +def calculate_scores(): + command = [ + "sbatch", + "scripts/sbatch/calculate_scores.sh" + ] + + # Print command to verify + subprocess.run(command) def run_consensus(par): @@ -55,7 +63,7 @@ def run_consensus(par): # Print command to verify subprocess.run(command) -def run_grn_inference_seqera(): +def run_grn_seqera(): # if False: # submit the job # !bash src/methods/multi_omics/celloracle/run.sh # if False: # get the results @@ -74,13 +82,22 @@ def run_grn_inference_seqera(): raise ValueError('define first') def run_grn_inference(): + # par = { + # 'methods': ['scglue'], + # 'models_dir': 'resources/grn_models/', + # 'multiomics_rna': 'resources/grn-benchmark/multiomics_rna.h5ad', + # 'multiomics_atac': 'resources/grn-benchmark/multiomics_atac.h5ad', + # 'num_workers': 20, + # 'mem': "120GB", + # 'time': "48:00:00" + # } par = { - 'methods': ['scglue'], + 'methods': ['scenicplus'], 'models_dir': 'resources/grn_models/', 'multiomics_rna': 'resources/grn-benchmark/multiomics_rna.h5ad', 'multiomics_atac': 'resources/grn-benchmark/multiomics_atac.h5ad', 'num_workers': 20, - 'mem': "120GB", + 'mem': "250GB", 'time': "48:00:00" } @@ -126,6 +143,8 @@ def run_grn_inference(): # Run sbatch command try: result = subprocess.run(['sbatch'] + full_tag + ['scripts/sbatch/grn_inference.sh', command], check=True, capture_output=True, text=True) + # result = subprocess.run(['bash'] + ['scripts/sbatch/grn_inference.sh', command], check=True, capture_output=True, text=True) + print(f"Job {method} submitted successfully.") print(result.stdout) # Print the standard output except subprocess.CalledProcessError as e: @@ -297,4 +316,9 @@ def elapsed_to_hours(elapsed_str): # Convert MaxRSS and MaxVMSize from KB to GB df_local['MaxRSS'] = df_local['MaxRSS'] / (1024 ** 2) # Convert KB to GB df_local['MaxVMSize'] = df_local['MaxVMSize'] / (1024 ** 2) # Convert KB to GB - return df_local \ No newline at end of file + return df_local + + +if __name__ == '__main__': + # run_grn_inference() + calculate_scores() \ No newline at end of file diff --git a/src/methods/multi_omics/scenicplus/main.py b/src/methods/multi_omics/scenicplus/main.py index a14b26651..e5d4153cd 100644 --- a/src/methods/multi_omics/scenicplus/main.py +++ b/src/methods/multi_omics/scenicplus/main.py @@ -703,17 +703,20 @@ def preprocess_rna(par): def snakemake_pipeline(par): import os snakemake_dir = os.path.join(par['temp_dir'], 'scplus_pipeline', 'Snakemake') - if os.path.exists(snakemake_dir): - import shutil - shutil.rmtree(snakemake_dir) # Replace 'snakemake_dir' with your directory path + # if os.path.exists(snakemake_dir): + # import shutil + # shutil.rmtree(snakemake_dir) os.makedirs(os.path.join(par['temp_dir'], 'scplus_pipeline'), exist_ok=True) os.makedirs(os.path.join(par['temp_dir'], 'scplus_pipeline', 'temp'), exist_ok=True) - subprocess.run(['scenicplus', 'init_snakemake', '--out_dir', os.path.join(par['temp_dir'], 'scplus_pipeline')]) - print('snake make initialized', flush=True) + + pipeline_dir = os.path.join(par['temp_dir'], 'scplus_pipeline') + if not os.path.exists(pipeline_dir): + subprocess.run(['scenicplus', 'init_snakemake', '--out_dir', pipeline_dir]) + print('snake make initialized', flush=True) # Load pipeline settings - with open(os.path.join(par['temp_dir'], 'scplus_pipeline', 'Snakemake', 'config', 'config.yaml'), 'r') as f: + with open(os.path.join(snakemake_dir, 'config', 'config.yaml'), 'r') as f: settings = yaml.safe_load(f) print('output_data:', settings['output_data'], flush=True) @@ -776,13 +779,28 @@ def snakemake_pipeline(par): # Run pipeline print('run snakemake ', flush=True) - + with contextlib.chdir(snakemake_dir): + # this fails to download atumatically so we do manually + if True: + url = "https://hgdownload.cse.ucsc.edu/goldenpath/hg38/bigZips/hg38.chrom.sizes" + response = requests.get(url) + with open('hg38.chrom.sizes', 'wb') as file: + file.write(response.content) + subprocess.run([ f'snakemake', - '--cores', str(par['num_workers']), - #'--unlock' + 'touch' ]) + + subprocess.run([ + f'snakemake', + '--unlock' + ]) + + subprocess.run([ + f'snakemake', + '--cores', str(par['num_workers'])]) def post_process(par): scplus_mdata = mudata.read(par['scplus_mdata']) diff --git a/src/methods/multi_omics/scenicplus/script.py b/src/methods/multi_omics/scenicplus/script.py index 876ecbe7b..5d6fca634 100644 --- a/src/methods/multi_omics/scenicplus/script.py +++ b/src/methods/multi_omics/scenicplus/script.py @@ -69,22 +69,22 @@ def main(par): par['MALLET_PATH'] = os.path.join(par['temp_dir'], 'Mallet-202108', 'bin', 'mallet') os.makedirs(par['atac_dir'], exist_ok=True) - print('------- download_databases -------') - download_databases(par) - print_memory_usage() - print('------- process_peak -------') - process_peak(par) - print_memory_usage() - print('------- run_cistopic -------') - run_cistopic(par) - print_memory_usage() - print('------- process_topics -------') - process_topics(par) - print_memory_usage() - print('------- preprocess_rna -------') - preprocess_rna(par) - print_memory_usage() - print('------- snakemake_pipeline -------') + # print('------- download_databases -------') + # download_databases(par) + # print_memory_usage() + # print('------- process_peak -------') + # process_peak(par) + # print_memory_usage() + # print('------- run_cistopic -------') + # run_cistopic(par) + # print_memory_usage() + # print('------- process_topics -------') + # process_topics(par) + # print_memory_usage() + # print('------- preprocess_rna -------') + # preprocess_rna(par) + # print_memory_usage() + # print('------- snakemake_pipeline -------') snakemake_pipeline(par) print_memory_usage() print('------- post_process -------') diff --git a/src/methods/multi_omics/scglue/main.py b/src/methods/multi_omics/scglue/main.py index 80180269a..598788822 100644 --- a/src/methods/multi_omics/scglue/main.py +++ b/src/methods/multi_omics/scglue/main.py @@ -264,12 +264,12 @@ def prune_grn(par): "--output", f"{par['temp_dir']}/pruned_grn.csv", "--top_n_targets", str(par['top_n_targets']), "--rank_threshold", str(par['rank_threshold']), - "--auc_threshold", "0", + "--auc_threshold", ".5", "--nes_threshold", "0", "--min_genes", "1", "--thresholds", "0.5", "0.7", "--top_n_regulators", "10", "50", "100", - "--max_similarity_fdr", "0.1", + "--max_similarity_fdr", "0.2", "--num_workers", f"{par['num_workers']}", "--cell_id_attribute", "obs_id", # be sure that obs_id is in obs and name is in var "--gene_attribute", "name" diff --git a/src/methods/single_omics/ppcor/script.R b/src/methods/single_omics/ppcor/script.R index e8c7e3a68..426f641ca 100644 --- a/src/methods/single_omics/ppcor/script.R +++ b/src/methods/single_omics/ppcor/script.R @@ -4,14 +4,16 @@ library(dplyr) ## VIASH START par <- list( - multiomics_rna = 'resources/grn-benchmark/multiomics_rna.h5ad', + multiomics_rna = 'resources/grn-benchmark/multiomics_rna_d0_hvg.h5ad', prediction = 'resources/grn_models/ppcor.csv', + tf_all = 'resources/prior/tf_all.csv', max_n_links = 50000 ) ## VIASH END print(par) -print(dim(par)) # input expression data +tf_names <- scan(par$tf_all, what = "", sep = "\n") + ad <- anndata::read_h5ad(par$multiomics_rna) inputExpr <- ad$X @@ -34,6 +36,9 @@ df <- df[order(df$weight, decreasing=TRUE),] # Remove self-interactions df_filtered <- df[df$source != df$target,] +# Filter to keep only connections where the source is a TF +df_filtered <- df_filtered %>% filter(source %in% tf_names) + # Reset index df_filtered$index = 1:nrow(df_filtered) diff --git a/src/metrics/regression_2/main.py b/src/metrics/regression_2/main.py index 82b3718d2..2aecca0b8 100644 --- a/src/metrics/regression_2/main.py +++ b/src/metrics/regression_2/main.py @@ -252,7 +252,7 @@ def static_approach( if n_features[n_features_dict[gene_name]] != 0: score = res['scores'][i]['avg-r2'] if clip_scores: - score = np.clip(score, 0, 1) + score = np.clip(score, 0, 1) r2.append(score) if len(r2)==0: raise ValueError('R2 score is empty') diff --git a/src/metrics/script_all.py b/src/metrics/script_all.py index 2c98506c1..dbcd95834 100644 --- a/src/metrics/script_all.py +++ b/src/metrics/script_all.py @@ -60,29 +60,33 @@ os.makedirs(par['scores_dir'], exist_ok=True) -for dataset_type in ['full', 'hvg']: - for apply_skeleton in [False, True]: - os.makedirs(f"{par['scores_dir']}/{dataset_type}/skeleton_{apply_skeleton}", exist_ok=True) - for layer in par['layers']: - par['layer'] = layer - i = 0 - for method in par['methods']: - print(method) - par['prediction'] = f"{par['models_dir']}/{dataset_type}/{method}.csv" - if not os.path.exists(par['prediction']): - print(f"{par['prediction']} doesnt exist. Skipped.") - continue - from regression_1.main import main - reg1 = main(par) - from regression_2.main import main - reg2 = main(par) - score = pd.concat([reg1, reg2], axis=1) - score.index = [method] - if i==0: - df_all = score - else: - df_all = pd.concat([df_all, score]) - df_all.to_csv(f"{par['scores_dir']}/{dataset_type}/skeleton_{apply_skeleton}/{layer}-{par['reg_type']}.csv") - print(df_all) - i+=1 +for max_n_links in [10000, 50000]: + par['max_n_links'] = max_n_links + for dataset_type in ['full', 'hvg']: + for apply_skeleton in [False, True]: + par['apply_skeleton'] = apply_skeleton + save_dir = f"{par['scores_dir']}/{dataset_type}/{max_n_links}/skeleton_{apply_skeleton}" + os.makedirs(save_dir, exist_ok=True) + for layer in par['layers']: + par['layer'] = layer + i = 0 + for method in par['methods']: + print(method) + par['prediction'] = f"{par['models_dir']}/{dataset_type}/{method}.csv" + if not os.path.exists(par['prediction']): + print(f"{par['prediction']} doesnt exist. Skipped.") + continue + from regression_1.main import main + reg1 = main(par) + from regression_2.main import main + reg2 = main(par) + score = pd.concat([reg1, reg2], axis=1) + score.index = [method] + if i==0: + df_all = score + else: + df_all = pd.concat([df_all, score]) + df_all.to_csv(f"{save_dir}/{layer}-{par['reg_type']}.csv") + print(df_all) + i+=1