ContextQMD
Back to Claude Scientific Skills

Whole slide and cohort runs

skills/deepspot-m/references/whole_slide.md

2.63.05.3 KB
Original Source

Whole slide and cohort runs

DeepSpot-M predicts per tile. A slide-scale virtual spatial transcriptomics map is a tiling step, a batched prediction loop, and an assembly step that puts the values back on the slide grid.

1. Pick the level that gives about 20x

Tiles must be 224x224 at roughly 20x, near 0.5 microns per pixel. Read the resolution off the slide rather than assuming level 0 is 20x, since many scanners write level 0 at 40x:

python
import openslide

slide = openslide.open_slide("slide.svs")
mpp_x = float(slide.properties.get(openslide.PROPERTY_NAME_MPP_X))
downsamples = slide.level_downsamples

level = min(
    range(slide.level_count),
    key=lambda i: abs(mpp_x * downsamples[i] - 0.5),
)

Tile at that level. A slide already scanned at 20x gives level 0; a 40x slide usually gives level 1.

2. Extract a tile grid

Use the histolab skill for tiling. A grid tiler at 224x224 with a tissue check covers the section and skips background:

python
from histolab.slide import Slide
from histolab.tiler import GridTiler

slide = Slide("slide.svs", processed_path="tiles/")

tiler = GridTiler(
    tile_size=(224, 224),
    level=level,
    check_tissue=True,
    tissue_percent=80.0,
    pixel_overlap=0,
)
tiler.extract(slide)

Keep each tile's coordinates. ScoreTiler.extract(slide, report_path="tiles_report.csv") writes a CSV with tile_name,x_coord,y_coord,level,..., which is the least fragile way to carry them. See the histolab skill for tissue masks, filters and the other tilers.

3. Predict in batches

Load the model once, then stream tiles through it. Reloading per batch redownloads nothing but rebuilds the model each time, which dominates the runtime of a slide:

python
from pathlib import Path

import torch
from PIL import Image
from deepspotm import DeepSpotM

TILE_PX = 224

def require_tile(tile):
    if tile.size != (TILE_PX, TILE_PX):
        raise ValueError(
            f"DeepSpot-M expects a {TILE_PX}x{TILE_PX} tile at about 20x "
            f"(~0.5 microns per pixel); got {tile.size[0]}x{tile.size[1]}."
        )
    return tile.convert("RGB")

def batched(items, size):
    for start in range(0, len(items), size):
        yield items[start : start + size]

device = "cuda" if torch.cuda.is_available() else "cpu"
model, image_processor = DeepSpotM.from_pretrained(
    "ratschlab/DeepSpotM", source="scgpt", device=device
)

genes = ["EPCAM", "CD3D", "PTPRC", "MKI67"]
tile_paths = sorted(Path("tiles/").glob("*.png"))

chunks = []
for paths in batched(tile_paths, 32):
    tiles = [require_tile(Image.open(p)) for p in paths]
    batch = torch.stack([image_processor(t) for t in tiles]).to(device)
    chunks.append(model.predict_genes(batch, genes).cpu())

expression = torch.cat(chunks).numpy()  # tiles by genes, log1p-CPM

Batch sizing: start at 32 tiles on a GPU and 8 on CPU. Memory grows with both the batch size and the number of genes requested in one call, so lower one when the other is large. Ask for the full gene list in each call rather than looping gene by gene, since the tile tokens are computed once per batch and reused across gene queries.

4. Assemble the slide map

Pair the matrix with the tile coordinates and the gene list. AnnData is the natural container, and it is what spatial analysis tools read:

python
import anndata as ad
import numpy as np
import pandas as pd

report = pd.read_csv("tiles_report.csv")
coords = report[["x_coord", "y_coord"]].to_numpy(dtype=float)

adata = ad.AnnData(
    X=expression,
    obs=pd.DataFrame({"tile_name": report["tile_name"]}).set_index("tile_name"),
    var=pd.DataFrame(index=pd.Index(genes, name="gene")),
)
adata.obsm["spatial"] = coords
adata.uns["deepspotm"] = {
    "source": "scgpt",
    "units": "log1p-CPM",
    "tile_px": 224,
    "level": int(level),
}
adata.write_h5ad("slide.h5ad")

Recording source, units and level in uns keeps the run readable later, and makes it obvious when two slides were produced under different settings.

5. Plot a gene

python
import matplotlib.pyplot as plt

values = adata[:, "EPCAM"].X.ravel()
plt.scatter(coords[:, 0], -coords[:, 1], c=values, s=6, cmap="viridis")
plt.gca().set_aspect("equal")
plt.colorbar(label="EPCAM (log1p-CPM)")

Negating the y coordinate puts the map in slide orientation, since slide coordinates grow downward.

6. Cohort scale

For many slides, run one slide per process and write one .h5ad per slide rather than holding a cohort in memory:

python
for svs in sorted(Path("cohort/").glob("*.svs")):
    out = Path("out") / f"{svs.stem}.h5ad"
    if out.exists():
        continue          # resume without recomputing finished slides
    run_slide(svs, out)   # steps 1 to 4 above

Points worth fixing across a cohort:

  • One source for every slide, so values stay comparable.
  • One gene list, stored in a file and read by every run.
  • The same target resolution, chosen per slide from its own metadata.
  • A skip-if-exists guard, so an interrupted cohort resumes where it stopped.

Concatenate afterwards with ad.concat(slides, label="slide_id") when a cohort-level matrix is needed. This is the shape of the run that produced the TCGA atlas of 28,664 slides across 32 cancer types.

Primary sources