scCODA / tascCODA¶
细胞组成的主模型。scCODA 看大类，tascCODA 看亚型（树：大类→亚型），对照为 CTRL。每个点是一例的占比。模型处理「占比加总为 100%」。看可信区间是否离开无变化；树图上先看中性粒、Pre-Neut、GMP。
In [1]:
import os
from pathlib import Path
os.environ.setdefault("QT_QPA_PLATFORM", "offscreen")
os.environ.setdefault("JAX_PLATFORMS", "cpu")
import anndata as ad
import jax.numpy as jnp
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from scipy.stats import mannwhitneyu
from statsmodels.stats.multitest import multipletests
from session_info2 import session_info
import numpyro
import pertpy as pt
from pertpy.tools._coda._base_coda import CompositionalModel2
from pertpy.tools._coda._base_coda import MCMC as _CodaMCMC
numpyro.enable_validation(False)
# pertpy writes n_total into AnnData.obsm as (n, 1); zeros become 0.5.
if not getattr(CompositionalModel2, "_n_total_squeezed", False):
_orig_run_mcmc = CompositionalModel2._CompositionalModel2__run_mcmc
def _run_mcmc_squeeze_n_total(
self, sample_adata, kernel, rng_key, copy=False, *args, **kwargs
):
_mcmc_run = _CodaMCMC.run
def _run(
mcmc_self,
rng_key,
counts,
covariates,
n_total,
ref_index,
sample_adata,
extra_fields=(),
):
# obsm forces n_total (n,1); 0.5 pseudocounts are not integers.
# numpyro DirichletMultinomial then rejects total_count.
counts_np = np.asarray(counts, dtype=np.float64)
counts_np = np.where(np.isclose(counts_np, 0.5), 1.0, counts_np)
counts_np = np.rint(counts_np).astype(np.int32)
n_total = jnp.asarray(counts_np.sum(axis=1).astype(np.int32))
counts = jnp.asarray(counts_np)
ref = np.asarray(ref_index).reshape(-1)
ref_index = int(ref[0]) if ref.size == 1 else jnp.asarray(ref)
return _mcmc_run(
mcmc_self,
rng_key,
counts,
covariates,
n_total,
ref_index,
sample_adata,
extra_fields=extra_fields,
)
_CodaMCMC.run = _run
try:
return _orig_run_mcmc(
self, sample_adata, kernel, rng_key, copy, *args, **kwargs
)
finally:
_CodaMCMC.run = _mcmc_run
CompositionalModel2._CompositionalModel2__run_mcmc = _run_mcmc_squeeze_n_total
CompositionalModel2._n_total_squeezed = True
# Same order as codes/plot_colors.R (annot_main_col / annot_sub_col).
ANNOT_MAIN_ORDER = [
"HSPC", "B", "Plasma", "T", "NK",
"Mono/DC", "pDC", "Pre-Neut", "Neut", "Ery",
]
ANNOT_SUB_ORDER = [
"HSC", "MEP", "GMP", "MDP",
"Pro-B", "Large Pre-B", "Small Pre-B", "Trans B", "Naive B", "Memory B",
"Plasma", "pDC",
"CD4 Tnaive", "CD4 Tcm", "CD4 Treg",
"CD8 Tnaive", "CD8 Tem", "CD8 CTL", "CD8 NKT", "MAIT", "gdT",
"CD56bright NK", "CD56dim NK",
"CD14 Mono", "CD16 Mono", "cDC2",
"Promyelocyte", "MMP8 Pre-Neut", "MMP9 Pre-Neut",
"S100A12 Neut", "SLC8A1 Neut", "ISG15 Neut",
"Early Ery", "Late Ery",
]
SUB_TO_MAIN = {
"HSC": "HSPC", "MEP": "HSPC", "GMP": "HSPC", "MDP": "HSPC",
"Pro-B": "B", "Large Pre-B": "B", "Small Pre-B": "B",
"Trans B": "B", "Naive B": "B", "Memory B": "B",
"Plasma": "Plasma",
"CD4 Tnaive": "T", "CD4 Tcm": "T", "CD4 Treg": "T",
"CD8 Tnaive": "T", "CD8 Tem": "T", "CD8 CTL": "T", "CD8 NKT": "T",
"MAIT": "T", "gdT": "T",
"CD56bright NK": "NK", "CD56dim NK": "NK",
"CD14 Mono": "Mono/DC", "CD16 Mono": "Mono/DC", "cDC2": "Mono/DC",
"pDC": "pDC",
"Promyelocyte": "Pre-Neut", "MMP8 Pre-Neut": "Pre-Neut", "MMP9 Pre-Neut": "Pre-Neut",
"S100A12 Neut": "Neut", "SLC8A1 Neut": "Neut", "ISG15 Neut": "Neut",
"Early Ery": "Ery", "Late Ery": "Ery",
}
def _node_label(n):
return (n.name or "").removeprefix("annot_main_")
def _node_rank(n):
name = _node_label(n)
main = name if name in ANNOT_MAIN_ORDER else SUB_TO_MAIN.get(name)
main_i = ANNOT_MAIN_ORDER.index(main) if main in ANNOT_MAIN_ORDER else 999
sub_i = ANNOT_SUB_ORDER.index(name) if name in ANNOT_SUB_ORDER else -1
return (main_i, sub_i, name)
def sort_coda_tree(tree):
"""Reorder children to match plot_colors.R lineage order (HSPC at top)."""
for n in tree.traverse("postorder"):
if n.children:
n.children.sort(key=_node_rank)
return tree
def plot_coda_tree(
tree,
ax=None,
node_value=None,
figsize=(8, 9),
pos_color="#F07A7A",
neg_color="#7EB8E8",
):
"""Draw an ete4 tree with matplotlib. Avoids ete4.treeview / Qt (hangs Jupyter)."""
if ax is None:
_, ax = plt.subplots(figsize=figsize)
sort_coda_tree(tree)
leaves = [n for n in tree.traverse("preorder") if n.is_leaf]
nleaf = len(leaves)
ymap = {id(leaf): float(nleaf - 1 - i) for i, leaf in enumerate(leaves)}
def y_of(node):
k = id(node)
if k not in ymap:
ymap[k] = float(np.mean([y_of(c) for c in node.children]))
return ymap[k]
def depth_of(node):
d = 0
cur = node
while cur.up is not None:
d += 1
cur = cur.up
return d
y_of(tree)
xs = {id(n): depth_of(n) for n in tree.traverse()}
vmax = 1.0
if node_value is not None and len(node_value):
vmax = max(float(np.nanmax(np.abs(np.asarray(node_value, dtype=float)))), 1e-8)
for n in tree.traverse():
xn, yn = xs[id(n)], y_of(n)
for c in n.children:
xc, yc = xs[id(c)], y_of(c)
ax.plot([xn, xn], [yn, yc], color="#8A8A8A", lw=0.8)
ax.plot([xn, xc], [yc, yc], color="#8A8A8A", lw=0.8)
name = _node_label(n)
if n.is_leaf:
ax.text(xn + 0.08, yn, name, va="center", ha="left", fontsize=8)
elif name:
ax.text(
xn - 0.06, yn + 0.15, name,
va="bottom", ha="right", fontsize=7, color="#666666",
)
if node_value is not None and n.name in getattr(node_value, "index", []):
v = float(node_value.loc[n.name])
if v != 0:
ax.scatter(
[xn], [yn],
s=25 + 200 * abs(v) / vmax,
c=pos_color if v > 0 else neg_color,
zorder=3, edgecolors="white", linewidths=0.4,
)
ax.set_axis_off()
ax.set_xlim(-0.8, max(xs.values()) + 2.8)
ax.set_ylim(-1, nleaf)
return ax
GROUP_COL = {"CTRL": "#7EB8E8", "ACI": "#F07A7A"}
GROUP_ORDER = ["CTRL", "ACI"]
def _canon_cell(name):
s = str(name)
for prefix in ("annot_main_", "annot_sub_"):
if s.startswith(prefix):
return s[len(prefix) :]
return s
def ordered_cell_types(names, kind="sub"):
"""Order cell types as names() of annot_main_col / annot_sub_col."""
order = ANNOT_MAIN_ORDER if kind == "main" else ANNOT_SUB_ORDER
rank = {n: i for i, n in enumerate(order)}
uniq = list(dict.fromkeys(map(str, names)))
return sorted(uniq, key=lambda n: (rank.get(_canon_cell(n), 999), n))
def reorder_coda_modality(mdata, kind="main", modality_key="coda"):
"""Put coda columns in plot_colors.R order (also fixes pertpy.plot_boxplots)."""
coda = mdata[modality_key]
order = ordered_cell_types(list(coda.var_names), kind=kind)
mdata.mod[modality_key] = coda[:, order].copy()
return order
def _wilcox_stars(p):
if p < 0.001:
return "***"
if p < 0.01:
return "**"
if p < 0.05:
return "*"
return "ns"
def plot_coda_boxplots(
data,
modality_key="coda",
cell_kind="sub",
figsize=(18, 6),
add_stats=True,
):
"""Grouped boxplots in plot_colors.R order, with Wilcoxon CTRL vs ACI."""
coda = data[modality_key]
x = np.asarray(coda.X, dtype=float)
props = x / np.clip(x.sum(axis=1, keepdims=True), 1e-12, None)
df = pd.DataFrame(props, columns=[str(c) for c in coda.var_names], index=coda.obs_names)
df["group"] = pd.Categorical(
coda.obs["group"].astype(str), categories=GROUP_ORDER, ordered=True
)
plot_df = df.melt(id_vars="group", var_name="Cell type", value_name="Proportion")
order = ordered_cell_types(plot_df["Cell type"].unique(), kind=cell_kind)
plot_df["Cell type"] = pd.Categorical(
plot_df["Cell type"], categories=order, ordered=True
)
fig, ax = plt.subplots(figsize=figsize)
sns.boxplot(
data=plot_df,
x="Cell type",
y="Proportion",
hue="group",
order=order,
hue_order=GROUP_ORDER,
palette=GROUP_COL,
fliersize=0,
width=0.7,
ax=ax,
)
sns.stripplot(
data=plot_df,
x="Cell type",
y="Proportion",
hue="group",
order=order,
hue_order=GROUP_ORDER,
palette=["#444444", "#444444"],
dodge=True,
size=3.5,
alpha=0.85,
legend=False,
ax=ax,
)
stats_rows = []
n_hue = 2
box_width = 0.7
step = box_width / n_hue
offsets = [-box_width / 2 + step / 2 + j * step for j in range(n_hue)]
yspan = max(float(plot_df["Proportion"].max()), 1e-6)
for i, ct in enumerate(order):
ctrl = plot_df.loc[
(plot_df["Cell type"] == ct) & (plot_df["group"] == "CTRL"), "Proportion"
]
aci = plot_df.loc[
(plot_df["Cell type"] == ct) & (plot_df["group"] == "ACI"), "Proportion"
]
p = float(
mannwhitneyu(
aci.to_numpy(), ctrl.to_numpy(), alternative="two-sided"
).pvalue
)
stats_rows.append(
{
"cell": _canon_cell(ct),
"n_ctrl": int(ctrl.size),
"n_aci": int(aci.size),
"median_ctrl": float(ctrl.median()) if len(ctrl) else np.nan,
"median_aci": float(aci.median()) if len(aci) else np.nan,
"p": p,
}
)
if not add_stats or p >= 0.05:
continue
star = _wilcox_stars(p)
y = max(ctrl.max(), aci.max()) + 0.04 * yspan
x0, x1 = i + offsets[0], i + offsets[1]
h = 0.015 * yspan
ax.plot([x0, x0, x1, x1], [y, y + h, y + h, y], color="0.35", lw=0.8)
ax.text(
(x0 + x1) / 2, y + h, star,
ha="center", va="bottom", fontsize=8, color="0.15",
)
stats = pd.DataFrame(stats_rows)
stats["p_adj"] = multipletests(stats["p"], method="fdr_bh")[1]
handles, labels = ax.get_legend_handles_labels()
ax.legend(
handles[:2], labels[:2], title="group",
loc="upper left", bbox_to_anchor=(1, 1),
)
ax.set_xlabel(None)
ax.tick_params(axis="x", rotation=90)
ax.set_ylim(0, plot_df["Proportion"].max() * 1.28)
ax.set_title(
"Wilcoxon rank-sum CTRL vs ACI; stars: uncorrected p "
"(* <0.05, ** <0.01, *** <0.001). BH FDR in the CSV.",
loc="left", fontsize=9, pad=8,
)
fig.tight_layout()
return ax, stats
ROOT = Path("..").resolve()
META = ROOT / "data/rdata/composition/cell_meta.csv"
RES = ROOT / "results/composition"
RES.mkdir(parents=True, exist_ok=True)
cell_meta = pd.read_csv(META)
cell_meta["group"] = pd.Categorical(
cell_meta["group"], categories=["CTRL", "ACI"], ordered=True
)
cell_meta.head()
/home/han/anaconda3/envs/pertpy/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html from .autonotebook import tqdm as notebook_tqdm
Out[1]:
| barcode | sample | group | annot_main | annot_sub | |
|---|---|---|---|---|---|
| 0 | ctrl_1_AAACCTGCAAGCGCTC-1 | CTRL-1 | CTRL | HSPC | HSC |
| 1 | ctrl_1_AAACCTGCAGACGCAA-1 | CTRL-1 | CTRL | B | Memory B |
| 2 | ctrl_1_AAACGGGCAGATGGCA-1 | CTRL-1 | CTRL | T | CD4 Tnaive |
| 3 | ctrl_1_AAACGGGCAGGTCCAC-1 | CTRL-1 | CTRL | Mono/DC | CD14 Mono |
| 4 | ctrl_1_AAACGGGGTACCATCA-1 | CTRL-1 | CTRL | B | Naive B |
In [2]:
adata = ad.AnnData(
X=np.zeros((cell_meta.shape[0], 1), dtype=np.float32),
obs=cell_meta.copy(),
)
adata.obs_names = cell_meta["barcode"].astype(str).values
adata
Out[2]:
AnnData object with n_obs × n_vars = 61190 × 1
obs: 'barcode', 'sample', 'group', 'annot_main', 'annot_sub'
1. scCODA (annot_main)¶
reference_cell_type="automatic" 选相对稳定的一类作参照。成分变化都是相对该参照。
大类水平。优先看中性粒这一枝的效应和可信区间。
In [3]:
sccoda = pt.tl.Sccoda()
sccoda_data = sccoda.load(
adata,
type="cell_level",
generate_sample_level=True,
cell_type_identifier="annot_main",
sample_identifier="sample",
covariate_obs=["group"],
)
reorder_coda_modality(sccoda_data, kind="main")
print(list(sccoda_data["coda"].var_names))
sccoda_data
['HSPC', 'B', 'Plasma', 'T', 'NK', 'Mono/DC', 'pDC', 'Pre-Neut', 'Neut', 'Ery']
Out[3]:
MuData object with n_obs × n_vars = 61201 × 11
2 modalities
rna: 61190 × 1
obs: 'barcode', 'sample', 'group', 'annot_main', 'annot_sub', 'scCODA_sample_id'
coda: 11 × 10
obs: 'sample', 'group'
var: 'n_cells'
In [4]:
reorder_coda_modality(sccoda_data, kind="main")
ax, sccoda_box_stats = plot_coda_boxplots(
sccoda_data, modality_key="coda", cell_kind="main", figsize=(12, 5)
)
sccoda_box_stats.to_csv(RES / "03_sccoda_box_wilcox.csv", index=False)
plt.savefig(RES / "03_sccoda_box.pdf", bbox_inches="tight")
plt.show()
sccoda_box_stats
Out[4]:
| cell | n_ctrl | n_aci | median_ctrl | median_aci | p | p_adj | |
|---|---|---|---|---|---|---|---|
| 0 | HSPC | 6 | 5 | 0.064219 | 0.055247 | 0.329004 | 0.329004 |
| 1 | B | 6 | 5 | 0.128233 | 0.019295 | 0.004329 | 0.007215 |
| 2 | Plasma | 6 | 5 | 0.021307 | 0.001527 | 0.008658 | 0.010823 |
| 3 | T | 6 | 5 | 0.311607 | 0.099641 | 0.004329 | 0.007215 |
| 4 | NK | 6 | 5 | 0.045914 | 0.015226 | 0.082251 | 0.091390 |
| 5 | Mono/DC | 6 | 5 | 0.097562 | 0.043791 | 0.008658 | 0.010823 |
| 6 | pDC | 6 | 5 | 0.015436 | 0.002356 | 0.004329 | 0.007215 |
| 7 | Pre-Neut | 6 | 5 | 0.048508 | 0.279479 | 0.004329 | 0.007215 |
| 8 | Neut | 6 | 5 | 0.030536 | 0.499515 | 0.004329 | 0.007215 |
| 9 | Ery | 6 | 5 | 0.181554 | 0.003193 | 0.004329 | 0.007215 |
In [5]:
sccoda_data = sccoda.prepare(
sccoda_data,
modality_key="coda",
formula="C(group, Treatment('CTRL'))",
reference_cell_type="automatic",
)
sccoda.run_nuts(sccoda_data, modality_key="coda")
sccoda.summary(sccoda_data, modality_key="coda")
sample: 100%|██████████| 11000/11000 [00:42<00:00, 258.61it/s, 255 steps of size 1.86e-02. acc. prob=0.75]
Compositional Analysis summary ┌────────────────────────────────────────────┬────────────────────────────────────────────────────────────────────┐ │ Name │ Value │ ├────────────────────────────────────────────┼────────────────────────────────────────────────────────────────────┤ │ Data │ Data: 11 samples, 10 cell types │ │ Reference cell type │ HSPC │ │ Formula │ C(group, Treatment('CTRL')) │ └────────────────────────────────────────────┴────────────────────────────────────────────────────────────────────┘
┌─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐ │ Intercepts │ ├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤ │ Final Parameter Expected Sample │ │ annot_main │ │ HSPC 1.045 526.577 │ │ B 1.498 828.319 │ │ Plasma -0.304 136.647 │ │ T 2.162 1609.048 │ │ NK 0.434 285.830 │ │ Mono/DC 1.184 605.102 │ │ pDC -0.295 137.882 │ │ Pre-Neut 0.452 291.022 │ │ Neut -0.125 163.432 │ │ Ery 1.665 978.869 │ └─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
┌─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐ │ Effects │ ├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤ │ Final Parameter Expected Sample log2-fold change │ │ Covariate Cell Type │ │ C(group, Treatment('CTRL'))T.ACI HSPC 0.000 237.568 -1.148 │ │ B -1.167 116.321 -2.832 │ │ Plasma 0.000 61.649 -1.148 │ │ T 0.000 725.930 -1.148 │ │ NK 0.000 128.954 -1.148 │ │ Mono/DC 0.000 272.995 -1.148 │ │ pDC 0.000 62.206 -1.148 │ │ Pre-Neut 2.545 1673.599 2.524 │ │ Neut 3.409 2228.287 3.769 │ │ Ery -2.079 55.219 -4.148 │ └─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
In [6]:
sccoda.plot_effects_barplot(
sccoda_data, modality_key="coda", parameter="Final Parameter"
)
plt.savefig(RES / "03_sccoda_effects.pdf", bbox_inches="tight")
plt.show()
cred = sccoda.credible_effects(sccoda_data, modality_key="coda")
cred.to_csv(RES / "03_sccoda_credible.csv")
cred
/home/han/anaconda3/envs/pertpy/lib/python3.11/site-packages/seaborn/axisgrid.py:854: FutureWarning: Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `x` variable to `hue` and set `legend=False` for the same effect. func(*plot_args, **plot_kwargs)
<Figure size 640x480 with 0 Axes>
Out[6]:
Covariate Cell Type
C(group, Treatment('CTRL'))[T.ACI] HSPC False
B True
Plasma False
T False
NK False
Mono/DC False
pDC False
Pre-Neut True
Neut True
Ery True
Name: Final Parameter, dtype: bool
2. tascCODA (annot_sub tree)¶
pen_args={"phi": 0} 不偏向聚到上层还是叶子。11 个样本后验会比较宽。
树图用 matplotlib（plot_coda_tree），不要调用 plot_draw_tree / plot_draw_effects：它们走 ete4+Qt 的 %%inline，无显示器时会卡死或把 kernel 打挂。
亚型挂在大类下面。先看 GMP、Pre-Neut、成熟粒所在的枝。
In [7]:
tasccoda = pt.tl.Tasccoda()
tasc_data = tasccoda.load(
adata,
type="cell_level",
cell_type_identifier="annot_sub",
sample_identifier="sample",
covariate_obs=["group"],
levels_orig=["annot_main", "annot_sub"],
add_level_name=True,
)
reorder_coda_modality(tasc_data, kind="sub")
# Do not display MuData/ete Tree — IPython may call Qt treeview and crash the kernel.
print(tasc_data)
print(tasc_data["coda"])
print(list(tasc_data["coda"].var_names))
print(tasc_data["coda"].uns["tree"])
MuData object with n_obs × n_vars = 61201 × 35
2 modalities
rna: 61190 × 1
obs: 'barcode', 'sample', 'group', 'annot_main', 'annot_sub', 'scCODA_sample_id'
coda: 11 × 34
obs: 'sample', 'group'
var: 'n_cells'
uns: 'tree'
AnnData object with n_obs × n_vars = 11 × 34
obs: 'sample', 'group'
var: 'n_cells'
uns: 'tree'
['HSC', 'MEP', 'GMP', 'MDP', 'Pro-B', 'Large Pre-B', 'Small Pre-B', 'Trans B', 'Naive B', 'Memory B', 'Plasma', 'pDC', 'CD4 Tnaive', 'CD4 Tcm', 'CD4 Treg', 'CD8 Tnaive', 'CD8 Tem', 'CD8 CTL', 'CD8 NKT', 'MAIT', 'gdT', 'CD56bright NK', 'CD56dim NK', 'CD14 Mono', 'CD16 Mono', 'cDC2', 'Promyelocyte', 'MMP8 Pre-Neut', 'MMP9 Pre-Neut', 'S100A12 Neut', 'SLC8A1 Neut', 'ISG15 Neut', 'Early Ery', 'Late Ery']
╭╴HSC
╭─┼╴MDP
│ ├╴MEP
│ ╰╴GMP
│ ╭╴Memory B
│ ├╴Naive B
├─┼╴Small Pre-B
│ ├╴Large Pre-B
│ ├╴Pro-B
│ ╰╴Trans B
│ ╭╴CD4 Tnaive
│ ├╴CD4 Tcm
│ ├╴gdT
│ ├╴CD8 NKT
├─┼╴CD8 Tem
│ ├╴CD8 CTL
─┤ ├╴MAIT
│ ├╴CD8 Tnaive
│ ╰╴CD4 Treg
│ ╭╴CD14 Mono
├─┼╴CD16 Mono
│ ╰╴cDC2
│ ╭╴MMP9 Pre-Neut
├─┼╴Promyelocyte
│ ╰╴MMP8 Pre-Neut
├─┬╴Early Ery
│ ╰╴Late Ery
├─┬╴CD56bright NK
│ ╰╴CD56dim NK
├─╌╴pDC
├─╌╴Plasma
│ ╭╴S100A12 Neut
╰─┼╴SLC8A1 Neut
╰╴ISG15 Neut
In [8]:
fig, ax = plt.subplots(figsize=(8, 9))
plot_coda_tree(tasc_data["coda"].uns["tree"], ax=ax)
fig.savefig(RES / "03_tasccoda_tree.pdf", bbox_inches="tight")
plt.show()
ax, tasc_box_stats = plot_coda_boxplots(
tasc_data, modality_key="coda", cell_kind="sub", figsize=(18, 6)
)
tasc_box_stats.to_csv(RES / "03_tasccoda_box_wilcox.csv", index=False)
plt.savefig(RES / "03_tasccoda_box.pdf", bbox_inches="tight")
plt.show()
tasc_box_stats
Out[8]:
| cell | n_ctrl | n_aci | median_ctrl | median_aci | p | p_adj | |
|---|---|---|---|---|---|---|---|
| 0 | HSC | 6 | 5 | 0.010360 | 0.002078 | 0.004329 | 0.007747 |
| 1 | MEP | 6 | 5 | 0.006436 | 0.000556 | 0.004329 | 0.007747 |
| 2 | GMP | 6 | 5 | 0.007718 | 0.038173 | 0.004329 | 0.007747 |
| 3 | MDP | 6 | 5 | 0.030616 | 0.011799 | 0.008658 | 0.011775 |
| 4 | Pro-B | 6 | 5 | 0.010138 | 0.000831 | 0.004329 | 0.007747 |
| 5 | Large Pre-B | 6 | 5 | 0.014560 | 0.001156 | 0.004329 | 0.007747 |
| 6 | Small Pre-B | 6 | 5 | 0.026845 | 0.001669 | 0.004329 | 0.007747 |
| 7 | Trans B | 6 | 5 | 0.013884 | 0.001527 | 0.017316 | 0.019625 |
| 8 | Naive B | 6 | 5 | 0.038385 | 0.004850 | 0.004329 | 0.007747 |
| 9 | Memory B | 6 | 5 | 0.023412 | 0.005112 | 0.004329 | 0.007747 |
| 10 | Plasma | 6 | 5 | 0.021307 | 0.001527 | 0.008658 | 0.011775 |
| 11 | pDC | 6 | 5 | 0.015436 | 0.002356 | 0.004329 | 0.007747 |
| 12 | CD4 Tnaive | 6 | 5 | 0.050911 | 0.005001 | 0.008658 | 0.011775 |
| 13 | CD4 Tcm | 6 | 5 | 0.036778 | 0.014225 | 0.004329 | 0.007747 |
| 14 | CD4 Treg | 6 | 5 | 0.008849 | 0.002622 | 0.008658 | 0.011775 |
| 15 | CD8 Tnaive | 6 | 5 | 0.024592 | 0.002208 | 0.008658 | 0.011775 |
| 16 | CD8 Tem | 6 | 5 | 0.064117 | 0.015130 | 0.004329 | 0.007747 |
| 17 | CD8 CTL | 6 | 5 | 0.048812 | 0.018115 | 0.125541 | 0.133387 |
| 18 | CD8 NKT | 6 | 5 | 0.032628 | 0.014714 | 0.246753 | 0.246753 |
| 19 | MAIT | 6 | 5 | 0.008898 | 0.004554 | 0.017316 | 0.019625 |
| 20 | gdT | 6 | 5 | 0.012588 | 0.003187 | 0.008658 | 0.011775 |
| 21 | CD56bright NK | 6 | 5 | 0.011407 | 0.001104 | 0.017316 | 0.019625 |
| 22 | CD56dim NK | 6 | 5 | 0.036080 | 0.014337 | 0.246753 | 0.246753 |
| 23 | CD14 Mono | 6 | 5 | 0.081585 | 0.041223 | 0.017316 | 0.019625 |
| 24 | CD16 Mono | 6 | 5 | 0.006582 | 0.002078 | 0.004329 | 0.007747 |
| 25 | cDC2 | 6 | 5 | 0.008561 | 0.001656 | 0.004329 | 0.007747 |
| 26 | Promyelocyte | 6 | 5 | 0.006574 | 0.047250 | 0.004329 | 0.007747 |
| 27 | MMP8 Pre-Neut | 6 | 5 | 0.007998 | 0.130525 | 0.004329 | 0.007747 |
| 28 | MMP9 Pre-Neut | 6 | 5 | 0.033103 | 0.128034 | 0.004329 | 0.007747 |
| 29 | S100A12 Neut | 6 | 5 | 0.013892 | 0.193920 | 0.004329 | 0.007747 |
| 30 | SLC8A1 Neut | 6 | 5 | 0.009561 | 0.212148 | 0.004329 | 0.007747 |
| 31 | ISG15 Neut | 6 | 5 | 0.005018 | 0.016699 | 0.013497 | 0.017650 |
| 32 | Early Ery | 6 | 5 | 0.115489 | 0.003193 | 0.004329 | 0.007747 |
| 33 | Late Ery | 6 | 5 | 0.067271 | 0.000514 | 0.053027 | 0.058159 |
In [9]:
tasc_data = tasccoda.prepare(
tasc_data,
modality_key="coda",
tree_key="tree",
formula="C(group, Treatment('CTRL'))",
reference_cell_type="automatic",
pen_args={"phi": 0},
)
# pertpy default: num_warmup=1000, num_samples=10000, rng_key=0
tasccoda.run_nuts(tasc_data, modality_key="coda")
tasccoda.summary(tasc_data, modality_key="coda")
sample: 100%|██████████| 11000/11000 [00:27<00:00, 405.39it/s, 63 steps of size 7.42e-02. acc. prob=0.93]
Compositional Analysis summary ┌────────────────────────────────────────────┬────────────────────────────────────────────────────────────────────┐ │ Name │ Value │ ├────────────────────────────────────────────┼────────────────────────────────────────────────────────────────────┤ │ Data │ Data: 11 samples, 34 cell types │ │ Reference cell type │ CD16 Mono │ │ Formula │ C(group, Treatment('CTRL')) │ └────────────────────────────────────────────┴────────────────────────────────────────────────────────────────────┘
┌─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐ │ Intercepts │ ├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤ │ Final Parameter Expected Sample │ │ annot_sub │ │ HSC 0.001 91.615 │ │ MEP -0.606 49.929 │ │ GMP 0.216 113.590 │ │ MDP 0.828 209.473 │ │ Pro-B 0.003 91.799 │ │ Large Pre-B 0.172 108.701 │ │ Small Pre-B 0.504 151.502 │ │ Trans B 0.070 98.160 │ │ Naive B 0.831 210.103 │ │ Memory B 0.670 178.859 │ │ Plasma 0.017 93.093 │ │ pDC 0.087 99.843 │ │ CD4 Tnaive 0.866 217.587 │ │ CD4 Tcm 1.060 264.171 │ │ CD4 Treg -0.127 80.608 │ │ CD8 Tnaive 0.259 118.581 │ │ CD8 Tem 1.406 373.380 │ │ CD8 CTL 1.068 266.293 │ │ CD8 NKT 0.816 206.975 │ │ MAIT 0.020 93.372 │ │ gdT 0.116 102.781 │ │ CD56bright NK -0.304 67.532 │ │ CD56dim NK 0.884 221.539 │ │ CD14 Mono 1.754 528.793 │ │ CD16 Mono -0.390 61.967 │ │ cDC2 -0.279 69.241 │ │ Promyelocyte -0.550 52.805 │ │ MMP8 Pre-Neut 0.373 132.900 │ │ MMP9 Pre-Neut 0.470 146.437 │ │ S100A12 Neut 0.165 107.942 │ │ SLC8A1 Neut 0.107 101.860 │ │ ISG15 Neut -1.766 15.652 │ │ Early Ery 1.899 611.305 │ │ Late Ery 0.898 224.662 │ └─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
┌─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐ │ Effects │ ├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤ │ Effect Expected Sample log2-fold change │ │ Covariate Cell Type │ │ C(group, Treatment('CTRL'))T.ACI HSC 0.000 36.844 -1.314 │ │ MEP 0.000 20.080 -1.314 │ │ GMP 1.221 154.888 0.447 │ │ MDP 0.000 84.243 -1.314 │ │ Pro-B -0.598 20.302 -2.177 │ │ Large Pre-B -0.598 24.040 -2.177 │ │ Small Pre-B -0.598 33.505 -2.177 │ │ Trans B -0.598 21.708 -2.177 │ │ Naive B -0.598 46.465 -2.177 │ │ Memory B -0.598 39.555 -2.177 │ │ Plasma 0.000 37.439 -1.314 │ │ pDC 0.000 40.153 -1.314 │ │ CD4 Tnaive 0.000 87.506 -1.314 │ │ CD4 Tcm 0.000 106.240 -1.314 │ │ CD4 Treg 0.000 32.418 -1.314 │ │ CD8 Tnaive 0.000 47.689 -1.314 │ │ CD8 Tem 0.000 150.160 -1.314 │ │ CD8 CTL 0.000 107.094 -1.314 │ │ CD8 NKT 0.000 83.238 -1.314 │ │ MAIT 0.000 37.551 -1.314 │ │ gdT 0.000 41.335 -1.314 │ │ CD56bright NK 0.000 27.159 -1.314 │ │ CD56dim NK 0.000 89.095 -1.314 │ │ CD14 Mono 0.000 212.662 -1.314 │ │ CD16 Mono 0.000 24.921 -1.314 │ │ cDC2 0.000 27.846 -1.314 │ │ Promyelocyte 2.513 262.094 2.311 │ │ MMP8 Pre-Neut 2.513 659.646 2.311 │ │ MMP9 Pre-Neut 2.513 726.838 2.311 │ │ S100A12 Neut 3.195 1059.660 3.295 │ │ SLC8A1 Neut 3.195 999.948 3.295 │ │ ISG15 Neut 3.195 153.654 3.295 │ │ Early Ery -1.612 49.043 -3.640 │ │ Late Ery -1.612 18.024 -3.640 │ └─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
┌─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐ │ Nodes │ ├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤ │ Covariate=C(group, Treatment('CTRL'))[T.ACI]_node │ ├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤ │ Final Parameter Is credible │ │ Node │ │ annot_main_HSPC 0.00 False │ │ annot_main_B -0.60 True │ │ annot_main_T 0.00 False │ │ annot_main_NK 0.00 False │ │ annot_main_Mono/DC 0.00 False │ │ annot_main_Pre-Neut 2.51 True │ │ annot_main_Neut 3.19 True │ │ annot_main_Ery -1.61 True │ │ Plasma 0.00 False │ │ pDC 0.00 False │ │ HSC 0.00 False │ │ MEP 0.00 False │ │ GMP 1.22 True │ │ MDP 0.00 False │ │ Pro-B 0.00 False │ │ Large Pre-B 0.00 False │ │ Small Pre-B 0.00 False │ │ Trans B 0.00 False │ │ Naive B 0.00 False │ │ Memory B 0.00 False │ │ CD4 Tnaive 0.00 False │ │ CD4 Tcm 0.00 False │ │ CD4 Treg 0.00 False │ │ CD8 Tnaive 0.00 False │ │ CD8 Tem 0.00 False │ │ CD8 CTL 0.00 False │ │ CD8 NKT 0.00 False │ │ MAIT 0.00 False │ │ gdT 0.00 False │ │ CD56bright NK 0.00 False │ │ CD56dim NK 0.00 False │ │ CD14 Mono 0.00 False │ │ CD16 Mono 0.00 False │ │ cDC2 0.00 False │ │ Promyelocyte 0.00 False │ │ MMP8 Pre-Neut 0.00 False │ │ MMP9 Pre-Neut 0.00 False │ │ S100A12 Neut 0.00 False │ │ SLC8A1 Neut 0.00 False │ │ ISG15 Neut 0.00 False │ │ Early Ery 0.00 False │ │ Late Ery 0.00 False │ └─────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
In [10]:
cov = tasc_data["coda"].uns["scCODA_params"]["covariate_names"][0]
node_df = tasc_data["coda"].uns["scCODA_params"]["node_df"]
node_key = cov + "_node"
if node_key in node_df.index.get_level_values(0):
node_effs = node_df.loc[(node_key,)].copy()
node_effs.index = node_effs.index.get_level_values("Node")
node_vals = node_effs["Final Parameter"]
else:
node_vals = None
print("node effects not found for", node_key, node_df.index[:5])
fig, ax = plt.subplots(figsize=(8, 9))
plot_coda_tree(tasc_data["coda"].uns["tree"], ax=ax, node_value=node_vals)
fig.savefig(RES / "03_tasccoda_effects_tree.pdf", bbox_inches="tight")
plt.show()
cred_t = tasccoda.credible_effects(tasc_data, modality_key="coda")
cred_t.to_csv(RES / "03_tasccoda_credible.csv")
cred_t
Out[10]:
Covariate Node
C(group, Treatment('CTRL'))[T.ACI]_node annot_main_HSPC False
annot_main_B True
annot_main_T False
annot_main_NK False
annot_main_Mono/DC False
annot_main_Pre-Neut True
annot_main_Neut True
annot_main_Ery True
Plasma False
pDC False
HSC False
MEP False
GMP True
MDP False
Pro-B False
Large Pre-B False
Small Pre-B False
Trans B False
Naive B False
Memory B False
CD4 Tnaive False
CD4 Tcm False
CD4 Treg False
CD8 Tnaive False
CD8 Tem False
CD8 CTL False
CD8 NKT False
MAIT False
gdT False
CD56bright NK False
CD56dim NK False
CD14 Mono False
CD16 Mono False
cDC2 False
Promyelocyte False
MMP8 Pre-Neut False
MMP9 Pre-Neut False
S100A12 Neut False
SLC8A1 Neut False
ISG15 Neut False
Early Ery False
Late Ery False
Name: Final Parameter, dtype: bool
3. Session info¶
软件版本，便于复现。
In [11]:
print(session_info())
pandas 2.3.3 anndata 0.12.19 jax 0.10.2 matplotlib 3.11.2 numpy 2.4.6 seaborn 0.13.2 scipy 1.17.1 statsmodels 0.15.0 numpyro 0.22.0 pertpy 1.0.3 mudata 0.3.10 ---- ---- Python 3.11.16 (main, Sep 2 2026, 23:39:52) [GCC 15.3.0] OS Linux-6.11.0-26-generic-x86_64-with-glibc2.39 Updated 2026-09-26 23:48