Reproducing PBMC analysis

[Sandeep-S-K]

Hello, really interesting work! Thanks for the tool.

I am trying to reproduce the plots from Figure 1 of the paper regarding the PBMC analysis using Spectra.

I annotated the data using the process from the notebooks provided.

I am using the following code to reproduce the analysis.

annotations =  #gene sets from supplementary table 1 of paper

adata  =  # annotated data

# gene sets for cell type allocation from the notebook Kartha_PBMC_annotation.ipynb

M_path_names = ['M_IL17A_response', 'M_angiogenic-effectors', 'M_mac_CSF1_response', 'M_gran_CSF2_response', 'mast_granule-exocytosis', 'DC_CD40L_response', 'DC_antigen-crosspresentation', 'DC_LPS_response', 'Mac_LPS_response', 'Mac_CSF1_response', 'Mac_IL4-IL13_response', 'Mac_IFNG_response', 'p-DC_CpG-TLR9_response']
BP_path_names = ['B_Breg_UP', 'B_effector-1_UP', 'B_CD40L_response', 'B_IL2_response', 'B_effector-2_UP', 'B_IgM-ligation_response']
global_path_names = ['all_pyrimidine_metabolism', 'all_ros_response', 'all_autophagic-cell-death', 'all_DNA-methylation', 'all_eicosanoid_metabolism', 'all_MHC-I-presentation', 'all_bile-acid_synthesis', 'all_MHC-II-presentation', 'all_autophagy-selective', 'all_heparan-sulfate_degradation', 'all_galactose_metabolism', 'all_transmembrane-transport-lysosome', 'all_steroid_metabolism', 'all_mitophagy', 'all_complement_production', 'all_mitophagy_regulation_positive', 'all_autophagy-peroxisome', 'all_DNA-demethylation', 'all_coagulation-factor_production', 'all_GABA-shunt', 'all_selenoamino-acid_metabolism', 'all_pyruvate_metabolism', 'all_creatinine_metabolism', 'all_peroxisome-component', 'all_glyoxylate-dicarboxylate_metabolism', 'all_DNA-repair', 'all_oxidative-phosphorylation', 'all_n-glycan_degradation', 'all_PI3K-AKT-mTOR_signaling', 'all_glutathione_metabolism', 'all_chondroitine-sulfate_degradation', 'all_o-glycan_synthesis', 'all_glycogenesis', 'all_pyroptosis', 'all_CoA_synthesis', 'all_retinol_metabolism', 'all_transmembrane-transport-golgi', 'all_transmembrane-transport-cellmembrane', 'all_protein-degradation-proteasome', 'all_circadian-rhythm', 'all_type-I-ifn-response', 'all_LYS_metabolism', 'all_multidrug-resistance', 'all_ROS-detoxification', 'all_phosphoinositide_signaling', 'all_IL6-JAK-STAT3_signaling', 'all_Beta-Ala_metabolism', 'all_thrombolysis-factor_production', 'all_TNF-via-NFkB_signaling', 'all_cyclic-nucleotide_metabolism', 'all_SASP', 'all_fatty-acid_synthesis', 'all_carnitine-shuttle', 'all_sphingolipid_metabolism', 'all_biotin_metabolism', 'all_posttranslation-modification', 'all_apoptosis', 'all_hedgehog_signaling', 'all_p53-signaling', 'all_TYR_metabolism', 'all_porphyrine-heme_metabolism', 'all_MET_metabolism', 'all_ARG-PRO_metabolism', 'all_lipophagy', 'all_TCA-cycle', 'all_fatty-acid-metabolism', 'all_CYS_metabolism', 'all_wnt-beta-catenin-signaling', 'all_glycerin-SER-THR_metabolism', 'all_autophagy_regulation_positive', 'all_pterin_synthesis', 'all_platelet-activation-factor_production', 'all_purine_metabolism', 'all_xenobiotics_metabolism', 'all_NOD-like-receptor_signaling', 'all_autophagy-chaperone-mediated', 'all_cholesterol_metabolism', 'all_macroautophagy', 'all_mTORC1_signaling', 'all_citric-acid-cycle', 'all_autophagy-of-mitochondria_regulation_positive', 'all_cholesterol-homeostasis', 'all_keratan-sulfate_degradation', 'all_iron-uptake-and-storage', 'all_GPI-anchor_synthesis', 'all_glycolysis', 'all_hyaluronan_metabolism', 'all_ALA-ASP_metabolism', 'all_purine_synthesis', 'all_G2M-transition', 'all_osmotic-stress-response', 'all_ethanol_metabolism', 'all_keratan-sulfate_synthesis', 'all_HIS_metabolism', 'all_ketone-body_metabolism', 'all_glycogenolysis', 'all_transmembrane-transport-mitochondrial', 'all_RIG-I-like-receptor_signaling', 'all_propanoate_metabolism', 'all_pyrimidine_synthesis', 'all_riboflavin_metabolism', 'all_type-II-ifn-response', 'all_actin-cytoskeleton_regulation', 'all_nucleotide_metabolism', 'all_VAL-LEU-ILE_metabolism', 'all_microautophagy-lysosomal', 'all_hypoxia-response', 'all_urea-cycle', 'all_n-glycan_synthesis', 'all_CYP_metabolism', 'all_fatty-acid-beta-oxidation-mitochondrial', 'all_JAK-STAT_signaling', 'all_reticulophagy', 'all_amino-sugar-nucleotide-sugar_metabolism', 'all_pentose-phosphate-pathway', 'all_lactate_production', 'all_autophagy-nucleus', 'all_MYC_targets', 'all_NOTCH_signaling', 'all_ascorbate-uptake', 'all_autophagy-of-mitochondria', 'all_fatty-acid-beta-oxidation-peroxisomal', 'all_TLR_signaling', 'all_chondroitine-and-heparan-sulfate_synthesis', 'all_inositol-phosphate_metabolism', 'all_ubiquinone_synthesis', 'all_PHE_metabolism', 'all_unfolded-protein-response', 'all_polyamines_metabolism', 'all_transmembrane-transport-ER', 'all_taurine-hypotaurine_metabolism', 'all_exocytosis', 'all_histone-methylation', 'all_fructose-mannose_metabolism', 'all_cytosolic-DNA-sensing_signaling', 'all_triacylglycerol_synthesis', 'all_mitotic-spindle-component', 'all_GLU_metabolism', 'all_glycerophospholipid_metabolism', 'all_thiamin_metabolism', 'all_TRP_metabolism', 'all_macroautophagy_regulation_positive', 'all_G1S-transition', 'all_DNA_synthesis', 'all_folate_metabolism', 'all_nucleophagy-late', 'all_NAD_metabolism', 'leuko_ROS_production', 'leuko_transendothelial-migration']
TNK_path_names = ['TNK_cytotoxicity-effectors', 'TNK_IL2-STAT5-signaling', 'TNK_IL2_response', 'TNK_PD-1_signaling', 'T_IL4_response', 'T_IL21_response', 'T_tcr-activation', 'NK_IL15_response', 'CD4-T_TH2_UP', 'CD4-T_TH17_UP', 'CD4-T_TH1_UP', 'CD4-T_TFH_UP', 'CD4-T_IL4_response', 'CD4-T_TH22_UP', 'CD4-T_IL12_response', 'CD4-T_TH9_UP', 'CD8-T_tumor-reactive-like_UP', 'CD8-T_terminal-exhaustion', 'CD8-T_IL12_response', 'Treg_FoxP3-stabilization']


annotations_per_cell_type ={'M':{},'BP':{},'TNK':{},'global':{}}

print('Gene sets excluded (not found in any of the types)')
for path in  annotations.keys():
    if path in M_path_names:
        annotations_per_cell_type['M'][path] = annotations[path]
    elif path in BP_path_names:
        annotations_per_cell_type['BP'][path] = annotations[path]
    elif path in TNK_path_names:
        annotations_per_cell_type['TNK'][path] = annotations[path]
    elif path in global_path_names:
        annotations_per_cell_type['global'][path] = annotations[path]
    else:
        print(path)

#λ of 0.01, δ of 0.001 and ρ of 0.001. - From the paper



model = spc.est_spectra(adata=adata, 
    gene_set_dictionary=annotations_per_cell_type , 
    use_highly_variable=True,
    cell_type_key= 'Cell_type', 
    use_cell_types = True,
    use_weights=True,
    lam=0.01, 
    delta=0.001, 
    kappa=None,
    rho=0.001, 
    n_top_vals=50,
    label_factors=True, 
    overlap_threshold=0.2,
    clean_gs = True, 
    min_gs_num = 3,
    verbose = True,
    num_epochs=10000 
)


len(adata.var[adata.var['spectra_vocab']]) # = 11,884  #which is different from the 11,840 mentioned in the paper 

The gene sets I am using are from the supplementary table 1. But it does not agree completely from the gene sets used in the notebooks, some seem to be missing. I'm not sure if I'm missing any steps.

It would be great if you could help me with this. Thanks!