zhao 2021 hierarchical

import%20marimo%0A%0A__generated_with%20%3D%20%220.23.9%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%0A%20%20%20%20return%20(mo%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Multi-donor%20differential%20expression%20with%20a%20*crossed%20hierarchical*%20%60scribe%60%20model%0A%0A%20%20%20%20Most%20differential-expression%20tutorials%20compare%20two%20conditions%20in%20a%20single%20sample.%20Real%20perturbation%20experiments%20are%20rarely%20that%20tidy%3A%20the%20same%20drug%20is%20applied%20to%20**several%20donors**%2C%20each%20with%20its%20own%20baseline%2C%20and%20we%20want%20the%20**population-level%20treatment%20effect**%20%E2%80%94%20not%20an%20artifact%20of%20which%20donor%20happened%20to%20dominate%20one%20arm.%0A%0A%20%20%20%20In%20this%20tutorial%20we%20fit%20**one%20joint%20model**%20to%20a%20crossed%20*donor%20%C3%97%20condition*%20design%20and%20read%20off%20the%20treatment%20effect%20while%20explicitly%20accounting%20for%20donor-to-donor%20heterogeneity.%20The%20running%20example%20is%20the%20**Zhao%20et%20al.%20(2021)**%20dataset%2C%20where%20peripheral%20cells%20from%20**seven%20donors**%20were%20profiled%20at%20baseline%20(%60control%60)%20and%20after%20treatment%20with%20**panobinostat**%2C%20a%20pan-**HDAC%20inhibitor**.%20Because%20HDAC%20inhibitors%20have%20a%20well-characterized%20transcriptional%20fingerprint%20(histone%20genes%20up%2C%20heat-shock%20and%20p53-stress%20programs%20up%2C%20MYC%20and%20the%20cell%20cycle%20down)%2C%20we%20can%20check%20the%20recovered%20signature%20against%20textbook%20biology.%0A%0A%20%20%20%20We%20will%20not%20be%20shy%20about%20the%20math%3A%20the%20whole%20point%20of%20a%20generative%20model%20is%20that%20*every*%20quantity%20%E2%80%94%20the%20treatment%20effect%2C%20the%20donor%20deviations%2C%20the%20differential-expression%20call%20%E2%80%94%20is%20a%20parameter%20with%20a%20posterior%2C%20so%20knowing%20what%20is%20being%20computed%20is%20what%20lets%20you%20trust%20(or%20question)%20the%20answer.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20The%20design%2C%20and%20why%20a%20*joint*%20model%0A%0A%20%20%20%20Seven%20donors%2C%20two%20conditions%20each%2C%20gives%20a%20%247%20%5Ctimes%202%24%20grid%20of%20**groups**%0A%20%20%20%20(we%20call%20each%20present%20combination%20a%20*leaf*).%20Three%20ways%20to%20analyze%20it%3A%0A%0A%20%20%20%201.%20**Pool%20everything**%2C%20ignore%20the%20donor%20labels%2C%20compare%20control%20vs.%20%20panobinostat.%20Fast%2C%20but%20a%20single%20donor%20with%20many%20cells%20or%20an%20extreme%20baseline%20can%20drive%20the%20result%2C%20and%20the%20donor%20structure%20is%20thrown%20away.%0A%20%20%20%202.%20**Fit%20each%20donor%20separately**%2C%20do%20seven%20pairwise%20comparisons%2C%20average%20the%20answers.%20Honest%20about%20pairing%2C%20but%20every%20fit%20sees%20only%20its%20own%20cells%20%E2%80%94%20no%20sharing%20of%20statistical%20strength%2C%20and%20no%20coherent%20population%20estimand.%0A%20%20%20%203.%20**One%20joint%20model**%20in%20which%20a%20cell's%20mean%20expression%20is%20its%20donor's%20baseline%20*plus*%20a%20shared%20treatment%20effect%20*plus*%20that%20donor's%20own%20deviation.%20This%20is%20what%20we%20do%20here.%20It%20keeps%20the%20pairing%20(control%20and%20treated%20cells%20from%20the%20same%20donor%20are%20linked%20through%20that%20donor's%20parameters)%2C%20shares%20strength%20across%20donors%2C%20and%20%E2%80%94%20crucially%20%E2%80%94%20lets%20us%20*separate*%20the%20treatment%20effect%20we%20care%20about%20from%20the%20donor%20heterogeneity%20we%20want%20to%20average%20over.%0A%0A%20%20%20%20The%20rest%20of%20the%20notebook%20builds%20option%203%20and%20shows%20how%20to%20interrogate%20it.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20%23%20Import%20basic%20packages%0A%20%20%20%20from%20pathlib%20import%20Path%0A%20%20%20%20import%20json%0A%20%20%20%20import%20pickle%0A%0A%20%20%20%20%23%20Import%20our%20main%20package%0A%20%20%20%20import%20scribe%0A%0A%20%20%20%20%23%20Import%20data%20tooling%0A%20%20%20%20import%20pertpy%0A%20%20%20%20import%20scanpy%20as%20sc%0A%0A%20%20%20%20%23%20Numerics%20and%20plotting%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20pandas%20as%20pd%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20seaborn%20as%20sns%0A%0A%20%20%20%20%23%20Set%20our%20plotting%20style%20(totally%20optional)%0A%20%20%20%20scribe.viz.matplotlib_style()%0A%20%20%20%20return%20Path%2C%20json%2C%20np%2C%20pd%2C%20pertpy%2C%20pickle%2C%20plt%2C%20sc%2C%20scribe%2C%20sns%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Loading%20and%20filtering%20the%20data%0A%0A%20%20%20%20We%20load%20the%20dataset%20through%20%60pertpy%60%2C%20drop%20genes%20that%20are%20never%20observed%2C%20and%20remove%20low-quality%20cells%20with%20fewer%20than%201%2C000%20UMIs.%20These%20are%20ordinary%20quality-control%20steps%3B%20nothing%20about%20normalization%20happens%20here%2C%20because%20in%20%60scribe%60%20normalization%20is%20part%20of%20the%20*model*%20(more%20on%20that%20below).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(pertpy%2C%20sc)%3A%0A%20%20%20%20%23%20Load%20the%20Zhao%202021%20perturbation%20dataset.%0A%20%20%20%20adata%20%3D%20pertpy.data.zhao_2021()%0A%0A%20%20%20%20%23%20Drop%20genes%20with%20zero%20total%20counts%20(does%20not%20change%20per-cell%20UMI%20totals).%0A%20%20%20%20sc.pp.filter_genes(adata%2C%20min_counts%3D1)%0A%0A%20%20%20%20%23%20Remove%20low-UMI%20cells.%0A%20%20%20%20_umi_thresh%20%3D%201_000%0A%20%20%20%20sc.pp.filter_cells(adata%2C%20min_counts%3D_umi_thresh)%0A%0A%20%20%20%20print(f%22%7Badata.n_obs%3A%2C%7D%20cells%20x%20%7Badata.n_vars%3A%2C%7D%20genes%20after%20QC%22)%0A%20%20%20%20adata%0A%20%20%20%20return%20(adata%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20We%20focus%20on%20the%20donors%20that%20received%20**panobinostat**%20and%20pull%20both%20their%20treated%20cells%20and%20their%20matched%20%60control%60%20cells%20into%20a%20single%20object%2C%20%60adata_joint%60.%20This%20is%20the%20crossed%20*donor%20%C3%97%20condition*%20table%20the%20model%20will%20see.%20The%20two%20grouping%20columns%20are%3A%0A%0A%20%20%20%20-%20%60perturbation%60%20%E2%80%94%20the%20**contrast%20of%20interest**%20(%60control%60%20vs%20%60panobinostat%60)%2C%0A%20%20%20%20-%20%60sample%60%20%E2%80%94%20the%20**donor**%2C%20which%20we%20want%20to%20account%20for%20but%20not%20test.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata%2C%20pd)%3A%0A%20%20%20%20%23%20Treat%20the%20control%20arm%20as%20dose%200%20(keep%20the%20column's%20categorical%20dtype%20happy).%0A%20%20%20%20adata.obs%5B%22dose_value%22%5D%20%3D%20adata.obs%5B%22dose_value%22%5D.cat.add_categories(%5B%220%22%5D)%0A%20%20%20%20adata.obs.loc%5Badata.obs%5B%22perturbation%22%5D%20%3D%3D%20%22control%22%2C%20%22dose_value%22%5D%20%3D%20%220%22%0A%0A%20%20%20%20%23%20Donors%20that%20received%20panobinostat...%0A%20%20%20%20_perturbation%20%3D%20%22panobinostat%22%0A%20%20%20%20_samples%20%3D%20adata.obs.loc%5B%0A%20%20%20%20%20%20%20%20adata.obs%5B%22perturbation%22%5D%20%3D%3D%20_perturbation%2C%20%22sample%22%0A%20%20%20%20%5D.unique()%0A%0A%20%20%20%20%23%20...and%20the%20matched%20control%20%2B%20treated%20cells%20from%20exactly%20those%20donors.%0A%20%20%20%20_mask%20%3D%20adata.obs%5B%22sample%22%5D.isin(_samples)%20%26%20adata.obs%5B%22perturbation%22%5D.isin(%0A%20%20%20%20%20%20%20%20%5B%22control%22%2C%20_perturbation%5D%0A%20%20%20%20)%0A%20%20%20%20adata_joint%20%3D%20adata%5B_mask%5D.copy()%0A%20%20%20%20perturbation%20%3D%20_perturbation%0A%0A%20%20%20%20%23%20The%20crossed%20donor%20x%20condition%20table%3A%20this%20is%20our%207%20x%202%20%3D%2014%20leaves.%0A%20%20%20%20crosstab%20%3D%20pd.crosstab(%0A%20%20%20%20%20%20%20%20adata_joint.obs%5B%22sample%22%5D%2C%20adata_joint.obs%5B%22perturbation%22%5D%0A%20%20%20%20)%0A%20%20%20%20crosstab%0A%20%20%20%20return%20adata_joint%2C%20perturbation%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Every%20donor%20appears%20in%20**both**%20columns%20%E2%80%94%20a%20fully%20crossed%2C%20balanced-ish%20design.%20Those%2014%20non-empty%20cells%20are%20the%20model's%20*leaves*.%20Reading%20**down%20a%20column**%20later%20will%20compare%20a%20donor%20against%20itself%20across%20the%20two%20conditions%20(the%20paired%20contrast)%3B%20reading%20**across%20a%20row**%20compares%20donors%20within%20a%20condition%20(the%20heterogeneity%20we%20average%20over).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Before%20modeling%2C%20the%20one%20plot%20we%20always%20look%20at%3A%20the%20distribution%20of%20total%20UMIs%20per%20cell.%20If%20library%20size%20varies%20by%20close%20to%20an%20order%20of%20magnitude%2C%20each%20cell%20is%20effectively%20sequenced%20at%20its%20own%20depth%2C%20and%20the%20model%20must%20be%20told%20so%20explicitly.%20In%20%60scribe%60%20that%20is%20the%20**per-cell%20capture%20probability**%20(%60variable_capture%3DTrue%60)%2C%20the%20principled%20replacement%20for%20%22divide%20by%20total%20counts%22.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata_joint%2C%20np%2C%20plt%2C%20sns)%3A%0A%20%20%20%20%23%20Plot%20total%20UMI%20count%20distribution%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(5%2C%204))%0A%20%20%20%20sns.ecdfplot(np.asarray(adata_joint.X.sum(axis%3D1)).ravel()%2C%20ax%3D_ax)%0A%20%20%20%20_ax.set_xscale(%22log%22)%0A%20%20%20%20_ax.set_xlabel(%22total%20UMI%20count%20per%20cell%22)%0A%20%20%20%20_ax.set_ylabel(%22ECDF%22)%0A%20%20%20%20_ax.set_ylim(-0.01%2C%201.01)%0A%20%20%20%20if%20_ax.legend_%20is%20not%20None%3A%0A%20%20%20%20%20%20%20%20_ax.legend_.remove()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20The%20model%2C%20written%20out%0A%0A%20%20%20%20%23%23%23%20Counts%3A%20a%20Negative%20Binomial%20with%20per-cell%20capture%0A%0A%20%20%20%20At%20its%20core%20%60scribe%60%20models%20each%20UMI%20count%20%24u_%7Bcg%7D%24%20(cell%20%24c%24%2C%20gene%20%24g%24)%20with%20a%20**Negative%20Binomial**%2C%20the%20distribution%20that%20already%20underlies%20DESeq2%2FedgeR%20and%20that%20drops%20out%20of%20a%20two-state%20promoter%20model%20of%20transcription.%20With%20%60parameterization%3D%22mean_odds%22%60%20we%20work%20in%20terms%20of%20a%20gene%20**mean%20expression**%20%24%5Cmu_g%24%20and%20a%20shared%20**odds**%20%24%5Cphi%20%3D%20p%2F(1-p)%24%2C%20related%20to%20the%20raw%20NB%20parameters%20by%20%24r_g%20%3D%20%5Cmu_g%5C%2C%5Cphi%24%20and%20%24p%20%3D%201%2F(1%2B%5Cphi)%24.%20The%20mean%E2%80%93odds%20view%20breaks%20the%20banana-shaped%20%24(r%2Cp)%24%20degeneracy%20and%20makes%20%24%5Cmu_g%24%20%E2%80%94%20the%20thing%20the%20data%20pins%20down%20directly%20%E2%80%94%20the%20primary%20object.%0A%0A%20%20%20%20Because%20library%20sizes%20vary%2C%20each%20cell%20gets%20its%20own%20**capture%20probability**%20%24%5Cnu_c%20%5Cin%20(0%2C1)%24%2C%20producing%20a%20cell-specific%20effective%20success%20probability%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Chat%7Bp%7D_c%20%5C%3B%3D%5C%3B%20%5Cfrac%7Bp%7D%7B%5Cnu_c%20%2B%20p%5C%2C(1-%5Cnu_c)%7D%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20So%20far%20this%20is%20the%20standard%20%60scribe%60%20NB-with-variable-capture%20model.%20The%20new%20part%20is%20what%20sits%20**on%20top%20of%20%24%5Cmu_g%24**.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Mean%20expression%3A%20an%20additive%20*crossed*%20hierarchy%0A%0A%20%20%20%20Each%20leaf%20%24%5Cell%24%20is%20a%20(condition%2C%20donor)%20pair.%20Instead%20of%20giving%20every%20leaf%20a%20free%20mean%2C%20we%20**decompose**%20the%20log-mean%20additively%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Clog%20%5Cmu_g%5E%7B(%5Cell)%7D%0A%20%20%20%20%5C%3B%3D%5C%3B%0A%20%20%20%20%5Cunderbrace%7B%5Clog%20%5Cmu_g%5E%7B%5Cmathrm%7Bpop%7D%7D%7D_%7B%5Ctext%7Bpopulation%20baseline%7D%7D%0A%20%20%20%20%5C%3B%2B%5C%3B%0A%20%20%20%20%5Cunderbrace%7B%5Calpha_g%5Cbig%5B%5Cmathrm%7Bcond%7D(%5Cell)%5Cbig%5D%7D_%7B%5Ctext%7Btreatment%20effect%20(fixed)%7D%7D%0A%20%20%20%20%5C%3B%2B%5C%3B%0A%20%20%20%20%5Cunderbrace%7B%5Cbeta_g%5Cbig%5B%5Cmathrm%7Bdonor%7D(%5Cell)%5Cbig%5D%7D_%7B%5Ctext%7Bdonor%20effect%20(random)%7D%7D%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Two%20factors%2C%20two%20very%20different%20roles%3A%0A%0A%20%20%20%20-%20**%24%5Calpha_g%24%20%E2%80%94%20the%20treatment%20effect%20%E2%80%94%20is%20a%20*fixed*%20effect.**%20It%20is%20the%20contrast%20we%20are%20actually%20asking%20about%2C%20so%20we%20give%20it%20a%20**weakly-informative%20Gaussian**%20prior%20with%20a%20*fixed*%20scale%20and%20**no%20learned%20shrinkage**%3A%20%24%5Calpha_%7Bg%2Cc%7D%20%3D%20%5Csigma_%5Calpha%5C%2C%20z_%7Bg%2Cc%7D%5E%7B%5Calpha%7D%24%2C%20%24z%5Csim%5Cmathcal%20N(0%2C1)%24.%20%20Why%20no%20adaptive%20shrinkage%3F%20A%20two-level%20factor%20has%20almost%20no%20information%20to%20estimate%20its%20own%20variance%2C%20and%20a%20learned%20scale%20would%20happily%20pull%20the%20contrast%20toward%20zero%20%E2%80%94%20exactly%20the%20effect%20we%20are%20trying%20to%20measure.%20The%20identified%20quantity%20is%20the%20**difference**%20%24%5Calpha_g%5B%5Ctext%7Bpanobinostat%7D%5D%20-%20%5Calpha_g%5B%5Ctext%7Bcontrol%7D%5D%24.%0A%0A%20%20%20%20-%20**%24%5Cbeta_g%24%20%E2%80%94%20the%20donor%20effect%20%E2%80%94%20is%20a%20*random*%20effect**%20with%20a%20**regularized%20horseshoe**%20prior.%20With%20seven%20donors%20there%20*is*%20information%20to%20learn%20how%20much%20donors%20vary%2C%20and%20the%20horseshoe%20adaptively%20shrinks%3A%20donors%20that%20look%20alike%20on%20a%20gene%20are%20pulled%20together%2C%20while%20a%20genuinely%20deviant%20donor%20is%20left%20alone.%20Each%20%24%5Cbeta_g%5Bd%5D%24%20is%20**zero-mean**%2C%20so%20it%20captures%20deviations%20*from*%20the%20population%20baseline%20rather%20than%20competing%20with%20it.%0A%0A%20%20%20%20The%20population%20baseline%20%24%5Clog%5Cmu_g%5E%7B%5Cmathrm%7Bpop%7D%7D%24%20is%20the%20only%20free%20intercept%2C%20which%20is%20what%20makes%20the%20decomposition%20identifiable.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Why%20this%20is%20the%20right%20shape%20for%20the%20question%0A%0A%20%20%20%20Splitting%20%24%5Clog%5Cmu_g%5E%7B(%5Cell)%7D%24%20into%20a%20fixed%20treatment%20term%20and%20a%20random%20donor%20term%20is%20the%20formal%20statement%20of%20%22estimate%20the%20average%20drug%20effect%20while%20controlling%20for%20the%20donor.%22%20The%20treatment%20effect%20is%20**shared**%20across%20donors%20(so%20all%2014%20leaves%20inform%20it)%2C%20the%20donor%20deviations%20soak%20up%20the%20baseline-and-response%20variability%20we%20do%20not%20want%20to%20mistake%20for%20a%20drug%20effect%2C%20and%20because%20control%20and%20treated%20cells%20of%20one%20donor%20pass%20through%20the%20*same*%20%24%5Cbeta_g%5Bd%5D%24%2C%20the%20design%20stays%20**paired**.%0A%0A%20%20%20%20Only%20the%20**expression**%20target%20(%24%5Cmu%24)%20gets%20this%20additive%20decomposition%3B%20the%20technical%20odds%20parameter%20is%20left%20free%20per%20leaf.%20That%20keeps%20the%20per-cell%20likelihood%20%E2%80%94%20the%20computational%20hot%20path%20%E2%80%94%20completely%20unchanged%3A%20the%20leaves%20are%20still%20an%20ordinary%20%22dataset%22%20axis%2C%20with%20structure%20layered%20*above*%20them.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Fitting%20the%20model%0A%0A%20%20%20%20The%20call%20below%20is%20the%20whole%20model.%20%60hierarchy%3D%5B...%5D%60%20declares%20the%20two%20grouping%20factors%20and%20their%20effect%20types%3B%20%60expression_dataset_prior%60%20picks%20the%20prior%20family%20per%20factor%20(a%20fixed-scale%20Gaussian%20for%20the%20treatment%20contrast%2C%20a%20horseshoe%20for%20the%20donors).%20Crossing%20is%20implicit%20%E2%80%94%20listing%20two%20factors%20with%20no%20nesting%20means%20%22donor%20crossed%20with%20condition.%22%0A%0A%20%20%20%20Fitting%20at%20this%20scale%20takes%20a%20while%2C%20so%20we%20**load%20a%20cached%20fit**%20when%20one%20is%20present%20and%20only%20run%20the%20optimizer%20otherwise.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Path%2C%20adata_joint%2C%20json%2C%20perturbation%2C%20pickle%2C%20scribe)%3A%0A%20%20%20%20%23%20Resolve%20the%20dataset%20directory%20from%20the%20local%20tutorial%20path%20config.%0A%20%20%20%20with%20(%0A%20%20%20%20%20%20%20%20Path(__file__)%0A%20%20%20%20%20%20%20%20.with_name(%22tutorial_paths.local.json%22)%0A%20%20%20%20%20%20%20%20.open(%22r%22%2C%20encoding%3D%22utf-8%22)%20as%20_f%0A%20%20%20%20)%3A%0A%20%20%20%20%20%20%20%20_data_root%20%3D%20json.load(_f)%5B%22SCRIBE_TUTORIAL_DATA_ROOT%22%5D%0A%20%20%20%20data_dir%20%3D%20Path(_data_root).expanduser()%20%2F%20%22zhao_2021%22%0A%20%20%20%20(data_dir%20%2F%20%22scribe_results%22).mkdir(parents%3DTrue%2C%20exist_ok%3DTrue)%0A%0A%20%20%20%20%23%20Define%20output%20path%0A%20%20%20%20_out_path%20%3D%20(%0A%20%20%20%20%20%20%20%20data_dir%0A%20%20%20%20%20%20%20%20%2F%20%22scribe_results%22%0A%20%20%20%20%20%20%20%20%2F%20f%22scribe_hierarchical_%7Bperturbation%7D_joint.pkl%22%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20Fit%2Fload%20pre-fit%20model%0A%20%20%20%20if%20_out_path.exists()%3A%0A%20%20%20%20%20%20%20%20with%20open(_out_path%2C%20%22rb%22)%20as%20_f%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20results_joint%20%3D%20pickle.load(_f)%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20results_joint%20%3D%20scribe.fit(%0A%20%20%20%20%20%20%20%20%20%20%20%20adata_joint%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20parameterization%3D%22mean_odds%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20variable_capture%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20unconstrained%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20positive_transform%3D%22exp%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Donor%20(sample)%20CROSSED%20with%20condition%20(perturbation)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20hierarchy%3D%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%202-level%20contrast%20of%20interest%20-%3E%20fixed%2C%20weakly-informative.%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20scribe.GroupLevel(%22perturbation%22%2C%20effect_type%3D%22fixed%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%207%20donors%20-%3E%20random%20effect%20with%20adaptive%20shrinkage.%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20scribe.GroupLevel(%22sample%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20expression_dataset_prior%3D%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22perturbation%22%3A%20%22gaussian%22%2C%20%20%23%20fixed-scale%20contrast%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22sample%22%3A%20%22horseshoe%22%2C%20%20%23%20shrink%20across%20donors%0A%20%20%20%20%20%20%20%20%20%20%20%20%7D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20prob_dataset_prior%3D%22gaussian%22%2C%20%20%23%20technical%20odds%3A%20leaf-exchangeable%0A%20%20%20%20%20%20%20%20%20%20%20%20early_stopping%3D%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22enabled%22%3A%20True%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22patience%22%3A%201000%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22restore_best%22%3A%20True%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%7D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20n_steps%3D100_000%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20batch_size%3D4096%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20gene_coverage%3D0.99%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20with%20open(_out_path%2C%20%22wb%22)%20as%20_f%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20pickle.dump(results_joint%2C%20_f)%0A%0A%20%20%20%20results_joint%0A%20%20%20%20return%20(results_joint%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Did%20the%20fit%20work%3F%20Three%20diagnostics%0A%0A%20%20%20%20%23%23%23%201.%20The%20ELBO%20loss%0A%0A%20%20%20%20Variational%20inference%20maximizes%20the%20**evidence%20lower%20bound%20(ELBO)**%3B%20the%20loss%20we%20plot%20is%20its%20negative.%20We%20do%20not%20read%20it%20quantitatively%20%E2%80%94%20we%20just%20want%20the%20textbook%20shape%3A%20a%20fast%20initial%20drop%20that%20flattens%20into%20a%20plateau.%20Spikes%20or%20upward%20drift%20would%20be%20warnings.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(results_joint%2C%20scribe)%3A%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_loss(results_joint%2C%20figsize%3D(7%2C%203))%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%202.%20Mean%20calibration%20%E2%80%94%20*per%20leaf*%0A%0A%20%20%20%20A%20converged%20loss%20does%20not%20prove%20the%20model%20describes%20the%20**counts**.%20The%20mean-calibration%20plot%20compares%2C%20for%20each%20gene%2C%20the%20**observed**%20mean%20count%20to%20the%20model's%20**predicted**%20mean%3B%20points%20on%20the%20diagonal%20mean%20the%20fit%20reproduces%20the%20first%20moment%20of%20the%20data.%0A%0A%20%20%20%20Because%20we%20have%20a%20crossed%20hierarchy%2C%20%60scribe%60%20lays%20the%20panels%20out%20on%20a%20**condition%20%C3%97%20donor%20grid**%3A%20rows%20are%20the%20two%20conditions%2C%20columns%20are%20the%20seven%20donors.%20This%20is%20the%20multi-factor%20layout%20doing%20its%20job%20%E2%80%94%20reading%20**down%20a%20column**%20shows%20control%20vs%20panobinostat%20*for%20one%20donor*%2C%20exactly%20the%20paired%20comparison%20the%20model%20encodes.%20A%20leaf%20that%20fell%20off%20the%20diagonal%20would%20point%20to%20a%20donor%20or%20condition%20the%20additive%20structure%20is%20failing%20to%20capture.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata_joint%2C%20results_joint%2C%20scribe)%3A%0A%20%20%20%20%23%20Per-leaf%20observed-vs-predicted%20means.%20The%20panels%20auto-arrange%20as%20a%0A%20%20%20%20%23%20condition%20(rows)%20x%20donor%20(cols)%20grid%20from%20the%20model's%20grouping%20spec.%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_mean_calibration(results_joint%2C%20counts%3Dadata_joint)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20For%20the%20most%20part%2C%20we%20have%20a%20very%20good%20fit%20as%20most%20points%20are%20close%20to%20the%20diagonal.%20The%20exceptions%20are%20samples%20%60PW051%60%20and%20%60PW053%60.%20However%2C%20because%20of%20this%20diagnostic%2C%20if%20wanted%2C%20we%20could%20look%20deeper%20into%20what%20could%20have%20caused%20this%20mismatch.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%203.%20Posterior%20predictive%20check%0A%0A%20%20%20%20Finally%2C%20a%20**posterior%20predictive%20check%20(PPC)**%3A%20simulate%20fresh%20UMI%20counts%20from%20the%20fitted%20model%20and%20overlay%20their%20distribution%20on%20the%20real%20histograms%2C%20for%20a%20spread%20of%20genes%20across%20the%20expression%20range.%20If%20the%20generative%20story%20is%20adequate%2C%20the%20simulated%20bands%20should%20track%20the%20observed%20counts.%0A%0A%20%20%20%20Because%20this%20is%20a%20multi-dataset%20fit%2C%20each%20simulated%20cell%20is%20drawn%20from%20**its%20own%20leaf's**%20parameters%2C%20so%20the%20plot%20below%20is%20an%20honest%20aggregate%20over%20all%20fourteen%20donor%20%C3%97%20condition%20leaves.%20That%20makes%20it%20a%20good%20global%20check%20%E2%80%94%20but%20it%20cannot%20tell%20us%20whether%20any%20*individual*%20fit%20is%20good.%20We%20start%20global%2C%20then%20zoom%20in.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata_joint%2C%20results_joint%2C%20scribe)%3A%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_ppc(%0A%20%20%20%20%20%20%20%20results_joint%2C%0A%20%20%20%20%20%20%20%20adata_joint%2C%0A%20%20%20%20%20%20%20%20n_genes%3D16%2C%0A%20%20%20%20%20%20%20%20n_rows%3D4%2C%0A%20%20%20%20%20%20%20%20n_samples%3D256%2C%0A%20%20%20%20%20%20%20%20figsize%3D(7%2C%207)%2C%0A%20%20%20%20)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%23%20Per-dataset%20PPC%3A%20the%20quality%20of%20an%20*individual*%20fit%0A%0A%20%20%20%20The%20mean-calibration%20grid%20flagged%20%60PW051%60%20and%20%60PW053%60.%20The%20PPC%20can%20zoom%20into%20exactly%20those%20leaves%3A%20%60plot_ppc(...%2C%20dataset%3D...)%60%20swaps%20in%20a%20single%20leaf's%20parameter%20view%20(%60results.get_dataset(leaf)%60%20under%20the%20hood)%20and%20the%20**observed%20cells%20for%20that%20leaf%20only**%2C%20so%20we%20judge%20one%20donor%20%C3%97%20condition%20fit%20at%20a%20time.%20Address%20a%20leaf%20either%20by%20integer%20index%20or%2C%20more%20readably%2C%20by%20a%20%60%7Bfactor%3A%20level%7D%60%20dict.%0A%0A%20%20%20%20The%20two%20donors%20below%20sit%20at%20opposite%20ends%20of%20the%20data-richness%20scale%2C%20which%20is%20the%20structural%20reason%20their%20fits%20differ%3A%20**%60PW030%60**%20contributes%20%E2%89%8819%2C000%20control%20cells%2C%20while%20**%60PW051%60**%20%E2%80%94%20one%20of%20the%20donors%20the%20calibration%20plot%20flagged%20%E2%80%94%20contributes%20only%20%E2%89%881%2C200.%20A%20data-rich%20leaf's%20predictive%20bands%20should%20hug%20the%20observed%20histograms%3B%20a%20sparse%20leaf's%20bands%20are%20wider%20and%20the%20fit%20is%20more%20stretched%2C%20which%20is%20what%20the%20calibration%20plot%20was%20picking%20up.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata_joint%2C%20results_joint%2C%20scribe)%3A%0A%20%20%20%20%23%20A%20data-rich%20donor%3A%20its%20individual%20fit%20should%20track%20the%20data%20tightly.%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_ppc(%0A%20%20%20%20%20%20%20%20results_joint%2C%0A%20%20%20%20%20%20%20%20adata_joint%2C%0A%20%20%20%20%20%20%20%20dataset%3D%7B%22sample%22%3A%20%22PW030%22%2C%20%22perturbation%22%3A%20%22control%22%7D%2C%0A%20%20%20%20%20%20%20%20n_genes%3D16%2C%0A%20%20%20%20%20%20%20%20n_rows%3D4%2C%0A%20%20%20%20%20%20%20%20n_samples%3D256%2C%0A%20%20%20%20%20%20%20%20figsize%3D(7%2C%207)%2C%0A%20%20%20%20)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(adata_joint%2C%20results_joint%2C%20scribe)%3A%0A%20%20%20%20%23%20A%20sparse%20donor%20the%20mean-calibration%20flagged%3A%20the%20same%20diagnostic%2C%20where%0A%20%20%20%20%23%20the%20individual%20fit%20is%20most%20stretched%20(%E2%89%881%2C200%20cells%20vs%20%E2%89%8819%2C000%20above).%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_ppc(%0A%20%20%20%20%20%20%20%20results_joint%2C%0A%20%20%20%20%20%20%20%20adata_joint%2C%0A%20%20%20%20%20%20%20%20dataset%3D%7B%22sample%22%3A%20%22PW051%22%2C%20%22perturbation%22%3A%20%22control%22%7D%2C%0A%20%20%20%20%20%20%20%20n_genes%3D16%2C%0A%20%20%20%20%20%20%20%20n_rows%3D4%2C%0A%20%20%20%20%20%20%20%20n_samples%3D256%2C%0A%20%20%20%20%20%20%20%20figsize%3D(7%2C%207)%2C%0A%20%20%20%20)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Reading%20the%20model's%20structure%0A%0A%20%20%20%20Before%20any%20differential-expression%20machinery%2C%20the%20fitted%20hierarchy%20is%20already%20interpretable.%20The%20additive%20effects%20are%20stored%20as%20parameters%20with%20posteriors%2C%20and%20%60scribe%60%20exposes%20them%20directly.%0A%0A%20%20%20%20%23%23%23%20The%20treatment%20effect%0A%0A%20%20%20%20%60results.get_factor_effect(%22perturbation%22)%60%20returns%20the%20fitted%20%24%5Calpha_g%24%20effect.%20Since%20%60perturbation%60%20is%20the%20fixed%20contrast%20factor%2C%20the%20identified%20quantity%20is%20the%20**contrast**%20%E2%80%94%20the%20per-gene%20log-mean%20shift%20from%20control%20to%20panobinostat%20%E2%80%94%20which%20we%20summarize%20by%20its%20posterior%20mean.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt%2C%20results_joint)%3A%0A%20%20%20%20%23%20The%20fixed%20treatment%20effect%3A%20alpha%5Bpanobinostat%5D%20-%20alpha%5Bcontrol%5D%2C%20in%20log-mu.%0A%20%20%20%20%23%20(Gene%20*names*%20live%20in%20the%20DE%20table%20below%3B%20here%20we%20work%20with%20the%20array.)%0A%20%20%20%20_tx%20%3D%20results_joint.get_factor_effect(%22perturbation%22)%0A%20%20%20%20treatment_effect%20%3D%20np.asarray(_tx.map_contrast(%22panobinostat%22%2C%20%22control%22))%0A%20%20%20%20print(f%22treatment-effect%20factor%3A%20%7B_tx!r%7D%22)%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20f%22%7Bint((np.abs(treatment_effect)%20%3E%20np.log(2)).sum())%3A%2C%7D%20genes%20shift%20by%20%22%0A%20%20%20%20%20%20%20%20f%22more%20than%202-fold%20(%7Clog-mu%20effect%7C%20%3E%20log%202)%22%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20Plot%20distribution%20of%20treatment%20effects%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(5%2C%203.2))%0A%20%20%20%20_ax.hist(treatment_effect%2C%20bins%3D80%2C%20color%3D%22steelblue%22)%0A%20%20%20%20_ax.axvline(0%2C%20ls%3D%22--%22%2C%20c%3D%22k%22%2C%20lw%3D1)%0A%20%20%20%20_ax.set_xlabel(r%22per-gene%20treatment%20effect%20%24%5Cbar%5Calpha_g%24%20(log-mean)%22)%0A%20%20%20%20_ax.set_ylabel(%22genes%22)%0A%20%20%20%20_fig%0A%20%20%20%20return%20(treatment_effect%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Donor%20heterogeneity%20%E2%80%94%20the%20reason%20we%20went%20joint%0A%0A%20%20%20%20%60results.get_factor_effect(%22sample%22)%60%20returns%20the%20**random**%20donor%20effects%20%24%5Cbeta_g%5Bd%5D%24%3A%20a%20%247%20%5Ctimes%20G%24%20matrix%20of%20zero-mean%2C%20log-mean%20deviations.%20Their%20spread%20across%20donors%20is%20the%20heterogeneity%20the%20model%20absorbed%20so%20that%20it%20did%20*not*%20leak%20into%20the%20treatment%20effect.%20If%20this%20spread%20were%20negligible%2C%20simple%20pooling%20would%20have%20been%20fine%3B%20if%20it%20is%20large%2C%20the%20joint%20model%20earned%20its%20keep.%0A%0A%20%20%20%20Below%20we%20compare%2C%20gene%20by%20gene%2C%20the%20**magnitude%20of%20the%20treatment%20effect**%20against%20the%20**donor-to-donor%20spread**.%20Genes%20where%20the%20donor%20spread%20rivals%20or%20exceeds%20the%20treatment%20effect%20are%20exactly%20the%20ones%20a%20naive%20pooled%20analysis%20would%20get%20wrong.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt%2C%20results_joint%2C%20treatment_effect)%3A%0A%20%20%20%20%23%20Extract%20random%20donor%20effect%0A%20%20%20%20_donor%20%3D%20results_joint.get_factor_effect(%22sample%22)%0A%20%20%20%20_eff%20%3D%20np.asarray(%0A%20%20%20%20%20%20%20%20_donor.effects()%0A%20%20%20%20)%20%20%23%20(n_donors%2C%20G)%20posterior-mean%20log-mu%20devs%0A%20%20%20%20_donor_spread%20%3D%20_eff.std(axis%3D0)%20%20%23%20per-gene%20SD%20across%20donors%0A%20%20%20%20_tx_abs%20%3D%20np.abs(treatment_effect)%0A%0A%20%20%20%20%23%20Plot%20donor-to-donor%20spread%20vs%20treatment%20effect%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(5%2C%205))%0A%20%20%20%20_ax.scatter(_donor_spread%2C%20_tx_abs%2C%20s%3D4%2C%20alpha%3D0.3)%0A%20%20%20%20_lim%20%3D%20float(np.nanmax(%5B_donor_spread.max()%2C%20_tx_abs.max()%5D))%20*%201.05%0A%20%20%20%20_ax.plot(%5B0%2C%20_lim%5D%2C%20%5B0%2C%20_lim%5D%2C%20ls%3D%22--%22%2C%20c%3D%22k%22%2C%20lw%3D1)%0A%20%20%20%20_ax.set_xlabel(%22donor-to-donor%20spread%20%20(SD%20of%20%24%5C%5Cbeta_g%5Bd%5D%24)%22)%0A%20%20%20%20_ax.set_ylabel(%22%7Ctreatment%20effect%7C%20%20%24%7C%5C%5Cbar%5C%5Calpha_g%7C%24%22)%0A%20%20%20%20_ax.set_title(f%22%7B_donor.n_levels%7D%20donors%3A%20heterogeneity%20vs%20treatment%22)%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Recovering%20the%20paired%20groups%0A%0A%20%20%20%20%60scribe%60%20also%20lets%20us%20pull%20back%20the%20leaf%20structure%20for%20any%20downstream%20slicing.%20%20%60results.get_group(sample%3D%22...%22)%60%20returns%20the%20leaves%20for%20one%20donor%2C%20indexed%20by%20condition%3B%20%60iter_groups(%22sample%22)%60%20walks%20every%20donor.%20Each%20value%20is%20an%20ordinary%20single-dataset%20results%20view%2C%20so%20you%20can%20drop%20it%20into%20any%20per-leaf%20analysis.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(results_joint)%3A%0A%20%20%20%20%23%20One%20donor's%20paired%20(control%2C%20panobinostat)%20views%2C%20keyed%20by%20condition.%0A%20%20%20%20_g%20%3D%20results_joint.get_group(sample%3Dresults_joint.group_levels(%22sample%22)%5B0%5D)%0A%20%20%20%20print(%22group%3A%22%2C%20repr(_g))%0A%20%20%20%20print(%22leaf%20labels%3A%22%2C%20_g.labels)%0A%20%20%20%20%23%20Enumerate%20all%20seven%20paired%20donor%20groups.%0A%20%20%20%20print(%22%5Cnall%20donor%20groups%3A%22)%0A%20%20%20%20for%20_level%2C%20_grp%20in%20results_joint.iter_groups(%22sample%22)%3A%0A%20%20%20%20%20%20%20%20print(f%22%20%20%7B_level%7D%3A%20conditions%20%3D%20%7B_grp.keys()%7D%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Differential%20expression%3A%20the%20donor-averaged%20treatment%20effect%0A%0A%20%20%20%20%23%23%23%20The%20estimand%0A%0A%20%20%20%20We%20want%20the%20population%20treatment%20effect%2C%20*paired%20within%20donor*.%20%60scribe%60's%20compositional%20DE%20works%20in%20**centered-log-ratio%20(CLR)**%20space%2C%20where%20each%20gene%20is%20compared%20to%20the%20geometric%20mean%20across%20genes%20%E2%80%94%20a%20**reference-free**%20coordinate%20in%20which%20differences%20are%20additive.%20For%20each%20donor%20%24d%24%20present%20in%20both%20arms%20and%20each%20posterior%20draw%20%24s%24%2C%20we%20form%20the%20within-donor%20CLR%20difference%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5CDelta_g%5E%7B(d%2Cs)%7D%0A%20%20%20%20%3D%20%5Cmathrm%7Bclr%7D%5C!%5Cbig(%5Crho_g%5E%7B(%5Ctext%7Bcontrol%7D%2C%5C%2Cd%2C%5C%2Cs)%7D%5Cbig)%0A%20%20%20%20-%20%5Cmathrm%7Bclr%7D%5C!%5Cbig(%5Crho_g%5E%7B(%5Ctext%7Bpanobinostat%7D%2C%5C%2Cd%2C%5C%2Cs)%7D%5Cbig)%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20and%20then%20**average%20over%20donors**%20to%20get%20the%20paired%20main%20effect%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbar%5CDelta_g%5E%7B(s)%7D%20%3D%20%5Csum_d%20w_d%5C%2C%20%5CDelta_g%5E%7B(d%2Cs)%7D%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20donor%20effect%20%24%5Cbeta_g%5Bd%5D%24%20cancels%20inside%20each%20within-donor%20difference%2C%20so%20%24%5Cbar%5CDelta_g%24%20is%20a%20clean%20draw%20from%20the%20posterior%20of%20the%20*average%20treatment%20effect*.%20The%20whole%20population%20contrast%20is%20one%20call%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(results_joint%2C%20scribe)%3A%0A%20%20%20%20%23%20Paired%20main%20effect%20across%20donors.%20delta%20%3D%20CLR(control)%20-%20CLR(panobinostat)%2C%0A%20%20%20%20%23%20so%20a%20gene%20UP%20in%20panobinostat%20has%20a%20NEGATIVE%20delta.%0A%20%20%20%20%23%20n_samples%20sets%20N%20in%20delta_samples%20(more%20draws%20-%3E%20smoother%20lfsr%2FPEFP).%0A%20%20%20%20%23%20batch_size%20chunks%20the%20draw%20and%20offloads%20it%20to%20host%20RAM%20(automatic%20for%0A%20%20%20%20%23%20n_samples%20%3E%20500)%2C%20keeping%20GPU%20memory%20free%20for%20the%20composition%20sampling.%0A%20%20%20%20results_de_raw%20%3D%20scribe.compare_groups(%0A%20%20%20%20%20%20%20%20results_joint%2C%0A%20%20%20%20%20%20%20%20%22perturbation%22%2C%0A%20%20%20%20%20%20%20%20%22control%22%2C%0A%20%20%20%20%20%20%20%20%22panobinostat%22%2C%0A%20%20%20%20%20%20%20%20n_samples%3D5_000%2C%0A%20%20%20%20%20%20%20%20batch_size%3D500%2C%0A%20%20%20%20)%0A%20%20%20%20results_de_raw%0A%20%20%20%20return%20(results_de_raw%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20A%20cautionary%20detour%3A%20CLR%20needs%20a%20stable%20reference%0A%0A%20%20%20%20CLR%20is%20principled%2C%20but%20it%20has%20a%20sharp%20edge.%20The%20reference%20it%20divides%20by%20is%20the%20**geometric%20mean%20of%20the%20log-fractions%20across%20*all*%20genes**%20%E2%80%94%20and%20a%20typical%20single-cell%20matrix%20has%20thousands%20of%20genes%20sitting%20at%20essentially%20zero%2C%20whose%20log-fractions%20are%20hugely%20negative%20and%20jitter%20from%20draw%20to%20draw.%20That%20jitter%20contaminates%20the%20reference%20and%20blows%20up%20the%20CLR%20contrast%20of%20the%20genes%20we%20actually%20care%20about.%0A%0A%20%20%20%20Look%20at%20the%20top%20%22hits%22%20from%20the%20unfiltered%20run%3A%20CLR%20shifts%20of%20%24%5Cpm%2040%24%20(an%20%24e%5E%7B40%7D%24-fold%20change%20is%20not%20biology)%20on%20near-zero%20pseudogenes%2C%20mitochondrial%20tRNAs%2C%20and%20stray%20neuronal%20markers.%20This%20is%20the%20reference%20exploding%2C%20not%20a%20drug%20effect.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20results_de_raw)%3A%0A%20%20%20%20_df_raw%20%3D%20results_de_raw.to_dataframe(%0A%20%20%20%20%20%20%20%20tau%3Dnp.log(1.1)%2C%0A%20%20%20%20%20%20%20%20target_pefp%3D0.05%2C%0A%20%20%20%20%20%20%20%20metrics%3D%22clr%22%2C%0A%20%20%20%20%20%20%20%20column_naming%3D%22prefixed%22%2C%0A%20%20%20%20)%0A%20%20%20%20%23%20Most%20extreme%20%7CCLR%20shift%7C%20--%20the%20unstable%2C%20near-zero%20genes.%0A%20%20%20%20_df_raw.reindex(%0A%20%20%20%20%20%20%20%20_df_raw%5B%22clr_delta_mean%22%5D.abs().sort_values(ascending%3DFalse).index%0A%20%20%20%20)%5B%5B%22gene%22%2C%20%22clr_delta_mean%22%2C%20%22clr_lfsr%22%5D%5D.head(12)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20The%20principled%20fix%3A%20keep%20the%20genes%20that%20carry%20the%20mass%2C%20pool%20the%20rest%20into%20%22other%22%0A%0A%20%20%20%20The%20cure%20is%20*not*%20to%20delete%20the%20quiet%20genes%20%E2%80%94%20that%20would%20break%20the%20simplex.%20Instead%20we%20**pool**%20the%20long%20tail%20of%20near-silent%20genes%20into%20a%20single%20%22other%22%20pseudo-gene%20before%20forming%20CLR%20coordinates.%20By%20the%20**Dirichlet%20closure%20property**%20(summing%20a%20subset%20of%20Dirichlet%20components%20gives%20another%20Dirichlet)%20this%20aggregation%20is%20**exact**%20%E2%80%94%20no%20information%20is%20lost%20%E2%80%94%20and%20it%20gives%20the%20geometric-mean%20reference%20something%20stable%20to%20stand%20on.%20CLR%20stays%20exactly%20as%20principled%3B%20we%20have%20only%20stopped%20asking%20it%20to%20resolve%20genes%20that%20carry%20no%20signal.%0A%0A%20%20%20%20How%20aggressively%20should%20we%20pool%3F%20Rather%20than%20an%20arbitrary%20UMI%20floor%2C%20we%20use%20a%20**compositional-coverage**%20rule%20%E2%80%94%20the%20natural%20choice%20*inside*%20a%20compositional%20framework%3A%20keep%20the%20smallest%20set%20of%20genes%20whose%20MAP%20composition%20reaches%20a%20target%20fraction%20of%20the%20total%20expressed%20mass%20(set%20below)%2C%20in%20*either*%20arm%2C%20and%20pool%20the%20rest.%20It%20is%20**scale-free**%20(invariant%20to%20sequencing%20depth)%20and%20keeps%20exactly%20the%20genes%20that%20actually%20make%20up%20the%20transcriptome%2C%20rather%20than%20a%20hand-picked%20count%20cutoff.%0A%0A%20%20%20%20This%20nests%20cleanly%20with%20a%20filter%20we%20already%20applied%20at%20*fit*%20time%3A%20%60gene_coverage%3D0.99%60%20pooled%20the%20un-modeled%20tail%20into%20the%20model's%20own%20%60_other%60%20pseudo-gene%2C%20so%20the%20fitted%20composition%20already%20sums%20to%20one%20over%20the%20full%20transcriptome.%20The%20DE%20composition%20samples%20that%20whole%20simplex%20%E2%80%94%20%60_other%60%20included%20%E2%80%94%20and%20our%20coverage%20mask%20simply%20folds%20%60_other%60%20back%20into%20the%20DE%20%22other%22%20tail%20(it%20is%20the%20compositional%20anchor%2C%20never%20a%20reportable%20gene).%20So%20the%20two%20filtering%20stages%20compose%3A%20closure%20is%20preserved%20end%20to%20end%2C%20and%20the%20CLR%20reference%20is%20always%20the%20full%20transcriptome.%0A%0A%20%20%20%20One%20API%20note%3A%20%60compare_groups%60%20does%20not%20retain%20the%20per-pair%20simplex%20(averaging%20simplices%20would%20not%20equal%20the%20paired%20CLR%20average)%2C%20so%20the%20in-place%20%60set_composition_coverage%60%20does%20not%20apply%20to%20it.%20Instead%20we%20**build**%20the%20mask%20from%20the%20model's%20mean%20expression%20%E2%80%94%20%60results_de_raw.composition_coverage_mask(...)%60%2C%20which%20reads%20only%20%60mu_map%60%20%E2%80%94%20and%20pass%20it%20up%20front%20as%20%60gene_mask%3D%60%2C%20so%20the%20pooling%20happens%20**inside%20each%20donor%20pair%2C%20before%20CLR**.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(results_de_raw%2C%20results_joint%2C%20scribe)%3A%0A%20%20%20%20COMPOSITION_COVERAGE%20%3D%200.95%0A%0A%20%20%20%20%23%20Keep%20the%20genes%20making%20up%20this%20fraction%20of%20the%20expressed%20mass%20in%20either%0A%20%20%20%20%23%20arm%3B%20pool%20the%20rest%20into%20%22other%22.%20The%20builder%20reads%20only%20the%20model's%20mean%0A%20%20%20%20%23%20expression%2C%20so%20it%20works%20even%20though%20compare_groups%20stores%20no%20simplex.%0A%20%20%20%20gene_mask%20%3D%20results_de_raw.composition_coverage_mask(COMPOSITION_COVERAGE)%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20f%22keeping%20%7Bint(gene_mask.sum())%3A%2C%7D%20%2F%20%7Bgene_mask.size%3A%2C%7D%20genes%20%22%0A%20%20%20%20%20%20%20%20f%22(%7Bgene_mask.mean()%3A.0%25%7D)%3B%20the%20rest%20pool%20into%20'other'%22%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20Re-run%20the%20paired%20contrast%20with%20the%20mask%20applied%20within%20each%20pair.%20The%20same%0A%20%20%20%20%23%20n_samples%20sets%20the%20posterior%20draw%20count%20behind%20the%20lfsr%2FPEFP%20estimates%3B%0A%20%20%20%20%23%20batch_size%20keeps%20the%20large%20draw%20within%20memory.%0A%20%20%20%20results_de%20%3D%20scribe.compare_groups(%0A%20%20%20%20%20%20%20%20results_joint%2C%0A%20%20%20%20%20%20%20%20%22perturbation%22%2C%0A%20%20%20%20%20%20%20%20%22control%22%2C%0A%20%20%20%20%20%20%20%20%22panobinostat%22%2C%0A%20%20%20%20%20%20%20%20gene_mask%3Dgene_mask%2C%0A%20%20%20%20%20%20%20%20n_samples%3D5_000%2C%0A%20%20%20%20%20%20%20%20batch_size%3D500%2C%0A%20%20%20%20)%0A%20%20%20%20results_de%0A%20%20%20%20return%20gene_mask%2C%20results_de%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20The%20compositional%20view%0A%0A%20%20%20%20With%20a%20stable%20reference%2C%20the%20standard%20%60scribe%60%20DE%20plots%20are%20well-defined.%20The%20**volcano**%20puts%20the%20CLR%20contrast%20on%20the%20%24x%24-axis%20and%20%24-%5Clog_%7B10%7D(%5Ctext%7Blfsr%7D)%24%20on%20the%20%24y%24-axis%2C%20where%20**lfsr**%20(local%20false%20sign%20rate)%20is%20the%20posterior%20probability%20that%20we%20have%20the%20*direction*%20of%20the%20effect%20wrong%20%E2%80%94%20the%20Bayesian%20replacement%20for%20a%20%24p%24-value.%20Genes%20called%20differentially%20expressed%20(combining%20an%20effect-size%20threshold%20%24%5Ctau%24%20and%20a%20target%20false-sign%20rate)%20are%20highlighted%3B%20the%20mean-expression%20panel%20shows%20each%20gene's%20CLR%20mean%20in%20the%20two%20arms.%0A%0A%20%20%20%20A%20note%20worth%20internalizing%20%E2%80%94%20and%20a%20real%20result%20on%20this%20dataset.%20Under%20a%20**broad**%20perturbation%20like%20an%20HDAC%20inhibitor%20*many*%20genes%20move%2C%20so%20the%20compositional%20reference%20itself%20drifts%2C%20and%20CLR's%20*confident*%20calls%20become%20the%20genes%20whose%20**proportion**%20swings%20hardest%20rather%20than%20the%20coherent%20expression%20program.%20Here%20those%20turn%20out%20to%20be%20**mitochondrial%20genes%20and%20MT-pseudogenes**%20(%60MT-ND6%60%2C%20%60MT-TD%60%2C%20%60MTRNR2L1%60%2C%20%E2%80%A6)%20%E2%80%94%20mitochondrial%20fraction%20is%20a%20classic%20single-cell%20confounder%20%E2%80%94%20while%20the%20canonical%20HDAC%20markers%20carry%20**high%20CLR%20lfsr**%3A%20once%20the%20whole%20composition%20is%20shifting%2C%20the%20sign%20of%20any%20single%20gene's%20*relative*%20change%20is%20genuinely%20uncertain.%20So%20CLR%20is%20behaving%20exactly%20as%20designed%20%E2%80%94%20guarding%20against%20being%20fooled%20by%20composition%20%E2%80%94%20but%20for%20this%20question%20it%20is%20the%20**more%20conservative**%20lens%3A%20at%20a%20strict%20error%20threshold%20its%20only%20confident%20calls%20are%20these%20compositional%20outliers.%20That%20is%20**not**%20the%20same%20as%20CLR%20contradicting%20the%20biology.%20As%20we%20show%20at%20the%20end%20of%20this%20section%2C%20CLR%20and%20the%20biological%20view%20agree%20in%20*direction*%20%E2%80%94%20CLR%20simply%20holds%20the%20coherent%20program%20at%20**lower%20confidence**%3A%20it%20is%20**under-confident%20here%2C%20not%20incoherent**.%20(Routinely%20dropping%20mitochondrial%20genes%20at%20QC%20would%20clean%20it%20up%3B%20we%20keep%20them%20in%20the%20headline%20run%20to%20make%20the%20point%20%E2%80%94%20and%20then%20drop%20them%2C%20by%20*name*%2C%20to%20see%20what's%20left.)%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20results_de%2C%20scribe)%3A%0A%20%20%20%20TAU%20%3D%20np.log(1.1)%0A%20%20%20%20TARGET_PEFP%20%3D%200.05%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_de_volcano(%0A%20%20%20%20%20%20%20%20results_de%2C%20mode%3D%22clr%22%2C%20tau%3DTAU%2C%20target_pefp%3DTARGET_PEFP%2C%20figsize%3D(5%2C%204)%0A%20%20%20%20)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%20TARGET_PEFP%2C%20TAU%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Cleaning%20up%20by%20identity%20%E2%80%94%20and%20what%20it%20reveals%0A%0A%20%20%20%20The%20mitochondrial%20dominance%20is%20a%20*known%20technical%20confounder*%2C%20so%20the%20principled%20fix%20is%20to%20drop%20those%20genes%20by%20**name**%20%E2%80%94%20mitochondrial%2C%20ribosomal%2C%20and%20hemoglobin%20families%20%E2%80%94%20a%20decision%20made%20from%20gene%20*identity*%2C%20a%20priori%20and%20**independent%20of%20the%20contrast**.%20(The%20tempting%20alternative%20%E2%80%94%20thresholding%20on%20the%20large%20CLR%20values%20themselves%20and%20re-running%20%E2%80%94%20is%20*circular*%3A%20it%20selects%20genes%20by%20the%20very%20effect%20we%20are%20testing%2C%20defining%20away%20real%20signal%20along%20with%20artifacts.%20Filtering%20on%20identity%2C%20or%20on%20expression%2Fprecision%2C%20avoids%20that.)%20%60exclude_gene_name_mask%60%20builds%20a%20keep-mask%20from%20name%20prefixes%20and%20regexes%3B%20we%20%60%26%60%20it%20with%20the%20composition-coverage%20mask%20and%20re-run%20the%20same%20CLR%20contrast.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(gene_mask%2C%20results_de_raw%2C%20results_joint%2C%20scribe)%3A%0A%20%20%20%20%23%20Drop%20mitochondrial%20(incl.%20MT-pseudogenes)%2C%20ribosomal%2C%20and%20hemoglobin%0A%20%20%20%20%23%20genes%20by%20name%20%E2%80%94%20known%20confounders%2C%20chosen%20by%20identity%2C%20not%20by%20effect%20size.%0A%20%20%20%20_nuisance_keep%20%3D%20results_de_raw.exclude_gene_name_mask(%0A%20%20%20%20%20%20%20%20prefixes%3D(%22MT-%22%2C%20%22RPL%22%2C%20%22RPS%22%2C%20%22MRPL%22%2C%20%22MRPS%22)%2C%0A%20%20%20%20%20%20%20%20patterns%3D(r%22%5EMT(RNR%7CCO%5Cd%7CND%5Cd%7CATP%5Cd%7CCYB)%22%2C%20r%22%5EHB%5BABDEGMQZ%5D%5Cd%3F%24%22)%2C%0A%20%20%20%20)%0A%20%20%20%20_gene_mask_clean%20%3D%20gene_mask%20%26%20_nuisance_keep%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20f%22dropping%20%7Bint((gene_mask%20%26%20~_nuisance_keep).sum())%7D%20mito%2Fribo%2FHB%20genes%3B%20%22%0A%20%20%20%20%20%20%20%20f%22keeping%20%7Bint(_gene_mask_clean.sum())%7D%20of%20%7Bgene_mask.size%3A%2C%7D%22%0A%20%20%20%20)%0A%20%20%20%20results_de_clean%20%3D%20scribe.compare_groups(%0A%20%20%20%20%20%20%20%20results_joint%2C%0A%20%20%20%20%20%20%20%20%22perturbation%22%2C%0A%20%20%20%20%20%20%20%20%22control%22%2C%0A%20%20%20%20%20%20%20%20%22panobinostat%22%2C%0A%20%20%20%20%20%20%20%20gene_mask%3D_gene_mask_clean%2C%0A%20%20%20%20%20%20%20%20n_samples%3D5_000%2C%0A%20%20%20%20%20%20%20%20batch_size%3D500%2C%0A%20%20%20%20)%0A%20%20%20%20return%20(results_de_clean%2C)%0A%0A%0A%40app.cell%0Adef%20_(TARGET_PEFP%2C%20TAU%2C%20results_de_clean%2C%20scribe)%3A%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_de_volcano(%0A%20%20%20%20%20%20%20%20results_de_clean%2C%0A%20%20%20%20%20%20%20%20mode%3D%22clr%22%2C%0A%20%20%20%20%20%20%20%20tau%3DTAU%2C%0A%20%20%20%20%20%20%20%20target_pefp%3DTARGET_PEFP%2C%0A%20%20%20%20%20%20%20%20figsize%3D(5%2C%204)%2C%0A%20%20%20%20)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20results_de%2C%20results_de_clean)%3A%0A%20%20%20%20%23%20With%20the%20confounders%20gone%2C%20what%20confident%20CLR%20calls%20remain%20%E2%80%94%20and%20do%20the%0A%20%20%20%20%23%20canonical%20markers%20recover%3F%20Compare%20the%20cleaned%20CLR%20lfsr%20to%20the%20biological%0A%20%20%20%20%23%20lfsr%20for%20the%20same%20markers%2C%20and%20the%20CLR%3C-%3Ebio%20correlation.%0A%20%20%20%20_clean%20%3D%20results_de_clean.to_dataframe(%0A%20%20%20%20%20%20%20%20tau%3Dnp.log(1.1)%2C%0A%20%20%20%20%20%20%20%20target_pefp%3D0.05%2C%0A%20%20%20%20%20%20%20%20metrics%3D%22clr%22%2C%0A%20%20%20%20%20%20%20%20column_naming%3D%22prefixed%22%2C%0A%20%20%20%20)%0A%20%20%20%20_n_de%20%3D%20int(_clean%5B%22clr_is_de%22%5D.sum())%0A%20%20%20%20_top%20%3D%20(%0A%20%20%20%20%20%20%20%20_clean%5B_clean%5B%22clr_lfsr%22%5D%20%3C%3D%200.05%5D%0A%20%20%20%20%20%20%20%20.reindex(%0A%20%20%20%20%20%20%20%20%20%20%20%20_clean%5B_clean%5B%22clr_lfsr%22%5D%20%3C%3D%200.05%5D%5B%22clr_delta_mean%22%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20.abs()%0A%20%20%20%20%20%20%20%20%20%20%20%20.sort_values(ascending%3DFalse)%0A%20%20%20%20%20%20%20%20%20%20%20%20.index%0A%20%20%20%20%20%20%20%20)%5B%5B%22gene%22%2C%20%22clr_delta_mean%22%2C%20%22clr_lfsr%22%5D%5D%0A%20%20%20%20%20%20%20%20.head(10)%0A%20%20%20%20)%0A%20%20%20%20print(f%22confident%20CLR%20calls%20after%20cleaning%20(lfsr%20%3C%3D%200.05)%3A%20%7B_n_de%7D%22)%0A%20%20%20%20_bio%20%3D%20results_de.to_dataframe(%0A%20%20%20%20%20%20%20%20tau%3Dnp.log(1.1)%2C%20metrics%3D%22all%22%2C%20column_naming%3D%22prefixed%22%0A%20%20%20%20)%5B%5B%22gene%22%2C%20%22bio_lfc_mean%22%5D%5D%0A%20%20%20%20_merged%20%3D%20_clean%5B%5B%22gene%22%2C%20%22clr_delta_mean%22%5D%5D.merge(_bio%2C%20on%3D%22gene%22)%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20%22cleaned-CLR%20vs%20biological%20log-fold-change%3A%20Pearson%20%22%0A%20%20%20%20%20%20%20%20f%22%7B_merged%5B'clr_delta_mean'%5D.corr(_merged%5B'bio_lfc_mean'%5D)%3A.3f%7D%2C%20%22%0A%20%20%20%20%20%20%20%20f%22Spearman%20%7B_merged%5B'clr_delta_mean'%5D.corr(_merged%5B'bio_lfc_mean'%5D%2C%20method%3D'spearman')%3A.3f%7D%22%0A%20%20%20%20)%0A%20%20%20%20_top%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20result%20is%20the%20cleanest%20statement%20of%20the%20whole%20section.%20Dropping%20the%20confounders%20does%20**not**%20surface%20the%20coherent%20program%20%E2%80%94%20it%20*empties%20the%20CLR%20hit%20list%20almost%20entirely*.%20The%20mitochondrial%20genes%20are%20gone%2C%20and%20what%20remains%20are%20at%20most%20a%20handful%20of%20confident%20calls%3B%20the%20canonical%20HDAC%20markers%20**keep%20their%20high%20CLR%20lfsr**%20(essentially%20unchanged)%2C%20and%20the%20correlation%20with%20the%20biological%20view%20barely%20moves.%20In%20other%20words%2C%20CLR's%20confidence%20here%20was%20concentrated%20*in%20the%20confounders%20themselves*%20%E2%80%94%20strip%20them%20and%20there%20is%20little%20confident%20*compositional*%20signal%20left%2C%20because%20the%20markers'%20limitation%20was%20never%20the%20reference%20or%20the%20confounders%20but%20their%20genuine%20per-draw%20uncertainty%20under%20a%20shifting%20composition.%0A%0A%20%20%20%20So%20name-based%20filtering%20is%20the%20right%2C%20non-circular%20hygiene%20%E2%80%94%20and%20it%20*confirms*%2C%20rather%20than%20fixes%2C%20the%20conclusion%3A%20for%20a%20broad%20perturbation%2C%20the%20coherent%20program%20is%20not%20hiding%20in%20the%20composition%20behind%20a%20few%20loud%20genes.%20It%20lives%20in%20the%20**biological**%20view%2C%20which%20we%20read%20next.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20The%20biological%20view%2C%20and%20the%20panobinostat%20signature%0A%0A%20%20%20%20The%20complementary%20lens%20is%20the%20**biological%20log-fold-change**%20(%60mode%3D%22bio%22%60)%3A%20the%20change%20in%20each%20gene's%20underlying%20NB%20**mean%20expression**%2C%20with%20its%20own%20lfsr.%20Where%20CLR%20asks%20%22how%20did%20the%20composition%20shift%2C%22%20the%20biological%20view%20asks%20%22how%20did%20this%20gene's%20expression%20change%22%20%E2%80%94%20the%20directly%20interpretable%20question%20for%20a%20global%20drug%20response.%20It%20is%20unaffected%20by%20the%20reference%20drift%20that%20made%20CLR%20uncertain%2C%20so%20the%20canonical%20markers%20that%20carried%20high%20*CLR*%20lfsr%20come%20back%20here%20with%20**biological%20lfsr%20%E2%89%88%200**%20(the%20posterior%20is%20certain%20of%20their%20direction).%20For%20a%20broad%20perturbation%20like%20this%2C%20the%20biological%20view%20is%20the%20**primary**%20readout%20and%20CLR%20is%20the%20compositional%20guardrail%3B%20the%20textbook%20signature%20lives%20here.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(TARGET_PEFP%2C%20TAU%2C%20results_de%2C%20scribe)%3A%0A%20%20%20%20_fig%20%3D%20scribe.viz.plot_de_volcano(%0A%20%20%20%20%20%20%20%20results_de%2C%20mode%3D%22bio%22%2C%20tau%3DTAU%2C%20target_pefp%3DTARGET_PEFP%2C%20figsize%3D(5%2C%204)%0A%20%20%20%20)%0A%20%20%20%20_fig.fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Let's%20read%20the%20actual%20hit%20list.%20We%20export%20the%20full%20table%20and%20orient%20the%20biological%20log-fold-change%20as%20**panobinostat%20vs%20control**%20(positive%20%3D%20up%20after%20treatment).%20One%20honest%20caveat%3A%20a%20fold-change%20ranking%20rewards%20the%20largest%20*relative*%20swings%2C%20and%20for%20the%20least-expressed%20retained%20genes%20those%20are%20the%20noisiest.%20So%20we%20read%20the%20table%20with%20**mean%20expression%20in%20view**%20and%20focus%20the%20headline%20on%20**well-expressed**%20hits%20(mean%20%24%5Cgeq%205%24%20UMIs%20in%20an%20arm)%2C%20then%20validate%20against%20known%20markers%20below.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20results_de)%3A%0A%20%20%20%20_df%20%3D%20results_de.to_dataframe(%0A%20%20%20%20%20%20%20%20tau%3Dnp.log(2.0)%2C%0A%20%20%20%20%20%20%20%20target_pefp%3D0.05%2C%0A%20%20%20%20%20%20%20%20metrics%3D%22all%22%2C%0A%20%20%20%20%20%20%20%20column_naming%3D%22prefixed%22%2C%0A%20%20%20%20)%0A%20%20%20%20%23%20bio_lfc_mean%20is%20log(mu_control%20%2F%20mu_panobinostat)%3B%20flip%20%2B%20rescale%20to%20log2(pano%2Fctrl).%0A%20%20%20%20de_table%20%3D%20_df.assign(%0A%20%20%20%20%20%20%20%20log2fc_pano%3D-_df%5B%22bio_lfc_mean%22%5D%20%2F%20np.log(2.0)%2C%0A%20%20%20%20%20%20%20%20lfsr%3D_df%5B%22bio_lfc_lfsr%22%5D%2C%0A%20%20%20%20%20%20%20%20mu_control%3D_df%5B%22bio_mu_A_mean%22%5D%2C%0A%20%20%20%20%20%20%20%20mu_pano%3D_df%5B%22bio_mu_B_mean%22%5D%2C%0A%20%20%20%20)%5B%5B%22gene%22%2C%20%22log2fc_pano%22%2C%20%22lfsr%22%2C%20%22mu_control%22%2C%20%22mu_pano%22%5D%5D%0A%0A%20%20%20%20%23%20DE%20call%3A%20confident%20direction%20(lfsr%20%3C%3D%200.05)%20and%20at%20least%20a%202-fold%20change.%0A%20%20%20%20_sig%20%3D%20de_table%5B%0A%20%20%20%20%20%20%20%20(de_table%5B%22lfsr%22%5D%20%3C%3D%200.05)%20%26%20(de_table%5B%22log2fc_pano%22%5D.abs()%20%3E%3D%201.0)%0A%20%20%20%20%5D%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20f%22%7Blen(_sig)%7D%20genes%20called%20DE%20(%7Clog2FC%7C%20%3E%3D%201%2C%20lfsr%20%3C%3D%200.05)%3A%20%22%0A%20%20%20%20%20%20%20%20f%22%7Bint((_sig%5B'log2fc_pano'%5D%20%3E%200).sum())%7D%20up%2C%20%22%0A%20%20%20%20%20%20%20%20f%22%7Bint((_sig%5B'log2fc_pano'%5D%20%3C%200).sum())%7D%20down%20in%20panobinostat%22%0A%20%20%20%20)%0A%20%20%20%20%23%20Headline%3A%20the%20well-expressed%20strong%20hits%2C%20split%20by%20direction%20so%20the%20up%20and%0A%20%20%20%20%23%20down%20tables%20don't%20overlap%20(positive%20log2fc%20%3D%20up%20in%20panobinostat).%0A%20%20%20%20_well%20%3D%20_sig%5B(_sig%5B%22mu_control%22%5D%20%3E%3D%205)%20%7C%20(_sig%5B%22mu_pano%22%5D%20%3E%3D%205)%5D%0A%20%20%20%20up_in_pano%20%3D%20(%0A%20%20%20%20%20%20%20%20_well%5B_well%5B%22log2fc_pano%22%5D%20%3E%200%5D%0A%20%20%20%20%20%20%20%20.sort_values(%22log2fc_pano%22%2C%20ascending%3DFalse)%0A%20%20%20%20%20%20%20%20.head(20)%0A%20%20%20%20%20%20%20%20.round(2)%0A%20%20%20%20)%0A%20%20%20%20down_in_pano%20%3D%20(%0A%20%20%20%20%20%20%20%20_well%5B_well%5B%22log2fc_pano%22%5D%20%3C%200%5D%0A%20%20%20%20%20%20%20%20.sort_values(%22log2fc_pano%22)%0A%20%20%20%20%20%20%20%20.head(20)%0A%20%20%20%20%20%20%20%20.round(2)%0A%20%20%20%20)%0A%20%20%20%20up_in_pano%0A%20%20%20%20return%20de_table%2C%20down_in_pano%0A%0A%0A%40app.cell%0Adef%20_(down_in_pano)%3A%0A%20%20%20%20down_in_pano%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Validating%20against%20known%20markers%0A%0A%20%20%20%20The%20cleanest%20sanity%20check%20is%20a%20targeted%20one%3A%20do%20the%20**canonical%20HDAC-inhibitor%20%2F%20panobinostat%20markers**%20come%20back%20with%20the%20sign%20biology%20predicts%2C%20and%20with%20the%20posterior%20confident%20about%20it%3F%20We%20pull%20a%20curated%20panel%20straight%20from%20the%20DE%20table.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(de_table)%3A%0A%20%20%20%20_markers%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%23%20expected%20UP%20after%20panobinostat%0A%20%20%20%20%20%20%20%20%22MT2A%22%2C%0A%20%20%20%20%20%20%20%20%22MT1X%22%2C%0A%20%20%20%20%20%20%20%20%22H1F0%22%2C%0A%20%20%20%20%20%20%20%20%22HIST1H1C%22%2C%0A%20%20%20%20%20%20%20%20%22HSPA1A%22%2C%0A%20%20%20%20%20%20%20%20%22HSPB1%22%2C%0A%20%20%20%20%20%20%20%20%22GADD45A%22%2C%0A%20%20%20%20%20%20%20%20%22DDIT3%22%2C%0A%20%20%20%20%20%20%20%20%22CDKN1A%22%2C%0A%20%20%20%20%20%20%20%20%22NEAT1%22%2C%0A%20%20%20%20%20%20%20%20%23%20expected%20DOWN%20after%20panobinostat%0A%20%20%20%20%20%20%20%20%22MYC%22%2C%0A%20%20%20%20%20%20%20%20%22CCND1%22%2C%0A%20%20%20%20%20%20%20%20%22MDM2%22%2C%0A%20%20%20%20%20%20%20%20%22CD163%22%2C%0A%20%20%20%20%20%20%20%20%22CD14%22%2C%0A%20%20%20%20%20%20%20%20%22MSR1%22%2C%0A%20%20%20%20%20%20%20%20%22FCGR3A%22%2C%0A%20%20%20%20%5D%0A%20%20%20%20marker_panel%20%3D%20(%0A%20%20%20%20%20%20%20%20de_table%5Bde_table%5B%22gene%22%5D.isin(_markers)%5D%0A%20%20%20%20%20%20%20%20.set_index(%22gene%22)%0A%20%20%20%20%20%20%20%20.reindex(_markers)%0A%20%20%20%20%20%20%20%20.round(3)%0A%20%20%20%20)%0A%20%20%20%20marker_panel%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Does%20it%20make%20biological%20sense%3F%0A%0A%20%20%20%20Yes%20%E2%80%94%20emphatically.%20Every%20canonical%20marker%20returns%20with%20the%20expected%20direction%20and%20**lfsr%20%24%5Capprox%200%24**%20(the%20posterior%20is%20essentially%20certain%20of%20the%20sign)%2C%20and%20the%20well-expressed%20hit%20list%20is%20a%20coherent%20program%3A%0A%0A%20%20%20%20**Up%20after%20treatment.**%0A%0A%20%20%20%20-%20**Metallothioneins**%20(%60MT2A%60%3A%20mean%2026%20%E2%86%92%20108%20UMIs%3B%20%60MT1X%60%2C%20%60MT1G%60%2C%20%60MT1H%60%2C%20%60MT1E%60)%20%E2%80%94%20a%20hallmark%20metal%2Foxidative-stress%20response%2C%20a%20dominant%20high-confidence%20block.%0A%20%20%20%20-%20**Linker%20histones**%20(%60H1F0%60%2C%20%60HIST1H1C%60)%20%E2%80%94%20the%20chromatin-remodeling%20response%20that%20is%20the%20defining%20footprint%20of%20HDAC%20inhibition.%0A%20%20%20%20-%20**Heat-shock%20and%20the%20integrated%20stress%20%2F%20p53%20program**%20(%60HSPA1A%60%2C%20%60HSPB1%60%2C%20%60GADD45A%60%2C%20%60DDIT3%60%2C%20%60NEAT1%60%2C%20and%20the%20cell-cycle%20brake%20**%60CDKN1A%60**%2Fp21)%2C%20plus%20acute-phase%20genes%20(%60SAA1%60%2C%20%60PI3%60).%0A%0A%20%20%20%20**Down%20after%20treatment.**%0A%0A%20%20%20%20-%20**%60MYC%60**%20and%20cyclins%20(%60CCND1%60%2C%20%60CCNE2%60)%20%E2%80%94%20proliferation%20shutting%20down%2C%20the%20expected%20cytostatic%20effect.%0A%20%20%20%20-%20**%60MDM2%60**%20%E2%80%94%20the%20p53%20E3%20ligase%20falling%2C%20consistent%20with%20**p53%20activation**.%0A%20%20%20%20-%20A%20coherent%20**myeloid%20%2F%20macrophage%20program**%20(%60CD163%60%2C%20%60CD14%60%2C%20%60MSR1%60%2C%20%60FCGR3A%60%2C%20%60FCER1G%60%2C%20%60SRGN%60)%20%E2%80%94%20the%20drug%20suppressing%2Fdepleting%20the%20monocyte%E2%80%93macrophage%20compartment.%0A%0A%20%20%20%20Crucially%2C%20this%20is%20the%20effect%20**averaged%20over%20seven%20donors%20with%20their%20heterogeneity%20modeled%20away**%20%E2%80%94%20not%20an%20artifact%20of%20one%20dominant%20sample.%20The%20hierarchy%20estimated%20the%20contrast%20we%20asked%20about%20while%20absorbing%20the%20donor%20variation%20we%20did%20not.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20A%20shortcut%3A%20differential%20expression%20straight%20from%20the%20fixed%20effect%0A%0A%20%20%20%20Everything%20above%20used%20%60estimand%3D%22paired_main_effect%22%60%20%E2%80%94%20it%20samples%20per-donor%20compositions%20and%20contrasts%20them%20in%20CLR%20space.%20But%20recall%20what%20the%20model%20encodes%3A%20the%20treatment%20is%20a%20**fixed%2C%20shared%20effect**%20%24%5Calpha%24%2C%20and%20within%20any%20donor%20the%20baseline%20%24%5Cbeta%24%20cancels%2C%20so%20the%20biological%20log-fold-change%20*is*%20%24%5Calpha%5B%5Ctext%7Bpanobinostat%7D%5D%20-%20%5Calpha%5B%5Ctext%7Bcontrol%7D%5D%24%20%E2%80%94%20the%20same%20for%20every%20donor.%20We%20can%20therefore%20read%20the%20treatment%20DEG%20**straight%20from%20the%20fitted%20%24%5Calpha%24**%2C%20with%20no%20composition%20sampling%20at%20all.%0A%0A%20%20%20%20%60scribe.compare_groups(...%2C%20estimand%3D%22effect%22)%60%20does%20exactly%20that%3A%20it%20pulls%20the%20%24%5Calpha%24%20contrast%20from%20the%20hierarchy%20(%60get_factor_effect%60%20under%20the%20hood)%2C%20so%20it%20is%20essentially%20**instant%20and%20memory-light**%20%E2%80%94%20the%20practical%20choice%20when%20the%20full%20paired-CLR%20pass%20is%20too%20heavy%2C%20and%20the%20*coherent*%20view%20for%20a%20broad%20perturbation%20like%20this.%20The%20returned%20%60delta%60%20is%20now%20a%20**log-mean**%20effect%20(not%20a%20CLR%20contrast)%2C%20carrying%20its%20own%20lfsr%3B%20the%20%60_other%60%20pseudo-gene%20is%20dropped%20automatically.%20Because%20it%20is%20the%20same%20quantity%20as%20the%20biological%20log-fold-change%20computed%20the%20long%20way%2C%20the%20two%20agree%20to%20numerical%20precision%20(we%20check%20the%20correlation%20below).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(gene_mask%2C%20np%2C%20results_de%2C%20results_joint%2C%20scribe)%3A%0A%20%20%20%20%23%20The%20treatment%20effect%2C%20read%20directly%20from%20the%20fitted%20alpha%20%E2%80%94%20no%20composition%0A%20%20%20%20%23%20sampling.%20Reuses%20the%20already-drawn%20posterior%2C%20so%20this%20is%20near-instant.%0A%20%20%20%20results_effect%20%3D%20scribe.compare_groups(%0A%20%20%20%20%20%20%20%20results_joint%2C%0A%20%20%20%20%20%20%20%20%22perturbation%22%2C%0A%20%20%20%20%20%20%20%20%22control%22%2C%0A%20%20%20%20%20%20%20%20%22panobinostat%22%2C%0A%20%20%20%20%20%20%20%20estimand%3D%22effect%22%2C%0A%20%20%20%20%20%20%20%20gene_mask%3Dgene_mask%2C%0A%20%20%20%20%20%20%20%20n_samples%3D5_000%2C%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20It%20is%20the%20same%20quantity%20as%20the%20biological%20log-fold-change%20(clr_*%20columns%0A%20%20%20%20%23%20here%20hold%20the%20log-mean%20effect%2C%20not%20a%20CLR%20contrast)%2C%20so%20they%20agree%3A%0A%20%20%20%20_eff%20%3D%20results_effect.to_dataframe(%0A%20%20%20%20%20%20%20%20tau%3Dnp.log(2.0)%2C%0A%20%20%20%20%20%20%20%20target_pefp%3D0.05%2C%0A%20%20%20%20%20%20%20%20metrics%3D%22clr%22%2C%0A%20%20%20%20%20%20%20%20column_naming%3D%22prefixed%22%2C%0A%20%20%20%20)%0A%20%20%20%20_bio%20%3D%20results_de.to_dataframe(%0A%20%20%20%20%20%20%20%20tau%3Dnp.log(2.0)%2C%20metrics%3D%22all%22%2C%20column_naming%3D%22prefixed%22%0A%20%20%20%20)%5B%5B%22gene%22%2C%20%22bio_lfc_mean%22%5D%5D%0A%20%20%20%20_merged%20%3D%20_eff.merge(_bio%2C%20on%3D%22gene%22)%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20%22effect-estimand%20delta%20vs%20biological%20log-fold-change%3A%20corr%20%3D%20%22%0A%20%20%20%20%20%20%20%20f%22%7B_merged%5B'clr_delta_mean'%5D.corr(_merged%5B'bio_lfc_mean'%5D)%3A.4f%7D%20%20%22%0A%20%20%20%20%20%20%20%20f%22(%7Bint(_eff%5B'clr_is_de'%5D.sum())%7D%20genes%20called%20DE)%22%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20The%20canonical%20markers%2C%20recovered%20by%20the%20fast%20alpha%20path%20(positive%0A%20%20%20%20%23%20log2FC%20%3D%20up%20in%20panobinostat)%20%E2%80%94%20the%20same%20signature%20as%20the%20bio%20view%20above.%0A%20%20%20%20effect_de%20%3D%20_eff.assign(log2fc_pano%3D-_eff%5B%22clr_delta_mean%22%5D%20%2F%20np.log(2.0))%5B%0A%20%20%20%20%20%20%20%20%5B%22gene%22%2C%20%22log2fc_pano%22%2C%20%22clr_lfsr%22%5D%0A%20%20%20%20%5D%0A%20%20%20%20_markers%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%22MT2A%22%2C%0A%20%20%20%20%20%20%20%20%22MT1X%22%2C%0A%20%20%20%20%20%20%20%20%22H1F0%22%2C%0A%20%20%20%20%20%20%20%20%22HIST1H1C%22%2C%0A%20%20%20%20%20%20%20%20%22CDKN1A%22%2C%0A%20%20%20%20%20%20%20%20%22MYC%22%2C%0A%20%20%20%20%20%20%20%20%22MDM2%22%2C%0A%20%20%20%20%20%20%20%20%22CD163%22%2C%0A%20%20%20%20%20%20%20%20%22CD14%22%2C%0A%20%20%20%20%20%20%20%20%22MSR1%22%2C%0A%20%20%20%20%20%20%20%20%22FCGR3A%22%2C%0A%20%20%20%20%5D%0A%20%20%20%20effect_de%5Beffect_de%5B%22gene%22%5D.isin(_markers)%5D.set_index(%22gene%22).reindex(%0A%20%20%20%20%20%20%20%20%5Bm%20for%20m%20in%20_markers%20if%20m%20in%20set(effect_de%5B%22gene%22%5D)%5D%0A%20%20%20%20).round(3)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20CLR%20and%20the%20biological%20view%20agree%20%E2%80%94%20conservative%2C%20not%20incoherent%0A%0A%20%20%20%20We%20have%20treated%20the%20compositional%20(CLR)%20and%20biological%20views%20as%20answering%20different%20questions%20%E2%80%94%20and%20they%20do%20%E2%80%94%20but%20it%20is%20worth%20asking%20directly%3A%20do%20they%20*disagree*%3F%20They%20do%20not.%20Plotting%20every%20gene's%20CLR%20contrast%20against%20its%20biological%20log-fold-change%20(left)%20shows%20a%20clear%20positive%20correlation%3A%20the%20canonical%20markers%20fall%20in%20the%20expected%20quadrants%20(metallothioneins%20%2F%20histones%20%2F%20p53-stress%20up%2C%20the%20myeloid%20program%20down)%2C%20and%20the%20two%20lenses%20rank%20the%20response%20the%20same%20way.%20The%20scatter%20also%20makes%20**compositional%20closure**%20visible%20%E2%80%94%20the%20cloud%20of%20genes%20that%20rise%20in%20absolute%20expression%20yet%20*fall*%20in%20proportion%2C%20squeezed%20as%20the%20high-mass%20drivers%20claim%20more%20of%20the%20simplex.%20That%20is%20real%20compositional%20structure%2C%20not%20error.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20pd%2C%20plt%2C%20results_de)%3A%0A%20%20%20%20_df%20%3D%20results_de.to_dataframe(%0A%20%20%20%20%20%20%20%20tau%3Dnp.log(2.0)%2C%20metrics%3D%22all%22%2C%20column_naming%3D%22prefixed%22%0A%20%20%20%20)%0A%20%20%20%20%23%20Orient%20both%20contrasts%20as%20panobinostat%20vs%20control%20(positive%20%3D%20up%20in%20pano).%0A%20%20%20%20_bio%20%3D%20-_df%5B%22bio_lfc_mean%22%5D.to_numpy()%20%2F%20np.log(2.0)%0A%20%20%20%20_clr%20%3D%20-_df%5B%22clr_delta_mean%22%5D.to_numpy()%20%2F%20np.log(2.0)%0A%20%20%20%20_lfsr%20%3D%20_df%5B%22clr_lfsr%22%5D.to_numpy()%0A%20%20%20%20_rho%20%3D%20pd.Series(_clr).corr(pd.Series(_bio)%2C%20method%3D%22spearman%22)%0A%20%20%20%20_g2i%20%3D%20%7Bg%3A%20i%20for%20i%2C%20g%20in%20enumerate(_df%5B%22gene%22%5D.astype(str))%7D%0A%0A%20%20%20%20%23%20As%20we%20relax%20the%20tolerated%20false%20proportion%20(PEFP)%2C%20how%20much%20of%20the%0A%20%20%20%20%23%20canonical%20signature%20does%20CLR%20recover%2C%20and%20do%20the%20calls%20keep%20the%20right%20sign%3F%0A%20%20%20%20_markers%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%22MT2A%22%2C%0A%20%20%20%20%20%20%20%20%22MT1X%22%2C%0A%20%20%20%20%20%20%20%20%22H1F0%22%2C%0A%20%20%20%20%20%20%20%20%22HIST1H1C%22%2C%0A%20%20%20%20%20%20%20%20%22CDKN1A%22%2C%0A%20%20%20%20%20%20%20%20%22NEAT1%22%2C%0A%20%20%20%20%20%20%20%20%22MYC%22%2C%0A%20%20%20%20%20%20%20%20%22MDM2%22%2C%0A%20%20%20%20%20%20%20%20%22CD163%22%2C%0A%20%20%20%20%20%20%20%20%22CD14%22%2C%0A%20%20%20%20%20%20%20%20%22FCGR3A%22%2C%0A%20%20%20%20%20%20%20%20%22MSR1%22%2C%0A%20%20%20%20%5D%0A%20%20%20%20_mset%20%3D%20set(_markers)%0A%20%20%20%20_pefps%20%3D%20%5B0.05%2C%200.10%2C%200.15%2C%200.20%2C%200.25%2C%200.30%2C%200.35%2C%200.40%5D%0A%20%20%20%20_nmark%2C%20_sign%20%3D%20%5B%5D%2C%20%5B%5D%0A%20%20%20%20for%20_pe%20in%20_pefps%3A%0A%20%20%20%20%20%20%20%20_d%20%3D%20results_de.to_dataframe(%0A%20%20%20%20%20%20%20%20%20%20%20%20tau%3Dnp.log(2.0)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20target_pefp%3D_pe%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20metrics%3D%22all%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20column_naming%3D%22prefixed%22%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20_c%20%3D%20_d%5B_d%5B%22clr_is_de%22%5D%5D%0A%20%20%20%20%20%20%20%20_nmark.append(len(set(_c%5B%22gene%22%5D.astype(str))%20%26%20_mset))%0A%20%20%20%20%20%20%20%20_s%20%3D%20np.sign(_c%5B%22clr_delta_mean%22%5D)%20%3D%3D%20np.sign(_c%5B%22bio_lfc_mean%22%5D)%0A%20%20%20%20%20%20%20%20_sign.append(float(_s.mean())%20if%20len(_c)%20else%20np.nan)%0A%0A%20%20%20%20_fig%2C%20(_axA%2C%20_axB)%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(11%2C%204.4))%0A%20%20%20%20_sc%20%3D%20_axA.scatter(_bio%2C%20_clr%2C%20c%3D_lfsr%2C%20s%3D6%2C%20cmap%3D%22viridis%22%2C%20alpha%3D0.5)%0A%20%20%20%20_fig.colorbar(_sc%2C%20ax%3D_axA%2C%20label%3D%22CLR%20lfsr%20(sign%20uncertainty)%22)%0A%20%20%20%20for%20_g%20in%20_markers%3A%0A%20%20%20%20%20%20%20%20if%20_g%20in%20_g2i%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_i%20%3D%20_g2i%5B_g%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20_axA.scatter(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_bio%5B_i%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_clr%5B_i%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20c%3D%22crimson%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20s%3D22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20edgecolor%3D%22k%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20lw%3D0.4%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20zorder%3D5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%20%20%20%20_axA.annotate(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_g%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20(_bio%5B_i%5D%2C%20_clr%5B_i%5D)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fontsize%3D6.5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20color%3D%22crimson%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20xytext%3D(2%2C%202)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20textcoords%3D%22offset%20points%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20_axA.axhline(0%2C%20c%3D%22k%22%2C%20lw%3D0.5)%0A%20%20%20%20_axA.axvline(0%2C%20c%3D%22k%22%2C%20lw%3D0.5)%0A%20%20%20%20_axA.set_xlabel(%22biological%20log2%20fold-change%20(pano%20%2F%20control)%22)%0A%20%20%20%20_axA.set_ylabel(%22CLR%20contrast%20(log2%20units%2C%20pano%20%2F%20control)%22)%0A%20%20%20%20_axA.set_title(f%22CLR%20vs%20biological%20DE%20%E2%80%94%20Spearman%20%7B_rho%3A.2f%7D%22)%0A%0A%20%20%20%20_axB.plot(_pefps%2C%20_nmark%2C%20%22o-%22%2C%20color%3D%22crimson%22)%0A%20%20%20%20_axB.set_xlabel(%22target%20PEFP%20(tolerated%20false%20proportion)%22)%0A%20%20%20%20_axB.set_ylabel(%22%23%20canonical%20markers%20called%22%2C%20color%3D%22crimson%22)%0A%20%20%20%20_axB.tick_params(axis%3D%22y%22%2C%20labelcolor%3D%22crimson%22)%0A%20%20%20%20_axB2%20%3D%20_axB.twinx()%0A%20%20%20%20_axB2.plot(_pefps%2C%20_sign%2C%20%22s--%22%2C%20color%3D%22steelblue%22)%0A%20%20%20%20_axB2.set_ylabel(%0A%20%20%20%20%20%20%20%20%22fraction%20of%20CLR%20calls%20agreeing%20with%20bio%20sign%22%2C%20color%3D%22steelblue%22%0A%20%20%20%20)%0A%20%20%20%20_axB2.set_ylim(0.7%2C%201.005)%0A%20%20%20%20_axB.set_title(%0A%20%20%20%20%20%20%20%20%22Conservative%2C%20not%20incoherent%3A%5Cnthe%20signature%20emerges%20as%20PEFP%20relaxes%22%0A%20%20%20%20)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20What%20separates%20the%20two%20views%20is%20**confidence%2C%20not%20direction.**%20CLR%20assigns%20the%20markers%20higher%20lfsr%2C%20so%20at%20a%20strict%20PEFP%20only%20the%20few%20ultra-confident%20driver%20genes%20clear%20the%20bar.%20But%20that%20is%20a%20property%20of%20the%20*threshold*%2C%20not%20the%20signal%3A%20relax%20the%20tolerated%20false%20proportion%20(right%20panel)%20and%20the%20full%20canonical%20signature%20emerges%20from%20CLR%2C%20the%20calls%20still%20agreeing%20with%20biology%20in%20sign.%0A%0A%20%20%20%20And%20CLR%20is%20the%20more%20conservative%20lens%20for%20a%20principled%20reason%20%E2%80%94%20it%20**assumes%20less%20and%20accounts%20for%20more.**%20It%20is%20provably%20invariant%20to%20capture-probability%20misspecification%20(the%20biological%20mean%20is%20*not*%20%E2%80%94%20it%20must%20trust%20the%20capture%20%2F%20library%20model)%3B%20it%20references%20each%20gene%20against%20the%20whole%20composition%3B%20and%20under%20the%20gene-specific%20%24p_g%24%20model%20it%20carries%20*shape*%20%2F%20overdispersion%20uncertainty%20into%20its%20lfsr%2C%20a%20term%20that%20simply%20cancels%20in%20the%20biological%20log-fold-change.%20Its%20wider%20intervals%20are%20**earned%20honesty**%20about%20what%20compositional%20data%20alone%20can%20claim%20%E2%80%94%20not%20noise.%0A%0A%20%20%20%20The%20scientific%20reading%2C%20then%2C%20is%20not%20%22believe%20whichever%20lens%20clears%20a%20%24p%24-value-like%20cutoff.%22%20It%20is%20that%20a%20**robust%2C%20assumption-light%20estimator%20(CLR)%20and%20a%20confident%2C%20model-dependent%20one%20(the%20biological%20mean)%20converge**%20on%20the%20same%20panobinostat%20signature.%20That%20convergence%20%E2%80%94%20read%20with%20knowledge%20of%20the%20marker%20biology%20%E2%80%94%20is%20what%20earns%20confidence%20in%20the%20result%2C%20far%20more%20than%20any%20single%20PEFP%20threshold%20could.%20The%20generative%20model%20hands%20us%20*both*%20lenses%20from%20one%20posterior%3B%20the%20scientific%20judgment%20is%20ours%20to%20bring.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Recap%0A%0A%20%20%20%20We%20built%20**one**%20joint%20generative%20model%20of%20a%20crossed%20*donor%20%C3%97%20condition*%20experiment%20and%2C%20from%20a%20single%20posterior%2C%20read%20off%3A%0A%0A%20%20%20%201.%20a%20**fit**%20that%20calibrates%20per%20leaf%20and%20passes%20posterior%20predictive%20checks%3B%0A%20%20%20%202.%20the%20**fitted%20structure**%20%E2%80%94%20a%20fixed%20treatment%20effect%20and%20the%20random%20donor%20deviations%20it%20was%20protected%20from%2C%20both%20as%20inspectable%20posteriors%3B%0A%20%20%20%203.%20a%20**donor-averaged%2C%20paired%20differential-expression**%20result%20whose%20**biological**%20log-fold-changes%20reproduce%20the%20canonical%20panobinostat%20%2F%20HDAC-inhibitor%20signature%20with%20high%20sign-confidence%20(lfsr%20%E2%89%88%200%20on%20every%20marker)%20%E2%80%94%20and%20the%20same%20calls%2C%20essentially%20for%20free%2C%20by%20reading%20the%20treatment%20effect%20%24%5Calpha%24%20directly%20(%60estimand%3D%22effect%22%60).%0A%0A%20%20%20%20Three%20lessons%20worth%20carrying%20forward.%20First%2C%20**CLR%20is%20principled%20but%20needs%20a%20stable%20reference**%3A%20pooling%20near-silent%20genes%20into%20%22other%22%20(an%20exact%20Dirichlet%20operation)%20is%20what%20lets%20the%20compositional%20contrast%20see%20signal%20instead%20of%20noise.%20Second%2C%20**CLR%20and%20the%20biological%20view%20answer%20different%20questions%2C%20and%20which%20one%20is%20the%20*confident*%20headline%20depends%20on%20the%20perturbation.**%20CLR%20asks%20%22did%20the%20*relative*%20composition%20shift%3F%22%20%E2%80%94%20invaluable%20as%20a%20guard%20against%20compositional%20artifacts%2C%20but%20under%20a%20broad%2C%20genome-wide%20response%20(like%20this%20one)%20the%20only%20calls%20it%20is%20*confident*%20about%20are%20the%20genes%20whose%20proportion%20swings%20hardest%20(here%2C%20mitochondrial%20genes).%20The%20biological%20log-fold-change%20asks%20%22did%20the%20*mean%20expression*%20change%3F%22%20%E2%80%94%20unaffected%20by%20the%20reference%20drift%2C%20and%20where%20the%20textbook%20signature%20emerged%20with%20high%20confidence.%20Third%20%E2%80%94%20and%20this%20is%20the%20subtle%20one%20%E2%80%94%20**conservative%20is%20not%20the%20same%20as%20wrong.**%20CLR%20is%20*under-confident%20here%2C%20not%20incoherent*%3A%20it%20agrees%20with%20the%20biological%20view%20in%20direction%2C%20recovers%20the%20full%20signature%20as%20the%20error%20tolerance%20relaxes%2C%20and%20is%20conservative%20precisely%20because%20it%20assumes%20less%20(capture-invariant)%20and%20accounts%20for%20more%20(global%20structure%2C%20gene-specific%20shape).%20The%20trustworthy%20conclusion%20comes%20from%20the%20**convergence**%20of%20the%20robust%20and%20the%20confident%20lens%20%E2%80%94%20read%20with%20scientific%20judgment%20%E2%80%94%20not%20from%20whichever%20one%20clears%20a%20threshold.%20The%20hierarchy%20ties%20it%20together%3A%20it%20estimates%20the%20effect%20we%20care%20about%20while%20explicitly%20accounting%20for%20the%20donor%20variation%20we%20do%20not.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A