Skip to content

Commit 88bce0a

Browse files
jthortonjthorton
committed
fix error estimation, clean imports
1 parent 6fb8bd6 commit 88bce0a

1 file changed

Lines changed: 66 additions & 59 deletions

File tree

‎openmm_septop/septop_analysis.ipynb‎

Lines changed: 66 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -23,12 +23,9 @@
2323
"outputs": [],
2424
"source": [
2525
"import numpy as np\n",
26-
"import glob\n",
27-
"import json\n",
28-
"import csv\n",
2926
"import os\n",
3027
"import pathlib\n",
31-
"from typing import Literal, List\n",
28+
"from typing import List\n",
3229
"from gufe.tokenization import JSON_HANDLER\n",
3330
"import pandas as pd\n",
3431
"from openff.units import unit\n",
@@ -219,6 +216,30 @@
219216
{
220217
"cell_type": "code",
221218
"execution_count": 6,
219+
"id": "30faf828",
220+
"metadata": {},
221+
"outputs": [],
222+
"source": [
223+
"def error_std(r):\n",
224+
" \"\"\"\n",
225+
" Calculate the error of the estimate as the std of the repeats\n",
226+
" \"\"\"\n",
227+
" return np.std([v[0].m for v in r[\"overall\"]])\n",
228+
"\n",
229+
"def error_mbar(r):\n",
230+
" \"\"\"\n",
231+
" Calculate the error of the estimate using the reported MBAR errors.\n",
232+
"\n",
233+
" This also takes into account that repeats may have been run for this edge by using the average MBAR error\n",
234+
" \"\"\"\n",
235+
" complex_errors = [x[1].m for x in r[\"complex\"]]\n",
236+
" solvent_errors = [x[1].m for x in r[\"solvent\"]]\n",
237+
" return math.sqrt(np.mean(complex_errors)**2 + np.mean(solvent_errors)**2)"
238+
]
239+
},
240+
{
241+
"cell_type": "code",
242+
"execution_count": 7,
222243
"id": "411fe035-2ae2-4f98-9bab-19764af724ce",
223244
"metadata": {},
224245
"outputs": [],
@@ -232,15 +253,12 @@
232253
" Dict of legs to write out.\n",
233254
" \"\"\"\n",
234255
" data = []\n",
256+
" # check the type of error which should be used based on the number of repeats\n",
257+
" repeats = {len(v['overall']) for v in legs.values()}\n",
258+
" error_func = error_mbar if 1 in repeats else error_std\n",
235259
" for ligpair, results in sorted(legs.items()):\n",
236260
" ddg = np.mean([v[0].m for v in results[\"overall\"]])\n",
237-
" # Use standard deviation as error when more than 1 repeat\n",
238-
" if len(results[\"overall\"]) > 1:\n",
239-
" error = np.std([v[0].m for v in results[\"overall\"]])\n",
240-
" if len(results[\"overall\"]) == 1:\n",
241-
" complex_error = results[\"complex\"][0][1].m\n",
242-
" solvent_error = results[\"solvent\"][0][1].m\n",
243-
" error = math.sqrt(complex_error**2 + solvent_error**2)\n",
261+
" error = error_func(results)\n",
244262
" m, u = format_estimate_uncertainty(ddg, error, unc_prec=2)\n",
245263
" data.append((ligpair[0], ligpair[1], m, u))\n",
246264
"\n",
@@ -250,7 +268,7 @@
250268
},
251269
{
252270
"cell_type": "code",
253-
"execution_count": 7,
271+
"execution_count": 8,
254272
"id": "213f1c7b-185f-403b-a98c-22d31a5e60e5",
255273
"metadata": {},
256274
"outputs": [],
@@ -263,9 +281,6 @@
263281
" legs : dict\n",
264282
" Dict of legs to write out.\n",
265283
" \"\"\"\n",
266-
" import networkx as nx\n",
267-
" import numpy as np\n",
268-
" from cinnabar.stats import mle\n",
269284
"\n",
270285
" DDGs = _generate_ddg(legs)\n",
271286
" fe_results = []\n",
@@ -296,7 +311,7 @@
296311
},
297312
{
298313
"cell_type": "code",
299-
"execution_count": 8,
314+
"execution_count": 9,
300315
"id": "24dfaef3-92ac-4602-85d4-400b8e6bd4cf",
301316
"metadata": {},
302317
"outputs": [],
@@ -322,15 +337,15 @@
322337
},
323338
{
324339
"cell_type": "code",
325-
"execution_count": 9,
340+
"execution_count": 10,
326341
"id": "3b81c26b-4e3d-4bca-bfb1-4f72c75ec1bf",
327342
"metadata": {},
328343
"outputs": [
329344
{
330345
"name": "stderr",
331346
"output_type": "stream",
332347
"text": [
333-
"/Users/hannahbaumann/miniforge3/envs/septop/lib/python3.12/site-packages/Bio/Application/__init__.py:39: BiopythonDeprecationWarning: The Bio.Application modules and modules relying on it have been deprecated.\n",
348+
"/Users/joshua/mambaforge/envs/septop/lib/python3.12/site-packages/Bio/Application/__init__.py:39: BiopythonDeprecationWarning: The Bio.Application modules and modules relying on it have been deprecated.\n",
334349
"\n",
335350
"Due to the on going maintenance burden of keeping command line application\n",
336351
"wrappers up to date, we have decided to deprecate and eventually remove these\n",
@@ -360,7 +375,7 @@
360375
},
361376
{
362377
"cell_type": "code",
363-
"execution_count": 10,
378+
"execution_count": 11,
364379
"id": "46996a74-709c-41f2-ac39-0f77fb33371e",
365380
"metadata": {},
366381
"outputs": [],
@@ -371,8 +386,8 @@
371386
},
372387
{
373388
"cell_type": "code",
374-
"execution_count": 11,
375-
"id": "d1a6ad61-1e5a-4d8a-9067-9ed428ef145c",
389+
"execution_count": 12,
390+
"id": "b2b6fdf0",
376391
"metadata": {},
377392
"outputs": [
378393
{
@@ -414,22 +429,22 @@
414429
" <th>1</th>\n",
415430
" <td>1</td>\n",
416431
" <td>7a</td>\n",
417-
" <td>0.65</td>\n",
418-
" <td>0.86</td>\n",
432+
" <td>0.6</td>\n",
433+
" <td>1.5</td>\n",
419434
" </tr>\n",
420435
" <tr>\n",
421436
" <th>2</th>\n",
422437
" <td>1</td>\n",
423438
" <td>7b</td>\n",
424-
" <td>0.15</td>\n",
425-
" <td>0.42</td>\n",
439+
" <td>0.1</td>\n",
440+
" <td>1.5</td>\n",
426441
" </tr>\n",
427442
" <tr>\n",
428443
" <th>3</th>\n",
429444
" <td>7a</td>\n",
430445
" <td>7b</td>\n",
431446
" <td>1.9</td>\n",
432-
" <td>2.1</td>\n",
447+
" <td>1.5</td>\n",
433448
" </tr>\n",
434449
" </tbody>\n",
435450
"</table>\n",
@@ -438,12 +453,12 @@
438453
"text/plain": [
439454
" ligand_i ligand_j DDG(i->j) (kcal/mol) uncertainty (kcal/mol)\n",
440455
"0 1 25 2.0 1.6\n",
441-
"1 1 7a 0.65 0.86\n",
442-
"2 1 7b 0.15 0.42\n",
443-
"3 7a 7b 1.9 2.1"
456+
"1 1 7a 0.6 1.5\n",
457+
"2 1 7b 0.1 1.5\n",
458+
"3 7a 7b 1.9 1.5"
444459
]
445460
},
446-
"execution_count": 11,
461+
"execution_count": 12,
447462
"metadata": {},
448463
"output_type": "execute_result"
449464
}
@@ -464,7 +479,7 @@
464479
},
465480
{
466481
"cell_type": "code",
467-
"execution_count": 12,
482+
"execution_count": 13,
468483
"id": "22e25226-0073-40a4-9cbc-da802a29fc25",
469484
"metadata": {},
470485
"outputs": [],
@@ -475,8 +490,8 @@
475490
},
476491
{
477492
"cell_type": "code",
478-
"execution_count": 13,
479-
"id": "519a7e3a-f092-4dc3-ab61-5697091954cd",
493+
"execution_count": 14,
494+
"id": "cda55931",
480495
"metadata": {},
481496
"outputs": [
482497
{
@@ -509,40 +524,40 @@
509524
" <tr>\n",
510525
" <th>0</th>\n",
511526
" <td>1</td>\n",
512-
" <td>-0.636553</td>\n",
513-
" <td>0.462091</td>\n",
527+
" <td>-0.675</td>\n",
528+
" <td>0.664267</td>\n",
514529
" </tr>\n",
515530
" <tr>\n",
516531
" <th>1</th>\n",
517532
" <td>25</td>\n",
518-
" <td>1.363447</td>\n",
519-
" <td>1.222100</td>\n",
533+
" <td>1.325</td>\n",
534+
" <td>1.311964</td>\n",
520535
" </tr>\n",
521536
" <tr>\n",
522537
" <th>2</th>\n",
523538
" <td>7a</td>\n",
524-
" <td>-0.319831</td>\n",
525-
" <td>0.720921</td>\n",
539+
" <td>-0.875</td>\n",
540+
" <td>0.903466</td>\n",
526541
" </tr>\n",
527542
" <tr>\n",
528543
" <th>3</th>\n",
529544
" <td>7b</td>\n",
530-
" <td>-0.407063</td>\n",
531-
" <td>0.535312</td>\n",
545+
" <td>0.225</td>\n",
546+
" <td>0.903466</td>\n",
532547
" </tr>\n",
533548
" </tbody>\n",
534549
"</table>\n",
535550
"</div>"
536551
],
537552
"text/plain": [
538553
" ligand DG (kcal/mol) uncertainty (kcal/mol)\n",
539-
"0 1 -0.636553 0.462091\n",
540-
"1 25 1.363447 1.222100\n",
541-
"2 7a -0.319831 0.720921\n",
542-
"3 7b -0.407063 0.535312"
554+
"0 1 -0.675 0.664267\n",
555+
"1 25 1.325 1.311964\n",
556+
"2 7a -0.875 0.903466\n",
557+
"3 7b 0.225 0.903466"
543558
]
544559
},
545-
"execution_count": 13,
560+
"execution_count": 14,
546561
"metadata": {},
547562
"output_type": "execute_result"
548563
}
@@ -561,7 +576,7 @@
561576
},
562577
{
563578
"cell_type": "code",
564-
"execution_count": 14,
579+
"execution_count": 15,
565580
"id": "8b2c1dd8-ffa3-4585-94a7-4a1ed1454f30",
566581
"metadata": {},
567582
"outputs": [],
@@ -572,7 +587,7 @@
572587
},
573588
{
574589
"cell_type": "code",
575-
"execution_count": 15,
590+
"execution_count": 16,
576591
"id": "08b72901-9c71-460b-b5da-9fb4a35e07f7",
577592
"metadata": {},
578593
"outputs": [
@@ -975,27 +990,19 @@
975990
"35 0.0 "
976991
]
977992
},
978-
"execution_count": 15,
993+
"execution_count": 16,
979994
"metadata": {},
980995
"output_type": "execute_result"
981996
}
982997
],
983998
"source": [
984999
"df_raw"
9851000
]
986-
},
987-
{
988-
"cell_type": "code",
989-
"execution_count": null,
990-
"id": "ae092ba5-3263-4e39-bc7d-82e8f1e1f58c",
991-
"metadata": {},
992-
"outputs": [],
993-
"source": []
9941001
}
9951002
],
9961003
"metadata": {
9971004
"kernelspec": {
998-
"display_name": "Python 3 (ipykernel)",
1005+
"display_name": "septop",
9991006
"language": "python",
10001007
"name": "python3"
10011008
},
@@ -1009,7 +1016,7 @@
10091016
"name": "python",
10101017
"nbconvert_exporter": "python",
10111018
"pygments_lexer": "ipython3",
1012-
"version": "3.12.10"
1019+
"version": "3.12.11"
10131020
}
10141021
},
10151022
"nbformat": 4,

0 commit comments

Comments
 (0)