Skip to content

Commit a0c8c39

Browse files
kennypavankennypavan
authored
Merge pull request #18 from ArpiarSaundersLab/robustness-validation
Example notebook added
2 parents 3940569 + 8ed8943 commit a0c8c39

1 file changed

Lines changed: 223 additions & 0 deletions

File tree

‎examples/recursieve_example.ipynb‎

Lines changed: 223 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,223 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# RecurSieve example\n",
8+
"\n",
9+
"Simulate a two-group dataset (2000 cells, 2000 genes) and run RecurSieve with default parameters."
10+
]
11+
},
12+
{
13+
"cell_type": "code",
14+
"execution_count": 1,
15+
"metadata": {},
16+
"outputs": [],
17+
"source": [
18+
"import numpy as np\n",
19+
"import pandas as pd\n",
20+
"import anndata as ad\n",
21+
"from recursieve import recursieve"
22+
]
23+
},
24+
{
25+
"cell_type": "markdown",
26+
"metadata": {},
27+
"source": [
28+
"## Simulate data\n",
29+
"\n",
30+
"Background counts are negative binomial with a per-cell library size. The first 10 genes are staged like the bimodal scenario in `manuscript/feature_comparision.py`: `control` is unimodal, `disease` is bimodal around the same center, so group means match and differential expression has little to find. In `disease`, cells switch between two states that raise one half of the staged genes and lower the other, so the genes move together. Each cell keeps its staged-block total, so library size is preserved."
31+
]
32+
},
33+
{
34+
"cell_type": "code",
35+
"execution_count": 2,
36+
"metadata": {},
37+
"outputs": [
38+
{
39+
"data": {
40+
"text/plain": [
41+
"AnnData object with n_obs × n_vars = 2000 × 2000\n",
42+
" obs: 'condition'"
43+
]
44+
},
45+
"execution_count": 2,
46+
"metadata": {},
47+
"output_type": "execute_result"
48+
}
49+
],
50+
"source": [
51+
"rng = np.random.default_rng(0)\n",
52+
"n_cells, n_genes, n_staged = 2000, 2000, 10\n",
53+
"center, sd, sep = 8.0, 1.5, 4.0\n",
54+
"\n",
55+
"groups = np.repeat([\"control\", \"disease\"], n_cells // 2)\n",
56+
"disease = groups == \"disease\"\n",
57+
"\n",
58+
"# background\n",
59+
"gene_means = rng.gamma(shape=0.6, scale=1.0, size=n_genes)\n",
60+
"gene_means[:n_staged] = 5.0\n",
61+
"library = rng.lognormal(mean=0.0, sigma=0.3, size=n_cells)\n",
62+
"mu = library[:, None] * gene_means[None, :]\n",
63+
"X = rng.poisson(rng.gamma(shape=10.0, scale=mu / 10.0))\n",
64+
"\n",
65+
"# staged pattern: unimodal in control, coordinated bimodal in disease\n",
66+
"pattern = rng.normal(center, sd, size=(n_cells, n_staged))\n",
67+
"state = rng.random(n_cells) < 0.5\n",
68+
"half = np.arange(n_staged) < n_staged // 2\n",
69+
"shift = np.where(state[:, None] == half[None, :], sep, -sep)\n",
70+
"pattern[disease] += shift[disease]\n",
71+
"\n",
72+
"# redistribute each cell's staged-block total by the pattern\n",
73+
"pattern = np.maximum(pattern, 0.0)\n",
74+
"prob = pattern / pattern.sum(axis=1, keepdims=True)\n",
75+
"totals = X[:, :n_staged].sum(axis=1)\n",
76+
"X[:, :n_staged] = [rng.multinomial(t, p) for t, p in zip(totals, prob)]\n",
77+
"\n",
78+
"adata = ad.AnnData(\n",
79+
" X=X.astype(np.float32),\n",
80+
" obs=pd.DataFrame({\"condition\": groups}, index=[f\"cell_{i}\" for i in range(n_cells)]),\n",
81+
" var=pd.DataFrame(index=[f\"gene_{i}\" for i in range(n_genes)]),\n",
82+
")\n",
83+
"staged = list(adata.var_names[:n_staged])\n",
84+
"adata"
85+
]
86+
},
87+
{
88+
"cell_type": "markdown",
89+
"metadata": {},
90+
"source": [
91+
"## Run RecurSieve"
92+
]
93+
},
94+
{
95+
"cell_type": "code",
96+
"execution_count": 3,
97+
"metadata": {},
98+
"outputs": [
99+
{
100+
"name": "stdout",
101+
"output_type": "stream",
102+
"text": [
103+
"Iteration 1\n",
104+
"Iteration 2\n",
105+
"flip_rate=0.1485\n",
106+
"Iteration 3\n",
107+
"flip_rate=0.0975\n",
108+
"Iteration 4\n",
109+
"flip_rate=0.0725\n",
110+
"Iteration 5\n",
111+
"flip_rate=0.0605\n",
112+
"Iteration 6\n",
113+
"flip_rate=0.0675\n",
114+
"Iteration 7\n",
115+
"flip_rate=0.1115\n",
116+
"Iteration 8\n",
117+
"flip_rate=0.1010\n",
118+
"Iteration 9\n",
119+
"flip_rate=0.1075\n",
120+
"Iteration 10\n",
121+
"flip_rate=0.0770\n",
122+
"Iteration 11\n",
123+
"flip_rate=0.0820\n",
124+
"Iteration 12\n",
125+
"flip_rate=0.1060\n",
126+
"Iteration 13\n",
127+
"flip_rate=0.0640\n",
128+
"Iteration 14\n",
129+
"flip_rate=0.0685\n",
130+
"Iteration 15\n",
131+
"flip_rate=0.0510\n",
132+
"Iteration 16\n",
133+
"flip_rate=0.0600\n",
134+
"Iteration 17\n",
135+
"flip_rate=0.0540\n",
136+
"Iteration 18\n",
137+
"flip_rate=0.0330\n",
138+
"Iteration 19\n",
139+
"flip_rate=0.0390\n",
140+
"Iteration 20\n",
141+
"flip_rate=0.0375\n",
142+
"Iteration 21\n",
143+
"flip_rate=0.0390\n",
144+
"Iteration 22\n",
145+
"flip_rate=0.0435\n",
146+
"Iteration 23\n",
147+
"flip_rate=0.0465\n",
148+
"Iteration 24\n",
149+
"flip_rate=0.0510\n",
150+
"Iteration 25\n",
151+
"flip_rate=0.0375\n",
152+
"Iteration 26\n",
153+
"flip_rate=0.0505\n",
154+
"Iteration 27\n",
155+
"flip_rate=0.0360\n",
156+
"Iteration 28\n",
157+
"flip_rate=0.0375\n",
158+
"stopping: no new flip-rate minimum in 10 iterations\n"
159+
]
160+
}
161+
],
162+
"source": [
163+
"model = recursieve(\n",
164+
" adata=adata,\n",
165+
" group1=\"control\",\n",
166+
" group2=\"disease\",\n",
167+
" field_name=\"condition\",\n",
168+
")"
169+
]
170+
},
171+
{
172+
"cell_type": "markdown",
173+
"metadata": {},
174+
"source": [
175+
"## Results\n",
176+
"\n",
177+
"`genes` is the full selected panel; `unique_gene_panel` is the subset not called by differential expression."
178+
]
179+
},
180+
{
181+
"cell_type": "code",
182+
"execution_count": 4,
183+
"metadata": {},
184+
"outputs": [
185+
{
186+
"name": "stdout",
187+
"output_type": "stream",
188+
"text": [
189+
"Panel: ['gene_4', 'gene_3', 'gene_0', 'gene_1', 'gene_2', 'gene_8', 'gene_5', 'gene_6', 'gene_7', 'gene_869', 'gene_56', 'gene_627', 'gene_1744', 'gene_1077', 'gene_317', 'gene_767', 'gene_1631', 'gene_914', 'gene_198', 'gene_1233', 'gene_1534', 'gene_553', 'gene_1779', 'gene_690', 'gene_933', 'gene_1105', 'gene_457', 'gene_578']\n",
190+
"Unique to RecurSieve: ['gene_3', 'gene_1', 'gene_2', 'gene_869', 'gene_56', 'gene_627', 'gene_1744', 'gene_1077', 'gene_317', 'gene_767', 'gene_1631', 'gene_914', 'gene_198', 'gene_1233', 'gene_1534', 'gene_1779', 'gene_690', 'gene_933', 'gene_1105', 'gene_457', 'gene_578']\n",
191+
"Staged genes recovered: ['gene_1', 'gene_2', 'gene_3']\n"
192+
]
193+
}
194+
],
195+
"source": [
196+
"print(\"Panel:\", model.genes)\n",
197+
"print(\"Unique to RecurSieve:\", model.unique_gene_panel)\n",
198+
"print(\"Staged genes recovered:\", sorted(set(model.unique_gene_panel) & set(staged)))"
199+
]
200+
}
201+
],
202+
"metadata": {
203+
"kernelspec": {
204+
"display_name": "env (3.12.3)",
205+
"language": "python",
206+
"name": "python3"
207+
},
208+
"language_info": {
209+
"codemirror_mode": {
210+
"name": "ipython",
211+
"version": 3
212+
},
213+
"file_extension": ".py",
214+
"mimetype": "text/x-python",
215+
"name": "python",
216+
"nbconvert_exporter": "python",
217+
"pygments_lexer": "ipython3",
218+
"version": "3.12.3"
219+
}
220+
},
221+
"nbformat": 4,
222+
"nbformat_minor": 4
223+
}

0 commit comments

Comments
 (0)