Import a saved Seurat object

# Work with local file paths.
from pathlib import Path

# Keep each import in a fresh directory.
from tempfile import TemporaryDirectory

# Download the example RDS file.
from huggingface_hub import download_bucket_files

# Display inspection results as tables.
import pandas as pd

# Inspect and import Seurat objects.
import cytearc

# Hide routine progress messages.
cytearc.configure_output(level="WARNING", progress=False)

# Show full import notices in the tables below.
pd.set_option("display.max_colwidth", None)

# Choose where to save the downloaded file.
data_directory = Path("cytearc_datasets")

# Create the download directory if needed.
data_directory.mkdir(exist_ok=True)

# Name the local RDS file.
rds_path = data_directory / "seurat_v4_1_3_pbmc_mye.rds"

# Download the public PBMC example.
download_bucket_files(
    "Nygen/cytebase",
    files=[("cytearc_tests/seurat_v4_1_3_pbmc_mye.rds", rds_path)],
    token=False,
    raise_on_missing_files=True,
)

# Keep the converted stores in a temporary directory.
work = TemporaryDirectory()

# Use this directory for each import below.
output_directory = Path(work.name)

# Inspect the assays, metadata and saved coordinates.
inspection = cytearc.inspect_seurat(rds_path)

# Check that the downloaded file matches the published example.
assert inspection.sourceDigest == "f84adf523a78aeb6e6681cf09e06a2a2fcd4e3fe857fdd89b17e90a1782fac3d"

# Show the active assay and cell count.
print("Active assay:", inspection.activeAssay, "| Cells:", inspection.nCells)

# List the assays and their feature-by-cell dimensions.
pd.DataFrame(inspection.assays, columns=["name", "dimensions", "importable"])

# List the saved coordinates, their cell-by-dimension shapes and their roles.
pd.DataFrame(inspection.reductions, columns=["name", "dimensions", "role", "assayUsed"])

# Open the RDS to read its metadata column names.
with cytearc.SeuratReader(rds_path) as reader:

    # Keep the names before closing the reader.
    metadata_names = reader.cellMetadata.columnNames

# Show the columns available for import.
metadata_names

# Choose the output store.
pbmc_store = output_directory / "pbmc.zarr"

# Select the counts, metadata and saved coordinates.
with cytearc.SeuratReader(
    rds_path,
    assays=["RNA"],
    reductions=["pca", "umap"],
    cell_metadata=["orig.ident", "nCount_RNA", "nFeature_RNA", "percent.mt"],
    reduction_roles={"pca": "coordinates", "umap": "umap"},
) as reader:

    # Write the store with RNA as its default assay.
    imported = cytearc.SeuratToZarr(
        reader, zarr_loc=str(pbmc_store), default_assay="RNA",
        nthreads=4, mem_budget="4G",
    ).dump()

# Show the imported assays and default.
print("Imported assays:", imported.assayNames, "| Default:", imported.defaultAssay)

# Show what the import left out.
pd.DataFrame(imported.notices, columns=["code", "objectPath", "message"])

# Open the imported store without filtering its cells.
ds = cytearc.DataStore(str(pbmc_store), min_features_per_cell=-1, nthreads=4)

# Inspect the imported cell and feature counts.
ds

# Plot the saved Seurat UMAP using the saved active identities.
ds.plots.embedding(
    layout=imported.reductionArtifacts["umap"],
    color_by=imported.activeIdentity,
)

# Inspect the small file before importing it.
partial_inspection = cytearc.inspect_seurat(partial_path)

# Show its active assay.
print("Active assay:", partial_inspection.activeAssay)

# Show the role inferred for the coordinates called layout.
print("Inferred role:", partial_inspection.reductions[0].role)

# Show which metadata could not be read.
pd.DataFrame(
    partial_inspection.cellMetadata.notices, columns=["code", "objectPath", "message"]
)

# Set the custom coordinates to a display role and keep all usable data.
with cytearc.SeuratReader(
    partial_path, reduction_roles={"layout": "umap"}
) as reader:

    # Import the data and let the writer choose the default assay.
    partial_import = cytearc.SeuratToZarr(
        reader, zarr_loc=str(output_directory / "partial.zarr"),
        nthreads=4, mem_budget="1G",
    ).dump()

# Show which assays were imported and which became the default.
print("Imported assays:", partial_import.assayNames, "| Default:", partial_import.defaultAssay)

# Show the role assigned to the custom layout.
print("Layout kind:", partial_import.reductionArtifacts["layout"].kind)

# Read the skipped-column and SCT notices.
pd.DataFrame(partial_import.notices, columns=["code", "objectPath", "message"])

# Select RNA counts, two metadata columns and the display layout.
with cytearc.SeuratReader(
    partial_path,
    assays=["RNA"],
    cell_metadata=["sample", "nCount_RNA"],
    reductions=["layout"],
    reduction_roles={"layout": "umap"},
) as reader:

    # Write a store with RNA as its explicit default.
    raw_import = cytearc.SeuratToZarr(
        reader, zarr_loc=str(output_directory / "raw.zarr"),
        default_assay="RNA", nthreads=4, mem_budget="1G",
    ).dump()

# Show the selected assay, its default and any import notices.
{
    "assays": raw_import.assayNames,
    "default": raw_import.defaultAssay,
    "notices": [item.code for item in raw_import.notices],
}