scverse datastructure for AIRR data

[grst]

Now that scirpy is part of scverse, we could think of an improved data structure for scAIRR data. See also the discussion at scverse/scanpy#1387.

The challenge with scAIRR data is that

• 1 cell can have n chains. Up to four of them are biologically meaningful but there could be more for technical reasons.
• Each chain has a lot of fields. See the AIRR rearrangement standard.

The current pragmatic solution is to store all fields in adata.obs.

• All columns from the airr rearrangement schema are repeated four times
• Excess chains are serialized into JSON and stored in an extra column. These chains are not used by scirpy, but enable lossless conversions.
• The downside is that there can easily be 100+ columns in adata.obs. Also serializing excess chains is not really elegant.
• The advantage is that it works really well with scanpy, i.e. any AIRR variable can immediately be used for grouping, plotting etc.

New options are

• mudata. AIRR data could be saved as a separate modality. Even if we keep the current reprepsentation of a wide data frame, it would at least not clutter the rest of adata.obs.
• awkward array support. Allows storing an arbitrary number of values per row. See first attempt to support awkward arrays anndata#647

The new representation should also aim at being a community standard for the scverse ecosystem and should build upon the AIRR rearrangement standard. Ideally, we could get additional stakeholders onboard, including conga, dandelion, tcrdist3 and possibly members of the AIRR community.

• what's the state of the AIRR single-cell schema? And what are its advantates over the rearrangement schema.

[bcorrie]

Hi All, I am active in the AIRR Standards group, can try and answer questions. We are looking at integrating scirpy into the iReceptor Gateway (gateway.ireceptor.org) to do analysis on AIRR Cell data - so have a vested interest... 8-)

The AIRR Cell schema was just announce at the AIRR meeting this week (https://www.antibodysociety.org/the-airr-community/meetings/airr-community-meeting-vi-exploring-new-frontiers/), with a tagged (v1.4) release expected on github in the next few weeks. This is an "Experimental" schema, meaning that we think it is pretty stable, but will probably need some minor changes before it is finalized, probably in an AIRR Standard v2.0 release sometime in the "near" future.

[bcorrie]

The Cell schema is intended to represent a Cell as an object. It is a light weight object, and contains a cell_id that is used to cross reference information in other collections such as Rearrangement and CellExpression. So if you have a cell_id from a Cell you can find all of the Rearrangements for the paired chains associated with the cell in an AIRR TSV file by searching for that same cell_id. Similarly you can find all of the CellExpression information by searching a CellExpression file for that same cell_id.

[bcorrie]

Here are the current docs on the master branch, which will be tagged with a v1.4 release soon.

https://docs.airr-community.org/en/latest/datarep/cell.html

[bcorrie]

One of the outstanding things we need to resolve in the AIRR Standard is what if any "efficient" file format should be supported. Currently the default standard for all AIRR objects is a JSON representation of the standard objects. For Rearrangement there is a defined TSV file. Such a file format does not yet exist for Cell and more importantly CellExpression. We could probably use some advice from this community...

We have an about to be release version of our AIRR repository (https://github.com/sfu-ireceptor/turnkey-service-php/tree/production-v4) that can store Cell and CellExpression data, but currently what you get out is a very big, verbose JSON file. It works, but...

[bcorrie]

BTW, I think it would be interesting to decorate a cell's adata.obs with AIRR Repertoire metadata. That is, it would be quite easy (I think) to decorate the adata.obs with all sorts of subject/sample metadata for each cell such as subject.disease.disease_daignosis, subject.age_min, subject.sex, sample.tissue etc so that one can use these criteria for coloring, group, etc as well. It seems to me that these fields in adata.obs would be more valuable than the most of the Rearrangement fields. I could see wanting some Rearrangement fields, maybe VDJ calls, CDR3 length, locus for the paired chains, but it seems like most of the Rearrangement fields would not be very interesting from a slicing/coloring/visualization point of view.

[grst]

Hi @bcorrie,

thanks for reaching out, I'm definitely interested in working with the AIRR community to make this as interoperable as possible.
Your description of the Cell schema is really helpful, before I wasn't sure it was to replace the Rearrangement schema, but now I understood it's built on top of it.

Actually scirpy internally uses something very very similar to the Cell schema already. Every file format read into scirpy is first represented as a list of AirrCell objects that have a cell_id, a list of Rearrangement entries and arbitrary cell-level attributes that are added as columns to adata.obs. Once the Cell schema is final, I could simply adjust my AirrCell implementation to match the cell-schema, which would make reading and writing AIRR-compliant data trivial.

The internal representation is one side of the medal, the other is representing the AIRR data as part of either an AnnData or MuData container, which allows to link the data to other modalities, allows efficient file storage in h5ad/h5mu and enables the interaction with scanpy or other scverse tools. The current way is to dump everything in .obs which is not ideal.

---

Such a file format does not yet exist for Cell and more importantly CellExpression. We could probably use some advice from this community...

Well, how about AnnData? If you are looking for a lower level of abstraction, I think hdf5, zarr or tiledb can handle these sparse data efficiently, but @ivirshup or @gtca are the experts here. Is it really necessary to define yet another format for single-cell expression?

Note that reading/writing gene expression formats would be out-of-scope for scirpy, this should be handled by scanpy or AnnData, or possibly by a future scverse IO package.

It's also worth noting here that the CZI is currently looking into standardizing Feature Observation Matrices (FOM) for single-cell data. Maybe it would be worth for the AIRR community and them getting in touch for talking about a AIRR-specific extension. Their public presence is still very sparse, but maybe @ivirshup can tell more about this initiative and make the contact.

Best,
Gregor

[bcorrie]

Such a file format does not yet exist for Cell and more importantly CellExpression. We could probably use some advice from this community...

Well, how about AnnData? If you are looking for a lower level of abstraction, I think hdf5, zarr or tiledb can handle these sparse data efficiently, but @ivirshup or @gtca are the experts here. Is it really necessary to define yet another format for single-cell expression?

Yes, AnnData seems logical and along the lines of what I was thinking... Definitely do not want YAF (yet another format) 8-)

Like you, this is likely out of scope for the AIRR Community - at least in the sense that the obvious thing for us to do would be to have the AIRR python library (https://github.com/airr-community/airr-standards/tree/master/lang/python) use scanpy to read and write these formats.

The CZI FOM is interesting as well...

[bcorrie]

Actually scirpy internally uses something very very similar to the Cell schema already. Every file format read into scirpy is first represented as a list of AirrCell objects that have a cell_id, a list of Rearrangement entries and arbitrary cell-level attributes that are added as columns to adata.obs. Once the Cell schema is final, I could simply adjust my AirrCell implementation to match the cell-schema, which would make reading and writing AIRR-compliant data trivial.

Very nice - one would hope that these would be closely aligned, and hopefully we can now make sure they are completely aligned 8-)

[grst]

Draft proposal for an awkward-array-based data structure

I've been playing around with the draft implementation of awkward array support in AnnData (scverse/anndata#647). This allows to store in .obsm (or in the future maybe also in .obs, see the comments there) ragged arrays of mixed types for each cell -- exactely what we need to represent the 1:n relationship of cells with chains.

An "awkward array" can be directly created from a list of AirrRearrangement dictionaries, e.g.

import awkward as ak
import scirpy as ir
airr_cells = ir.io.to_airr_cells(adata)
chains_as_awkarr = ak.Array([cell.chains for cell in airr_cells])

# an array of length 3000 (= number of cells), variable width (= number of chains) and 
# AIRR Rearrangement fields on the 3rd dimension
> chains_as_awkarr
<Array [[{c_call: 'TRBC2', ... v_cigar: None}]] type='3000 * var * {"c_call": op...'>

To get, e.g. all junction_aa sequences of the first chain of each cell, this awkward array can be sliced:

> chains_as_awkarr[:, 0, "junction_aa"]
<Array ['CASSLMRLAGDTQYF', ... 'CATEEYGNKLVF'] type='3000 * string'>

It is therefore straightforward and computationally efficient to retrieve individual elements.

In summary: We simply take AIRR compliant data, put it in an efficient data structure and store it an AnnData. Not a lot new to define from our side

Caveats and possible solutions

1. Primary and Secondary chains

The list of chains does not know the concept of primary and secondary chains. This can be solved by adding an index array:

# the index arrays point to the element in the list of chains
primary_vj = [0, 0, 1, 2, 0, 1, ...]
primary_vdj = [1, 3, 5, 3, 2, ...]

The awkward array can then be sliced

> chains_as_awkarr[:, primary_vj, "junction_aa"]

to retrieve e,g. the junction AA sequences for the primary VJ chains.

2. Getting values for plotting

The current solution to put everything in obs was amongst others a pragmatic solution to make all AIRR data immediately available to scanpy plotting functions. This is not possible anymore with the proposed new implementation. Arguably, most of those columns are used very rarely for those purposes anyway.

If they are, here are some solutions to do so

• implement a ir.get function to retrieve data from the awkward array (see also Easy-access getter functions #184). The result can be re-assigned to adata.obs for plotting and/or grouping
• re-implement all scanpy plotting functions as e.g. ir.pl.umap, where color is retrieved from the awkward array instead of .obs.

3. Where to store the data

The current draft implementation only supports awkward arrays in .obsm. IMO it would be completely fine to have an .obsm["airr"]. In the future, it may be possible to have a column .obs["airr"].

It was also previously suggested to use mudata as AIRR data can be considered a modality. Since we don't have any data for .X and .var, I don't see a clear advantage of using mudata, though.

---

I'm looking forward to feedback, especially from people that develop ecosystem packages for
AIRR analysis (@zktuong, ..., ?)

[grst]

@bcorrie, this should make it possible that we merely define a mapping how the AIRR standard is represented as AnnData object rather than converting between a scirpy and an AIRR representation. Need to hash that out in more detail, but roughly

• Cell metadata -> adata.obs
• Rearrangement chains -> adata.obsm["airr"] (as awkward array)
• CellExpression -> adata.X or mudata (in case of CITE-seq)

[ivirshup]

The CZI FOM is interesting as well...

@bcorrie from discussion with the SOMA team, the intent was to keep anndata and SOMA (formerly FOM) pretty consistent with each other. Ideally they are the same thing eventually.

Though I haven't heard much on SOMA development recently, so hopefully we are still on the same page here.

Where to store the data

@grst, I've recently heard people wanting to treat their re-arrangment data as a separate modality. E.g. have an AnnData with just re-arrangement info which could then be merged with other modalities. I think the main point here was to avoid having to arbitrarily associate AIRR data with another modality.

I think @bio-la and @crichgriffin could elaborate on this.

[grst]

I've recently heard people wanting to treat their re-arrangment data as a separate modality.

Wouldn't it anyway work with both anndata and mudata to use .obsm["airr"].

# anndata
ir.tl.something(adata)
# anndata, override obsm_key
ir.tl.something(adata, airr_key="ir")
# mudata
ir.tl.something(mdata.mod["AIRR"])

Or would you prefer something like

# mudata only, use `.mod["AIRR"]` by default
ir.tl.something(mdata)
# ... or override default modality
ir.tl.something(mdata, mod="myairr")

[ivirshup]

I don't remember the details, but one of the issues was with subsetting. Wanting to filter AIRR obs, but not the RNA.

I think it would work with AnnData right now, I think the issue was more about want to use scirpy with a MuData where RNA and AIRR were in separate modalities? I think this sorta worked, but there were issues with mudata operations.

[crichgriffin]

In the context of multimodal data and using MuData, I like the idea of keeping the AIRR data separate from any specific modality. My current workflow is to store the scirpy anndata object as a separate modality in muon, where I also have rna and adt data, as separate modalities. But I agree that it doesn't make a lot of sense, since the scirpy anndata object is mostly empty.

I agree with @grst that storing airr data could potentially work in both AnnData and MuData by using .obsm["airr"] except (iirc) in AnnData .obsm row dimensions match .obs row dimensions, and I assume the same applies to MuData.
What would you do in the instances of cells not having AIRR data? E.g. in PBMC data, monocytes would not have airr data, but you would want to keep them in your AnnData/MuData.

[grst]

What would you do in the instances of cells not having AIRR data?

fill the corresponding rows with NaNs? This is how scirpy also handles this currently.
Although this is a good point in favor of mudata. It would also make ir.pp.merge_with_ir superfluous.

I think it should be easy to support both. Just need to think what to advertise as default in the tutorial.