Mapping between features for MS-proteomics support

[lucas-diedrich]

Hi, great package!

---

I’m a member of the scverse proteomics interest group (led by @mikelkou) and, together with @vbrennsteiner co-developer of alphapepttools that uses anndata as central data structure to perform exploratory downstream analysis of LC/MS-proteomics data.

Within the proteomics interest group we observed that many properties of MuData would be very useful for the analysis of MS proteomics data. However, MuData currently lacks a first-class way to represent and query explicit relationships between feature indices across modalities.

Why?

In LC/MS-proteomics, we have a naturally arising hierarchical feature structure

• At the lowest level, we detect + quantify fragments from charged peptides (precursors) in the mass spectrometry (MS) instruments. The precursor-level data is relatively large (N samples x ~100 000 features)
• Proteomics search engines identify the precursor sequences and match them to their corresponding proteins. Ultimately, search engines derive protein-specific intensities (N samples x 3000-10000 features). Note that there is an N:M relationship between precursors and proteins, i.e. many precursors can map to one protein and sometimes a precursor could potentially be derived from different (homologous) proteins.

While the protein intensities are most commonly used for the biological interpretation of the data, many analyses depend on the knowledge of the precursor or peptide-level intensities. This includes, for example,

• software packages for differential expression testing, e.g. (alphaquant from our lab)
• peptide-level assays
• but also basic quality control/data inspection tasks.

Currently, we are able to read any of the desired feature levels (precursor < peptides < proteins < genes) into anndata objects. And we can obviously also bind them to mudata objects. However, at the current state, mudata does not support the efficient mapping/querying of features between these different layers of evidence to answer questions like:

• "Which protein belongs to the peptide with the sequence "PETPIDE"
• "Plot a boxplot of the intensity distributions of precursors that map to the protein MAPK"

Suggested feature

What would simplify the development of proteomics tools in the scverse ecosystem is an explicit and queryable way to link features across modalities within a MuData object — for example, a directed acyclic graph (DAG) or similar abstraction that formalizes relationships between variable indices across containers.

Current state

There are already ongoing efforts/suggestions into this direction

• @mffrank implemented a proof-of-principle DAG mapping for mudata
• @ilan-gold suggested to use pandas extension index arrays to achieve this behavior.

We are opening this issue

1. Because we wanted to ask whether such functionality would be considered in scope for mudata itself (or better suited for an extension)
2. If so, discuss possible implementation approaches (e.g. those suggested by @mffrank or @ilan-gold)
3. Avoid redundant or fragmented efforts across projects.

While we focus here on LC/MS proteomics, similar hierarchical or many-to-many feature relationships also arise in other multi-modal settings (e.g. ATAC-RNA, Protein-RNA, ...) so we think that this would generally be a super useful feature.

We would be really happy to help in any way to enable this feature, please let us know what the best way to contribute is.

---

@mikelkou @josenimo @LucaMarconato @melonora Theodoros Visvikis @mffrank @ilan-gold @soroorh

[ilan-gold]

@mffrank implemented a proof-of-principle DAG mapping for mudata
@ilan-gold suggested to use pandas extension index arrays to achieve this behavior.

Just want to clarify, I am proposing implementing a DAG within a pandas index via extension arrays (so just building off hte DAG idea I heard on the call, it's a great one). I am not sure concretely, but conceptually, I think this will work. MuData would then maintain an internally "joined" var represnetation that encodes the DAG structure in its index, and its columns would then apply to each member of the index accordingly (so if a precursor has a label, so does everything downstream of it in the DAG etc.)

This is highly conceptual but I think the "cleanest" route worth pursuing

[ilan-gold]

Concretely, a first pass could be just trying out MultiIndex and seeing where it fails: https://pandas.pydata.org/docs/user_guide/advanced.html#multiindex-advanced-indexing

On the face of it, a DAG seems mappable onto this data structure. But I dont know the performance at the scale we will be doing things.

[LucaMarconato]

Thanks @lucas-diedrich for starting the discussion. Supporting queries of data based on relationships between the various entities is something that has been asked multiple times also from the users of spatialdata. I could see value in exploring a general approach applicable to anndata, mudata and spatialdata.

If you have already a concrete implementation path in mind, I'd go with that (to keep momentum going). But in parallel to that, I would find it useful to create a shared design document where we could draft in code or pseudocode, how the syntax for making such queries would look like (or is already implemented), across several use cases. This would help with choosing a consistent design and guide implementation.

If you like the idea, would you like to start such document and add the examples from @mffrank you mentioned? I'll be happy to fill in examples from the spatialdata side.

[ilan-gold]

@LucaMarconato Maybe you could share the examples here? I could imagine havin a DAG pandas index package for example. But I'm also not 100% convinced that we even need that - I would personally like to understand what MultiIndex can't do for this use-case (or yours).

[LucaMarconato]

Some examples:

1. map indices of cells across multiple tables (or tables/other elements, like a geopandas dataframe). Common queries are all the possible SQL-like joins, such as "I have the cells from this segmentation, get me all the cells that are present also in this other segmentation". Or "I have a table for 5 sample, get me only the cells corresponding to table (representing a single sample)".
2. map transcripts belonging to a particular cell. Example query: "get me all the transcripts belonging to this cell"
3. support for N-1 relationships, such as for representing nuclei belonging to a cell. In Xenium data we can have polynucleated cells, or cells without nuclei.

We have some APIs to simplify some of the queries in 1 from the above, such as scverse/spatialdata#894, but we don't have a general approach.

In particular:

• we never offer a graph representation for the relationships above
• we don't support most of the relationships above, letting the users use custom columns if they want
• databases offer a solution for the above, but they are overkill. So we implement some of the relationships above (such as to map table annotations to external elements). Still experimenting more with some relational databases may bring value, especially if we want to offer a solution valid for several use cases.

[ilan-gold]

This all sounds like the obs version of this issue, which I am in favor of bringing into mudata as well. I think there is already support for this because different AnnData objects can have different obs spaces in mudata. Maybe we need a separate issue for that even if the underlying implementation will be similar (some sort of DAG index).

[jan-engelmann]

Interesting discussion. Thanks @ilan-gold

We have implemented a similar hierarchy on obs level for combining genetic/SNP data with single-cell data by synchronizing on donor/sample level.
See the DonorData structure in https://github.com/theislab/cellink/tree/main
So e.g. filtering for a rare cell type drops donor-level data of donors which don't have that cell type.

For genetics it's also interesting to link rare variants to genes. Or annotations per variant which are different depending on the gene. E.g. overlapping genes / isoforms

cc @LArnoldt @HolEv

[ilan-gold]

@jan-engelmann Could you point to how that is implemented? Did you guys consider a multi-index? Are you using a DAG?

It sounds like what you have is along the lines of what @LucaMarconato wants. Very cool!

[jan-engelmann]

Sure! The code is here: https://github.com/theislab/cellink/blob/main/src/cellink/_core/donordata.py

At the moment it only implements the donor-sync functionality. DonorData uses composition and syncs it's two slots .G and .C which are either Anndata or Mudata objects.

The other feature level links/hierarchies are sth we need in the future. (Sorry if that wasnt clear).
We have not yet looked into MultiIndex. At the moment we have to store some variant annotations in the uns

Maybe @HolEv or @LArnoldt can comment with more detail