---
jupyter:
  jupytext:
    text_representation:
      extension: .Rmd
      format_name: rmarkdown
      format_version: '1.2'
  kernelspec:
    display_name: R
    language: R
    name: ir
---

# pRoloc integration: for pRoloc users

The [grassp data portal](https://grassp.apps.czbiohub.org/datasets) 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:**
{doc}`pRoloc integration: for grassp users <notebooks/proloc_tutorial>`. 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

```r
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.

```{r}
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)
```

## Getting a dataset

Browse [the portal](https://grassp.apps.czbiohub.org/datasets), 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:

```{r}
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
```

```{r}
x <- grassp_as_msnset(h5ad_path)
x
```

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:

```{r}
dim(x)
head(fvarLabels(x), 12)
assayDataElementNames(x)
```

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:

```{r}
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")]
```

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`:

```{r}
uns <- experimentData(x)@other$grassp_uns
unlist(uns[c("title", "publication_journal", "publication_doi")])
```

### 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:

```{r}
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
)
```

We will use `marker_lilley`, the set `pRoloc::pRolocmarkers()` also ships:

```{r}
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))
```

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

```{r}
plot2D(x, fcol = FCOL, main = "pRoloc default PCA (scaled)")
addLegend(x, fcol = FCOL, where = "topright", cex = 0.6, ncol = 2)
```

```{r}
options(repr.plot.height = 4)
plot2D(x, method = "scree")
options(repr.plot.height = 6)
```

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**:

```{r}
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]))
```

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:

```{r}
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))
```

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:

```{r}
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)
```

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()`.

```{r}
c(max_abs_difference = max(abs(aligned - stored)))
stopifnot(all(abs(diag(cor(recomputed, stored))) > 0.9999))
```

### 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:

```{r}
umap <- fData(x)$X_umap
colnames(umap) <- c("UMAP1", "UMAP2")
par(mfrow = c(1, 2))
plot2D(stored[, 1:2], method = "none", methargs = list(x), fcol = FCOL, main = "portal X_pca")
plot2D(umap, method = "none", methargs = list(x), fcol = FCOL, main = "portal X_umap")
par(mfrow = c(1, 1))
```

## 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:

```{r}
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))
```

```{r}
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.

```{r}
summary(QSep(x, fcol = FCOL))
sapply(marker_cols, function(f) median(summary(QSep(x, fcol = f), verbose = FALSE)))
```

## 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:

```r
params <- svmOptimisation(x, fcol = FCOL, times = 100,
                          class.weights = classWeights(x, fcol = FCOL))
plot(params)
getParams(params)
```

```{r}
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)
```

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.

```{r}
# `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")
```

```{r}
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?

```{r}
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)))
```

## 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:

```{r}
suppressPackageStartupMessages(library(pRolocdata))
data(dunkley2006)
dunkley2006
```

```{r}
out <- file.path(tempdir(), "dunkley2006.h5ad")
grassp_write_msnset(dunkley2006, out, overwrite = TRUE)
```

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:

```{r}
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")
```

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 {doc}`../api/io`.

## Session info

```{r}
sessionInfo()
```