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
Load a dataset and pick a marker set.
Write it out as a plain h5ad.
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).
Read the results back with
anndata.read_h5adand plot them.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
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 .obsp — eSet 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")
<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()
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/.varpare the only slots that cannot cross:eSethas 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, sogr.pp.neighborsrebuilds it.Multi-localisation must cross as pRoloc’s binary
Markersmatrix — one-hot the labels into.obsmyourself (as a DataFrame, so the class names come along) and it crosses like any other matrix column.pRoloc::mrkMatToVeccollapses 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 inuns[f"{key}_categories"]instead, and the R side says so when it happens.Optimisation and MCMC side objects (
GenRegRes,MAPParams,bandleParams) live outsidefData, so only the summaries pRoloc writes intofDatacome across. Save those withsaveRDS()if you need them — re-running MCMC to recover one is expensive.Nothing is renamed, so a column called
svmin R is a column calledsvmhere. If you want grassp’s own naming (svm_probabilitiesand friends), rename it yourself — the bridge deliberately does not guess.pRoloc::tagmMapTrainfails 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 ast(A) %*% Awith 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 sametagm.map.allocation,.probabilityand.outliercolumns and is a practical alternative.