Tutorial: SOMA shapes

As of TileDB-SOMA 1.15 we’re proud to support a more intutive and extensible notion of shape.

In this notebook, we’ll go through how you use shapes for the dataframes and arrays within your SOMA experiments, when and how you can resize, and options for experiments created before TileDB-SOMA 1.15.

The dataset used is from Peripheral Blood Mononuclear Cells (PBMC), which is freely available from 10X Genomics.

(Please also see the Academy tutorial.)

We’ll start by importing tiledbsoma.

import tiledbsoma

The shape feature

As we’ve seen in other tutorials in this series, the SOMA data model brings across many familiar concepts from AnnData. This includes the ability to ask component dataframes and arrays what their shapes are.

First, let’s unpack and open an experiment.

import tarfile
import tempfile

uri = tempfile.mktemp()
with tarfile.open("data/pbmc3k-sparse.tgz") as handle:
    handle.extractall(uri)
exp = tiledbsoma.Experiment.open(uri)

The obs dataframe has a domain, which is a soft limit on what values can be written to it. You’ll get an exception if you try to read or write soma_joinid values outside this range, which is an important data-integrity reassurance.

The domain we see here matches with the data populated inside of it.

(This will usually be the case. It might not, if you’ve created the dataframe but not written any data to it yet – at that point it’s empty but it still has a shape.)

If you have more data – more cells – to add to the experiment later, you will be able resize the obs, up to the maxdomain which is a hard limit.

exp.obs.domain

((0, 2637),)

exp.obs.maxdomain

((0, 9223372036854773758),)

exp.obs.read().concat().to_pandas()

	soma_joinid	obs_id	n_genes	percent_mito	n_counts	louvain
0	0	AAACATACAACCAC-1	781	0.030178	2419.0	CD4 T cells
1	1	AAACATTGAGCTAC-1	1352	0.037936	4903.0	B cells
2	2	AAACATTGATCAGC-1	1131	0.008897	3147.0	CD4 T cells
3	3	AAACCGTGCTTCCG-1	960	0.017431	2639.0	CD14+ Monocytes
4	4	AAACCGTGTATGCG-1	522	0.012245	980.0	NK cells
...	...	...	...	...	...	...
2633	2633	TTTCGAACTCTCAT-1	1155	0.021104	3459.0	CD14+ Monocytes
2634	2634	TTTCTACTGAGGCA-1	1227	0.009294	3443.0	B cells
2635	2635	TTTCTACTTCCTCG-1	622	0.021971	1684.0	B cells
2636	2636	TTTGCATGAGAGGC-1	454	0.020548	1022.0	B cells
2637	2637	TTTGCATGCCTCAC-1	724	0.008065	1984.0	CD4 T cells

2638 rows × 6 columns

We’ll see more about this on experiment-level resizes below, as well as in the tutorial on TileDB-SOMA’s append mode.

The var dataframe’s domain is similar:

var = exp.ms["RNA"].var
var.domain

((0, 1837),)

var.maxdomain

((0, 9223372036854773968),)

Likewise, the N-dimensional arrays within the experiment have their shapes as well.

There’s an important difference: while the dataframe domain gives you the inclusive lower and upper bounds for soma_joinid writes, the shape for the N-dimensional arrays is the upper bound plus 1.

Since there are 2638 cells and 1838 genes here, X’s shape reflects that.

exp.obs.domain

((0, 2637),)

exp.ms["RNA"].var.domain

((0, 1837),)

exp.ms["RNA"].X["data"].shape

(2638, 1838)

exp.ms["RNA"].X["data"].maxshape

(9223372036854773759, 9223372036854773759)

The other N-dimensional arrays are similar:

obsm = exp.ms["RNA"].obsm
list(obsm.keys())

['X_draw_graph_fr', 'X_pca', 'X_tsne', 'X_umap']

obsp = exp.ms["RNA"].obsp
list(obsp.keys())

['connectivities', 'distances']

[
    obsm["X_pca"].shape,
    obsm["X_pca"].maxshape,
]

[(2638, 50), (9223372036854773759, 9223372036854773759)]

[
    obsp["distances"].shape,
    obsp["distances"].maxshape,
]

[(2638, 2638), (9223372036854773759, 9223372036854773759)]

In particular, the X array in this experiment – and in most experiments – is sparse. That means there needn’t be a number in every row or cell of the matrix. Nonetheless, the shape serves as a soft limit for reads and writes: you’ll get an exception trying to read or write outside of these.

As a convenience, you can see all the experiment’s objects’ shapes at once as follows:

import tiledbsoma.io
tiledbsoma.io.show_experiment_shapes(exp.uri)

As with AnnData, as a general rule you’ll see the following:

• An X array’s shape is nobs x nvar

• An obsm array’s shape is nobs x some number, maybe 50

• An obsp array’s shape is nobs x nobs

• A varm array’s shape is var x some number, maybe 50

• A varp array’s shape is nvar x nvar

When and how to resize at the experiment level

The primary reason you’d resize a dataframe or an array within an experiment is to append more data. For example, say you have an experiment with the results of Monday’s lab run on a sample of 100,000 cells. Then maybe on Tuesday you’ll want to add that day’s lab run of an additional 70,000 cells to the same experiment, for a new total of 170,000 cells. It’s also possible that Tuesday’s data might include some infrequently expressed genes that didn’t appear in Monday’s data.

Because the shapes are soft limits, reading or writing beyond which will result in an exception, you’d need to resize the experiment to accommodate new shapes for the dataframes and arrays in the experiment to allow for new nobs = 170,000.

Please see the append-mode tutorial for how to do that using tiledbsoma.io.register_anndatas and tiledbsoma.io.resize_experiment

While you can resize each dataframe and array in the experiment one at a time – see “Advanced usage”, below in this notebook – by var the most common case is tiledbsoma.io.resize_experiment, which exists to make this simple and convenient.

How to upgrade older experiments

Experiments created by TileDB-SOMA 1.15 and higher will look as shown above. Let’s take a look at an experiment from before TileDB-SOMA 1.15.

import tarfile
import tempfile

import tiledbsoma.io

uri = tempfile.mktemp()
with tarfile.open("data/pbmc3k-sparse-pre-1.15.tgz") as handle:
    handle.extractall(uri)
expold = tiledbsoma.Experiment.open(uri)

This is the same PBMC3K data as above. Compare the old and new shapes:

expold.obs.domain

((0, 9223372036854773758),)

expold.obs.maxdomain

((0, 9223372036854773758),)

expold.obs.tiledbsoma_has_upgraded_domain

False

[ expold.ms["RNA"].X["data"].shape, expold.ms["RNA"].X["data"].maxshape, expold.ms["RNA"].X["data"].tiledbsoma_has_upgraded_shape ]

[(9223372036854773759, 9223372036854773759),
 (9223372036854773759, 9223372036854773759),
 False]

Note that for the pre-1.15 experiment, the shape is huge – like the maxshape – and tiledbsoma_has_upgraded_domain is False.

To make the old experiment look like the new experiment, simply call upgrade_experiment_shapes, and re-open:

tiledbsoma.io.upgrade_experiment_shapes(expold.uri)

True

expold = tiledbsoma.open(expold.uri)

[ expold.ms["RNA"].X["data"].shape, expold.ms["RNA"].X["data"].maxshape, expold.ms["RNA"].X["data"].tiledbsoma_has_upgraded_shape ]

[(2638, 1838), (9223372036854773759, 9223372036854773759), True]

Additionally, you can call tiledbsoma.io.show_experiment_shapes(expold.uri) before and after doing the upgrade.

To run a pre-check, you can do

tiledbsoma.io.upgrade_experiment_shapes(expold.uri, check_only=True)

This won’t change anything – it’ll simply tell you if the operation will be possible.

Advanced usage: dataframes with non-standard index columns

In the SOMA data model, the SparseNDArray and DenseNDArray objects always have int64 dimensions named soma_dim_0, soma_dim_1, and up, and they have a numeric soma_data attribute for the contents of the array. Furthermore, this is always the case.

exp.ms["RNA"].X["data"].schema

soma_dim_0: int64 not null
soma_dim_1: int64 not null
soma_data: float not null

For dataframes, though, while there must be a soma_joinid column of type int64, you can have one or more other index columns in addtion – or, soma_joinid can be a non-index column.

This means that in the default, simplest, and most common case, you can think of a dataframe has having a shape just as the N-dimensional arrays do.

exp.obs.schema

soma_joinid: int64 not null
obs_id: large_string
n_genes: int64
percent_mito: float
n_counts: float
louvain: dictionary<values=string, indices=int32, ordered=0>

exp.obs.index_column_names

('soma_joinid',)

But really, dataframes are capable of more than that, via the index-column names you specify at creation time.

Let’s create a couple dataframes, with the same data, but different choices of index-column names.

sdfuri1 = tempfile.mktemp()
sdfuri2 = tempfile.mktemp()

import pyarrow as pa

schema = pa.schema([
    ("soma_joinid", pa.int64()),
    ("mystring", pa.string()),
    ("myint", pa.int32()),
    ("myfloat", pa.float32()),
])

data = pa.Table.from_pydict({
    "soma_joinid": [0, 1],
    "mystring": ["hello", "world"],
    "myint": [33, 44],
    "myfloat": [4.5, 5.5],
})

with tiledbsoma.DataFrame.create(
    sdfuri1,
    schema=schema,
    index_column_names=["soma_joinid", "mystring"],
    domain=[(0, 9), None],
) as sdf1:
        sdf1.write(data)

Now let’s look at the domain and maxdomain for these dataframes.

sdf1 = tiledbsoma.DataFrame.open(sdfuri1)

sdf1.index_column_names

('soma_joinid', 'mystring')

Here we see the soma_joinid slot of the dataframe’s domain is as requested.

Another point is that domain cannot be specified for string-type index columns.

You can set them at create one of two ways:

domain=[(0, 9), None],

or

domain=[(0, 9), ('', '')],

and in either case the domain slot for a string-typed index column will read back as ('', '').

sdf1.domain

((0, 9), ('', ''))

sdf1.maxdomain

((0, 9223372036854775796), ('', ''))

Now let’s look at our other dataframe. Here soma_joinid is not an index column at all. This is fine, as long as within the data you write to it, the index-column values uniquely identify each row.

with tiledbsoma.DataFrame.create(
    sdfuri2,
    schema=schema,
    index_column_names=["myfloat", "myint"],
    domain=[(0, 999), (-1000, 1000)],
) as sdf2:
        sdf2.write(data)

sdf2 = tiledbsoma.DataFrame.open(sdfuri2)

sdf2.index_column_names

('myfloat', 'myint')

The domain reads back as written.

sdf2.domain

((0.0, 999.0), (-1000, 1000))

sdf2.maxdomain

((-3.4028234663852886e+38, 3.4028234663852886e+38), (-2147483648, 2147481645))

Advanced usage: using resize at the dataframe/array level using the SOMA API

Above we saw a simple and convenient way to resize all the dataframes and arrays within an experiment.

However, should you choose to do so, you can apply these one dataframe or array at a time.

For N-dimensional arrays that have been upgraded, or that were created using TileDB-SOMA 1.15 or higher, simply do the following:

• If the array’s .tiledbsoma_has_upgraded_shape reports False, invoke the .tiledbsoma_upgrade_shape method.

• Otherwise invoke the .resize method.

Let’s do a fresh unpack of a pre-1.15 experiment:

import tarfile
import tempfile

uri = tempfile.mktemp()
with tarfile.open("data/pbmc3k-sparse-pre-1.15.tgz") as handle:
    handle.extractall(uri)
expold = tiledbsoma.Experiment.open(uri)
X = expold.ms["RNA"].X["data"]

Here we see that the X array has not been upgraded, and that its shape reports the same as maxshape:

X.tiledbsoma_has_upgraded_shape

False

X.shape

(9223372036854773759, 9223372036854773759)

Now let’s give the X array the new-style shape:

X.non_empty_domain()

((0, 2637), (0, 1837))

with tiledbsoma.Experiment.open(uri, "w") as exp:
    exp.ms["RNA"].X["data"].tiledbsoma_upgrade_shape([X.non_empty_domain()[0][1]+1, X.non_empty_domain()[1][1]+1])

Next let’s re-open and see what happened:

expold = tiledbsoma.Experiment.open(expold.uri)
X = expold.ms["RNA"].X["data"]

X.tiledbsoma_has_upgraded_shape

True

X.shape

(2638, 1838)

X.maxshape

(9223372036854773759, 9223372036854773759)

Furthermore, we can resize it even farther:

with tiledbsoma.Experiment.open(expold.uri, "w") as exp:
    exp.ms["RNA"].X["data"].resize([7200, 1848])

expold = tiledbsoma.Experiment.open(expold.uri)
X = expold.ms["RNA"].X["data"]

X.shape

(7200, 1848)

For dataframes, the process is similar. If you want to expand only the soft limits for soma_joinid, you can use some simpler methods:

• If the dataframe’s tiledbsoma_has_upgraded_domain reports False, invoke .tiledbsoma_upgrade_domain

• Otherwise invoke the .change_domain method.

expold.obs.tiledbsoma_has_upgraded_domain

False

expold.obs.domain

((0, 9223372036854773758),)

expold.obs.maxdomain

((0, 9223372036854773758),)

expold.obs.non_empty_domain()

((0, 2637),)

with tiledbsoma.Experiment.open(expold.uri, "w") as exp:
    exp.obs.tiledbsoma_upgrade_domain([[0, expold.obs.non_empty_domain()[0][1]+1]])

expold = tiledbsoma.Experiment.open(expold.uri)

expold.obs.tiledbsoma_has_upgraded_domain

True

expold.obs.domain

((0, 2638),)

expold.obs.maxdomain

((0, 9223372036854773758),)