| """per-class pathway module scoring cKO/mutant vs WT across pan_skin, hematopoiesis, pancreas.
|
| pancreas contrast is Veres stage 6 vs stage 5 (HUMAN gene symbols)."""
|
| from __future__ import annotations
|
|
|
| import argparse
|
| import warnings
|
| from pathlib import Path
|
|
|
| import anndata as ad
|
| import numpy as np
|
| import pandas as pd
|
| import scanpy as sc
|
| import scipy.sparse as sp
|
| from scipy.stats import mannwhitneyu
|
|
|
| warnings.filterwarnings("ignore")
|
| sc.settings.verbosity = 0
|
|
|
| import os as _os
|
| from pathlib import Path as _Path
|
| PANDA_ROOT = _Path(_os.environ.get("PANDA_ROOT", str(_Path(__file__).resolve().parents[2])))
|
| ROOT = Path(str(PANDA_ROOT))
|
| MIN_PER_GROUP = 15
|
|
|
|
|
| def _M(genes, type_, citation, direction=None):
|
| return {"genes": list(genes), "type": type_, "citation": citation,
|
| "direction": direction}
|
|
|
|
|
|
|
| SKIN_MODULES = {
|
| "MITF_regulon": _M(["Mitf","Dct","Tyr","Pmel","Mlana","Tyrp1","Slc24a5",
|
| "Slc45a2","Sox10","Pax3","Kit","Ednrb"],
|
| "lineage", "Steingrimsson-2004"),
|
| "Wnt_signaling": _M(["Wnt3","Wnt5a","Wnt7a","Wnt10b","Ctnnb1","Lef1",
|
| "Tcf4","Tcf7","Axin2","Dkk1","Sfrp1","Fzd7","Lrp5"],
|
| "signaling", "Nusse-2017"),
|
| "BMP_signaling": _M(["Bmp2","Bmp4","Bmp5","Bmp7","Bmpr1a","Bmpr1b",
|
| "Bmpr2","Smad1","Smad5","Id1","Id2","Id3"],
|
| "signaling", "Botchkarev-2003"),
|
| "TGFB_signaling": _M(["Tgfb1","Tgfb2","Tgfbr1","Tgfbr2","Smad3","Smad7"],
|
| "signaling", "Massague-2012"),
|
| "FGF_signaling": _M(["Fgf1","Fgf2","Fgf7","Fgf9","Fgf10","Fgfr1","Fgfr2",
|
| "Etv1","Etv4","Etv5","Spry2","Dusp6"],
|
| "signaling", "Ornitz-2015"),
|
| "Notch_signaling": _M(["Notch1","Notch2","Notch3","Jag1","Dll1","Hes1",
|
| "Hes5","Hey1","Hey2","Rbpj"],
|
| "signaling", "Andersson-2011"),
|
| "Hedgehog": _M(["Shh","Ptch1","Smo","Gli1","Gli2","Gli3"],
|
| "signaling", "St-Jacques-1998"),
|
| "Eda_ectodysplasin": _M(["Eda","Edar","Edaradd","Nfkb1","Nfkb2","Rela"],
|
| "signaling", "Mikkola-2009"),
|
| "EMT": _M(["Zeb1","Zeb2","Snai1","Snai2","Twist1","Twist2",
|
| "Vim","Cdh2","Fn1","Prrx1"],
|
| "lineage", "Thiery-2009"),
|
| "Cell_cycle": _M(["Ccnd1","Ccne1","Ccna2","Ccnb1","Cdk1","Cdk2",
|
| "Cdk4","Mki67","Top2a","Pcna","Mcm2","Mcm3"],
|
| "cycle", "Whitfield-2002"),
|
| "KC_differentiation":_M(["Krt1","Krt10","Ivl","Lor","Flg","Flg2","Klk5",
|
| "Klk7","Cdsn"],
|
| "lineage", "Fuchs-2007"),
|
| "Basal_keratinocyte":_M(["Krt5","Krt14","Krt15","Trp63","Itga6","Itgb1",
|
| "Itga3"],
|
| "lineage", "Blanpain-2007"),
|
| "Sweat_gland": _M(["Foxi3","Foxa1","En1","Krt8","Krt18","Krt19",
|
| "Muc5b","Aqp5","Cutl1"],
|
| "lineage", "Lu-2016"),
|
| "Hair_placode": _M(["Wnt10b","Shh","Lef1","Foxi3","Edar","Bmp4","Msx2"],
|
| "lineage", "Millar-2002"),
|
| "Neural_crest": _M(["Sox10","Sox9","Sox2","Pax3","Foxd3","Nes","Tfap2a"],
|
| "lineage", "Simoes-Costa-2015"),
|
| "Apoptosis": _M(["Bax","Bak1","Bad","Bcl2","Casp3","Casp9","Trp53",
|
| "Cdkn1a"],
|
| "stress", "Youle-2008"),
|
| "Melanogenesis_late":_M(["Tyrp1","Slc45a2","Oca2","Gpnmb","Pmel","Silv",
|
| "Mlph","Rab27a","Melana"],
|
| "lineage", "Raposo-2007"),
|
| "Sebogenesis": _M(["Elovl3","Awat2","Scd1","Scd3","Adipoq","Fasn",
|
| "Mgst1","Srebf1","Pparg"],
|
| "metabolism", "Zouboulis-2016"),
|
| "Immune_Th1_Th2_Th17":_M(["Tbx21","Gata3","Rorc","Ifng","Il4","Il13","Il17a",
|
| "Il17f","Il22","Foxp3"],
|
| "immune", "Zhu-2010"),
|
| "DNA_damage": _M(["Trp53","Cdkn1a","Atm","Atr","Brca1","Chek1",
|
| "Chek2","Rad51","Mre11a","H2ax","Nbn"],
|
| "stress", "Ciccia-2010"),
|
| "Autophagy": _M(["Atg5","Atg7","Atg12","Becn1","Map1lc3b","Sqstm1",
|
| "Ulk1","Atg3","Atg16l1"],
|
| "stress", "Mizushima-2011"),
|
| "Senescence": _M(["Cdkn2a","Cdkn2b","Cdkn1a","Il6","Cxcl1","Serpine1",
|
| "Glb1","Lmnb1"],
|
| "stress", "Coppe-2010"),
|
| "Epidermal_junction":_M(["Cdh1","Dsg1a","Dsg1b","Dsg2","Dsg3","Cldn1",
|
| "Cldn4","Ocln","Tjp1","Cldn23","Dsp","Pkp1"],
|
| "junction", "Green-2010"),
|
| "ECM_collagen": _M(["Col1a1","Col1a2","Col3a1","Col4a1","Col6a1",
|
| "Col17a1","Lum","Dcn","Fbn1","Postn","Fbln1"],
|
| "ecm", "Ricard-Blum-2011"),
|
| "Endothelial_tip_stalk":_M(["Dll4","Notch1","Hes1","Kdr","Kit","Cxcr4",
|
| "Angpt2","Nrp1","Flt1","Cdh5","Pecam1"],
|
| "lineage", "Blanco-2013"),
|
| "Fibroblast_wound": _M(["Postn","Fap","Aspn","Tnc","Acta2","Prrx1","Pdgfra",
|
| "Ly6a"],
|
| "lineage", "Rinkevich-2015"),
|
| "Fatty_acid_oxidation":_M(["Cpt1a","Acadm","Acadl","Acadvl","Hadha","Hadhb",
|
| "Ppara","Ppargc1a","Ucp2"],
|
| "metabolism", "Houten-2010"),
|
| "Nrf2_oxidative_stress":_M(["Nfe2l2","Nqo1","Gclc","Hmox1","Slc7a11",
|
| "Txnrd1","Gsta3","Gpx2","Keap1"],
|
| "stress", "Ma-2013"),
|
| "IFN_gamma": _M(["Ifng","Stat1","Ifit1","Ifit2","Ifit3","Isg15",
|
| "Irf1","Cxcl9","Cxcl10","Gbp2"],
|
| "immune", "Schoggins-2011"),
|
| "IL6_JAK_STAT": _M(["Il6","Stat3","Socs3","Jak1","Jak2","Il6ra",
|
| "Il6st","Osm"],
|
| "signaling", "Heinrich-2003"),
|
| "Pigment_regulation":_M(["Asip","Kitl","Kit","Bcl2","Mc1r","Pomc","Adcy8",
|
| "Ednrb","Edn3"],
|
| "signaling", "Slominski-2004"),
|
| }
|
|
|
|
|
|
|
| HSC_MODULES = {
|
| "Kit_signaling": _M(["Kit","Kitl","Sox4","Gata2","Runx1","Meis1"],
|
| "signaling", "Lennartsson-2012"),
|
| "Kit_ligand": _M(["Kit","Kitl"], "signaling", "Broudy-1997"),
|
| "MYC_targets": _M(["Myc","Nolc1","Nop58","Ncl","Npm1","Fbl","Eif4e",
|
| "Nop56","Ldha","Odc1"],
|
| "lineage", "Dang-2012"),
|
| "Cell_cycle": _M(["Ccnd1","Ccne1","Ccna2","Ccnb1","Cdk1","Cdk2",
|
| "Cdk4","Mki67","Top2a","Pcna","Mcm2","Mcm3","Mcm5"],
|
| "cycle", "Whitfield-2002"),
|
| "DNA_replication": _M(["Mcm2","Mcm3","Mcm4","Mcm5","Mcm6","Mcm7","Pcna",
|
| "Rfc4","Pola1","Pole","Rpa1","Rpa2"],
|
| "cycle", "Bell-2002"),
|
| "Integrated_stress": _M(["Atf4","Ddit3","Ppp1r15a","Ppp1r15b","Eif2ak3",
|
| "Eif2s1","Atf3","Gadd45a"],
|
| "stress", "Pakos-Zebrucka-2016"),
|
| "Apoptosis_pro": _M(["Bax","Bak1","Bid","Bad","Bim","Puma","Noxa",
|
| "Casp3","Casp9"],
|
| "stress", "Youle-2008"),
|
| "Apoptosis_anti": _M(["Bcl2","Bcl2l1","Mcl1","Bcl2l2","Bcl2l10","Xiap"],
|
| "stress", "Adams-2018"),
|
| "Erythropoiesis_early":_M(["Gata1","Klf1","Epo","Epor","Tal1","Zfpm1",
|
| "Gypa","Lmo2"],
|
| "lineage", "Palis-2014"),
|
| "Erythropoiesis_late":_M(["Alas2","Hba-a1","Hba-a2","Hbb-b1","Hbb-b2","Slc4a1",
|
| "Ank1","Blvrb","Car1","Car2"],
|
| "lineage", "Palis-2014"),
|
| "Granulopoiesis": _M(["Cebpa","Cebpe","Elane","Mpo","Prtn3","Csf3r",
|
| "S100a8","S100a9","Ctsg","Ltf","Lcn2","Mmp8"],
|
| "lineage", "Rosenbauer-2007"),
|
| "Lymphopoiesis_B": _M(["Rag1","Rag2","Dntt","Vpreb1","Vpreb3","Igll1",
|
| "Cd19","Pax5","Ebf1"],
|
| "lineage", "Nutt-2011"),
|
| "Lymphopoiesis_T": _M(["Il7r","Cd3d","Cd3e","Cd3g","Lck","Zap70","Gata3",
|
| "Tcf7","Runx3"],
|
| "lineage", "Rothenberg-2014"),
|
| "Megakaryopoiesis": _M(["Nfe2","Gata1","Fli1","Runx1","Itga2b","Pf4",
|
| "Gp1bb","Gp9","Mpl","Vwf"],
|
| "lineage", "Tijssen-2013"),
|
| "Basophil_mast": _M(["Cpa3","Ms4a2","Gata2","Hdc","Mcpt8","Prss34",
|
| "Fcer1a","Il4","Il6"],
|
| "lineage", "Voehringer-2013"),
|
| "Hemostasis": _M(["Vwf","F5","F13a1","Fga","Fgb","Fgg","Serpine1",
|
| "Plat","Plau","Plg"],
|
| "signaling", "Furie-2008"),
|
| "OXPHOS_ETC": _M(["Ndufa1","Ndufa2","Ndufb1","Ndufb2","Sdha","Sdhb",
|
| "Cox4i1","Cox5a","Cox6a1","Atp5a1","Atp5b","Uqcrq"],
|
| "metabolism", "Mishra-2016"),
|
| "Glycolysis": _M(["Hk1","Hk2","Pfkm","Pfkl","Aldoa","Gapdh","Pgk1",
|
| "Pkm","Ldha","Eno1","Tpi1","Pgam1"],
|
| "metabolism", "Vander-Heiden-2009"),
|
| "TCA": _M(["Cs","Aco2","Idh2","Idh3a","Sdha","Fh1","Mdh2",
|
| "Ogdh","Sucla2"],
|
| "metabolism", "Chandel-2015"),
|
| "Redox_glutathione": _M(["Gpx1","Gpx2","Gpx3","Gpx4","Gsr","Prdx1","Prdx2",
|
| "Prdx3","Prdx4","Prdx5","Prdx6","Sod1","Sod2","Cat"],
|
| "stress", "Ho-2007"),
|
| "Wnt_hemato": _M(["Wnt3a","Wnt5a","Ctnnb1","Lef1","Tcf7","Axin2",
|
| "Fzd4","Fzd7"],
|
| "signaling", "Reya-2003"),
|
| "Notch_hemato": _M(["Notch1","Notch2","Jag1","Hes1","Dll1","Dll4",
|
| "Rbpj","Hey1"],
|
| "signaling", "Bigas-2018"),
|
| "TGFb_hemato": _M(["Tgfb1","Tgfb2","Tgfbr1","Tgfbr2","Smad2","Smad3",
|
| "Smad4","Smad7"],
|
| "signaling", "Blank-2015"),
|
| "IFN_signaling": _M(["Ifnar1","Ifnar2","Stat1","Stat2","Ifit1","Ifit2",
|
| "Ifit3","Isg15","Irf7","Mx1"],
|
| "immune", "Essers-2009"),
|
| "Complement": _M(["C1qa","C1qb","C1qc","C3","C4b","Cfp","Cfh","Cfd"],
|
| "immune", "Ricklin-2016"),
|
| "NK_cytotoxicity": _M(["Ncr1","Klrk1","Prf1","Gzmb","Gzmk","Nkg7","Klrd1",
|
| "Klrb1c","Klra8"],
|
| "immune", "Vivier-2011"),
|
| "Mast_cell_degran": _M(["Ms4a2","Fcer1a","Cpa3","Kit","Hdc","Tpsb2",
|
| "Prss34","Mcpt4"],
|
| "immune", "Galli-2011"),
|
| "Autophagy": _M(["Atg5","Atg7","Atg12","Becn1","Map1lc3b","Sqstm1",
|
| "Ulk1","Atg3","Atg16l1"],
|
| "stress", "Warr-2013"),
|
| "Senescence": _M(["Cdkn2a","Cdkn2b","Cdkn1a","Il6","Cxcl1","Serpine1",
|
| "Glb1","Lmnb1"],
|
| "stress", "Chang-2016"),
|
| "LT_HSC_quiescence": _M(["Hlf","Meis1","Mecom","Procr","Fgd5","Mllt3","Egr1",
|
| "Rgs1","Cdkn1c","Ndn","Mpl"],
|
| "lineage", "Cabezas-Wallscheid-2017"),
|
| }
|
|
|
|
|
|
|
| PANCREAS_MODULES = {
|
| "Insulin_secretion": _M(["INS","IAPP","CHGA","CHGB","SCG5","ERO1B","PCSK1",
|
| "PCSK2","SLC30A8","G6PC2"],
|
| "hormone", "Rorsman-2013"),
|
| "Glucose_sensing": _M(["SLC2A2","GCK","KCNJ11","ABCC8","SIRT1","GLUT1",
|
| "SLC2A1"],
|
| "signaling", "Matschinsky-2013"),
|
| "Alpha_master_TF": _M(["ARX","IRX1","IRX2","MAFB","POU3F4","GCG","TTR"],
|
| "lineage", "Collombat-2003"),
|
| "Beta_master_TF_embryonic":_M(["NKX6-1","MNX1","NEUROD1","PDX1","NKX2-2",
|
| "HNF1B"],
|
| "lineage", "Gu-2004"),
|
| "Beta_master_TF_adult":_M(["MAFA","UCN3","SIX3","INS","IAPP","G6PC2"],
|
| "lineage", "Blum-2012"),
|
| "Neurog3_EP_cascade":_M(["NEUROG3","PAX4","FEV","INSM1","NEUROD1","SOX4",
|
| "CBFA2T3","BTBD17"],
|
| "lineage", "Gradwohl-2000"),
|
| "Endocrine_maturation":_M(["RFX3","RFX6","ISL1","FOXA2","PAX6","NKX2-2"],
|
| "lineage", "Piccand-2014"),
|
| "Exocrine_acinar": _M(["PRSS1","PRSS2","CEL","CPA1","CTRB1","AMY2A",
|
| "ELOVL5","PTF1A","CELA1"],
|
| "lineage", "Kawaguchi-2002"),
|
| "Ductal_epithelial": _M(["KRT19","KRT7","SOX9","MUC1","ONECUT1","HES1",
|
| "HNF1B","CFTR"],
|
| "lineage", "Solar-2009"),
|
| "Foregut_endoderm": _M(["SOX17","FOXA1","FOXA2","ONECUT1","PROX1","HNF1A",
|
| "HNF1B","GATA4","GATA6"],
|
| "lineage", "Zorn-2009"),
|
| "Cilium_Foxj1": _M(["FOXJ1","CFAP43","CFAP157","NPHP1","IFT88","DNAH5",
|
| "TEKT1","SPAG6"],
|
| "lineage", "Choksi-2014"),
|
| "Delta_master": _M(["SST","HHEX","LEPR","GHSR"],
|
| "hormone", "Rorsman-2018"),
|
| "Gamma_master": _M(["PPY","PYY","SLC38A4"],
|
| "hormone", "Wang-2016"),
|
| "Epsilon_ghrelin": _M(["GHRL","ACSL1"],
|
| "hormone", "Prado-2004"),
|
| "ER_stress_pancreas":_M(["ATF6","XBP1","ERN1","DDIT3","HSPA5","HSPA1A",
|
| "HSPA1B","EIF2AK3"],
|
| "stress", "Back-2012"),
|
| "Unfolded_protein_response":_M(["ATF4","ATF6","XBP1","HERPUD1","BAK1","BAX",
|
| "EDEM1","DERL1"],
|
| "stress", "Walter-2011"),
|
| "Hormone_processing":_M(["PCSK1","PCSK2","CPE","CHGA","CHGB","SCG2","SCG5",
|
| "PAM"],
|
| "hormone", "Docherty-1997"),
|
| "Insulin_receptor_signaling":_M(["INSR","IRS1","IRS2","AKT2","PDX1","FOXO1",
|
| "GSK3B","MTOR"],
|
| "signaling", "Kulkarni-1999"),
|
| "Mesenchyme_pancreatic":_M(["NKX3-2","BMP4","SOX9","FGF10","COL1A1","COL3A1",
|
| "DCN"],
|
| "lineage", "Landsman-2011"),
|
| "Fatty_acid_oxidation":_M(["CPT1A","ACADM","ACADL","HADHA","PPARA","PPARGC1A",
|
| "ACOX1"],
|
| "metabolism", "Houten-2010"),
|
| "Glycolysis": _M(["HK1","HK2","PFKM","PFKL","ALDOA","GAPDH","PGK1",
|
| "PKM","LDHA","ENO1","TPI1"],
|
| "metabolism", "Vander-Heiden-2009"),
|
| "TCA": _M(["CS","ACO2","IDH2","IDH3A","SDHA","FH","MDH2",
|
| "OGDH","SUCLA2"],
|
| "metabolism", "Chandel-2015"),
|
| "OXPHOS_ETC": _M(["NDUFA1","NDUFA2","NDUFB1","SDHA","SDHB","COX4I1",
|
| "COX5A","COX6A1","ATP5A1","ATP5B","UQCRQ"],
|
| "metabolism", "Mishra-2016"),
|
| "Redox_glutathione": _M(["GPX1","GPX2","GPX3","GPX4","GSR","PRDX1","PRDX2",
|
| "PRDX3","PRDX4","PRDX5","PRDX6","SOD1","SOD2","CAT"],
|
| "stress", "Ho-2007"),
|
| "Wnt_pancreas": _M(["WNT3A","WNT5A","CTNNB1","LEF1","TCF7","AXIN2",
|
| "FZD7"],
|
| "signaling", "Murtaugh-2008"),
|
| "Notch_pancreas": _M(["NOTCH1","NOTCH2","JAG1","HES1","DLL1","DLL4",
|
| "RBPJ","HEY1"],
|
| "signaling", "Apelqvist-1999"),
|
| "TGFb_pancreas": _M(["TGFB1","TGFB2","TGFBR1","TGFBR2","SMAD2","SMAD3",
|
| "SMAD4","SMAD7"],
|
| "signaling", "Sanvito-1994"),
|
| "Immune_pancreas": _M(["PTPRC","CD68","ADGRE1","CD3D","CD3E","CD4","CD8A",
|
| "CD19"],
|
| "immune", "Homo-2015"),
|
| "Endothelial_pancreas":_M(["PECAM1","CDH5","KDR","VWF","PLVAP","FLT1",
|
| "TEK","ENG"],
|
| "lineage", "Cleaver-2019"),
|
| "Cell_cycle": _M(["CCND1","CCNE1","CCNA2","CCNB1","CDK1","CDK2",
|
| "CDK4","MKI67","TOP2A","PCNA","MCM2","MCM3"],
|
| "cycle", "Whitfield-2002"),
|
| }
|
|
|
|
|
| def load_pan_skin():
|
| RAW = ROOT / "data/raw/GSE220977_combined.h5ad"
|
| PRED = ROOT / "discovery/pan_skin/marker/dingwall_predictions.csv"
|
|
|
| CKO_GSMS = {"GSM6833482", "GSM6833483"}
|
| WT_GSMS = {"GSM6833478", "GSM6833479", "GSM6833480", "GSM6833481"}
|
|
|
| a = ad.read_h5ad(RAW)
|
| pred = pd.read_csv(PRED)
|
| common = a.obs_names.intersection(pd.Index(pred["cell_id"].astype(str)))
|
| a = a[list(common)].copy()
|
| pred_map = dict(zip(pred["cell_id"].astype(str), pred["pred_label"]))
|
| a.obs["pred_label"] = pd.Categorical([pred_map.get(c,"unknown")
|
| for c in a.obs_names])
|
| samp = a.obs["sample"].astype(str)
|
| a.obs["group"] = np.where(samp.isin(list(CKO_GSMS)), "En1-cKO",
|
| np.where(samp.isin(list(WT_GSMS)), "WT", "other"))
|
| a = a[a.obs["group"].isin(["En1-cKO","WT"])].copy()
|
| return a, "En1-cKO", "WT", SKIN_MODULES
|
|
|
|
|
| def load_hematopoiesis():
|
| D_DIR = ROOT / "data/corpus/hematopoiesis/held_out_unlabeled/dahlin_extract"
|
| PRED = ROOT / "discovery/hematopoiesis/marker/nestorowa_anchor_predictions.csv"
|
| if not PRED.exists():
|
| PRED = ROOT / "discovery/hematopoiesis/marker/97_nestorowa_anchor_predictions.csv"
|
|
|
|
|
| DAHLIN_PRED_CANDIDATES = [
|
| ROOT / "discovery/hematopoiesis/marker/dahlin_predictions.csv",
|
| ROOT / "discovery/hematopoiesis/marker/92_dahlin_predictions.csv",
|
| ]
|
| for p in DAHLIN_PRED_CANDIDATES:
|
| if p.exists():
|
| PRED = p
|
| break
|
|
|
| GT = {"SIGAB1":"WT","SIGAC1":"WT","SIGAD1":"WT","SIGAF1":"WT","SIGAG1":"WT",
|
| "SIGAH1":"WT","SIGAG8":"Kit_W41","SIGAH8":"Kit_W41"}
|
|
|
| parts = []
|
| for f in sorted(D_DIR.glob("*.txt.gz")):
|
| sample = f.name.split("_")[1].split(".")[0]
|
| df = pd.read_csv(f, sep="\t", compression="gzip", index_col=0)
|
| X = sp.csr_matrix(df.values.T.astype(np.float32))
|
| obs = pd.DataFrame(index=[f"{sample}_{bc}" for bc in df.columns.astype(str)])
|
| obs["sample"] = sample
|
| obs["group"] = GT.get(sample, "unknown")
|
| var = pd.DataFrame(index=df.index.astype(str))
|
| parts.append(ad.AnnData(X=X, obs=obs, var=var))
|
| a = ad.concat(parts, join="outer", label="_batch")
|
|
|
| try:
|
| import mygene
|
| mg = mygene.MyGeneInfo()
|
| res = mg.querymany(a.var_names.astype(str).tolist(), scopes="ensembl.gene",
|
| fields="symbol", species="mouse", verbose=False)
|
| id2sym = {r["query"]: r["symbol"] for r in res if "symbol" in r}
|
| syms = pd.Series(a.var_names.astype(str)).map(id2sym).values
|
| keep = pd.notna(syms)
|
| a = a[:, keep].copy()
|
| a.var_names = syms[keep]
|
| a.var_names_make_unique()
|
| except Exception as e:
|
| print(f"[warn] mygene mapping failed: {e}")
|
|
|
| a = a[a.obs["group"].isin(["Kit_W41","WT"])].copy()
|
|
|
| if PRED.exists():
|
| pred = pd.read_csv(PRED)
|
| common = a.obs_names.intersection(pd.Index(pred["cell_id"].astype(str)))
|
| a = a[list(common)].copy()
|
| pred_map = dict(zip(pred["cell_id"].astype(str), pred["pred_label"]))
|
| a.obs["pred_label"] = pd.Categorical([pred_map.get(c,"unknown")
|
| for c in a.obs_names])
|
| else:
|
|
|
| a.obs["pred_label"] = pd.Categorical(["all"] * a.n_obs)
|
| print(f"[warn] no Dahlin prediction file found; using pred_label='all'")
|
|
|
| return a, "Kit_W41", "WT", HSC_MODULES
|
|
|
|
|
| def load_pancreas():
|
| SHARON_DIR = ROOT / "data/corpus/pancreas/held_out_unlabeled/sharon_extract"
|
| PRED = ROOT / "discovery/pancreas/marker/veres_predictions.csv"
|
|
|
| parts = []
|
| for meta_file in sorted(SHARON_DIR.glob("*.cell_metadata.tsv.gz")):
|
| counts_file = str(meta_file).replace("cell_metadata", "processed_counts")
|
| if not Path(counts_file).exists():
|
| continue
|
| meta = pd.read_csv(meta_file, sep="\t", compression="gzip")
|
| counts = pd.read_csv(counts_file, sep="\t", compression="gzip", index_col=0)
|
| obs = meta.set_index("library.barcode")
|
| obs = obs.loc[obs.index.intersection(counts.index)]
|
| counts_al = counts.loc[obs.index]
|
| X = sp.csr_matrix(counts_al.values.astype(np.float32))
|
| aa = ad.AnnData(X=X, obs=obs,
|
| var=pd.DataFrame(index=counts_al.columns))
|
| aa.var_names_make_unique()
|
| parts.append(aa)
|
| a = ad.concat(parts, join="outer")
|
|
|
| pred = pd.read_csv(PRED)
|
|
|
| pred["cell_id_stripped"] = pred["cell_id"].astype(str).str.replace(r"^veres_", "", regex=True)
|
| common = a.obs_names.intersection(pd.Index(pred["cell_id_stripped"]))
|
| a = a[list(common)].copy()
|
| pred_map = dict(zip(pred["cell_id_stripped"], pred["pred_label"]))
|
| a.obs["pred_label"] = pd.Categorical([pred_map.get(c,"unknown")
|
| for c in a.obs_names])
|
|
|
|
|
| stage = a.obs["Stage"].astype(str)
|
| a.obs["group"] = np.where(stage == "6", "Stage6",
|
| np.where(stage == "5", "Stage5", "other"))
|
| a = a[a.obs["group"].isin(["Stage6","Stage5"])].copy()
|
| return a, "Stage6", "Stage5", PANCREAS_MODULES
|
|
|
|
|
| LOADERS = {
|
| "pan_skin": load_pan_skin,
|
| "hematopoiesis": load_hematopoiesis,
|
| "pancreas": load_pancreas,
|
| }
|
|
|
|
|
| def score_modules(sub, modules):
|
| var_set = set(sub.var_names.astype(str))
|
| for name, spec in modules.items():
|
| present = [g for g in spec["genes"] if g in var_set]
|
| if not present:
|
| sub.obs[f"pw_{name}"] = 0.0
|
| continue
|
| try:
|
| sc.tl.score_genes(sub, gene_list=present, score_name=f"pw_{name}",
|
| random_state=0, use_raw=False)
|
| except Exception:
|
| sub.obs[f"pw_{name}"] = 0.0
|
| return sub
|
|
|
|
|
| def run_system(system_name):
|
| print(f"[load] {system_name}", flush=True)
|
| a, g1, g2, modules = LOADERS[system_name]()
|
| print(f"[load] {a.n_obs} cells, {sum(a.obs['group']==g1)} {g1}, "
|
| f"{sum(a.obs['group']==g2)} {g2}, {len(modules)} modules", flush=True)
|
|
|
| sc.pp.normalize_total(a, target_sum=1e4)
|
| sc.pp.log1p(a)
|
|
|
| OUT = ROOT / f"discovery/{system_name}/marker"
|
| OUT.mkdir(parents=True, exist_ok=True)
|
|
|
| classes = sorted(a.obs["pred_label"].astype(str).unique())
|
| rows = []
|
| for cls in classes:
|
| mask = (a.obs["pred_label"].astype(str) == cls).values
|
| n1 = int((mask & (a.obs["group"].values == g1)).sum())
|
| n2 = int((mask & (a.obs["group"].values == g2)).sum())
|
| if n1 < MIN_PER_GROUP or n2 < MIN_PER_GROUP:
|
| print(f"[pw] {cls}: skip (n_{g1}={n1}, n_{g2}={n2})")
|
| continue
|
| sub = a[mask].copy()
|
| sub = score_modules(sub, modules)
|
| grp = sub.obs["group"].values
|
| for mod_name, spec in modules.items():
|
| s = sub.obs[f"pw_{mod_name}"].astype(float).values
|
| v1 = s[grp == g1]; v2 = s[grp == g2]
|
| try:
|
| _, pval = mannwhitneyu(v1, v2, alternative="two-sided")
|
| except Exception:
|
| pval = 1.0
|
| delta = float(v1.mean() - v2.mean())
|
| rows.append({
|
| "class": cls,
|
| "module_name": mod_name,
|
| "module_type": spec["type"],
|
| "citation": spec["citation"],
|
| "direction": spec["direction"],
|
| "n_g1": n1,
|
| "n_g2": n2,
|
| "group_g1": g1,
|
| "group_g2": g2,
|
| "delta": round(delta, 4),
|
| "mannu_p": float(pval),
|
| })
|
| print(f"[pw] {cls}: {n1} {g1}, {n2} {g2} — scored")
|
|
|
| df = pd.DataFrame(rows)
|
| if df.empty:
|
| print("[pw] no eligible classes — done")
|
| return
|
|
|
| n_tests = len(df)
|
| df["mannu_p_adj_bonferroni"] = np.minimum(df["mannu_p"] * n_tests, 1.0)
|
|
|
| csv_path = OUT / "57_pathway_analysis.csv"
|
| df.to_csv(csv_path, index=False)
|
| print(f"[pw] wrote {csv_path} ({len(df)} rows, n_tests={n_tests})",
|
| flush=True)
|
|
|
| pivot_delta = df.pivot(index="module_name", columns="class", values="delta")
|
| pivot_padj = df.pivot(index="module_name", columns="class",
|
| values="mannu_p_adj_bonferroni")
|
| pivot_delta.to_csv(OUT / "57_pathway_class_by_module_delta.tsv", sep="\t")
|
| pivot_padj.to_csv(OUT / "57_pathway_class_by_module_padj.tsv", sep="\t")
|
| print(f"[pw] wrote heatmap TSVs to {OUT}")
|
|
|
| sig = df[df["mannu_p_adj_bonferroni"] < 0.01].sort_values(
|
| "mannu_p_adj_bonferroni")
|
| print(f"\n[pw] top Bonferroni-significant shifts (padj<0.01, "
|
| f"n={len(sig)}):")
|
| if len(sig):
|
| print(sig[["class","module_name","module_type","delta",
|
| "mannu_p_adj_bonferroni"]].head(30).to_string(index=False))
|
|
|
|
|
| def main():
|
| p = argparse.ArgumentParser()
|
| p.add_argument("--system", required=True,
|
| choices=list(LOADERS.keys()),
|
| help="pan_skin | hematopoiesis | pancreas")
|
| args = p.parse_args()
|
| run_system(args.system)
|
|
|
|
|
| if __name__ == "__main__":
|
| main()
|
|
|
