pRoloc integration: for pRoloc users#
The grassp data portal publishes around a hundred
curated subcellular fractionation datasets, uniformly processed, each with several published
marker sets already mapped onto it. They are stored as .h5ad, which sounds like a Python
format — but grasspio reads it natively in R.
So this tutorial needs no Python at all. You download a file, get an
MSnSet, and from there it is ordinary pRoloc. It also goes the other way: the last section
takes a pRolocdata object and writes it out as h5ad.
Note
This tutorial is for someone who works primarily in R and wants to use grassp’s data or a grassp method. If you work primarily in Python, read the companion tutorial instead: pRoloc integration: for grassp users. It is the same bridge seen from the other end — you preprocess and plot in Python, hand the object to R for a classifier pRoloc has and grassp does not, and read the results back.
Requirements#
install.packages(c("remotes", "BiocManager"))
# We need rdf5 to read the on disk h5ad format that grassp uses, and pRoloc/pRolocdata for the processing.
BiocManager::install(c("rhdf5", "pRoloc", "pRolocdata"))
# grasspio is a small R package that grassp provides to read/write h5ad files. It uses scVerse's anndataR package under the hood.
remotes::install_github("czbiohub-sf/grassp", subdir = "r/grasspio")
Note that grasspio requires R >=4.5.
suppressPackageStartupMessages({
library(grasspio)
library(pRoloc)
library(MSnbase)
})
# IRkernel decides figure size from these rather than from knitr chunk options.
options(repr.plot.width = 9, repr.plot.height = 6, width = 100)
# Print results as an R console would, rather than as IRkernel's rich HTML. A named vector
# otherwise renders as a <dl> definition list instead of aligned monospace, and the theme paints
# a light box behind any non-table text/html output in dark mode. Plots are unaffected.
options(jupyter.rich_display = FALSE)
Warning message in fun(libname, pkgname):
“mzR has been built against a different Rcpp version (1.1.1)
than is installed on your system (1.1.2). This might lead to errors
when loading mzR. If you encounter such issues, please send a report,
including the output of sessionInfo() to the Bioc support forum at
https://support.bioconductor.org/. For details see also
https://github.com/sneumann/mzR/wiki/mzR-Rcpp-compiler-linker-issue.”
Getting a dataset#
Browse the portal, pick a dataset, and take its
.h5ad. The files sit at a predictable URL, so you can download one in your browser or fetch it
from R:
https://public.czbiohub.org/proteinxlocation/datasets/<name>.h5ad # processed
https://public.czbiohub.org/proteinxlocation/datasets_raw/<name>.h5ad # pre-enrichment
The datasets/ objects are already preprocessed: enriched, replicates collapsed, and each
protein’s profile already normalised to sum to 1. The datasets_raw/ objects are raw intensities
with replicates preserved — useful, but they need normalising before pRoloc’s distance-based
methods mean anything.
Here we take Currie_2024_AC16_Control, a LOPIT-DC experiment on AC16 cardiomyocytes:
options(timeout = 600) # the default 60 s is not enough for the larger objects
h5ad_path <- file.path(tempdir(), "Currie_2024_AC16_Control.h5ad")
download.file(
"https://public.czbiohub.org/proteinxlocation/datasets/Currie_2024_AC16_Control.h5ad",
h5ad_path,
mode = "wb"
)
round(file.size(h5ad_path) / 1024^2, 1) # MiB
[1] 3.8
x <- grassp_as_msnset(h5ad_path)
x
MSnSet (storageMode: lockedEnvironment)
assayData: 2538 features, 10 samples
element names: exprs, log_intensities, original_intensities, pvals
protocolData: none
phenoData
sampleNames: F1 F2 ... F10 (10 total)
varLabels: development_stage tissue ... PCs (33 total)
varMetadata: labelDescription
featureData
featureNames: A0AVT1 A1L0T0 ... Q9Y6Y8 (2538 total)
fvarLabels: protein_name gene_symbol ... harmonized_annotation_propagated_probabilities
(41 total)
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:
- - - Processing information - - -
Imported from grassp h5ad [Currie_2024_AC16_Control.h5ad]: Wed Aug 12 16:34:32 2026
MSnbase version: 2.36.0
That is the whole conversion. There is nothing special about the file: a portal dataset is a plain
h5ad, and so is anything grassp writes — the objects describe themselves, so there is no version
block to check and no metadata to be missing. The one thing grassp_as_msnset() does change is the
nan_to_unknown default, which rewrites NA markers to the literal "unknown" that pRoloc
requires, because markerMSnSet() fails outright on NA.
Getting oriented#
Everything in the AnnData object has landed somewhere. Per-protein annotations are fData,
per-fraction ones are pData, and the extra matrices are additional assayData elements:
dim(x)
head(fvarLabels(x), 12)
assayDataElementNames(x)
[1] 2538 10
[1] "protein_name" "gene_symbol" "author_annotation" "author_markers"
[5] "bandle_probability" "dl_candidate" "author_markers_data10" "marker_lilley"
[9] "marker_christopher" "marker_geladaki" "marker_itzhak" "marker_villaneuva"
[1] "exprs" "log_intensities" "original_intensities" "pvals"
Some fData columns are themselves matrices — the per-protein × per-compartment arrays that
came from AnnData’s obsm. That is pRoloc’s own idiom for score matrices, so they arrive in the
shape pRoloc expects:
vapply(
fData(x),
function(col) if (is.matrix(col)) paste(dim(col), collapse = " x ") else "",
character(1)
)[c("X_pca", "X_umap", "harmonized_annotation_propagated_probabilities")]
X_pca X_umap
"2538 x 9" "2538 x 2"
harmonized_annotation_propagated_probabilities
"2538 x 16"
The publication metadata came too. An MSnSet has no slot for arbitrary key-value metadata, so
grasspio parks whatever it cannot map on experimentData(x)@other$grassp_uns:
uns <- experimentData(x)@other$grassp_uns
unlist(uns[c("title", "publication_journal", "publication_doi")])
title publication_journal
"Currie_2024_AC16_Control" "Nature Communications"
publication_doi
"https://doi.org/10.1038/s41467-024-46600-5"
Pick a marker column — there is no single one#
The portal maps several published marker sets onto every dataset. pRoloc’s fcol is a per-call
argument, so all of them are usable and you choose per analysis:
marker_cols <- grep("^marker_|^author_markers", fvarLabels(x), value = TRUE)
data.frame(
column = marker_cols,
classes = vapply(marker_cols, function(f) length(getMarkerClasses(x, fcol = f)), integer(1)),
labelled = vapply(marker_cols, function(f) sum(fData(x)[[f]] != "unknown"), integer(1)),
row.names = NULL
)
column classes labelled
1 author_markers 12 334
2 author_markers_data10 12 334
3 marker_lilley 12 403
4 marker_christopher 12 948
5 marker_geladaki 10 349
6 marker_itzhak 12 725
7 marker_villaneuva 12 426
8 marker_hein2025 20 2266
9 marker_hein2025_gt 16 756
We will use marker_lilley, the set pRoloc::pRolocmarkers() also ships:
FCOL <- "marker_lilley"
getMarkerClasses(x, fcol = FCOL)
dim(markerMSnSet(x, fcol = FCOL))
c(missing = sum(is.na(exprs(x))), rows_summing_to_1 = sum(abs(rowSums(exprs(x)) - 1) < 1e-6))
[1] "40S Ribosome" "60S Ribosome" "Actin Cytoskeleton" "Cytosol"
[5] "ER" "Golgi" "Lysosome" "Mitochondrion"
[9] "Nucleus" "Peroxisome" "PM" "Proteasome"
[1] 403 10
missing rows_summing_to_1
0 2538
No preprocessing is needed on an enriched object: filterNA() and normalise() would both be
no-ops here.
One trap before you plot. pRoloc defaults to fcol = "markers", and a portal object has no
such column — so a bare plot2D(x) errors with fcol %in% fvarLabels(object) is not TRUE.
Always name your column, or pass fcol = NULL.
Plotting#
plot2D(x, fcol = FCOL, main = "pRoloc default PCA (scaled)")
addLegend(x, fcol = FCOL, where = "topright", cex = 0.6, ncol = 2)
Note the tenth bar is empty. The profiles are normalised to sum to 1, so the centred matrix is rank 9 — which is also why the portal stores nine PCA components for ten fractions, not ten.
Comparing pRoloc’s PCA with the portal’s own#
The portal object already carries PCA coordinates in fData(x)$X_pca, computed by grassp when
the dataset was built. Recomputing them in pRoloc is a good way to check you are looking at the
same object — but only if you match the convention, and pRoloc’s default does not:
stored <- fData(x)$X_pca
colnames(stored) <- paste0("PC", seq_len(ncol(stored)))
default_pca <- plot2D(x, fcol = FCOL, plot = FALSE)
colnames(default_pca) # pRoloc prints the variance explained into the names
c(PC1 = cor(default_pca[, 1], stored[, 1]), PC2 = cor(default_pca[, 2], stored[, 2]))
[1] "PC1 (44.54%)" "PC2 (27.41%)"
PC1 PC2
0.9198377 0.2361334
PC2 correlates at about 0.24. Nothing is broken: plot2D() passes scale = TRUE to prcomp,
scaling every fraction to unit variance, while grassp centres without scaling. Scaling does not
merely stretch the axes — it changes the eigenvectors and their order, so comparing
component-by-component across the two conventions compares different directions.
Match the convention and they agree:
matched <- plot2D(x, fcol = FCOL, plot = FALSE, methargs = list(center = TRUE, scale = FALSE))
colnames(matched)
vapply(1:2, function(k) cor(matched[, k], stored[, k]), numeric(1))
[1] "PC1 (39.82%)" "PC2 (32.92%)"
[1] 1 -1
The variance-explained figures in the column names are the quickest tell that you have the right convention: 39.8% / 32.9% here versus 44.5% / 27.4% for the scaled default.
cor(PC2) = -1 is not a disagreement either. An eigenvector is only defined up to sign, so some
components come back negated; that is arithmetic, not a discrepancy. Fix it explicitly:
recomputed <- prcomp(exprs(x), center = TRUE, scale. = FALSE)$x[, 1:ncol(stored)]
flip <- sign(diag(cor(recomputed, stored)))
round(flip, 0)
aligned <- sweep(recomputed, 2, flip, `*`)
round(diag(cor(aligned, stored)), 6)
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9
1 -1 1 -1 1 -1 -1 1 1
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9
1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 0.999998 0.999996
Those correlations round to 1, but the coordinates themselves still differ by around 1e-4: the
portal ships the embedding its pipeline produced, not a bit-exact function of the matrix
alongside it. Treat it as a fixed reference embedding — compare with a correlation or a tolerance
of about 1e-3, never with all.equal().
c(max_abs_difference = max(abs(aligned - stored)))
stopifnot(all(abs(diag(cor(recomputed, stored))) > 0.9999))
max_abs_difference
0.000172202
Drawing grassp’s embeddings with pRoloc#
plot2D(method = "none") takes coordinates you already have and styles them like any other
pRoloc plot. That works for the stored PCA and, usefully, for the UMAP — which spares you
compiling a UMAP package just to look at it:
Profiles and marker quality#
plotDist() is most informative when you give it two things: the proteins assigned to a
compartment, and — via markers = — the subset of those that are actually markers for it. The
markers are drawn in colour over the rest in grey, so you can see at a glance whether the
assigned population really follows the marker profile.
The portal makes that easy because it carries both an author annotation and several published
marker sets. One catch: they do not share a vocabulary. marker_lilley says
"Mitochondrion" where the authors’ own column says "MITOCHONDRION", and "Golgi" is "GA".
Nothing in the bridge translates compartment names — that is deliberate, since silently renaming
someone’s labels is worse than a mismatch — so map them yourself:
to_author <- c(Mitochondrion = "MITOCHONDRION", ER = "ER")
par(mfrow = c(1, 2))
for (cl in names(to_author)) {
assigned <- featureNames(x)[fData(x)$author_annotation == to_author[[cl]]]
markers <- featureNames(x)[fData(x)[[FCOL]] == cl]
# plotDist() subsets the object it is *given* by `markers`, so a marker the authors placed
# somewhere else is out of bounds. Pass the intersection; the shortfall is itself worth
# seeing, since it counts markers the two annotations disagree about.
shown <- intersect(markers, assigned)
plotDist(x[assigned, ], markers = shown, pch = 1,
main = sprintf("%s\n%d assigned, %d of %d markers",
cl, length(assigned), length(shown), length(markers)))
}
par(mfrow = c(1, 1))
options(repr.plot.height = 5)
plotConsProfiles(mrkConsProfiles(x, fcol = FCOL))
options(repr.plot.height = 6)
QSep() quantifies how well separated the marker classes are — the single most useful sanity
check on a fractionation experiment, and something the portal makes unusually easy because you
can ask it of every marker set at once. Higher is better: it is the between-class distance in
units of within-class spread, so a value near 1 means a class is no further from its neighbours
than from itself.
summary(QSep(x, fcol = FCOL))
sapply(marker_cols, function(f) median(summary(QSep(x, fcol = f), verbose = FALSE)))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.8742 1.5359 2.3911 2.5345 3.1159 7.2937
author_markers author_markers_data10 marker_lilley marker_christopher
3.173759 3.173759 2.391133 2.894165
marker_geladaki marker_itzhak marker_villaneuva marker_hein2025
3.066839 1.703440 2.609419 1.634307
marker_hein2025_gt
1.719568
Classification#
svmClassification() needs hyperparameters. A real analysis gets them from svmOptimisation(),
which is a grid search over repeated stratified cross-validation and takes minutes, so it is
shown here but not run:
params <- svmOptimisation(x, fcol = FCOL, times = 100,
class.weights = classWeights(x, fcol = FCOL))
plot(params)
getParams(params)
set.seed(1) # e1071's probability scaling consumes the RNG; without this, runs differ
x <- svmClassification(x, fcol = FCOL, sigma = 0.1, cost = 16, scores = "all")
# `scores = "all"` writes the per-class matrix but NOT the scalar score that orgQuants() and
# getPredictions() look for, so derive it.
fData(x)$svm.scores <- apply(fData(x)$svm.all.scores, 1, max)
dim(fData(x)$svm.all.scores)
[1] "marker_lilley"
[1] 2538 12
The output columns are named svm, svm.scores, svm.all.scores regardless of which fcol you
trained on. The bare [1] "marker_lilley" above the result is a stray print() inside
MLInterfaces, not something you did.
# `mcol` is the second fcol-like trap: orgQuants() and getPredictions() default it to "markers"
# too, and use it to decide which proteins were training data.
ts <- orgQuants(x, fcol = "svm", scol = "svm.scores", mcol = FCOL, t = 0.75, verbose = FALSE)
ts[is.na(ts)] <- Inf # a class with too few markers has no quantile
x <- getPredictions(x, fcol = "svm", scol = "svm.scores", mcol = FCOL, t = ts, verbose = FALSE)
table(fData(x)$svm.pred == "unknown")
FALSE TRUE
941 1597
plot2D(x, fcol = "svm.pred", methargs = list(center = TRUE, scale = FALSE),
main = "SVM predictions, thresholded")
A second classifier is a cheap sanity check, and QSep closes the loop quantitatively — do the
predicted labels separate better than the markers you started from?
set.seed(1)
x <- knnClassification(x, fcol = FCOL, k = 5) # note: knnClassification rejects `verbose`
c(svm_vs_knn_agreement = mean(fData(x)$svm == fData(x)$knn))
c(markers = median(summary(QSep(x, fcol = FCOL), verbose = FALSE)),
svm.pred = median(summary(QSep(x, fcol = "svm.pred"), verbose = FALSE)))
[1] "marker_lilley"
svm_vs_knn_agreement
0.822695
markers svm.pred
2.391133 2.866938
Going the other way: a native pRoloc object out to h5ad#
The bridge is symmetric, and useful even if you never touch Python — neither pRoloc nor MSnbase
ships an h5ad exporter. Any MSnSet will do, including the ones pRolocdata distributes:
suppressPackageStartupMessages(library(pRolocdata))
data(dunkley2006)
dunkley2006
MSnSet (storageMode: lockedEnvironment)
assayData: 689 features, 16 samples
element names: exprs
protocolData: none
phenoData
sampleNames: M1F1A M1F4A ... M2F11B (16 total)
varLabels: membrane.prep fraction replicate
varMetadata: labelDescription
featureData
featureNames: AT1G09210 AT1G21750 ... AT4G39080 (689 total)
fvarLabels: assigned evidence ... markers (8 total)
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
pubMedIds: 16618929
Annotation:
- - - Processing information - - -
Loaded on Thu Jul 16 22:53:08 2015.
Normalised to sum of intensities.
Added markers from 'mrk' marker vector. Thu Jul 16 22:53:08 2015
MSnbase version: 1.17.12
out <- file.path(tempdir(), "dunkley2006.h5ad")
grassp_write_msnset(dunkley2006, out, overwrite = TRUE)
Wrote /var/folders/zk/39l3k6s15hz3z9k65hb14xpr0000gq/T//RtmpuUka6l/dunkley2006.h5ad (689 features x 16 fractions; 0 matrix column(s); 0 extra assay element(s)).
Back in Python:
adata = anndata.read_h5ad("dunkley2006.h5ad")
Scalar fData columns become obs, matrix-valued ones become obsm data frames that carry their
own class names, extra assayData elements become layers, and "unknown" becomes NA on
the way out because that is grassp’s encoding. Reading it back gives the same object:
back <- grassp_as_msnset(out)
c(
same_dim = identical(dim(back), dim(dunkley2006)),
exprs_equal = isTRUE(all.equal(unname(exprs(back)), unname(exprs(dunkley2006)),
tolerance = 1e-8)),
same_columns = identical(sort(fvarLabels(back)), sort(fvarLabels(dunkley2006)))
)
getMarkerClasses(back, fcol = "markers")
689 of 689 profiles do not sum to 1 (observed range 4-4). pRoloc's distance-based methods and its plots assume sum-normalised profiles; see normalise() or gr.pp.normalize_total().
same_dim exprs_equal same_columns
TRUE TRUE TRUE
[1] "ER lumen" "ER membrane" "Golgi" "Mitochondrion" "Plastid" "PM"
[7] "Ribosome" "TGN" "vacuole"
That file is what a grassp user would open with anndata.read_h5ad("dunkley2006.h5ad") — no
grassp function required, because it is an ordinary h5ad. The mapping table, and the few things
that cannot cross, are in vignette("grasspio") and in IO: io.
Session info#
sessionInfo()
R version 4.5.3 (2026-03-11)
Platform: x86_64-apple-darwin13.4.0
Running under: macOS Tahoe 26.6
Matrix products: default
BLAS/LAPACK: /opt/homebrew/Caskroom/miniconda/base/envs/grassp-r64/lib/libopenblasp-r0.3.34.dylib; LAPACK version 3.12.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Los_Angeles
tzcode source: system (macOS)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] pRolocdata_1.48.0 pRoloc_1.51.1 BiocParallel_1.44.0 MLInterfaces_1.90.0
[5] cluster_2.1.8.3 annotate_1.88.0 XML_3.99-0.23 AnnotationDbi_1.72.0
[9] IRanges_2.44.0 MSnbase_2.36.0 ProtGenerics_1.42.0 S4Vectors_0.48.0
[13] mzR_2.44.0 Rcpp_1.1.2 Biobase_2.70.0 BiocGenerics_0.56.0
[17] generics_0.1.4 grasspio_0.1.0
loaded via a namespace (and not attached):
[1] splines_4.5.3 pbdZMQ_0.3-14 filelock_1.0.3
[4] tibble_3.3.1 hardhat_1.4.3 preprocessCore_1.72.0
[7] pROC_1.19.0.1 rpart_4.1.27 lifecycle_1.0.5
[10] httr2_1.3.0 doParallel_1.0.17 globals_0.19.1
[13] lattice_0.22-9 MASS_7.3-66 MultiAssayExperiment_1.36.1
[16] dendextend_1.19.1 magrittr_2.0.5 limma_3.66.0
[19] plotly_4.12.1 otel_0.2.0 reticulate_1.46.0
[22] MsCoreUtils_1.22.1 DBI_1.3.0 RColorBrewer_1.1-3
[25] lubridate_1.9.5 abind_1.4-8 GenomicRanges_1.62.1
[28] purrr_1.2.2 mixtools_2.0.0.1 AnnotationFilter_1.34.0
[31] nnet_7.3-20 ipred_0.9-15 lava_1.9.2
[34] listenv_1.0.0 gdata_3.0.1 parallelly_1.48.0
[37] ncdf4_1.24 codetools_0.2-20 DelayedArray_0.36.0
[40] tidyselect_1.2.1 Spectra_1.20.1 farver_2.1.2
[43] viridis_0.6.5 matrixStats_1.5.0 BiocFileCache_3.0.0
[46] base64enc_0.1-6 Seqinfo_1.0.0 jsonlite_2.0.0
[49] caret_7.0-1 e1071_1.7-17 survival_3.8-9
[52] iterators_1.0.14 foreach_1.5.2 segmented_2.2-1
[55] tools_4.5.3 progress_1.2.3 glue_1.8.1
[58] prodlim_2026.03.11 gridExtra_2.3.1 SparseArray_1.10.8
[61] BiocBaseUtils_1.12.0 xfun_0.60 MatrixGenerics_1.22.0
[64] IRdisplay_1.1 dplyr_1.2.1 withr_3.0.3
[67] BiocManager_1.30.27 fastmap_1.2.0 rhdf5filters_1.22.0
[70] digest_0.6.39 timechange_0.4.0 R6_2.6.1
[73] colorspace_2.1-3 gtools_3.9.5 lpSolve_5.6.23
[76] biomaRt_2.66.1 RSQLite_3.53.3 tidyr_1.3.2
[79] hexbin_1.28.6 data.table_1.18.4 recipes_1.3.3
[82] FNN_1.1.4.1 class_7.3-23 prettyunits_1.2.0
[85] PSMatch_1.14.0 httr_1.4.8 htmlwidgets_1.6.4
[88] S4Arrays_1.10.1 ModelMetrics_1.2.2.2 pkgconfig_2.0.3
[91] gtable_0.3.6 timeDate_4052.112 blob_1.3.0
[94] S7_0.2.2 impute_1.84.0 XVector_0.50.0
[97] htmltools_0.5.9 MALDIquant_1.22.3 clue_0.3-68
[100] scales_1.4.0 png_0.1-9 gower_1.0.2
[103] knitr_1.51 MetaboCoreUtils_1.18.1 reshape2_1.4.5
[106] uuid_1.2-2 coda_0.19-4.1 nlme_3.1-170
[109] curl_7.1.0 anndataR_1.3.1 rhdf5_2.54.1
[112] repr_1.1.7 proxy_0.4-29 cachem_1.1.0
[115] stringr_1.6.0 parallel_4.5.3 mzID_1.48.0
[118] vsn_3.78.1 pillar_1.11.1 grid_4.5.3
[121] vctrs_0.7.3 pcaMethods_2.2.0 randomForest_4.7-1.2
[124] dbplyr_2.6.0 xtable_1.8-8 evaluate_1.0.5
[127] mvtnorm_1.4-2 cli_3.6.6 compiler_4.5.3
[130] rlang_1.3.0 crayon_1.5.3 future.apply_1.20.2
[133] labeling_0.4.3 LaplacesDemon_16.1.8 mclust_6.1.3
[136] QFeatures_1.20.0 affy_1.88.0 plyr_1.8.9
[139] fs_2.1.0 stringi_1.8.9 viridisLite_0.4.3
[142] Biostrings_2.78.0 lazyeval_0.2.3 Matrix_1.7-5
[145] IRkernel_1.3.2 hms_1.1.4 bit64_4.8.2
[148] future_1.75.0 Rhdf5lib_1.32.0 ggplot2_4.0.3
[151] KEGGREST_1.50.0 statmod_1.5.2 SummarizedExperiment_1.40.0
[154] kernlab_0.9-33 igraph_2.3.3 memoise_2.0.1
[157] affyio_1.80.0 sampling_2.11 bit_4.6.0