pRoloc integration: for grassp users#

pRoloc is an R/Bioconductor framework for spatial proteomics. It and grassp have overlapping functionality and hold nearly the same data model, but they share no file format, making it difficult to move data or compare results between the two.

grassp and pRoloc exchange objects as h5ad in both directions, to enable seamless integration. For example, you could preprocess and plot in Python, hand the object to R for pRoloc’s classifiers, and read it back to Python. What comes back is as close to what you sent as the two data models allow, including the results of pRolocs classification. Practically this opens the possibility to use functionality that is not shared between the two frameworks, such as BANDLE in pRoloc for differential localization detection, or the independent diffusion annotation approach in grassp (see the diffusion tutorial).

Note

This tutorial is for someone who works primarily in Python and wants to use a specific pRoloc functionality. If you work primarily in R, read the companion tutorial instead: pRoloc integration: for pRoloc users. It stays in R throughout — the ~100 datasets on the grassp portal are h5ad files, so grassp_as_msnset() reads them directly with no Python involved — and vignette("grasspio") is the reference for that side.

What this tutorial does

  1. Load a dataset and pick a marker set.

  2. Write it out as a plain h5ad.

  3. Run pRoloc’s SVM and k-NN in R (the R code is shown; its output ships with this tutorial so the notebook builds without R).

  4. Read the results back with anndata.read_h5ad and plot them.

  5. Merge them onto a session you already have instead.

Installation#

The Python side needs nothing extra — the exchange format is h5ad, which anndata already writes. The R side is a companion package, grasspio, that lives in the same repository:

install.packages(c("remotes", "BiocManager"))
BiocManager::install(c("pRoloc", "rhdf5"))
remotes::install_github("czbiohub-sf/grassp", subdir = "r/grasspio")

grasspio uses the scVerse package anndataR, to convert AnnData files into in-memory R objects under the hood. For now that arrives from a fork, pinned in grasspio’s own Remotes: field and installed automatically: matrix-valued fData columns are written as obsm data frames, and the fix that indexes such a frame by the parent’s obs_names is newer than the last anndataR release.

import warnings

import anndata
import matplotlib.pyplot as plt
import numpy as np
import scanpy as sc

import grassp as gr
/Users/mfrank/code/grassp/.venv-dev/lib/python3.12/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

Loading the data#

Currie_2024_AC16_Control is a LOPIT-DC experiment on AC16 cardiomyocytes: 2538 proteins across 10 fractions, with each protein’s profile already normalised to sum to 1. That normalisation matters here — pRoloc’s distance-based methods and its plots assume it.

adata = gr.ds.load_dataset("Currie_2024_AC16_Control")
adata
AnnData object with n_obs × n_vars = 2538 × 10
    obs: 'protein_name', 'gene_symbol', 'author_annotation', 'author_markers', 'bandle_probability', 'dl_candidate', 'author_markers_data10', 'marker_lilley', 'marker_christopher', 'marker_geladaki', 'marker_itzhak', 'marker_villaneuva', 'marker_hein2025', 'marker_hein2025_gt', 'n_detected', 'n_samples_by_intensity', 'mean_intensity', 'log1p_mean_intensity', 'pct_dropout_by_intensity', 'total_intensity', 'log1p_total_intensity', 'unsupervised_annotation', 'harmonized_annotation_propagated', 'harmonized_annotation_propagated_probability'
    var: 'development_stage', 'tissue', 'sex', 'organism', 'disease', 'cell_line', 'cell_type', 'fraction', 'fraction_order', 'subcellular_enrichment', 'n_proteins_by_intensity', 'log1p_n_proteins_by_intensity', 'pct_dropout_by_intensity', 'n_merged_samples', 'enriched_vs', 'enrichment_strategy', 'mass_spectrometry_method', 'cell_type_ontology_term_id', 'perturbation', 'disease_ontology_term_id', 'organism_ontology_term_id', 'sex_ontology_term_id', 'development_stage_ontology_term_id', 'tissue_ontology_term_id', 'cell_line_ontology_term_id', 'tissue_type', 'tissue_general_ontology_term_id', 'tissue_general', 'organ_ontology_term_id', 'organ', 'system_ontology_term_id', 'system'
    uns: 'author_annotation_colors', 'author_markers_data10_colors', 'dl_candidate_colors', 'harmonized_annotation_propagated_colors', 'marker_christopher_colors', 'marker_geladaki_colors', 'marker_hein2025_colors', 'marker_hein2025_gt_colors', 'marker_itzhak_colors', 'marker_lilley_colors', 'marker_villaneuva_colors', 'neighbors', 'pca', 'publication_authors', 'publication_date', 'publication_doi', 'publication_journal', 'publication_title', 'schema_reference', 'schema_version', 'serves_as_reference', 'title', 'umap', 'unsupervised_annotation_colors'
    obsm: 'X_pca', 'X_umap', 'X_umap3D', 'harmonized_annotation_propagated_one_hot_labels', 'harmonized_annotation_propagated_probabilities'
    varm: 'PCs'
    layers: 'log_intensities', 'original_intensities', 'pvals'
    obsp: 'connectivities', 'distances'

The dataset ships several published marker sets. We’ll use Lilley’s, the one pRoloc::pRolocmarkers() provides, and copy it to obs["markers"] — the column name every pRoloc function defaults to.

In grassp, unlabelled proteins are NaN; in pRoloc they are the literal string "unknown", and pRoloc genuinely needs it (markerMSnSet and unknownMSnSet fail outright on NA). That is the one real semantic difference between the frameworks, and it is handled entirely on the R side: grassp_as_msnset(nan_to_unknown = TRUE) fills the sentinel in as the MSnSet is built, and grassp_write_msnset(unknown_to_na = TRUE) strips it on the way back. So the sentinel never touches a file, and you never have to think about it here.

Both directions convert every text column, not a nominated one. That matters: fcol is a per-call argument in pRoloc, so one MSnSet can carry markers, markers.orig and pd.markers at once and point different functions at different ones, exactly as an AnnData can.

adata.obs["markers"] = adata.obs["marker_lilley"]

print(
    f"{adata.obs['markers'].notna().sum()} markers across "
    f"{adata.obs['markers'].nunique()} compartments\n"
)
print(adata.obs["markers"].value_counts().to_string())

sc.pl.umap(adata, color="markers")
403 markers across 12 compartments

markers
Mitochondrion         65
ER                    56
Nucleus               45
60S Ribosome          44
Cytosol               42
PM                    37
Proteasome            31
40S Ribosome          30
Actin Cytoskeleton    28
Lysosome              13
Peroxisome             8
Golgi                  4
../../_images/be80479ab5ee3259c3f845d182a7993f6246604930eb392e6d7ff65aba262996.png

Exporting for pRoloc#

There is no exporter to call: write the object as h5ad and hand over the file. The checks that used to live in a wrapper now live where the constraints actually are — grassp_as_msnset() refuses duplicate or blank protein IDs (R rownames must be unique) and says something about missing values and profiles that don’t sum to 1, because those are pRoloc’s requirements rather than h5ad’s.

Everything crosses, because the expected next step is to read the object back. An MSnSet’s assayData is a Biobase environment holding any number of equal-dimension matrices, so .layers become extra assay elements; pData is the same class as fData, so .varm becomes matrix-valued pData columns; and anything in .uns that has no MSnSet slot at all rides along on experimentData@other. All of them subset correctly in R alongside exprs, which is what makes them safe rather than merely stored.

The one thing that genuinely cannot cross is .obspeSet has no pairwise slot, and pRoloc has no graph structure to map onto. It is also the slot least worth carrying: a neighbour graph is derived from .X, so gr.pp.neighbors rebuilds it in one call.

adata.write_h5ad("experiment.h5ad")

print("wrote experiment.h5ad")
print("stays behind (eSet has no pairwise slot):", list(adata.obsp))
wrote experiment.h5ad
stays behind (eSet has no pairwise slot): ['connectivities', 'distances']

To send less — a big object with a slow classifier is a reason to — subset before writing, with the ordinary anndata tools: adata[:, keep].write_h5ad(...), or del sub.obsm["X_pca"] on a copy.

Two things worth knowing. Only .X becomes exprs(), the matrix pRoloc’s functions operate on, so if the one you want is a layer, make it .X first — this dataset ships log_intensities, original_intensities and pvals alongside .X, and they cross as extra assayData elements either way. And .uns is written wholesale, so an entry h5py cannot serialise will abort the write; if that happens, drop the offending key.

Over in R#

Reading the artifact takes no arguments — an h5ad object already describes itself:

library(grasspio)
library(pRoloc)

x <- grassp_as_msnset("experiment.h5ad")
x
#> MSnSet (storageMode: lockedEnvironment)
#> assayData: 2538 features, 10 samples

getMarkerClasses(x, fcol = "markers")
#>  [1] "40S Ribosome" "60S Ribosome" "Actin Cytoskeleton" "Cytosol" "ER" ...

From there it is an ordinary MSnSet, so the whole pRoloc workflow applies:

## Support vector machine. In real work, get the hyperparameters from
## svmOptimisation(x, fcol = "markers", times = 100, xval = 5,
##                 class.weights = classWeights(x, fcol = "markers"))
x <- svmClassification(x, fcol = "markers", sigma = 0.1, cost = 16, scores = "all")

## `scores = "all"` stores the per-class matrix but NOT the scalar winning score, while
## orgQuants()/getPredictions() look for <fcol>.scores -- so derive it from the matrix.
fData(x)$svm.scores <- apply(fData(x)$svm.all.scores, 1, max)
ts <- orgQuants(x, fcol = "svm", scol = "svm.scores", t = 0.75)
ts[is.na(ts)] <- Inf   # a class with too few markers has no quantile
x <- getPredictions(x, fcol = "svm", scol = "svm.scores", t = ts)

## k nearest neighbours
x <- knnClassification(x, fcol = "markers", k = 5, scores = "prediction")

## Several classes here have fewer markers than there are fractions; minMarkers()
## demotes those to "unknown" in a new `markers10` column.
x <- minMarkers(x, n = 10, fcol = "markers")

grassp_write_msnset(x, "proloc_tutorial_results.h5ad", overwrite = TRUE)

That exact script ships next to this notebook as proloc_tutorial.R. Its output is published alongside the portal datasets, and the next cell fetches it — so the rest of the notebook runs whether or not you have R, and the results below are genuine pRoloc 1.51.1 output rather than a simulation.

from pathlib import Path
from urllib.request import urlretrieve

# The R side's output, published next to the portal datasets rather than committed to the repo
# (it is 4 MB). Regenerate it with proloc_tutorial.R and re-upload if the R workflow changes.
RESULTS = Path("proloc_tutorial_results.h5ad")
if not RESULTS.exists():
    urlretrieve(
        "https://public.czbiohub.org/proteinxlocation/internal/proloc_tutorial_results.h5ad",
        RESULTS,
    )
print(f"{RESULTS}{RESULTS.stat().st_size / 1024**2:.1f} MiB")
proloc_tutorial_results.h5ad — 4.1 MiB

Reading the results back#

anndata.read_h5ad, and that is the whole of it. There is no importer, because there is nothing to import: the file the R side wrote is an AnnData, so the object prints its own inventory.

annotated = anndata.read_h5ad(RESULTS)
annotated
AnnData object with n_obs × n_vars = 2538 × 10
    obs: 'protein_name', 'gene_symbol', 'author_annotation', 'author_markers', 'bandle_probability', 'dl_candidate', 'author_markers_data10', 'marker_lilley', 'marker_christopher', 'marker_geladaki', 'marker_itzhak', 'marker_villaneuva', 'marker_hein2025', 'marker_hein2025_gt', 'n_detected', 'n_samples_by_intensity', 'mean_intensity', 'log1p_mean_intensity', 'pct_dropout_by_intensity', 'total_intensity', 'log1p_total_intensity', 'unsupervised_annotation', 'harmonized_annotation_propagated', 'harmonized_annotation_propagated_probability', 'markers', 'svm', 'svm.scores', 'svm.pred', 'knn', 'knn.scores', 'markers10'
    var: 'development_stage', 'tissue', 'sex', 'organism', 'disease', 'cell_line', 'cell_type', 'fraction', 'fraction_order', 'subcellular_enrichment', 'n_proteins_by_intensity', 'log1p_n_proteins_by_intensity', 'pct_dropout_by_intensity', 'n_merged_samples', 'enriched_vs', 'enrichment_strategy', 'mass_spectrometry_method', 'cell_type_ontology_term_id', 'perturbation', 'disease_ontology_term_id', 'organism_ontology_term_id', 'sex_ontology_term_id', 'development_stage_ontology_term_id', 'tissue_ontology_term_id', 'cell_line_ontology_term_id', 'tissue_type', 'tissue_general_ontology_term_id', 'tissue_general', 'organ_ontology_term_id', 'organ', 'system_ontology_term_id', 'system'
    uns: 'author_annotation_colors', 'author_markers_data10_colors', 'dl_candidate_colors', 'harmonized_annotation_propagated_colors', 'marker_christopher_colors', 'marker_geladaki_colors', 'marker_hein2025_colors', 'marker_hein2025_gt_colors', 'marker_itzhak_colors', 'marker_lilley_colors', 'marker_villaneuva_colors', 'neighbors', 'pca', 'processing', 'publication_authors', 'publication_date', 'publication_doi', 'publication_journal', 'publication_title', 'schema_reference', 'schema_version', 'serves_as_reference', 'title', 'umap', 'unsupervised_annotation_colors'
    obsm: 'X_pca', 'X_umap', 'X_umap3D', 'harmonized_annotation_propagated_one_hot_labels', 'harmonized_annotation_propagated_probabilities', 'svm.all.scores'
    varm: 'PCs'
    layers: 'log_intensities', 'original_intensities', 'pvals'

There is no mapping table and no list of supported methods in any of this. A scalar fData column becomes an .obs column and a matrix-valued one becomes an .obsm entry; names and values cross unchanged. So svm, knn, markers10 and svm.pred all arrive without the bridge knowing anything about them, and a classifier added to pRoloc tomorrow needs no change here.

Note that the label columns are already Categoricals with NaN, and svm.scores is already a float. That is not Python inspecting the data — the R side converted "unknown" back to NA on the way out and wrote those columns as R factors, which h5ad represents natively. anndataR maps types faithfully in both directions (numeric, integer and logical; Categorical ↔ factor, with the level order and the ordered flag intact), and the one gap — it has no nullable-string encoding, so a character NA would arrive as the literal string "NA" — is exactly what the write-as-a-factor step exists to avoid.

One name is worth knowing about, though nothing in the bridge treats it specially: svm.all.scores is also what ksvmClassification writes its scores under, so in pRoloc the column is genuinely ambiguous — read it as “whichever SVM you ran”.

The point of the exercise is what that object contains — not just pRoloc’s new columns, but everything you sent, down to the nested uns entries an MSnSet has no slot for.

One thing a wrapper used to do for you: compartment colours. gr.pp.set_sensible_compartment_colors decides for itself which of the new columns look like compartment annotations, so point it at them once and every grassp plot below picks the palette up.

gr.pp.set_sensible_compartment_colors(annotated)

print("new from pRoloc:", [c for c in annotated.obs.columns if c not in adata.obs.columns])
print()
for slot in ("obs", "var"):
    lost = set(getattr(adata, slot).columns) - set(getattr(annotated, slot).columns)
    print(f"{slot + ' columns lost:':22s} {sorted(lost) or 'none'}")
for slot in ("obsm", "varm", "layers", "uns", "obsp"):
    lost = set(getattr(adata, slot)) - set(getattr(annotated, slot))
    print(f"{slot + ' lost:':22s} {sorted(lost) or 'none'}")
print()
print("uns['neighbors'] params:", annotated.uns["neighbors"]["params"])
new from pRoloc: ['svm', 'svm.scores', 'svm.pred', 'knn', 'knn.scores', 'markers10']

obs columns lost:      none
var columns lost:      none
obsm lost:             none
varm lost:             none
layers lost:           none
uns lost:              none
obsp lost:             ['connectivities', 'distances']

uns['neighbors'] params: {'method': 'umap', 'metric': 'euclidean', 'n_neighbors': 20, 'random_state': 0, 'use_rep': 'X'}

pRoloc’s own column names come through verbatim, so a workflow you know from R reads the same way here:

annotated.obs[["markers", "svm", "svm.scores", "svm.pred", "knn", "markers10"]].head(8)
markers svm svm.scores svm.pred knn markers10
A0AVT1 Cytosol Cytosol 1.000000 Cytosol Cytosol Cytosol
A1L0T0 NaN ER 0.759524 NaN ER NaN
A2RRP1 NaN ER 0.883894 ER ER NaN
A5PLN9 NaN Actin Cytoskeleton 0.246245 NaN Actin Cytoskeleton NaN
A5YKK6 NaN ER 0.330107 NaN ER NaN
A6NDG6 Cytosol Cytosol 1.000000 Cytosol Cytosol Cytosol
A6NHT5 NaN 60S Ribosome 0.512753 60S Ribosome 60S Ribosome NaN
E9PRG8 NaN 40S Ribosome 0.371695 NaN 60S Ribosome NaN

"unknown" arrived as NaN, which is what grassp’s annotators expect — they select markers with .notna(). Note how getPredictions thresholding moved proteins into NaN in svm.pred relative to the unthresholded svm:

print("unlabelled proteins")
for column, note in [
    ("markers", "(input)"),
    ("svm", "(unthresholded)"),
    ("svm.pred", "(t = 0.75)"),
    ("markers10", "(minMarkers)"),
]:
    print(f"  {column:10s}{note:18s}{annotated.obs[column].isna().sum():>5}")
unlabelled proteins
  markers   (input)            2135
  svm       (unthresholded)       0
  svm.pred  (t = 0.75)         1597
  markers10 (minMarkers)       2147

The probability matrix#

The most interesting part of the mapping. pRoloc stores per-protein × per-compartment scores inside a single fData column — a matrix nested in a data frame. That is exactly what .obsm is for.

It comes back as a DataFrame, so the class names are attached to the numbers rather than recorded in a side table — with pRoloc’s own <class>.svm.scores decoration left intact, because renaming them would be a change we have no need to make. An embedding, which has no class names to carry, stays a plain array: X_umap goes out and comes back as one.

scores = annotated.obsm["svm.all.scores"]

print("obsm['svm.all.scores']:", type(scores).__name__, scores.shape)
print("classes (pRoloc's own names):")
print("   ", list(scores.columns)[:3], "...")
print()
print("rows sum to 1:      ", np.allclose(scores.to_numpy().sum(axis=1), 1))
print("X_umap stays an array:", type(annotated.obsm["X_umap"]).__name__)
obsm['svm.all.scores']: DataFrame (2538, 12)
classes (pRoloc's own names):
    ['Cytosol.svm.scores', 'Peroxisome.svm.scores', 'PM.svm.scores'] ...

rows sum to 1:       True
X_umap stays an array: ndarray

grassp’s plotting takes the column names as arguments, so pRoloc’s own names go straight in. The UMAP came back with the object, so there is nothing to recompute — here the SVM call with point transparency scaled by its confidence:

gr.pl.umap_prob(annotated, color="svm", color_prob="svm.scores")
../../_images/beb8c01d9c890653f8b5455ed336c938cc75bb30abd02adb3022409360362e00.png
<Axes: title={'center': 'svm (opacity ~ svm.scores)'}, xlabel='UMAP1', ylabel='UMAP2'>

And the raw class-probability matrix, as a per-compartment heatmap over the marker proteins:

probs = annotated.obsm["svm.all.scores"].rename(columns=lambda c: c.replace(".svm.scores", ""))

labelled = annotated.obs["markers"].notna()
ordered = probs[labelled.values].assign(_c=annotated.obs.loc[labelled, "markers"].astype(str))
ordered = ordered.sort_values("_c").drop(columns="_c")

fig, ax = plt.subplots(figsize=(6, 7))
im = ax.imshow(ordered.to_numpy(), aspect="auto", cmap="magma", vmin=0, vmax=1)
ax.set_xticks(range(ordered.shape[1]))
ax.set_xticklabels(ordered.columns, rotation=90)
ax.set_yticks([])
ax.set_ylabel(f"{int(labelled.sum())} marker proteins, grouped by compartment")
ax.set_title("pRoloc SVM class probabilities")
fig.colorbar(im, ax=ax, label="probability", shrink=0.6)
plt.tight_layout()
../../_images/22b80c0740f71b56085229f9a5d97dd442c939964b98f07793ce3eff9a4df839.png

Do the two frameworks agree?#

grassp has its own RBF-SVM annotator. Running it on the same data, the same markers and the same hyperparameters is a useful check that nothing was scrambled in transit. pRoloc’s sigma is the RBF width scikit-learn calls gamma, and its cost is C, so we can match them exactly — though the two wrap different libraries (e1071 vs scikit-learn), so expect close agreement rather than identity.

One subtlety that is easy to get wrong: gr.tl.svm_annotation defaults to min_probability=0.5, so it returns NaN for calls it is not confident about, while pRoloc’s unthresholded svm column is populated everywhere. Comparing the two naively therefore compares only the proteins grassp was sure about. Both numbers are worth seeing.

for key, min_probability in [("svm_confident", 0.5), ("svm_all", 0.0)]:
    gr.tl.svm_annotation(
        annotated,
        gt_col="markers",
        C=16,
        gamma=0.1,
        min_probability=min_probability,
        key_added=key,
    )
    both = annotated.obs[["svm", key]].dropna()
    agreement = (both["svm"].astype(str) == both[key].astype(str)).mean()
    print(
        f"min_probability={min_probability}: agree on {agreement:.1%} "
        f"of the {len(both)} proteins grassp labelled"
    )
min_probability=0.5: agree on 89.3% of the 1155 proteins grassp labelled
min_probability=0.0: agree on 68.0% of the 2538 proteins grassp labelled
/Users/mfrank/code/grassp/.venv-dev/lib/python3.12/site-packages/sklearn/svm/_base.py:239: FutureWarning: The `probability` parameter was deprecated in 1.9 and will be removed in version 1.11. Use `CalibratedClassifierCV(SVC(), ensemble=False)` instead of `SVC(probability=True)`
  warnings.warn(
/Users/mfrank/code/grassp/.venv-dev/lib/python3.12/site-packages/sklearn/svm/_base.py:239: FutureWarning: The `probability` parameter was deprecated in 1.9 and will be removed in version 1.11. Use `CalibratedClassifierCV(SVC(), ensemble=False)` instead of `SVC(probability=True)`
  warnings.warn(

The other way to come back#

If you still have the object in your session and would rather not rebuild it, copy over the columns you actually want. You know which ones those are, so there is no function for it:

adata.obs["svm.pred"] = res.obs["svm.pred"].reindex(adata.obs_names)
adata.obsm["svm.all.scores"] = res.obsm["svm.all.scores"].reindex(adata.obs_names)

The reindex is the load-bearing part. Real pRoloc workflows shrink objects routinely (filterNA, markerMSnSet, unknownMSnSet, plain x[i, ]), so the artifact often covers fewer proteins than the object you are merging onto — and a silently half-populated annotation reads exactly like a result. Aligning on obs_names makes the shortfall NaN, which is what every grassp annotator treats as unlabelled. The cell below trims the results object first, standing in for a filterNA() on the R side, so you can see that happen.

Note that assignment replaces, under whatever name you give it — so to keep pRoloc’s svm and grassp’s own side by side, rename one as you copy it.

# Standing in for a filterNA() in R: the results now cover 10 fewer proteins than `adata`.
res = annotated[10:].copy()

adata.obs["svm.pred"] = res.obs["svm.pred"].reindex(adata.obs_names)
adata.obsm["svm.all.scores"] = res.obsm["svm.all.scores"].reindex(adata.obs_names)

print(f"results cover {res.n_obs} of the {adata.n_obs} proteins in the session")
print("proteins the artifact never saw:", int(adata.obs["svm.pred"].isna().sum()), "NaN")
print("dtype preserved:", adata.obs["svm.pred"].dtype)
print("class names preserved:", list(adata.obsm["svm.all.scores"].columns)[:2], "...")
print()
print("untouched by the merge: X, layers, obsp, var — nothing was rebuilt")
print("  layers:", sorted(adata.layers))
print("  obsp:  ", sorted(adata.obsp))
results cover 2528 of the 2538 proteins in the session
proteins the artifact never saw: 1602 NaN
dtype preserved: category
class names preserved: ['Cytosol.svm.scores', 'Peroxisome.svm.scores'] ...

untouched by the merge: X, layers, obsp, var — nothing was rebuilt
  layers: ['log_intensities', 'original_intensities', 'pvals']
  obsp:   ['connectivities', 'distances']

Going the other way#

The bridge is symmetric, and useful R-first too — pRoloc and MSnbase ship no exporter of their own, so grassp_write_msnset is a way to hand a classified MSnSet to anything that reads h5ad. That direction is the subject of the companion vignette, vignette("grasspio") — which also notes that grassp portal datasets are plain h5ad, so R users can read them directly.

The one piece that belongs on this side: for the many pRolocdata datasets distributed as R .rda files, read_prolocdata() reads them without any R at all, via the pure-Python rdata parser (pip install grassp[proloc]).

Next: BANDLE#

The reason to build this bridge is the methods with no Python equivalent. BANDLE infers differential localisation between two conditions, and expects a list of MSnSets per condition with matching features and equal channel counts.

The same write_h5ad / read_h5ad pair carries that workflow — one file per condition. Whatever bandle.* columns it writes into fData arrive structurally, exactly like svm and knn above, with no special case anywhere in the bridge — which is the whole point of not having a list of method names.

Limitations worth knowing#

  • .obsp/.varp are the only slots that cannot cross: eSet has no pairwise slot, and pRoloc’s one neighbour-ish representation (nndist()) writes flat positional indices, which are silently wrong after any subsetting. A graph is derived from .X, so gr.pp.neighbors rebuilds it.

  • Multi-localisation must cross as pRoloc’s binary Markers matrix — one-hot the labels into .obsm yourself (as a DataFrame, so the class names come along) and it crosses like any other matrix column. pRoloc::mrkMatToVec collapses any protein with more than one label to "unknown", so don’t round-trip through the vector encoding.

  • A class name containing / cannot be a DataFrame column, because HDF5 reads it as a path separator. pRoloc produces them — hyperLOPIT’s classes include “Endoplasmic reticulum/Golgi apparatus” — so those matrices arrive as plain arrays with their names in uns[f"{key}_categories"] instead, and the R side says so when it happens.

  • Optimisation and MCMC side objects (GenRegRes, MAPParams, bandleParams) live outside fData, so only the summaries pRoloc writes into fData come across. Save those with saveRDS() if you need them — re-running MCMC to recover one is expensive.

  • Nothing is renamed, so a column called svm in R is a column called svm here. If you want grassp’s own naming (svm_probabilities and friends), rename it yourself — the bridge deliberately does not guess.

  • pRoloc::tagmMapTrain fails on this dataset with “x is not a symmetric matrix”. It is tempting to blame the 10-fraction design, and a co-linearity message appears just beforehand, but that message is a red herring — it fails the same way on full-rank data. The cause is upstream: LaplacesDemon::is.symmetric.matrix() tests a scatter matrix built as t(A) %*% A with exact ==, and BLAS returns it asymmetric at around 3e-16. grassp’s own port, gr.tl.tagm_map_train / gr.tl.tagm_map_predict, writes the same tagm.map.allocation, .probability and .outlier columns and is a practical alternative.