RFC

This page describes the design ideas of the mulink repository.

Terminology

Motivation

Modern -omics experiments often generate measurements across multiple biological layers (modalities), generating multiple views of a biological specimen. Often, the same biological concepts (e.g. a transcript and the corresponding protein) are measured across these modalities. It is often useful to keep track of these conceptual relationships.

mulink aims to formalize the relationship between features from different modalities by explicitly indicating their relationship.

Note that we interpret modalities very loosely in this context. They might represent distinct measurements (e.g. matched transcriptomics + proteomics measurements), derived measurements (e.g. protein quantification data which was aggregated from charged precursor intensities in MS-based proteomics), or even other loose relationships (e.g. receptor-ligand interactions for cell-cell-communication analysis). The common denominator is that features are in relation to one another.

Data representation

We represent the feature-feature relation (feature mapping) as directed acyclic graph. In this representation, each feature is a node and the relationship to another feature is represented as a directed graph edge (think: arrow).

The meaning/semantics of the edge might differ depending on the application.

Examples

MS-proteomics: In MS-proteomics biologically meaningful protein intensities represent a inferred measurement from the observed measurements of charged peptides (precursors) that were generated via tryptic digests from proteins. Multiple precursors can be derived from a single protein (N precursors map to 1 protein). Due to sequence homology, a precursor might also be derived from different proteins (1 precursor might map to M proteins). This corresponds to a N:M mapping and can be represented as graph:

        flowchart LR

    r1[precursor 1] -->|derived from| p1[protein 1]
    r2[precursor 2] -->|derived from| p1[protein 1]
    r3[precursor 3] -->|derived from| p1[protein 1]
    r3[precursor 4] -->|derived from| p2[protein 2]
    

Gene set enrichment analysis decoupler represents gene sets as bipartite graph. The features are gene set terms (e.g. GO ontology terms) that map to genes.

        flowchart LR

    gs1[Gene Set 1] -->|contains| g1[Gene 1]
    gs1[Gene Set 1] -->|contains| g2[Gene 2]
    gs1[Gene Set 1] -->|contains| g3[Gene 3]

    

Important concepts

We follow the networkx conventions for the representation of graphs.

Graph

A structure consisting of nodes (vertices) and edges that connect them.

Directed Graph

A graph in which the edges are directed (arrows).

Directed Acyclic Graph

A directed graph that does not contain cycles, i.e. trails where the first and last node are equal.

A minimal directed acyclic graph:

        flowchart LR

    A --> B
    

A minimal directed cyclic graph:

        flowchart LR

    A --> B
    B --> A
    

Adjacency matrix

A matrix representation of a graph (see also adjacency_matrix). If the graph consists of $p$ vertices, the adjacency matrix has a shape of $p \times p$. A non-zero value at row $u$ and column $v$ indicates that node $u$ has an edge pointing to node $v$.

![ Important ] Negative edge weights do not invert connections. A negative edge weight at position (u, v) still indicates a connection directed from node u to node v.

E.g. the corresponding adjacency matrix of this graph

        flowchart LR

    A --> B
    

would be

$$\begin{matrix} & A & B \\ \hline A & 0 & 1 \\ B & 0 & 0 \\ \end{matrix}$$

Descendant + Ancestor

A node $v$ is a descendant of node $u$ if there is a path (a sequence of directed edges) from $u$ to $v$

$$u \rightarrow ... \rightarrow v$$

The node $u$ is then the ancestor of node $v$.

Transitive closure

If a node $c$ is reachable from a node $a$ via node $b$, then there is also a direct edge between $a$ and $c$

A minimal directed graph:

        flowchart LR

    A --> B --> C
    

With transitive closure

        flowchart LR

    A --> B --> C
    A --> C
    

Related projects

QFeatures implements a similar approach to represent quantitative proteomics data.

Implementation

Interoperability with mudata

We extend the mudata framework that natively supports multimodal data with a custom namespace. We encode the desired graph structure with the existing varp attribute (see graph representation)

The general syntax is:

import mudata
import mulink

mdata = mudata.MuData(...)
mdata.link.<function>

Considered Alternatives

• Creating a custom class that inherits from mudata. This strongly couples the development to changes in the mudata class

• Creating a custom class via composition of mudata with the desired feature-mapping graph. This reduces the interoperability with existing tools and does not provide significant benefits as the feature is largely supported in within the existing framework. Further, serialization and on-disk format would have to be reimplemented for a custom class.

Graph representation

We store the adjacency matrix of the graph that encodes the feature mapping as compressed sparse row (csr) matrix in the varp attribute of the mudata.MuData object.

Advantages

• The varp attribute is automatically synchronized with the mudata index, i.e. it stays up to date upon querying, etc.

• The varp attribute can be serialized as is with the mudata framework.

• The sparse representation of the graph adjacency matrix is memory efficient.

Potential disadvantages

• The csr matrix format is optimized for row-wise lookups (i.e. bottom-up actions). This might lead to performance issues for inverse (top-down) actions.

Considered alternatives

• Storing the feature mapping in the mdata.uns attribute might also be feasible but would loose the advantage of the synchronization of the varp attribute with the var axis.

Multiple link types

Per default, we store the feature mapping matrix under the key mdata.varp["feature_mapping"]. Natively, multiple link types can be supported by using a custom key via the key parameter in the respective functions.

If only explicit mappings between a subset of modalities are desired, the rest of the adjacency matrix can be intentionally set to zeros/be left uninitialized.

Querying

Given the directed graph structure, we enable querying of direct descendant nodes and ancestor nodes.

This returns a mudata object containing all observations but only the relevant subset of features.

subset = mdata.link.query_descendants(["feature1", "feature2"], include_self=True)
# > returns mdata subset that includes feature1, feature2 and their descendants

assert mdata.n_obs == subset.n_obs
assert "feature1" in subset.var_names
assert "feature2" in subset.var_names

mdata.link.query_ancestors(["feature1", "feature2"], include_self=True)
# > returns mdata subset that includes feature1, feature2 and their ancestors

assert mdata.n_obs == subset.n_obs
assert mdata.n_obs == subset.n_obs
assert "feature1" in subset.var_names
assert "feature2" in subset.var_names

A specific modality can be extracted as anndata.AnnData object from these subsets with the existing mudata syntax:

subset.mod["mod"]

Self-behaviour

Depending on the application, it might be useful to keep the query. We allow to optionally keep or remove the query components from the subset with the include_self argument.

Considered alternatives

• Add self referring edges to features: This would remove the flexibility of either including or excluding the nodes or would require additional logic to enable this. Additionally, the resulting graph structure would contain cycles, which reduces the usability of the graph itself. The semantics of self-referring edges are not clear in many contexts (e.g. feature aggregation, …)

        flowchart TD
    A --> B
    A --> A
    B --> B
    

Multi-hop traversal

![ Note ] Work in progress

Currently, we only support querying for direct ancestors. These means that in hierarchical structures, the graph should be transitively closed.

Aggregation

![ Note ] Work in progress

Plotting

![ Note ] Work in progress

Constraints and validation

![ Note ] Work in progress

Currently, no validation or constraints are enfored for the adjacency matrix/graph, except for the shape (must be equivalent to the number of features in the mudata object $p \times p$)

Future constrains

Potential constraints that could be enforced:

• minimal connectivity: Every feature should contain at least one incoming or outgoing connection

• maximal connectivity: Feature might only be allowed a maximum number of incoming/outgoing edges.

• acyclicity: Graph should not contain cycles

• self-referential: Features should not map to features from the same modality

• hierarchy: Features from one modality might only be allowed to map to features from a subset of other modalities

As the package should support various use cases, the enforcement of these constraints will likely have to be context-depedent and configurable. The constraints would also have to be serializable, e.g. in the mdata.uns attribute or the mdata.varp['feature_mapping'].attrs attribute.

Encoding and serialization

![ Note ] Work in progress

Some package-specific functionalities might be facilitated by storing additional metadata with the feature mapping. This might include

• indicator that the adjacency matrix conforms with mulink conventions

• versioning

• encoding of additional constraints/validations

Non-goals

mulink is not a general graph analysis library and is meant to faciliate relating different modalities in mudata objects; for complex graph algorithms, export the adjacency matrix and use graph analysis libraries directly.