BCR Embedding with scibex — Stephenson 2021 COVID-19 Dataset¶
This tutorial walks through a complete scibex workflow on the Stephenson 2021 COVID-19 PBMC dataset (5 000 BCR-containing cells):
- Load a pre-processed scirpy
MuData - Embed heavy- and light-chain CDR3 sequences with scibex
- Visualise the embedding space
- Train a logistic-regression classifier to predict patient outcome (Home vs Death)
scibex wraps the R Ibex package and stores results directly in AnnData.obsm, following scverse conventions. It operates on the airr modality of a scirpy MuData — no manual sequence extraction needed.
Prerequisites¶
pip install scibex scirpy scanpy muon
scibex also requires R with the Ibex package installed. Run the setup cell below. It installs everything and is a no-op if already installed.
import os
os.chdir("../..")
import scibex as ib
# Install R dependencies (fast no-op if already installed)
ib.install_r_deps()
import warnings
import pandas as pd
import scanpy as sc
import scirpy as ir
warnings.filterwarnings("ignore")
sc.settings.verbosity = 1
sc.settings.set_figure_params(dpi=80, frameon=False)
1. Load the dataset¶
The dataset is a MuData with two modalities:
gex— gene expression (pre-processed; PCA/UMAP already computed)airr— BCR contigs in AIRR format, withobsm["chain_indices"]already populated by scirpy
mdata = ir.datasets.stephenson2021_5k()
mdata
mdata["gex"].obs[["Status_on_day_collection_summary", "Outcome"]].value_counts().head(10)
2. Embed BCR sequences with scibex¶
ib.tl.ibex reads CDR3 sequences directly from the scirpy chain_indices in the airr modality and writes an [N, D] embedding array to mdata["airr"].obsm[key_added]. Cells that lack a usable chain receive rows filled with zeros (the default fill_value=0.0; pass fill_value=float("nan") to get NaN instead). Pass verbose=True to see how many cells were affected.
Available methods:
method |
encoder_model |
Description |
|---|---|---|
"encoder" |
"CNN" |
CDR3 convolutional autoencoder |
"encoder" |
"VAE" |
CDR3 variational autoencoder |
"encoder" |
"CNN.EXP" |
CDR1+2+3 CNN (needs cdr1_aa/cdr2_aa AIRR fields) |
"encoder" |
"VAE.EXP" |
CDR1+2+3 VAE |
"geometric" |
— | Fast rule-based geometric transform (BLOSUM62-based) |
Available encoder_input options: "atchleyFactors", "kideraFactors", "crucianiProperties", "MSWHIM", "tScales", "OHE"
For EXP models, strategy controls how cells with partial CDR missingness (CDR1 or CDR2 absent but CDR3 present) are handled:
"lenient"(default) — substitute missing CDRs with"NA"(Ibex's R convention) and embed anyway"strict"— any missing CDR → zero row
# Heavy-chain CDR1+2+3 embedding (CNN.EXP with Atchley factors)
# verbose=True reports how many cells lacked complete CDR1/CDR2/CDR3 (filled with zeros)
ib.tl.ibex(
mdata,
chain="Heavy",
method="encoder",
encoder_model="CNN.EXP",
encoder_input="atchleyFactors",
key_added="X_ibex_heavy",
)
# Light-chain CDR1+2+3 embedding
ib.tl.ibex(
mdata,
chain="Light",
method="encoder",
encoder_model="CNN.EXP",
encoder_input="atchleyFactors",
key_added="X_ibex_light",
)
print("Heavy embedding:", mdata["airr"].obsm["X_ibex_heavy"].shape)
print("Light embedding:", mdata["airr"].obsm["X_ibex_light"].shape)
3. Visualise the heavy-chain embedding¶
We compute a UMAP from the heavy-chain Ibex embedding and colour cells by COVID-19 outcome and clinical status.
Note:
gexandairrare stored with different cell orderings in thisMuData. We reindex the embeddings to thegexcell ordering via barcode before handing off to scanpy.
gex = mdata["gex"]
# Align heavy embeddings to gex cell ordering via barcode index
heavy_df = pd.DataFrame(
mdata["airr"].obsm["X_ibex_heavy"],
index=mdata["airr"].obs_names,
)
gex.obsm["X_ibex_heavy"] = heavy_df.reindex(gex.obs_names).values
sc.pp.neighbors(gex, use_rep="X_ibex_heavy", n_neighbors=30, metric="cosine")
sc.tl.umap(gex)
sc.pl.umap(gex, color=["Outcome", "Status_on_day_collection_summary"], ncols=2, wspace=0.4)
4. Predict COVID-19 outcome from paired BCR embeddings¶
We concatenate the heavy- and light-chain embeddings into a single feature matrix, then train a regularised logistic regression. LogisticRegressionCV selects the regularisation strength via leave-one-out CV on the training data. We evaluate with stratified 5-fold cross-validation and report balanced accuracy (appropriate given the ~18:1 Home/Death imbalance).
from sklearn.linear_model import LogisticRegressionCV
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
# Align light embeddings to gex cell ordering
light_df = pd.DataFrame(
mdata["airr"].obsm["X_ibex_light"],
index=mdata["airr"].obs_names,
)
outcome = mdata["gex"].obs["Outcome"]
X = pd.concat(
[heavy_df.reindex(outcome.index), light_df.reindex(outcome.index)],
axis=1,
).values
y = outcome.values
# Keep only cells with a known outcome label
labeled = y != "unknown"
X_clf, y_clf = X[labeled], y[labeled]
print(f"Labeled cells: {labeled.sum()} (Home={(y_clf == 'Home').sum()}, Death={(y_clf == 'Death').sum()})")
pipe = make_pipeline(
StandardScaler(),
LogisticRegressionCV(
max_iter=500,
cv=3,
class_weight="balanced",
random_state=42,
),
)
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)
scores = cross_val_score(pipe, X_clf, y_clf, cv=cv, scoring="balanced_accuracy")
print(f"5-fold CV balanced accuracy: {scores.mean():.3f} ± {scores.std():.3f}")
print(f"Per-fold: {scores.round(3)}")
Session info¶
sc.logging.print_header()