use scarf for single-cell imaging dataset, trying to generate a zarr array from scratch

[FrenchSauce]

hi!

I would like to use scarf functionalities but with a single cell imaging dataset, so instead of genes my features are shape, intensity, textures. I used to do it with Seurat or Scanpy but this time I have a huge dataset (million of cells) and using scarf could be a very nice optimization, thanks for developing such a tool (and for the good documentation)!
so since I don't have a scRNA or scATAC dataset I am trying to generate a zarr array from scratch before converting it to DataStore object, this is what i have done ( inspired by this thread #54) :

import zarr
from scarf.writers import create_zarr_obj_array

ncells = len(df)
nfeatures = len(measurements)

z = zarr.group()
counts  = z.create_group('RNA')
# for now I am trying with an array filled with zeros
arr = counts.zeros('counts', shape=(ncells, nfeatures), chunks=(100000, nfeatures), dtype='i4')

z.RNA.attrs["is_assay"] = True
z.RNA.attrs["misc"] = {}

g = z.create_group("RNA/featureData")
create_zarr_obj_array(g, "ids", measurements,'str')
create_zarr_obj_array(g, "names", measurements,'str')
create_zarr_obj_array(g, "I", [True for _ in range(len(measurements))], "bool")
g = z.create_group("cellData")
create_zarr_obj_array(g, "ids",  df['id_roi'].to_numpy())
create_zarr_obj_array(g, "names", df['id_roi'].to_numpy())
create_zarr_obj_array(g, "I", [True for _ in range(ncells)], "bool") 

zarr.save('E:/example.zarr', z)

The current tree looks like:

/
 ├── RNA
 │   ├── counts (26718720, 304) int32
 │   └── featureData
 │       ├── I (304,) bool
 │       ├── ids (304,) object
 │       └── names (304,) object
 └── cellData
     ├── I (26718720,) bool
     ├── ids (26718720,) <U9
     └── names (26718720,) <U9

but when trying to convert it to a DataStore object:

import scarf
ds = scarf.DataStore(
    'E:/example.zarr',
    nthreads=4,
    default_assay='RNA'
)

I get this error:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Input In [15], in <cell line: 2>()
      1 import scarf
----> 2 ds = scarf.DataStore(
      3     'E:/Fabien/TMA CRC 2022 dataset2/example.zarr',
      4     assay_types = {'RNA':'RNA'},
      5     nthreads=4,
      6     default_assay='RNA'
      7 )

File C:\Hyapp\Anaconda3-2020.11\envs\tmaclust\lib\site-packages\scarf\datastore.py:3419, in DataStore.__init__(self, zarr_loc, assay_types, default_assay, min_features_per_cell, min_cells_per_feature, mito_pattern, ribo_pattern, nthreads, zarr_mode, synchronizer)
   3417 if synchronizer is None:
   3418     synchronizer = zarr.ThreadSynchronizer()
-> 3419 super().__init__(
   3420     zarr_loc=zarr_loc,
   3421     assay_types=assay_types,
   3422     default_assay=default_assay,
   3423     min_features_per_cell=min_features_per_cell,
   3424     min_cells_per_feature=min_cells_per_feature,
   3425     mito_pattern=mito_pattern,
   3426     ribo_pattern=ribo_pattern,
   3427     nthreads=nthreads,
   3428     zarr_mode=zarr_mode,
   3429     synchronizer=synchronizer,
   3430 )

File C:\Hyapp\Anaconda3-2020.11\envs\tmaclust\lib\site-packages\scarf\datastore.py:2526, in MappingDatastore.__init__(self, **kwargs)
   2525 def __init__(self, **kwargs):
-> 2526     super().__init__(**kwargs)

File C:\Hyapp\Anaconda3-2020.11\envs\tmaclust\lib\site-packages\scarf\datastore.py:603, in GraphDataStore.__init__(self, **kwargs)
    602 def __init__(self, **kwargs):
--> 603     super().__init__(**kwargs)

File C:\Hyapp\Anaconda3-2020.11\envs\tmaclust\lib\site-packages\scarf\datastore.py:102, in BaseDataStore.__init__(self, zarr_loc, assay_types, default_assay, min_features_per_cell, min_cells_per_feature, mito_pattern, ribo_pattern, nthreads, zarr_mode, synchronizer)
    100 self.nthreads = nthreads
    101 # The order is critical here:
--> 102 self.cells = self._load_cells()
    103 self._defaultAssay = self._load_default_assay(default_assay)
    104 self._load_assays(min_cells_per_feature, assay_types)

File C:\Hyapp\Anaconda3-2020.11\envs\tmaclust\lib\site-packages\scarf\datastore.py:123, in BaseDataStore._load_cells(self)
    121 if "cellData" not in self.z:
    122     raise KeyError("ERROR: cellData not found in zarr file")
--> 123 return MetaData(self.z["cellData"])

File C:\Hyapp\Anaconda3-2020.11\envs\tmaclust\lib\site-packages\zarr\core.py:788, in Array.__getitem__(self, selection)
    786     result = self.vindex[selection]
    787 else:
--> 788     result = self.get_basic_selection(pure_selection, fields=fields)
    789 return result

File C:\Hyapp\Anaconda3-2020.11\envs\tmaclust\lib\site-packages\zarr\core.py:907, in Array.get_basic_selection(self, selection, out, fields)
    904     self._load_metadata()
    906 # check args
--> 907 check_fields(fields, self._dtype)
    909 # handle zero-dimensional arrays
    910 if self._shape == ():

File C:\Hyapp\Anaconda3-2020.11\envs\tmaclust\lib\site-packages\zarr\indexing.py:849, in check_fields(fields, dtype)
    847 if fields:
    848     if dtype.names is None:
--> 849         raise IndexError("invalid 'fields' argument, array does not have any fields")
    850     try:
    851         if isinstance(fields, str):
    852             # single field selection

IndexError: invalid 'fields' argument, array does not have any fields

I am not sure what I am missing

[parashardhapola]

Hi!

Thanks for the appreciation! and I'm glad that you are trying this with Scarf.

I'm not sure what exactly is causing the problem.

Can you please send a portion of your data (df and measurements) to me. Maybe first 1000 lines of df and complete measurement to my email: parashar.dhapola@med.lu.se? This will help me reproduce the error.

[parashardhapola]

Hi!

I had a crack at this and it seems to be working. Try it out at your end. The issue was that somehow the count matrix was not created within the RNA group.

import zarr
from scarf.writers import create_zarr_obj_array
import pandas as pd
import scarf

df = pd.DataFrame({'feat_A': [1, 2, 3], 'feat_B': [4, 5, 6], 'id_roi': ['C1', 'C2', 'C3']})
measurements = ['feat_A', 'feat_B']

ncells = len(df)
nfeatures = len(measurements)

modality_name = 'RNA'

z = zarr.open('custom.zarr')

g = z.create_group(modality_name, overwrite=True)
g.create_dataset('counts', shape=(ncells, nfeatures), chunks=(100000, nfeatures), dtype='i4')
z[modality_name].attrs["is_assay"] = True
z[modality_name].attrs["misc"] = {}

g = z.create_group("RNA/featureData", overwrite=True)
create_zarr_obj_array(g, "ids", measurements,'str')
create_zarr_obj_array(g, "names", measurements,'str')
create_zarr_obj_array(g, "I", [True for _ in range(len(measurements))], "bool")

g = z.create_group("cellData", overwrite=True)
create_zarr_obj_array(g, "ids",  df['id_roi'].to_numpy())
create_zarr_obj_array(g, "names", df['id_roi'].to_numpy())
create_zarr_obj_array(g, "I", [True for _ in range(ncells)], "bool")

ds = scarf.DataStore(
    'custom.zarr',
    nthreads=4,
    default_assay='RNA',
    min_cells_per_feature=1,
    min_features_per_cell=1
)

[FrenchSauce]

thanks a lot for your fast answer! I missunderstood how zarr opening operate, it works perfectly now