Tumor region (core vs border) composition analysis with scCODA¶

Purpose¶

This notebook runs a Bayesian compositional differential abundance analysis (scCODA) to test whether cell-type composition differs between tumor border vs core, and whether those region effects differ by condition/day.

Concretely, it fits a scCODA model with an interaction:

• Condition_day * roi with:
• baseline roi = core
• baseline Condition_day = IL-15#4 (for the day 4 run) or IL-15#8 (for the day 8 run)

The main output is an effect table per cell type and condition level reporting:

• posterior mean/median of the interaction term (border vs core effect relative to baseline)
• HDI intervals
• inclusion probability–derived local FDR and q-values

Inputs¶

• data/sccoda_tumor_regions_input.csv

Expected columns in the input CSV:

• batch : sample identifier
• roi : tumor region label (must map cleanly to "core" or "border")
• cluster : cell type / cluster label
• Condition : condition label
• Harvest_Day : numeric day
• Condition_day : combined condition+day label used in the model (e.g., IL-15#4, Both#8, etc.)
• roi_area : (optional) used only if you extend/enable area normalization utilities

Outputs¶

This notebook currently creates pandas DataFrames in memory (does not save by default):

• df_sum_4, df_sum_8 : scCODA summary tables (effects for border vs core interaction)
• df_raw_4, df_raw_8 : raw per-(batch × celltype) counts with border/core totals and log-odds ratio estimate
• df_summary : concatenated summary across baselines (IL-15#4 and IL-15#8)

Recommended (optional) export step for reproducibility:

• Save df_summary and df_raw_* to results/tumor_regions_sccoda/ as CSV (or Excel) so downstream notebooks/figures don’t require rerunning HMC.

Dependencies¶

• sccoda (CompositionalAnalysis + HMC sampling)
• arviz (HDI computation)
• pandas, numpy
• matplotlib (only needed if you add plots)

Run order¶

Run top-to-bottom:

1. Imports + random seed
2. Load sccoda_tumor_regions_input.csv
3. Define scooda_region_composition helper class
4. Fit scCODA twice (baseline = IL-15#4 and IL-15#8)
5. Combine outputs into df_summary and (recommended) export

Notes / common pitfalls¶

• ROI labels: roi is normalized to lowercase and mapped to "core"/"border". Any unexpected ROI strings will break the “two ROI levels” assumption.
• Reference cell type: scCODA uses reference_cell_type='Tumor'. The cluster column must contain a "Tumor" level, or the model will fail.
• Runtime: HMC sampling can be slow (defaults: 2000 draws, 1000 tune, 2 chains). Consider lowering for testing, then restoring for final runs.
• scCODA version differences: the code includes logic to handle different sample_hmc argument names across versions.

In [ ]:

# mamba activate env_scoda 
import scanpy as sc
import numpy as np
import os
import matplotlib.pyplot as plt

import pandas as pd

import alphashape as alphashape

# mamba activate env_scoda import scanpy as sc import numpy as np import os import matplotlib.pyplot as plt import pandas as pd import alphashape as alphashape

In [ ]:

import arviz as az

try:
    from sccoda.util import comp_ana as scd
except ImportError as e:
    raise ImportError(
        "Missing dependency: sccoda. Install with: pip install sccoda"
    ) from e

# scCODA HMC backend deps (often required)
try:
    import jax
    import numpyro
except ImportError as e:
    raise ImportError(
        "Missing dependency: jax/numpyro. Install with: pip install 'jax[cpu]' numpyro"
    ) from e

import arviz as az try: from sccoda.util import comp_ana as scd except ImportError as e: raise ImportError( "Missing dependency: sccoda. Install with: pip install sccoda" ) from e # scCODA HMC backend deps (often required) try: import jax import numpyro except ImportError as e: raise ImportError( "Missing dependency: jax/numpyro. Install with: pip install 'jax[cpu]' numpyro" ) from e

Load input table for tumor region analysis¶

Set path_tables to the repository data/ directory (or wherever you placed the input CSV). This table should contain one row per cell with columns including batch, roi, cluster, and Condition_day.

In [ ]:

# Large input files for this notebook are hosted on Zenodo; download and unpack the Zenodo bundle locally and set the base data path below to your local folder.
DATA_ROOT = "data_external/nguyenlab-il15act-dataset-v1"
path_tables = os.path.join(DATA_ROOT, "02_tables_inputs")
df_in = pd.read_csv(os.path.join(path_tables, "sccoda_tumor_regions_input.csv"))

# Large input files for this notebook are hosted on Zenodo; download and unpack the Zenodo bundle locally and set the base data path below to your local folder. DATA_ROOT = "data_external/nguyenlab-il15act-dataset-v1" path_tables = os.path.join(DATA_ROOT, "02_tables_inputs") df_in = pd.read_csv(os.path.join(path_tables, "sccoda_tumor_regions_input.csv"))

scCODA wrapper class¶

scooda_region_composition:

• normalizes ROI labels
• constructs a sample-by-celltype count matrix (pd.crosstab)
• fits scCODA with a Condition_day * roi interaction
• returns (a) an effect summary table and (b) raw per-sample count summaries

In [21]:

class scooda_region_composition:
    ROI_BASE  = "core"        # compare border vs core (core is baseline)
    ROI_OTHER = "border"
    # --- choose reference cell type for scCODA log-ratio (must exist in cluster) ---
    REF_CELLTYPE = "Tumor"
    # --- HDI probability ---
    HDI_PROB = 0.94
    # --- HMC settings ---
    hmc_draws = 2000
    hmc_tune = 1000
    chains = 2
    target_accept = 0.9
    
    def __init__(self, df_input,BASE_COND = "IL-15#4"):
        self.df_input = df_input
        self.BASE_COND =BASE_COND    # compare everything relative to this
        self.ROI_BASE  = "core"        # compare border vs core (core is baseline)
        self.ROI_OTHER = "border"
        # --- choose reference cell type for scCODA log-ratio (must exist in cluster) ---
        self.REF_CELLTYPE = "Tumor"
        # --- HDI probability ---
        self.HDI_PROB = 0.94
        # --- HMC settings ---
        self.hmc_draws = 2000
        self.hmc_tune = 1000
        self.chains = 2
        self.target_accept = 0.9
        self.data,self.res = self.run_sccooda()
    
    @staticmethod
    def area_logOddsRatio(df_arr):
        df_border = df_arr.query('roi=="border"')
        df_core = df_arr.query('roi=="core"')
        frac =  df_border['roi_area'].iloc[0]/df_core['roi_area'].iloc[0]
        return np.log(frac) 

    def logit(frac):
        return np.log(frac/(1-frac))
    
    @staticmethod
    def compute_param(df_curr,df_immun,ct):
        def logit(frac):
            return np.log(frac/(1-frac))
        if df_curr.empty or (len(df_curr.roi.unique())!=2):
            cell = ct
            df_ct = df_immun.query(f'cluster==@cell')
            total = int(np.mean(df_ct.groupby(['batch']).size().values))
            border = int(np.mean(df_ct.query('roi=="border"').groupby(['batch']).size().values))
            core = int(np.mean(df_ct.query('roi=="core"').groupby(['batch']).size().values))
        else:
            total = len(df_curr)
            border = len(df_curr.query('roi == "border"'))
            core = len(df_curr.query('roi == "core"'))
        frac_border = border / total
        if  frac_border==0.0:
            frac_border = 1e-5
        logit_border = logit(frac_border)
        frac_core = core / total
        if  frac_core==0.0:
            frac_core = 1e-5
        logit_core = logit(frac_core)
        loddr_roi = logit_border - logit_core
        return total, border, core, loddr_roi
    
    
    
    @staticmethod
    def add_fdr_flag(df_effects, res_obj, fdr=0.05):
        """
        Merge scCODA effect_df (after set_fdr) into df_effects using
        (celltype, covariate). Returns df with final_parameter + credible flag.
        """
        res_obj.set_fdr(est_fdr=fdr)

        eff = res_obj.effect_df.copy().reset_index()

        # Normalize expected column names across versions
        rename_map = {}
        for col in eff.columns:
            if col.lower() == "covariate":
                rename_map[col] = "covariate"
            if col.lower() in ["cell type", "cell_type", "celltype"]:
                rename_map[col] = "celltype"
            if col.lower() in ["final parameter", "final_parameter"]:
                rename_map[col] = "final_parameter"

        eff = eff.rename(columns=rename_map)

        if not set(["covariate", "celltype", "final_parameter"]).issubset(set(eff.columns)):
            raise ValueError(
                "Could not find required columns in res.effect_df.\n"
                f"Columns present: {list(eff.columns)}"
            )

        out = df_effects.merge(
            eff[["covariate", "celltype", "final_parameter"]],
            on=["covariate", "celltype"],
            how="left"
        )

        out[f"credible_FDR{int(fdr*100):02d}"] = out["final_parameter"].fillna(0) != 0
        return out
    
    def run_sccooda(self):
        df = self.df_input.copy()

        # Normalize ROI labels to exactly "core" / "border"
        df["roi"] = (
            df["roi"]
            .astype(str)
            .str.strip()
            .str.lower()
        )

        # (optional) map common variants to core/border if you have them
        roi_map = {
            "Core": "core",
            "Border": "border",
            "CORE": "core",
            "BORDER": "border",
            "core ": "core",
            "border ": "border",
        }
        df["roi"] = df["roi"].replace(roi_map)

        # Make sample ID
        df["sample"] = df["batch"].astype(str) + "_" + df["roi"].astype(str)

        # Covariates (1 row per sample)
        cov = (
            df[["sample", "Condition_day", "roi"]]
            .drop_duplicates(subset=["sample"])
            .set_index("sample")
        )

        # Count table: sample x cluster
        Y = pd.crosstab(df["sample"], df["cluster"].astype(str))

        # Merge covariates + counts
        df_model = cov.join(Y, how="left").fillna(0)

        count_cols = Y.columns.tolist()
        df_model[count_cols] = df_model[count_cols].apply(pd.to_numeric, errors="coerce").fillna(0).astype(int)

        # Drop empty samples (safety)
        df_model = df_model.loc[df_model[count_cols].sum(axis=1) > 0].copy()


        # ------------------------------------------------------------
        # 2) Build scCODA data object
        # ------------------------------------------------------------
        data = ccd.from_pandas(df_model, covariate_columns=["Condition_day", "roi"])

        # ------------------------------------------------------------
        # 3) Fit scCODA model:
        # Condition_day * roi
        # with baseline = BASE_COND and roi baseline = ROI_BASE
        # ------------------------------------------------------------
        formula = (
            f"C(Condition_day, Treatment('{self.BASE_COND}')) * "
            f"C(roi, Treatment('{self.ROI_BASE}'))"
        )

        model = ca.CompositionalAnalysis(
            data,
            formula=formula,
            reference_cell_type='Tumor'
        )

        # --- HMC settings ---
        hmc_draws = 2000
        hmc_tune = 1000
        chains = 2
        target_accept = 0.9

        # Handle scCODA version differences in sample_hmc signature
        sig = inspect.signature(model.sample_hmc)
        kwargs = {}

        if "num_results" in sig.parameters:
            kwargs["num_results"] = self.hmc_draws
        if "num_burnin" in sig.parameters:
            kwargs["num_burnin"] = self.hmc_tune
        if "num_warmup" in sig.parameters:
            kwargs["num_warmup"] = self.hmc_tune
        if "target_accept_prob" in sig.parameters:
            kwargs["target_accept_prob"] = self.target_accept
        if "target_accept" in sig.parameters:
            kwargs["target_accept"] = self.target_accept
        if "chains" in sig.parameters:
            kwargs["chains"] = self.chains
        elif "n_chains" in sig.parameters:
            kwargs["n_chains"] = self.chains

        res = model.sample_hmc(**kwargs)
        return data,res
    
    def table_generate(self):
        res = self.res
        data = self.data
        beta = res.posterior["beta"]
        cov_names = [str(c) for c in beta.coords["covariate"].values]
        cell_types = list(beta.coords["cell_type"].values)

        def interaction_cov_name(level, base=self.BASE_COND, roi_ref=self.ROI_BASE, roi_other=self.ROI_OTHER):
            return (
                f"C(Condition_day, Treatment('{base}'))[T.{level}]:"
                f"C(roi, Treatment('{roi_ref}'))[T.{roi_other}]"
            )

        # pick which condition levels to report (exclude the baseline itself)
        all_levels = sorted(data.obs["Condition_day"].unique())
        levels = [x for x in all_levels if x != self.BASE_COND]

        rows = []
        for ct in cell_types:
            for level in levels:
                term = interaction_cov_name(level)

                # some terms might be missing if data is sparse; skip safely
                if term not in cov_names:
                    continue

                draws = beta.sel(cell_type=ct, covariate=term).values.ravel()
                hdi = az.hdi(draws, hdi_prob=self.HDI_PROB)

                rows.append({
                    "celltype": str(ct),
                    "level": str(level),
                    "contrast": f"{level} minus {self.BASE_COND}",
                    "covariate": term,
                    "delta_mean": float(draws.mean()),
                    "delta_median": float(np.median(draws)),
                    "hdi_lo": float(hdi[0]),
                    "hdi_hi": float(hdi[1]),
                    "P_gt0": float((draws > 0).mean()),
                    "P_lt0": float((draws < 0).mean()),
                })

        df_interaction = pd.DataFrame(rows).sort_values(["celltype", "level"]).reset_index(drop=True)
        df_final = self.add_fdr_flag(df_interaction, res, fdr=0.05)
        eff = res.effect_df.copy().reset_index()
        
        rename_map = {}

        for col in eff.columns:
            cl = col.lower()
            if cl == "covariate":
                rename_map[col] = "covariate"
            elif cl in ["cell type", "cell_type", "celltype"]:
                rename_map[col] = "celltype"
            elif ("inclusion" in cl and "prob" in cl) or ("inclusion probability" in cl):
                rename_map[col] = "inclusion_prob"
            elif "posterior" in cl and "inclusion" in cl:
                rename_map[col] = "inclusion_prob"

        eff = eff.rename(columns=rename_map)

        if "inclusion_prob" not in eff.columns:
            raise ValueError(
                "Could not find inclusion probability in res.effect_df.\n"
                f"Columns present: {list(eff.columns)}\n"
                "Try printing res.effect_df.head() and we will map the exact column name."
            )

        # Merge inclusion probability into df_final
        df_final = df_final.merge(
            eff[["covariate", "celltype", "inclusion_prob"]],
            on=["covariate", "celltype"],
            how="left"
        )

        # local FDR proxy (per-effect): 1 - inclusion_prob
        df_final["local_fdr"] = 1.0 - df_final["inclusion_prob"]
        # Bayesian q-value style number:
        # sort by strongest evidence first (largest inclusion_prob)
        tmp = df_final.sort_values("inclusion_prob", ascending=False).copy()
        tmp["fdr_at_threshold"] = tmp["local_fdr"].expanding().mean()
        tmp["q_value"] = tmp["fdr_at_threshold"][::-1].cummin()[::-1]

        # merge back
        df_final = df_final.merge(
            tmp[["covariate", "celltype", "q_value", "fdr_at_threshold"]],
            on=["covariate", "celltype"],
            how="left"
        )
        df_final = df_final[['celltype', 'level', 'delta_mean', 'hdi_lo', 'hdi_hi', 'local_fdr', 'q_value']]
        return df_final
    
    def raw_table(self):
        df_immun = self.df_input.copy()
        table_key = df_immun[['Condition','batch','Harvest_Day']]
        table_key.drop_duplicates(['batch'],inplace=True)
        b2cond = {b:cond for b,cond in zip(table_key['batch'],table_key['Condition'])}
        b2day = {b:day for b,day in zip(table_key['batch'],table_key['Harvest_Day'])} 
        clusters = df_immun['cluster'].unique()
        container = []  
        for batch in df_immun.batch.astype(str).unique() :
            for cells in clusters:
                df_curr = df_immun.query('batch==@batch and cluster==@cells')
                total, border, core, loddr_roi = self.compute_param(df_curr,df_immun,cells)
                container.append(pd.DataFrame({'cell':[cells],
                                        'batch':[batch],
                                        'total':[total],
                                        'border':[border],
                                        'core':[core],
                                        # 'logodds_area':[logOddsArea],
                                        'loddr_roi':[loddr_roi],}))
        df_values = pd.concat(container).reset_index()
        df_values['harvest_day'] = df_values['batch'].map(b2day)
        df_values['condition'] = df_values['batch'].map(b2cond)
        return df_values

class scooda_region_composition: ROI_BASE = "core" # compare border vs core (core is baseline) ROI_OTHER = "border" # --- choose reference cell type for scCODA log-ratio (must exist in cluster) --- REF_CELLTYPE = "Tumor" # --- HDI probability --- HDI_PROB = 0.94 # --- HMC settings --- hmc_draws = 2000 hmc_tune = 1000 chains = 2 target_accept = 0.9 def __init__(self, df_input,BASE_COND = "IL-15#4"): self.df_input = df_input self.BASE_COND =BASE_COND # compare everything relative to this self.ROI_BASE = "core" # compare border vs core (core is baseline) self.ROI_OTHER = "border" # --- choose reference cell type for scCODA log-ratio (must exist in cluster) --- self.REF_CELLTYPE = "Tumor" # --- HDI probability --- self.HDI_PROB = 0.94 # --- HMC settings --- self.hmc_draws = 2000 self.hmc_tune = 1000 self.chains = 2 self.target_accept = 0.9 self.data,self.res = self.run_sccooda() @staticmethod def area_logOddsRatio(df_arr): df_border = df_arr.query('roi=="border"') df_core = df_arr.query('roi=="core"') frac = df_border['roi_area'].iloc[0]/df_core['roi_area'].iloc[0] return np.log(frac) def logit(frac): return np.log(frac/(1-frac)) @staticmethod def compute_param(df_curr,df_immun,ct): def logit(frac): return np.log(frac/(1-frac)) if df_curr.empty or (len(df_curr.roi.unique())!=2): cell = ct df_ct = df_immun.query(f'cluster==@cell') total = int(np.mean(df_ct.groupby(['batch']).size().values)) border = int(np.mean(df_ct.query('roi=="border"').groupby(['batch']).size().values)) core = int(np.mean(df_ct.query('roi=="core"').groupby(['batch']).size().values)) else: total = len(df_curr) border = len(df_curr.query('roi == "border"')) core = len(df_curr.query('roi == "core"')) frac_border = border / total if frac_border==0.0: frac_border = 1e-5 logit_border = logit(frac_border) frac_core = core / total if frac_core==0.0: frac_core = 1e-5 logit_core = logit(frac_core) loddr_roi = logit_border - logit_core return total, border, core, loddr_roi @staticmethod def add_fdr_flag(df_effects, res_obj, fdr=0.05): """ Merge scCODA effect_df (after set_fdr) into df_effects using (celltype, covariate). Returns df with final_parameter + credible flag. """ res_obj.set_fdr(est_fdr=fdr) eff = res_obj.effect_df.copy().reset_index() # Normalize expected column names across versions rename_map = {} for col in eff.columns: if col.lower() == "covariate": rename_map[col] = "covariate" if col.lower() in ["cell type", "cell_type", "celltype"]: rename_map[col] = "celltype" if col.lower() in ["final parameter", "final_parameter"]: rename_map[col] = "final_parameter" eff = eff.rename(columns=rename_map) if not set(["covariate", "celltype", "final_parameter"]).issubset(set(eff.columns)): raise ValueError( "Could not find required columns in res.effect_df.\n" f"Columns present: {list(eff.columns)}" ) out = df_effects.merge( eff[["covariate", "celltype", "final_parameter"]], on=["covariate", "celltype"], how="left" ) out[f"credible_FDR{int(fdr*100):02d}"] = out["final_parameter"].fillna(0) != 0 return out def run_sccooda(self): df = self.df_input.copy() # Normalize ROI labels to exactly "core" / "border" df["roi"] = ( df["roi"] .astype(str) .str.strip() .str.lower() ) # (optional) map common variants to core/border if you have them roi_map = { "Core": "core", "Border": "border", "CORE": "core", "BORDER": "border", "core ": "core", "border ": "border", } df["roi"] = df["roi"].replace(roi_map) # Make sample ID df["sample"] = df["batch"].astype(str) + "_" + df["roi"].astype(str) # Covariates (1 row per sample) cov = ( df[["sample", "Condition_day", "roi"]] .drop_duplicates(subset=["sample"]) .set_index("sample") ) # Count table: sample x cluster Y = pd.crosstab(df["sample"], df["cluster"].astype(str)) # Merge covariates + counts df_model = cov.join(Y, how="left").fillna(0) count_cols = Y.columns.tolist() df_model[count_cols] = df_model[count_cols].apply(pd.to_numeric, errors="coerce").fillna(0).astype(int) # Drop empty samples (safety) df_model = df_model.loc[df_model[count_cols].sum(axis=1) > 0].copy() # ------------------------------------------------------------ # 2) Build scCODA data object # ------------------------------------------------------------ data = ccd.from_pandas(df_model, covariate_columns=["Condition_day", "roi"]) # ------------------------------------------------------------ # 3) Fit scCODA model: # Condition_day * roi # with baseline = BASE_COND and roi baseline = ROI_BASE # ------------------------------------------------------------ formula = ( f"C(Condition_day, Treatment('{self.BASE_COND}')) * " f"C(roi, Treatment('{self.ROI_BASE}'))" ) model = ca.CompositionalAnalysis( data, formula=formula, reference_cell_type='Tumor' ) # --- HMC settings --- hmc_draws = 2000 hmc_tune = 1000 chains = 2 target_accept = 0.9 # Handle scCODA version differences in sample_hmc signature sig = inspect.signature(model.sample_hmc) kwargs = {} if "num_results" in sig.parameters: kwargs["num_results"] = self.hmc_draws if "num_burnin" in sig.parameters: kwargs["num_burnin"] = self.hmc_tune if "num_warmup" in sig.parameters: kwargs["num_warmup"] = self.hmc_tune if "target_accept_prob" in sig.parameters: kwargs["target_accept_prob"] = self.target_accept if "target_accept" in sig.parameters: kwargs["target_accept"] = self.target_accept if "chains" in sig.parameters: kwargs["chains"] = self.chains elif "n_chains" in sig.parameters: kwargs["n_chains"] = self.chains res = model.sample_hmc(**kwargs) return data,res def table_generate(self): res = self.res data = self.data beta = res.posterior["beta"] cov_names = [str(c) for c in beta.coords["covariate"].values] cell_types = list(beta.coords["cell_type"].values) def interaction_cov_name(level, base=self.BASE_COND, roi_ref=self.ROI_BASE, roi_other=self.ROI_OTHER): return ( f"C(Condition_day, Treatment('{base}'))[T.{level}]:" f"C(roi, Treatment('{roi_ref}'))[T.{roi_other}]" ) # pick which condition levels to report (exclude the baseline itself) all_levels = sorted(data.obs["Condition_day"].unique()) levels = [x for x in all_levels if x != self.BASE_COND] rows = [] for ct in cell_types: for level in levels: term = interaction_cov_name(level) # some terms might be missing if data is sparse; skip safely if term not in cov_names: continue draws = beta.sel(cell_type=ct, covariate=term).values.ravel() hdi = az.hdi(draws, hdi_prob=self.HDI_PROB) rows.append({ "celltype": str(ct), "level": str(level), "contrast": f"{level} minus {self.BASE_COND}", "covariate": term, "delta_mean": float(draws.mean()), "delta_median": float(np.median(draws)), "hdi_lo": float(hdi[0]), "hdi_hi": float(hdi[1]), "P_gt0": float((draws > 0).mean()), "P_lt0": float((draws < 0).mean()), }) df_interaction = pd.DataFrame(rows).sort_values(["celltype", "level"]).reset_index(drop=True) df_final = self.add_fdr_flag(df_interaction, res, fdr=0.05) eff = res.effect_df.copy().reset_index() rename_map = {} for col in eff.columns: cl = col.lower() if cl == "covariate": rename_map[col] = "covariate" elif cl in ["cell type", "cell_type", "celltype"]: rename_map[col] = "celltype" elif ("inclusion" in cl and "prob" in cl) or ("inclusion probability" in cl): rename_map[col] = "inclusion_prob" elif "posterior" in cl and "inclusion" in cl: rename_map[col] = "inclusion_prob" eff = eff.rename(columns=rename_map) if "inclusion_prob" not in eff.columns: raise ValueError( "Could not find inclusion probability in res.effect_df.\n" f"Columns present: {list(eff.columns)}\n" "Try printing res.effect_df.head() and we will map the exact column name." ) # Merge inclusion probability into df_final df_final = df_final.merge( eff[["covariate", "celltype", "inclusion_prob"]], on=["covariate", "celltype"], how="left" ) # local FDR proxy (per-effect): 1 - inclusion_prob df_final["local_fdr"] = 1.0 - df_final["inclusion_prob"] # Bayesian q-value style number: # sort by strongest evidence first (largest inclusion_prob) tmp = df_final.sort_values("inclusion_prob", ascending=False).copy() tmp["fdr_at_threshold"] = tmp["local_fdr"].expanding().mean() tmp["q_value"] = tmp["fdr_at_threshold"][::-1].cummin()[::-1] # merge back df_final = df_final.merge( tmp[["covariate", "celltype", "q_value", "fdr_at_threshold"]], on=["covariate", "celltype"], how="left" ) df_final = df_final[['celltype', 'level', 'delta_mean', 'hdi_lo', 'hdi_hi', 'local_fdr', 'q_value']] return df_final def raw_table(self): df_immun = self.df_input.copy() table_key = df_immun[['Condition','batch','Harvest_Day']] table_key.drop_duplicates(['batch'],inplace=True) b2cond = {b:cond for b,cond in zip(table_key['batch'],table_key['Condition'])} b2day = {b:day for b,day in zip(table_key['batch'],table_key['Harvest_Day'])} clusters = df_immun['cluster'].unique() container = [] for batch in df_immun.batch.astype(str).unique() : for cells in clusters: df_curr = df_immun.query('batch==@batch and cluster==@cells') total, border, core, loddr_roi = self.compute_param(df_curr,df_immun,cells) container.append(pd.DataFrame({'cell':[cells], 'batch':[batch], 'total':[total], 'border':[border], 'core':[core], # 'logodds_area':[logOddsArea], 'loddr_roi':[loddr_roi],})) df_values = pd.concat(container).reset_index() df_values['harvest_day'] = df_values['batch'].map(b2day) df_values['condition'] = df_values['batch'].map(b2cond) return df_values

Fit scCODA models (day 4 baseline and day 8 baseline)¶

We run the model twice with different baselines so that contrasts at day 4 and day 8 can be interpreted cleanly:

• baseline IL-15#4 for day 4 comparisons
• baseline IL-15#8 for day 8 comparisons

In [22]:

scooda_4_obj = scooda_region_composition(df_input = df_in,BASE_COND = "IL-15#4")
df_sum_4 = scooda_4_obj.table_generate()
df_raw_4 = scooda_4_obj.raw_table()
scooda_8_obj = scooda_region_composition(df_input = df_in,BASE_COND = "IL-15#8")
df_sum_8 = scooda_8_obj .table_generate()
df_raw_8 = scooda_8_obj.raw_table()

scooda_4_obj = scooda_region_composition(df_input = df_in,BASE_COND = "IL-15#4") df_sum_4 = scooda_4_obj.table_generate() df_raw_4 = scooda_4_obj.raw_table() scooda_8_obj = scooda_region_composition(df_input = df_in,BASE_COND = "IL-15#8") df_sum_8 = scooda_8_obj .table_generate() df_raw_8 = scooda_8_obj.raw_table()

Zero counts encountered in data! Added a pseudocount of 0.5.

100%|██████████| 2000/2000 [00:17<00:00, 115.26it/s]

MCMC sampling finished. (26.639 sec)
Acceptance rate: 48.8%
Zero counts encountered in data! Added a pseudocount of 0.5.

100%|██████████| 2000/2000 [00:17<00:00, 113.98it/s]

MCMC sampling finished. (28.330 sec)
Acceptance rate: 60.2%

In [23]:

df_sum_4['baseline'] = 'IL-15#4'
df_sum_8['baseline'] = 'IL-15#8'
df_summary = pd.concat([df_sum_4, df_sum_8], axis=0)

df_sum_4['baseline'] = 'IL-15#4' df_sum_8['baseline'] = 'IL-15#8' df_summary = pd.concat([df_sum_4, df_sum_8], axis=0)

Combine and (recommended) export results¶

df_summary concatenates the scCODA interaction effects across the two baseline runs. For reproducibility, consider saving this table to disk so later steps don’t need to rerun HMC.

In [24]:

df_summary

df_summary

Out[24]:

	celltype	level	delta_mean	hdi_lo	hdi_hi	local_fdr	q_value	baseline
0	CAFs	Both#4	-0.210970	-9.777551e-01	8.945419e-02	0.126	0.050917	IL-15#4
1	CAFs	Both#8	-0.012006	-6.292036e-02	3.738368e-04	0.816	0.421984	IL-15#4
2	CAFs	IL-15#8	-0.022049	-2.031053e-01	3.349162e-01	0.112	0.044091	IL-15#4
3	CAFs	IL-21#4	-0.078417	-5.374433e-01	-3.897494e-52	0.820	0.440578	IL-15#4
4	CAFs	IL-21#8	0.076208	2.479009e-34	3.416394e-01	0.461	0.216647	IL-15#4
...	...	...	...	...	...	...	...	...
79	Tumor	IL-15#4	0.000000	0.000000e+00	0.000000e+00	1.000	0.572915	IL-15#8
80	Tumor	IL-21#4	0.000000	0.000000e+00	0.000000e+00	1.000	0.533053	IL-15#8
81	Tumor	IL-21#8	0.000000	0.000000e+00	0.000000e+00	1.000	0.539197	IL-15#8
82	Tumor	TCR#4	0.000000	0.000000e+00	0.000000e+00	1.000	0.545182	IL-15#8
83	Tumor	TCR#8	0.000000	0.000000e+00	0.000000e+00	1.000	0.583083	IL-15#8

168 rows × 8 columns