Spatial Proteomics: CODEX, IMC & Phenocycler Analysis

conda create -n spatprot python=3.10 -y
conda activate spatprot

# Core analysis
pip install scimap anndata scanpy napari

# Cell segmentation
pip install deepcell tensorflow

# Image processing
pip install tifffile imageio scikit-image

# Verify
python -c "import scimap; print(scimap.__version__)"

import tifffile
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors

# Load a multiplexed TIFF (channels × height × width)
# CODEX/Phenocycler output is typically a multi-channel TIFF
image = tifffile.imread("codex_tissue.tiff")
print(f"Image shape: {image.shape}")  # (C, H, W) — C channels

# Load channel names (antibody panel)
panel = pd.read_csv("panel.csv")
# panel columns: channel_number, marker_name, target
print(f"Markers: {panel['marker_name'].tolist()}")

# Example panel for a tumour immunology study:
# DAPI, CD3, CD4, CD8, CD20, CD56, CD68, CD163, CK, PD-L1,
# Ki67, FOXP3, CD45, CD11b, HLA-DR, VISTA, TIM3, LAG3...

# Extract specific channels by marker name
def get_channel(image, panel, marker):
    idx = panel[panel["marker_name"] == marker].index[0]
    return image[idx]

dapi = get_channel(image, panel, "DAPI")
cd8  = get_channel(image, panel, "CD8")
ck   = get_channel(image, panel, "CK")   # cytokeratin (tumour cells)

# Visualise composite (RGB overlay)
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
axes[0].imshow(dapi, cmap="Blues",  vmax=np.percentile(dapi, 99))
axes[0].set_title("DAPI (nuclei)", fontweight="bold")
axes[1].imshow(cd8,  cmap="Greens", vmax=np.percentile(cd8, 99))
axes[1].set_title("CD8 (cytotoxic T cells)", fontweight="bold")
axes[2].imshow(ck,   cmap="Reds",   vmax=np.percentile(ck, 99))
axes[2].set_title("CK (tumour cells)", fontweight="bold")
for ax in axes: ax.axis("off")
plt.tight_layout()
plt.savefig("channel_overview.pdf", dpi=300, bbox_inches="tight")
plt.show()

from deepcell.applications import Mesmer
import numpy as np

# Initialise Mesmer (downloads weights on first run ~500 MB)
app = Mesmer()

# Prepare input: Mesmer needs [batch, H, W, 2]
# Channel 0: nuclear marker (DAPI)
# Channel 1: membrane/cytoplasm marker (e.g. CD45 or CK)
nuclear = dapi
membrane = get_channel(image, panel, "CD45")

# Stack and add batch + channel dimensions
X = np.stack([nuclear, membrane], axis=-1)
X = np.expand_dims(X, axis=0)  # add batch dim -> (1, H, W, 2)

# Normalise
X = X.astype(np.float32)
X[..., 0] = X[..., 0] / np.percentile(X[..., 0], 99.9)
X[..., 1] = X[..., 1] / np.percentile(X[..., 1], 99.9)
X = np.clip(X, 0, 1)

# Run segmentation
# Returns (1, H, W, 2): channel 0 = whole cell, channel 1 = nuclear
segmentation_predictions = app.predict(
    X,
    image_mpp=0.5,           # microns per pixel (check your microscope settings)
    compartment="whole-cell" # or "nuclear" or "both"
)

cell_mask = segmentation_predictions[0, ..., 0]
print(f"Cells detected: {cell_mask.max()}")

# Visualise segmentation
fig, axes = plt.subplots(1, 2, figsize=(12, 6))
axes[0].imshow(dapi, cmap="Blues", vmax=np.percentile(dapi, 99))
axes[0].set_title("DAPI", fontweight="bold")
axes[1].imshow(cell_mask, cmap="tab20")
axes[1].set_title(f"Cell segmentation ({cell_mask.max()} cells)", fontweight="bold")
for ax in axes: ax.axis("off")
plt.savefig("segmentation.pdf", dpi=300, bbox_inches="tight")
plt.show()

from skimage.measure import regionprops_table

# Extract mean intensity per cell per channel
props = regionprops_table(
    cell_mask.astype(int),
    intensity_image=image.transpose(1, 2, 0),   # H × W × C
    properties=["label", "centroid", "area",
                 "intensity_mean"]
)

df = pd.DataFrame(props)

# Rename intensity columns to marker names
intensity_cols = [c for c in df.columns if "intensity_mean" in c]
rename_map = {intensity_cols[i]: panel["marker_name"].iloc[i]
              for i in range(min(len(intensity_cols), len(panel)))}
df = df.rename(columns=rename_map)
df = df.rename(columns={"centroid-0": "y", "centroid-1": "x"})

print(f"Cell table shape: {df.shape}")
print(df[["label","x","y","area","CD3","CD8","CK","CD68"]].head())

# Save cell table
df.to_csv("cell_table.csv", index=False)

import scimap as sm
import anndata as ad
import scanpy as sc

# Build AnnData from cell table
marker_cols = panel["marker_name"].tolist()
X = df[marker_cols].values

adata = ad.AnnData(X=X)
adata.var_names = marker_cols
adata.obs = df[["label","x","y","area"]].reset_index(drop=True)
adata.obs.index = adata.obs["label"].astype(str)
adata.obsm["spatial"] = df[["x","y"]].values

# Arcsinh normalisation (standard for mass cytometry data)
sm.pp.rescale(adata, gate=0.5)

# Automated phenotyping using marker thresholds
# Define phenotyping hierarchy
phenotype_rules = {
    "T cell":      {"CD3": "pos"},
    "CD4 T cell":  {"CD3": "pos", "CD4": "pos"},
    "CD8 T cell":  {"CD3": "pos", "CD8": "pos"},
    "Treg":        {"CD3": "pos", "FOXP3": "pos"},
    "B cell":      {"CD20": "pos"},
    "NK cell":     {"CD56": "pos", "CD3": "neg"},
    "Macrophage":  {"CD68": "pos"},
    "M2 Macro":    {"CD68": "pos", "CD163": "pos"},
    "Tumour":      {"CK": "pos"},
    "Endothelial": {"CD31": "pos"},
}

sm.tl.phenotype_cells(
    adata,
    phenotype=phenotype_rules,
    gate=0.5,
    label="phenotype"
)

# Check phenotype distribution
counts = adata.obs["phenotype"].value_counts()
print("Cell type counts:")
print(counts.to_string())

# Plot pie chart of cell composition
fig, ax = plt.subplots(figsize=(8, 8))
wedges, texts, autotexts = ax.pie(
    counts.values,
    labels=counts.index,
    autopct="%1.1f%%",
    startangle=90,
    pctdistance=0.82
)
ax.set_title("Cell type composition", fontsize=14, fontweight="bold")
plt.tight_layout()
plt.savefig("cell_composition_pie.pdf", dpi=300, bbox_inches="tight")
plt.show()

# Build spatial KNN graph
sm.tl.spatial_interaction(
    adata,
    phenotype="phenotype",
    method="knn",
    knn=10,
    label="spatial_interaction"
)

# Visualise interaction heatmap
sm.pl.spatial_interaction(
    adata,
    spatial_interaction="spatial_interaction",
    figsize=(10, 8),
    title="Cell-cell spatial interactions"
)
plt.savefig("spatial_interactions.pdf", dpi=300, bbox_inches="tight")

# Find spatial communities (tissue niches)
sm.tl.spatial_cluster(
    adata,
    df_name="spatial_interaction",
    cluster_method="kmeans",
    k=5,
    label="tissue_niche"
)

# Map niches back to tissue
sm.pl.image_viewer(
    adata,
    imageid="tissue",
    overlay="tissue_niche",
    point_size=5
)

import scipy.spatial as spatial
import seaborn as sns

# For each cell type pair: is the distance shorter than expected by chance?
# (permutation-based spatial enrichment test)

coords = adata.obsm["spatial"]
phenotypes = adata.obs["phenotype"].values

cell_types = list(counts.index)
n_types = len(cell_types)

# Build KDTree for efficient distance queries
tree = spatial.cKDTree(coords)

# Count observed neighbours within radius r
r = 50  # micrometers (adjust for your image scale)

interaction_matrix = np.zeros((n_types, n_types))

for i, ct_source in enumerate(cell_types):
    source_idx = np.where(phenotypes == ct_source)[0]
    source_coords = coords[source_idx]

    for j, ct_target in enumerate(cell_types):
        target_idx = np.where(phenotypes == ct_target)[0]
        target_coords = coords[target_idx]

        # Count target cells within radius r of each source cell
        target_tree = spatial.cKDTree(target_coords)
        neighbours = tree.query_ball_point(source_coords, r)
        avg_neighbours = np.mean([len(n) for n in neighbours])
        interaction_matrix[i, j] = avg_neighbours

# Plot interaction heatmap
fig, ax = plt.subplots(figsize=(10, 8))
sns.heatmap(
    pd.DataFrame(interaction_matrix, index=cell_types, columns=cell_types),
    cmap="Reds", annot=True, fmt=".1f",
    linewidths=0.5, ax=ax
)
ax.set_title(f"Mean neighbours within {r}μm radius", fontsize=13, fontweight="bold")
ax.set_xlabel("Target cell type")
ax.set_ylabel("Source cell type")
plt.tight_layout()
plt.savefig("interaction_proximity.pdf", dpi=300, bbox_inches="tight")
plt.show()

# Plot cell phenotypes overlaid on tissue
fig, axes = plt.subplots(1, 2, figsize=(16, 8))

# Left: DAPI background with phenotype overlay
axes[0].imshow(dapi, cmap="Greys_r", vmax=np.percentile(dapi, 99))

# Define colours per phenotype
colours = {
    "CD8 T cell":  "#2196F3",
    "CD4 T cell":  "#4CAF50",
    "Treg":        "#8BC34A",
    "B cell":      "#9C27B0",
    "NK cell":     "#FF9800",
    "Macrophage":  "#F44336",
    "M2 Macro":    "#E91E63",
    "Tumour":      "#FFEB3B",
    "Endothelial": "#00BCD4",
    "Other":       "#9E9E9E"
}

for ct, colour in colours.items():
    mask = adata.obs["phenotype"] == ct
    if mask.sum() > 0:
        xy = adata.obsm["spatial"][mask.values]
        axes[0].scatter(xy[:, 0], xy[:, 1], s=2, c=colour, label=ct, alpha=0.7)

axes[0].legend(loc="upper right", markerscale=5, fontsize=7,
               framealpha=0.8, ncol=2)
axes[0].set_title("Cell phenotypes on tissue", fontweight="bold")
axes[0].axis("off")

# Right: same plot coloured by tissue niche
niche_colours = plt.cm.Set2(np.linspace(0, 1, 5))
for i, niche in enumerate(adata.obs["tissue_niche"].unique()):
    mask = adata.obs["tissue_niche"] == niche
    xy = adata.obsm["spatial"][mask.values]
    axes[1].scatter(xy[:, 0], xy[:, 1], s=2,
                    c=[niche_colours[i]], label=f"Niche {niche}", alpha=0.7)
axes[1].set_title("Tissue niches", fontweight="bold")
axes[1].legend(markerscale=5, fontsize=8)
axes[1].axis("off")

plt.tight_layout()
plt.savefig("spatial_phenotype_map.pdf", dpi=300, bbox_inches="tight")
plt.show()