{ "nbformat": 4, "nbformat_minor": 5, "metadata": { "colab": {"name": "PX1 malaria susceptibility partial reanalysis", "provenance": []}, "kernelspec": {"display_name": "Python 3", "name": "python3"}, "language_info": {"name": "python"} }, "cells": [ { "cell_type": "markdown", "metadata": {"id": "intro"}, "source": [ "# PX1/PIN malaria susceptibility reanalysis\n", "\n", "**Status: original experiment partially reproduced.** This notebook downloads two small public XLSX inputs, validates them, reruns the unadjusted published-data comparison and saves aggregate results. It is not the article's exact depth-filtered or COI-adjusted model. No patient/sample identifiers are written to results.\n", "\n", "Run in Google Colab with Runtime → Run all. Data use is subject to the source terms recorded in `data_manifest.json`; verify reuse rights before redistribution." ] }, { "cell_type": "code", "execution_count": null, "metadata": {"id": "deps"}, "outputs": [], "source": ["%pip -q install pandas scipy openpyxl matplotlib requests"] }, { "cell_type": "code", "execution_count": null, "metadata": {"id": "fetch"}, "outputs": [], "source": [ "from pathlib import Path\n", "from hashlib import sha256\n", "from io import BytesIO\n", "import json, requests\n", "import numpy as np\n", "import pandas as pd\n", "from scipy.stats import mannwhitneyu\n", "import matplotlib.pyplot as plt\n", "\n", "PX1_URL = 'https://media.springernature.com/original/springer-static/esm/art:10.1038%2Fs41591-026-04590-5/MediaObjects/41591_2026_4590_MOESM5_ESM.xlsx'\n", "WHO_URL = 'https://cdn.who.int/media/docs/default-source/malaria/drug-resistance/compendium-of-molecular-markers-for-antimalarial-drug-resistance.xlsx?download=true&sfvrsn=ef5bc699_8'\n", "CACHE = Path('/content/px1-malaria-data'); CACHE.mkdir(exist_ok=True)\n", "def fetch(url, filename, cap=2_000_000):\n", " dest = CACHE / filename\n", " if not dest.exists():\n", " response = requests.get(url, timeout=60); response.raise_for_status()\n", " if len(response.content) > cap or not response.content.startswith(b'PK'):\n", " raise ValueError(f'Unexpected/oversize XLSX response for {filename}')\n", " dest.write_bytes(response.content)\n", " data = dest.read_bytes()\n", " if len(data) > cap or not data.startswith(b'PK'): raise ValueError(f'Invalid cached XLSX: {filename}')\n", " print(filename, len(data), sha256(data).hexdigest())\n", " return data\n", "px1_bytes = fetch(PX1_URL, 'px1-source-data.xlsx')\n", "who_bytes = fetch(WHO_URL, 'who-markers.xlsx')\n", "px1 = pd.read_excel(BytesIO(px1_bytes), sheet_name='Source_dataset1')\n", "who = pd.read_excel(BytesIO(who_bytes), sheet_name='Artemisinins (Pf)')\n", "required = {'sample_name','site','region','year','parasitemia','px1_haplotype','wsaf','k13_genotype','DHA','Lumefantrine','Mefloquine','RSA'}\n", "missing = required - set(px1.columns)\n", "if missing: raise ValueError(f'Missing source columns: {sorted(missing)}')\n", "if px1.sample_name.isna().any() or px1.sample_name.astype(str).str.strip().eq('').any(): raise ValueError('Blank sample IDs')\n", "if px1.sample_name.duplicated().any(): raise ValueError('Duplicate sample IDs')\n", "print('Source rows:', len(px1), '| unique IDs:', px1.sample_name.nunique())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {"id": "analysis"}, "outputs": [], "source": [ "ASSAYS = [('DHA','Dihydroartemisinin (DHA)','nM IC50'), ('Lumefantrine','Lumefantrine','nM IC50'), ('Mefloquine','Mefloquine','nM IC50'), ('RSA','Ring-stage survival assay','% survival')]\n", "allowed_k13 = {'WT','A675V','C469Y','C469F','R561H','T348I','A675V_T348I'}\n", "d = px1.copy()\n", "d['wsaf'] = pd.to_numeric(d.wsaf, errors='coerce')\n", "for assay, _, _ in ASSAYS: d[assay] = pd.to_numeric(d[assay], errors='coerce')\n", "d['k13_genotype'] = d.k13_genotype.fillna('').astype(str).str.strip()\n", "d['px1_haplotype'] = d.px1_haplotype.fillna('').astype(str).str.strip()\n", "eligible = d[d.px1_haplotype.isin(['PIN','LMD']) & (d.wsaf > .9) & d.k13_genotype.isin(allowed_k13)].copy()\n", "print({'raw':len(d), 'PIN_or_LMD':int(d.px1_haplotype.isin(['PIN','LMD']).sum()), 'WSAF_gt_0.9':int((d.px1_haplotype.isin(['PIN','LMD']) & (d.wsaf > .9)).sum()), 'K13_eligible':len(eligible)})\n", "\n", "def to_int32(value):\n", " value = int(value) & 0xffffffff\n", " return value - 0x100000000 if value >= 0x80000000 else value\n", "\n", "def cathedral_rng(seed):\n", " state = to_int32(seed)\n", " def random():\n", " nonlocal state\n", " unsigned = state & 0xffffffff\n", " state = to_int32(unsigned ^ ((unsigned << 13) & 0xffffffff))\n", " unsigned = state & 0xffffffff\n", " state = to_int32(unsigned ^ (unsigned >> 17))\n", " unsigned = state & 0xffffffff\n", " state = to_int32(unsigned ^ ((unsigned << 5) & 0xffffffff))\n", " return (state & 0xffffffff) / 4294967296.0\n", " return random\n", "\n", "def bootstrap_median_diff(a, b, seed, reps=4000):\n", " random = cathedral_rng(seed)\n", " a = np.asarray(a, dtype=float); b = np.asarray(b, dtype=float)\n", " draws = np.empty(reps)\n", " for i in range(reps):\n", " aa = np.array([a[int(random() * len(a))] for _ in range(len(a))])\n", " bb = np.array([b[int(random() * len(b))] for _ in range(len(b))])\n", " draws[i] = np.median(aa) - np.median(bb)\n", " return np.quantile(draws, [.025,.975]).tolist()\n", "\n", "def rank_sum_p(a, b):\n", " combined = sorted([(float(v), 0) for v in a] + [(float(v), 1) for v in b])\n", " rank_sum_a = 0.0; tie_correction = 0; i = 0\n", " while i < len(combined):\n", " j = i + 1\n", " while j < len(combined) and combined[j][0] == combined[i][0]: j += 1\n", " size = j - i; rank = ((i + 1) + j) / 2\n", " rank_sum_a += sum(rank for k in range(i, j) if combined[k][1] == 0)\n", " tie_correction += size ** 3 - size; i = j\n", " n1, n2 = len(a), len(b); n = n1 + n2\n", " u = rank_sum_a - n1 * (n1 + 1) / 2; mean = n1 * n2 / 2\n", " variance = (n1 * n2 / 12) * ((n + 1) - tie_correction / (n * (n - 1)))\n", " if variance <= 0: return float(u), 1.0\n", " z = (abs(u - mean) - 0.5) / np.sqrt(variance)\n", " x = abs(z) / np.sqrt(2); t = 1 / (1 + 0.5 * x); poly = 0.17087277\n", " for coefficient in [-0.82215223,1.48851587,-1.13520398,0.27886807,-0.18628806,0.09678418,0.37409196,1.00002368]: poly = coefficient + t * poly\n", " return float(u), float(t * np.exp(-x*x - 1.26551223 + t*poly))\n", "\n", "def analyze_subset(frame, subset, seed_offset=0):\n", " out=[]\n", " for idx,(assay,label,unit) in enumerate(ASSAYS):\n", " valid=frame[frame[assay].notna()]\n", " a=valid.loc[valid.px1_haplotype=='PIN',assay].to_numpy()\n", " b=valid.loc[valid.px1_haplotype=='LMD',assay].to_numpy()\n", " if len(a)>=2 and len(b)>=2:\n", " u,p=rank_sum_p(a,b)\n", " else: p=None\n", " out.append({'subset':subset,'assay':assay,'label':label,'unit':unit,'n':len(valid),'pin_n':len(a),'pin_median':float(np.median(a)) if len(a) else None,'pin_q1':float(np.quantile(a,.25)) if len(a) else None,'pin_q3':float(np.quantile(a,.75)) if len(a) else None,'lmd_n':len(b),'lmd_median':float(np.median(b)) if len(b) else None,'lmd_q1':float(np.quantile(b,.25)) if len(b) else None,'lmd_q3':float(np.quantile(b,.75)) if len(b) else None,'median_difference_pin_minus_lmd':float(np.median(a)-np.median(b)) if len(a) and len(b) else None,'bootstrap_95':bootstrap_median_diff(a,b,20261029+idx*101+seed_offset) if len(a) and len(b) else None,'p_value_approx':p,'u_statistic':u if len(a)>=2 and len(b)>=2 else None,'missing_outcome_excluded':len(frame)-len(valid)})\n", " validp=[(i,r['p_value_approx']) for i,r in enumerate(out) if r['p_value_approx'] is not None]\n", " ranked=sorted(validp,key=lambda x:x[1]); prev=1.0\n", " for rank in range(len(ranked)-1,-1,-1):\n", " i,p=ranked[rank]; prev=min(prev,p*len(ranked)/(rank+1)); out[i]['bh_q_value']=min(1.0,prev)\n", " for r in out: r.setdefault('bh_q_value',None)\n", " return out\n", "results = analyze_subset(eligible,'All eligible PIN/LMD records') + analyze_subset(eligible[eligible.k13_genotype=='WT'],'K13 wild type',1000)\n", "csv_rows=[]\n", "for r in results:\n", " ci=r['bootstrap_95'] or [None,None]\n", " csv_rows.append({'subset':r['subset'],'assay':r['assay'],'label':r['label'],'unit':r['unit'],'n':r['n'],'pin_n':r['pin_n'],'pin_median':r['pin_median'],'pin_q1':r['pin_q1'],'pin_q3':r['pin_q3'],'lmd_n':r['lmd_n'],'lmd_median':r['lmd_median'],'lmd_q1':r['lmd_q1'],'lmd_q3':r['lmd_q3'],'median_difference_pin_minus_lmd':r['median_difference_pin_minus_lmd'],'bootstrap_95_low':ci[0],'bootstrap_95_high':ci[1],'p_value':r['p_value_approx'],'bh_q_value':r['bh_q_value'],'u_statistic':r['u_statistic'],'missing_outcome_excluded':r['missing_outcome_excluded'],'test':'Mann-Whitney U, tie-corrected normal approximation with continuity correction'})\n", "csv_columns=['subset','assay','label','unit','n','pin_n','pin_median','pin_q1','pin_q3','lmd_n','lmd_median','lmd_q1','lmd_q3','median_difference_pin_minus_lmd','bootstrap_95_low','bootstrap_95_high','p_value','bh_q_value','u_statistic','missing_outcome_excluded','test']\n", "pd.DataFrame(csv_rows,columns=csv_columns).to_csv('phenotype-comparisons.csv',index=False)\n", "pd.DataFrame(results).round(4)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {"id": "who-plot"}, "outputs": [], "source": [ "markers = who[(who['Gene(s)'].astype(str).str.strip().str.lower()=='pfk13') & who['Alteration(s)'].astype(str).str.upper().isin(['C469Y','A675V'])][['Gene(s)','Alteration(s)','Classification']].drop_duplicates()\n", "if len(markers)!=2: raise ValueError('Expected WHO K13 C469Y and A675V classifications; check source version/schema')\n", "display(markers)\n", "\n", "core=[r for r in results if r['subset']=='All eligible PIN/LMD records' and r['assay'] in ['DHA','Lumefantrine']]\n", "fig,ax=plt.subplots(figsize=(8,3.5))\n", "for y,r in enumerate(core):\n", " lo,hi=r['bootstrap_95']; mid=r['median_difference_pin_minus_lmd']\n", " ax.errorbar(mid,y,xerr=[[mid-lo],[hi-mid]],fmt='o',capsize=4,label=r['label'])\n", "ax.axvline(0,color='gray',lw=1,ls='--'); ax.set_yticks(range(len(core)),[r['label'] for r in core]); ax.set_xlabel('PIN minus LMD median (nM IC50); bootstrap 95% interval'); ax.set_title('Unadjusted source-table comparisons')\n", "plt.tight_layout(); plt.savefig('core-phenotype-contrast-colab.png',dpi=160); plt.show()\n", "\n", "summary={'schemaVersion':'cathedral.malaria-resistance.result.v1','outcomeClass':'Original experiment partially reproduced','rawRows':len(px1),'eligibleRows':len(eligible),'articleReportedDepthQualifiedN':465,'readDepthAvailableInWorkbook':False,'results':results,'whoK13Markers':markers.to_dict(orient='records'),'limits':['The >=25x article inclusion threshold could not be applied because numeric row-level depth is absent.','The article exact cohort and COI-adjusted mixed-effects model are not reproduced.','Ex-vivo outcomes are not clinical treatment outcomes; association is not causation.','PX1/PIN is not WHO-classified as a validated marker; WHO classifications here apply to K13 variants.']}\n", "Path('px1-malaria-summary-colab.json').write_text(json.dumps(summary,indent=2,allow_nan=False))\n", "print('Saved aggregate CSV, figure and JSON summary in the Colab working directory.')" ] }, { "cell_type": "markdown", "metadata": {"id": "limitations"}, "source": [ "## Interpretation\n", "\n", "Higher IC50 is interpreted as lower ex-vivo susceptibility in the assay. This analysis tests distributional association only. It neither establishes PX1 causality nor predicts clinical treatment outcomes. The results differ in sample counts from the paper because the downloaded workbook does not expose all exact depth/COI joins used in the publication. Check [`../results/report.md`](../results/report.md) and [`../data_manifest.json`](../data_manifest.json) for the full local report and provenance." ] } ] }