Toxicology composition analysis (tissue-wise) + severity associations¶

Purpose¶

This notebook analyzes tissue-level shifts in cell-type composition in the toxicology arm of the study and quantifies how histopathology severity relates to cell-type abundance.

It contains two main analyses:

1. scCODA compositional modeling per tissue (condition-driven composition changes)
2. Severity association analysis per (tissue × cell type), using:
   • Spearman correlation between severity score and logit(cell-type proportion)
   • An OLS fit for a smoothed trend used in visualization

Inputs¶

Primary input¶

• tox_cytokines_andata.h5ad (AnnData), loaded via scanpy.read_h5ad(...).

Required AnnData fields¶

This notebook expects these columns in adata.obs:

• tissue (e.g., Liver/Lung/Eye/Spleen/Kidney)
• batch (sample identifier for compositional modeling / clustering)
• annotation (cell-type labels)
• condition (experimental condition)
• response (severity category; mapped to 0–5 below)
• animal (used to construct per-sample covariates for scCODA)

Outputs¶

This notebook currently builds outputs in-memory and prints tables. For reproducible runs, we recommend saving:

• scCODA results per tissue (model summary + credible effects)
• severity association table across all tissues:
  • recommended filename: toxicology_severity_correlation_all_tissues.csv

Tip: add an OUTDIR = Path("./results/tox") and write all tables/figures there.

Recommended run order¶

Run top-to-bottom:

1. Load AnnData
2. Run scCODA per tissue (reference cell type is tissue-specific)
3. Build per-(batch × cell type) proportions and logit-transform
4. Map severity categories to numeric scores (0–5)
5. Compute Spearman correlations per (tissue × cell type)
6. Fit OLS per stratum for smoothed trend (for plotting)

Notes / common pitfalls¶

• If a tissue has only one condition represented, scCODA is skipped for that tissue.
• Correlations are skipped for strata with too few observations or no variation.
• response must match the severity mapping dictionary (normal/minimal/.../severe).

In [1]:

# mamba activate env_scoda 
import sccoda
import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
from matplotlib.colors import ListedColormap
from matplotlib.colors import Normalize, LinearSegmentedColormap
from matplotlib.lines import Line2D
import seaborn as sns
import os
import gzip
import numpy as np
import statsmodels.api as sm
from scipy.special import softmax
import random
import re

from sccoda.util import cell_composition_data as ccd
from sccoda.util import comp_ana as ca

from typing import Optional, Dict
import inspect
import io
import contextlib
import pandas as pd

import math
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
import anndata as ad


from sccoda.util import cell_composition_data as ccd
from sccoda.util import comp_ana as ca

from typing import Optional, Dict
import inspect
import io
import contextlib
import pandas as pd



# Use this ONLY if you fully understand the risks and want to suppress the warning
pd.options.mode.chained_assignment = None # Setting to None suppresses the warning

# To set it back to the default (which shows the warning):
# pd.options.mode.chained_assignment = 'warn'

import anndata as ad

# mamba activate env_scoda import sccoda import scanpy as sc import pandas as pd import matplotlib.pyplot as plt from matplotlib.backends.backend_pdf import PdfPages from matplotlib.colors import ListedColormap from matplotlib.colors import Normalize, LinearSegmentedColormap from matplotlib.lines import Line2D import seaborn as sns import os import gzip import numpy as np import statsmodels.api as sm from scipy.special import softmax import random import re from sccoda.util import cell_composition_data as ccd from sccoda.util import comp_ana as ca from typing import Optional, Dict import inspect import io import contextlib import pandas as pd import math RANDOM_SEED = 8927 rng = np.random.default_rng(RANDOM_SEED) import anndata as ad from sccoda.util import cell_composition_data as ccd from sccoda.util import comp_ana as ca from typing import Optional, Dict import inspect import io import contextlib import pandas as pd # Use this ONLY if you fully understand the risks and want to suppress the warning pd.options.mode.chained_assignment = None # Setting to None suppresses the warning # To set it back to the default (which shows the warning): # pd.options.mode.chained_assignment = 'warn' import anndata as ad

2026-02-25 13:45:27.300524: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used.
2026-02-25 13:45:27.384017: I tensorflow/core/platform/cpu_feature_guard.cc:210] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
To enable the following instructions: AVX2 AVX512F FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
2026-02-25 13:45:30.626379: I external/local_xla/xla/tsl/cuda/cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used.

Manuscript linkage¶

This notebook supports the toxicology composition and severity-association analyses shown in the toxicology figure panels (e.g., “Figure 4a–d” in the manuscript numbering used for ms_8).

Outputs typically used downstream:

• scCODA tissue-wise condition effects
• Tissue × cell-type severity association table (Spearman ρ, p-value, and smoothed trend)

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"
pathout = os.path.join(DATA_ROOT, "01_processed_anndata")
adata_work = sc.read_h5ad(os.path.join(pathout, "tox_cytokines_andata.h5ad"))

# 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" pathout = os.path.join(DATA_ROOT, "01_processed_anndata") adata_work = sc.read_h5ad(os.path.join(pathout, "tox_cytokines_andata.h5ad"))

In [3]:

def pick_reference_celltype(Y: pd.DataFrame) -> str:
    present_all = (Y > 0).all(axis=0)
    if present_all.any():
        candidates = Y.columns[present_all]
        return Y[candidates].mean(axis=0).sort_values(ascending=False).index[0]
    return Y.mean(axis=0).sort_values(ascending=False).index[0]


def run_sccoda_one_tissue(
    adata,
    tissue: str,
    celltype_col: str = "annotation",
    sample_col: str = "batch",
    condition_col: str = "condition",
    response_col: str = "response",
    reference_cell_type: Optional[str] = None,
    ref_ct: Optional[Dict[str, str]] = None,   # <-- ADDED (optional dict)
    hmc_draws: int = 2000,
    hmc_tune: int = 1000,
    chains: int = 2,
    target_accept: float = 0.9,
):
    obs = adata.obs[adata.obs["tissue"].astype(str) == str(tissue)].copy()

    cov = (
        obs[[sample_col, condition_col, response_col, "animal"]]
        .drop_duplicates(subset=[sample_col])
        .set_index(sample_col)
    )

    Y = pd.crosstab(obs[sample_col].astype(str), obs[celltype_col].astype(str))
    Y = Y.loc[cov.index]
    Y = Y.loc[:, Y.sum(axis=0) > 0]

    if cov[condition_col].astype(str).nunique() < 2:
        print(f"[SKIP] {tissue}: only one condition present")
        return None

    # --- CHANGED: allow tissue-specific reference from ref_ct dict ---
    ref = reference_cell_type
    if ref is None and ref_ct is not None:
        ref = ref_ct.get(str(tissue))

    # If dict provided a ref that isn't present in this tissue's Y, fallback safely
    if ref is None or ref not in Y.columns:
        ref = pick_reference_celltype(Y)

    # --- CHANGED: harden df so scCODA sees only counts as cell types ---
    df = pd.concat([cov[[condition_col]].copy(), Y], axis=1)
    count_cols = [c for c in df.columns if c != condition_col]
    df[count_cols] = df[count_cols].apply(pd.to_numeric, errors="coerce").fillna(0).astype(int)
    df = df.loc[df[count_cols].sum(axis=1) > 0]  # drop empty samples if any

    data = ccd.from_pandas(df, covariate_columns=[condition_col])

    model = ca.CompositionalAnalysis(
        data,
        formula=condition_col,
        reference_cell_type=ref
    )

    # --- version tolerant sample_hmc ---
    sig = inspect.signature(model.sample_hmc)
    kwargs = {}

    if "num_results" in sig.parameters:
        kwargs["num_results"] = hmc_draws
    if "num_burnin" in sig.parameters:
        kwargs["num_burnin"] = hmc_tune
    if "num_warmup" in sig.parameters:
        kwargs["num_warmup"] = hmc_tune

    if "target_accept_prob" in sig.parameters:
        kwargs["target_accept_prob"] = target_accept
    if "target_accept" in sig.parameters:
        kwargs["target_accept"] = target_accept

    if "chains" in sig.parameters:
        kwargs["chains"] = chains
    elif "n_chains" in sig.parameters:
        kwargs["n_chains"] = chains

    res = model.sample_hmc(**kwargs)

    # --- capture printed summary as TEXT (because summary() returns None) ---
    buf = io.StringIO()
    with contextlib.redirect_stdout(buf):
        res.summary()
    summary_text = buf.getvalue()

    print(f"\n=== {tissue} ===")
    print("df rows (samples used):", df.shape[0], "| df cols (celltypes):", len(count_cols), "| ref celltype:", ref)

    return {"tissue": tissue, "ref_celltype": ref, "result": res, "summary_text": summary_text}

def pick_reference_celltype(Y: pd.DataFrame) -> str: present_all = (Y > 0).all(axis=0) if present_all.any(): candidates = Y.columns[present_all] return Y[candidates].mean(axis=0).sort_values(ascending=False).index[0] return Y.mean(axis=0).sort_values(ascending=False).index[0] def run_sccoda_one_tissue( adata, tissue: str, celltype_col: str = "annotation", sample_col: str = "batch", condition_col: str = "condition", response_col: str = "response", reference_cell_type: Optional[str] = None, ref_ct: Optional[Dict[str, str]] = None, # <-- ADDED (optional dict) hmc_draws: int = 2000, hmc_tune: int = 1000, chains: int = 2, target_accept: float = 0.9, ): obs = adata.obs[adata.obs["tissue"].astype(str) == str(tissue)].copy() cov = ( obs[[sample_col, condition_col, response_col, "animal"]] .drop_duplicates(subset=[sample_col]) .set_index(sample_col) ) Y = pd.crosstab(obs[sample_col].astype(str), obs[celltype_col].astype(str)) Y = Y.loc[cov.index] Y = Y.loc[:, Y.sum(axis=0) > 0] if cov[condition_col].astype(str).nunique() < 2: print(f"[SKIP] {tissue}: only one condition present") return None # --- CHANGED: allow tissue-specific reference from ref_ct dict --- ref = reference_cell_type if ref is None and ref_ct is not None: ref = ref_ct.get(str(tissue)) # If dict provided a ref that isn't present in this tissue's Y, fallback safely if ref is None or ref not in Y.columns: ref = pick_reference_celltype(Y) # --- CHANGED: harden df so scCODA sees only counts as cell types --- df = pd.concat([cov[[condition_col]].copy(), Y], axis=1) count_cols = [c for c in df.columns if c != condition_col] df[count_cols] = df[count_cols].apply(pd.to_numeric, errors="coerce").fillna(0).astype(int) df = df.loc[df[count_cols].sum(axis=1) > 0] # drop empty samples if any data = ccd.from_pandas(df, covariate_columns=[condition_col]) model = ca.CompositionalAnalysis( data, formula=condition_col, reference_cell_type=ref ) # --- version tolerant sample_hmc --- sig = inspect.signature(model.sample_hmc) kwargs = {} if "num_results" in sig.parameters: kwargs["num_results"] = hmc_draws if "num_burnin" in sig.parameters: kwargs["num_burnin"] = hmc_tune if "num_warmup" in sig.parameters: kwargs["num_warmup"] = hmc_tune if "target_accept_prob" in sig.parameters: kwargs["target_accept_prob"] = target_accept if "target_accept" in sig.parameters: kwargs["target_accept"] = target_accept if "chains" in sig.parameters: kwargs["chains"] = chains elif "n_chains" in sig.parameters: kwargs["n_chains"] = chains res = model.sample_hmc(**kwargs) # --- capture printed summary as TEXT (because summary() returns None) --- buf = io.StringIO() with contextlib.redirect_stdout(buf): res.summary() summary_text = buf.getvalue() print(f"\n=== {tissue} ===") print("df rows (samples used):", df.shape[0], "| df cols (celltypes):", len(count_cols), "| ref celltype:", ref) return {"tissue": tissue, "ref_celltype": ref, "result": res, "summary_text": summary_text}

In [4]:

ref_ct = {
    'Lung': 'Lymphatic endothelial cells',
    'Liver': 'Endothelial cells',
    'Kidney': 'Fibroblasts',
    'Eye': 'Muscle',
    'Spleen': 'Endothelial cells'
}

results = {}
for ts in ["Liver", "Lung", "Eye", "Spleen", "Kidney"]:
    out = run_sccoda_one_tissue(adata_work, ts, ref_ct=ref_ct)  # <-- CHANGED: pass dict
    if out is not None:
        results[ts] = out
        # with open(f"sccoda_{ts}_summary_select_ref.txt", "w") as f:
        #     f.write(out["summary_text"])

ref_ct = { 'Lung': 'Lymphatic endothelial cells', 'Liver': 'Endothelial cells', 'Kidney': 'Fibroblasts', 'Eye': 'Muscle', 'Spleen': 'Endothelial cells' } results = {} for ts in ["Liver", "Lung", "Eye", "Spleen", "Kidney"]: out = run_sccoda_one_tissue(adata_work, ts, ref_ct=ref_ct) # <-- CHANGED: pass dict if out is not None: results[ts] = out # with open(f"sccoda_{ts}_summary_select_ref.txt", "w") as f: # f.write(out["summary_text"])

2026-02-25 13:45:48.559053: E external/local_xla/xla/stream_executor/cuda/cuda_platform.cc:51] failed call to cuInit: INTERNAL: CUDA error: Failed call to cuInit: UNKNOWN ERROR (303)
2026-02-25 13:45:51.579484: I external/local_xla/xla/service/service.cc:163] XLA service 0x155088007300 initialized for platform Host (this does not guarantee that XLA will be used). Devices:
2026-02-25 13:45:51.579509: I external/local_xla/xla/service/service.cc:171]   StreamExecutor device (0): Host, Default Version
  0%|          | 0/2000 [00:00<?, ?it/s]2026-02-25 13:45:51.755182: I tensorflow/compiler/mlir/tensorflow/utils/dump_mlir_util.cc:269] disabling MLIR crash reproducer, set env var `MLIR_CRASH_REPRODUCER_DIRECTORY` to enable.
WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1772045153.006754 2130107 device_compiler.h:196] Compiled cluster using XLA!  This line is logged at most once for the lifetime of the process.
100%|██████████| 2000/2000 [00:12<00:00, 156.93it/s]

MCMC sampling finished. (21.595 sec)
Acceptance rate: 30.1%

=== Liver ===
df rows (samples used): 8 | df cols (celltypes): 10 | ref celltype: Endothelial cells
Zero counts encountered in data! Added a pseudocount of 0.5.

100%|██████████| 2000/2000 [00:13<00:00, 151.88it/s]

MCMC sampling finished. (22.120 sec)
Acceptance rate: 36.8%

=== Lung ===
df rows (samples used): 7 | df cols (celltypes): 16 | ref celltype: Lymphatic endothelial cells
Zero counts encountered in data! Added a pseudocount of 0.5.

100%|██████████| 2000/2000 [00:10<00:00, 184.32it/s]

MCMC sampling finished. (18.582 sec)
Acceptance rate: 70.2%

=== Eye ===
df rows (samples used): 6 | df cols (celltypes): 9 | ref celltype: Muscle

100%|██████████| 2000/2000 [00:12<00:00, 162.51it/s]

MCMC sampling finished. (20.558 sec)
Acceptance rate: 46.6%

=== Spleen ===
df rows (samples used): 7 | df cols (celltypes): 12 | ref celltype: Endothelial cells
Zero counts encountered in data! Added a pseudocount of 0.5.

100%|██████████| 2000/2000 [00:13<00:00, 151.48it/s]

WARNING:tensorflow:5 out of the last 5 calls to <function CompositionalModel.sampling.<locals>.sample_mcmc at 0x154da376d310> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has reduce_retracing=True option that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for  more details.
MCMC sampling finished. (22.297 sec)
Acceptance rate: 45.1%

=== Kidney ===
df rows (samples used): 9 | df cols (celltypes): 13 | ref celltype: Fibroblasts

In [5]:

def logit(p):
    eps = 1e-4
    p = np.clip(p, eps, 1 - eps)
    return np.log(p / (1 - p))

def campute_df(df_in,tissue):
    df_tiss = df_in[df_in['tissue']==tissue].copy()
    batches = df_tiss['batch'].unique().tolist()
    data_container = []
    for b in batches:
        df_tissu_curr = df_tiss[df_tiss['batch']==b].copy()
        total_cells = len(df_tissu_curr)
        condtion = df_tissu_curr['condition'].iloc[0]
        response = df_tissu_curr['response'].iloc[0]
        celltypes = df_tissu_curr['annotation'].unique().tolist() 
        for ct in celltypes:
            n_ct = len(df_tissu_curr[df_tissu_curr['annotation']==ct])
            prop = n_ct / total_cells
            logit_prop = logit(prop)
            data_container.append({
                'tissue': tissue,
                'condition': condtion,
                'batch': b,
                'annotation': ct,
                'n_cells': n_ct,
                'n_total': total_cells,
                'prop': prop,
                'logit_p': logit_prop,
                'response': response
            })
    return pd.DataFrame(data_container)    

def logit(p): eps = 1e-4 p = np.clip(p, eps, 1 - eps) return np.log(p / (1 - p)) def campute_df(df_in,tissue): df_tiss = df_in[df_in['tissue']==tissue].copy() batches = df_tiss['batch'].unique().tolist() data_container = [] for b in batches: df_tissu_curr = df_tiss[df_tiss['batch']==b].copy() total_cells = len(df_tissu_curr) condtion = df_tissu_curr['condition'].iloc[0] response = df_tissu_curr['response'].iloc[0] celltypes = df_tissu_curr['annotation'].unique().tolist() for ct in celltypes: n_ct = len(df_tissu_curr[df_tissu_curr['annotation']==ct]) prop = n_ct / total_cells logit_prop = logit(prop) data_container.append({ 'tissue': tissue, 'condition': condtion, 'batch': b, 'annotation': ct, 'n_cells': n_ct, 'n_total': total_cells, 'prop': prop, 'logit_p': logit_prop, 'response': response }) return pd.DataFrame(data_container)

In [6]:

from tqdm import tqdm
df_in = adata_work.obs[['tissue', 'batch','annotation','response','condition']].copy()
tissues = df_in['tissue'].unique().tolist()
table_correlation_container = []
for tissue in tqdm(tissues):
    table_correlation_container.append(
        campute_df(df_in,tissue)
    )
table_correlation_prep = pd.concat(table_correlation_container)
sev_map = {
    "normal": 0,
    "minimal": 1,
    "minimal-mild": 2,
    "mild": 3,
    "mild-moderate": 4,
    "severe": 5,
}
table_correlation_prep['sev_score'] = table_correlation_prep['response'].map(sev_map)
table_correlation_prep['celltype_batch_pair'] = table_correlation_prep['annotation'] + '_' + table_correlation_prep['batch']

from tqdm import tqdm df_in = adata_work.obs[['tissue', 'batch','annotation','response','condition']].copy() tissues = df_in['tissue'].unique().tolist() table_correlation_container = [] for tissue in tqdm(tissues): table_correlation_container.append( campute_df(df_in,tissue) ) table_correlation_prep = pd.concat(table_correlation_container) sev_map = { "normal": 0, "minimal": 1, "minimal-mild": 2, "mild": 3, "mild-moderate": 4, "severe": 5, } table_correlation_prep['sev_score'] = table_correlation_prep['response'].map(sev_map) table_correlation_prep['celltype_batch_pair'] = table_correlation_prep['annotation'] + '_' + table_correlation_prep['batch']

100%|██████████| 5/5 [00:00<00:00, 12.79it/s]

In [8]:

table_correlation_prep

table_correlation_prep

Out[8]:

	tissue	condition	batch	annotation	n_cells	n_total	prop	logit_p	response	sev_score	celltype_batch_pair
0	Liver	IL-15	0	Hepatocytes	3093	10278	0.300934	-0.842854	severe	5	Hepatocytes_0
1	Liver	IL-15	0	Endothelial cells	1505	10278	0.146429	-1.762886	severe	5	Endothelial cells_0
2	Liver	IL-15	0	T-cells	2231	10278	0.217066	-1.282849	severe	5	T-cells_0
3	Liver	IL-15	0	cDC1	1601	10278	0.155770	-1.690047	severe	5	cDC1_0
4	Liver	IL-15	0	Myeloid cells	961	10278	0.093501	-2.271622	severe	5	Myeloid cells_0
...	...	...	...	...	...	...	...	...	...	...	...
79	Spleen	IL-15	8	Stromal cells	1096	39403	0.027815	-3.553965	mild	3	Stromal cells_8
80	Spleen	IL-15	8	Endothelial cells	2828	39403	0.071771	-2.559795	mild	3	Endothelial cells_8
81	Spleen	IL-15	8	Progenitor	3471	39403	0.088090	-2.337186	mild	3	Progenitor_8
82	Spleen	IL-15	8	NK cells	215	39403	0.005456	-5.205488	mild	3	NK cells_8
83	Spleen	IL-15	8	Megakaryocytes	88	39403	0.002233	-6.102025	mild	3	Megakaryocytes_8

425 rows × 11 columns

Severity association model (used for plotting and summary tables)¶

For each (tissue × cell type) stratum we quantify the relationship between:

• severity score (ordinal; mapped to numeric 0–5), and
• cell-type abundance, summarized as the logit-transformed proportion of that cell type within each batch/sample.

What we compute¶

1. Spearman rank correlation (ρ) between severity and logit(proportion)

   • Reported with a two-sided p-value
   • Strata are skipped if:
     • fewer than 3 observations, or
     • no variation in severity, or
     • no variation in logit(proportion)

2. OLS regression of severity on logit(proportion) for a smoothed trend line

   • Used only for visualization (“smoothed severity score”)
   • If multiple batches exist, we use cluster-robust standard errors (cluster = batch)

In [11]:

from scipy.stats import spearmanr
def ols_prediction(dfm, batch_col="batch", x_col="logit_p", y_col="sev_score"):
    d = dfm.dropna(subset=[x_col, y_col, batch_col]).copy()

    # Need at least 3 rows for intercept+slope with df_resid > 0
    if len(d) < 3:
        return pd.Series([np.nan] * len(dfm), index=dfm.index)

    # Need variation in x
    if d[x_col].nunique() < 2:
        return pd.Series([np.nan] * len(dfm), index=dfm.index)

    X = sm.add_constant(d[x_col])
    y = d[y_col]

    n_groups = d[batch_col].nunique()

    # Cluster-robust only if >=2 groups and enough df
    if n_groups >= 2 and (len(d) - X.shape[1]) > 0:
        fit = sm.OLS(y, X).fit(cov_type="cluster", cov_kwds={"groups": d[batch_col]})
    else:
        fit = sm.OLS(y, X).fit()

    pred = fit.predict(X)
    # align back to original dfm index
    out = pd.Series(np.nan, index=dfm.index, dtype=float)
    out.loc[d.index] = pred
    return out

def compute_spearman_corr(df, tissue):
    df_tiss = df[df["tissue"] == tissue].copy()
    df_tiss = df_tiss.dropna(subset=["sev_score", "logit_p", "batch"]).copy()

    records = []
    for ct, df_ct_cur in df_tiss.groupby("annotation", sort=False):
        # Spearman requires at least 2 points with some variation
        if len(df_ct_cur) < 3 or df_ct_cur["sev_score"].nunique() < 2 or df_ct_cur["logit_p"].nunique() < 2:
            df_ct_cur = df_ct_cur.copy()
            df_ct_cur["score_smoothed"] = np.nan
            df_ct_cur["spearman_rho"] = np.nan
            df_ct_cur["spearman_pval"] = np.nan
            records.append(df_ct_cur)
            continue

        rho, p_rho = spearmanr(df_ct_cur["sev_score"], df_ct_cur["logit_p"])
        df_ct_cur = df_ct_cur.copy()
        df_ct_cur["score_smoothed"] = ols_prediction(df_ct_cur)
        df_ct_cur["spearman_rho"] = rho
        df_ct_cur["spearman_pval"] = p_rho
        records.append(df_ct_cur)

    return pd.concat(records, axis=0)

from scipy.stats import spearmanr def ols_prediction(dfm, batch_col="batch", x_col="logit_p", y_col="sev_score"): d = dfm.dropna(subset=[x_col, y_col, batch_col]).copy() # Need at least 3 rows for intercept+slope with df_resid > 0 if len(d) < 3: return pd.Series([np.nan] * len(dfm), index=dfm.index) # Need variation in x if d[x_col].nunique() < 2: return pd.Series([np.nan] * len(dfm), index=dfm.index) X = sm.add_constant(d[x_col]) y = d[y_col] n_groups = d[batch_col].nunique() # Cluster-robust only if >=2 groups and enough df if n_groups >= 2 and (len(d) - X.shape[1]) > 0: fit = sm.OLS(y, X).fit(cov_type="cluster", cov_kwds={"groups": d[batch_col]}) else: fit = sm.OLS(y, X).fit() pred = fit.predict(X) # align back to original dfm index out = pd.Series(np.nan, index=dfm.index, dtype=float) out.loc[d.index] = pred return out def compute_spearman_corr(df, tissue): df_tiss = df[df["tissue"] == tissue].copy() df_tiss = df_tiss.dropna(subset=["sev_score", "logit_p", "batch"]).copy() records = [] for ct, df_ct_cur in df_tiss.groupby("annotation", sort=False): # Spearman requires at least 2 points with some variation if len(df_ct_cur) < 3 or df_ct_cur["sev_score"].nunique() < 2 or df_ct_cur["logit_p"].nunique() < 2: df_ct_cur = df_ct_cur.copy() df_ct_cur["score_smoothed"] = np.nan df_ct_cur["spearman_rho"] = np.nan df_ct_cur["spearman_pval"] = np.nan records.append(df_ct_cur) continue rho, p_rho = spearmanr(df_ct_cur["sev_score"], df_ct_cur["logit_p"]) df_ct_cur = df_ct_cur.copy() df_ct_cur["score_smoothed"] = ols_prediction(df_ct_cur) df_ct_cur["spearman_rho"] = rho df_ct_cur["spearman_pval"] = p_rho records.append(df_ct_cur) return pd.concat(records, axis=0)

In [12]:

table_correlation = table_correlation_prep.copy()
tissues = table_correlation['tissue'].unique().tolist()
table_correlation_container = []
for tissue in tqdm(tissues):
    table_correlation_container.append(compute_spearman_corr(table_correlation_prep, tissue))

table_correlation = table_correlation_prep.copy() tissues = table_correlation['tissue'].unique().tolist() table_correlation_container = [] for tissue in tqdm(tissues): table_correlation_container.append(compute_spearman_corr(table_correlation_prep, tissue))

100%|██████████| 5/5 [00:00<00:00,  8.50it/s]

In [13]:

pd.concat(table_correlation_container, axis=0)

pd.concat(table_correlation_container, axis=0)

Out[13]:

	tissue	condition	batch	annotation	n_cells	n_total	prop	logit_p	response	sev_score	celltype_batch_pair	score_smoothed	spearman_rho	spearman_pval
0	Liver	IL-15	0	Hepatocytes	3093	10278	0.300934	-0.842854	severe	5	Hepatocytes_0	3.068054	-0.204834	0.626544
10	Liver	IL-15	10	Hepatocytes	3083	8458	0.364507	-0.555855	mild-moderate	4	Hepatocytes_10	2.742304	-0.204834	0.626544
20	Liver	IL-15	15	Hepatocytes	3393	7968	0.425828	-0.298892	mild-moderate	4	Hepatocytes_15	2.450645	-0.204834	0.626544
31	Liver	Both	25	Hepatocytes	3426	6310	0.542948	0.172215	minimal-mild	2	Hepatocytes_25	1.915927	-0.204834	0.626544
40	Liver	Both	30	Hepatocytes	3569	7247	0.492480	-0.030084	minimal	1	Hepatocytes_30	2.145541	-0.204834	0.626544
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
35	Spleen	Both	33	Megakaryocytes	26	41618	0.000625	-7.377567	minimal	1	Megakaryocytes_33	1.575991	0.327327	0.473597
47	Spleen	Both	38	Megakaryocytes	65	40038	0.001623	-6.421572	minimal	1	Megakaryocytes_38	2.064595	0.327327	0.473597
59	Spleen	IL-15	4	Megakaryocytes	45	40161	0.001120	-6.792868	severe	5	Megakaryocytes_4	1.874827	0.327327	0.473597
71	Spleen	TCR	43	Megakaryocytes	42	32194	0.001305	-6.640560	normal	0	Megakaryocytes_43	1.952671	0.327327	0.473597
83	Spleen	IL-15	8	Megakaryocytes	88	39403	0.002233	-6.102025	mild	3	Megakaryocytes_8	2.227914	0.327327	0.473597

425 rows × 14 columns