# Reading in a 1.1 million cell HDF5 dataset

**URL:** <https://discourse.scverse.org/t/reading-in-a-1-1-million-cell-hdf5-dataset/375>\
**Category:** scRNA-seq\
**Tags:** h5\
**Created:** [March 25, 2022, 12:08am UTC](https://discourse.scverse.org/t/reading-in-a-1-1-million-cell-hdf5-dataset/375 "2022-03-25T00:08:26Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![salwanbutrus](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/salwanbutrus/32/184_2.png) [@salwanbutrus](https://discourse.scverse.org/u/salwanbutrus)\
**Post date:** [March 25, 2022, 12:08am UTC](https://discourse.scverse.org/t/reading-in-a-1-1-million-cell-hdf5-dataset/375/1 "2022-03-25T00:08:26Z")

</div>

Hello,

I am trying to read in a huge dataset from the [Allen](https://portal.brain-map.org/atlases-and-data/rnaseq/mouse-whole-cortex-and-hippocampus-10x) atlas. File name is **Gene expression matrix (HDF5)** 5.3 GB. I tried the following script but even on UC Berkeley’s Savio computational cluster, I’m pausing it after 10 hours, stuck at the “counts =” step. Is there a more efficient way to read the data in? Should I just let it run for much longer than 10 hours?

Thank you!

````auto
from scipy.sparse import csr_matrix
import numpy as np
import h5py
import scanpy as sc
import anndata

hf = h5py.File('expression_matrix.hdf5', 'r')

counts = csr_matrix(np.array(hf.get('data').get('counts')))
genes = list(hf.get('data').get('gene'))
cells = list(hf.get('data').get('samples'))

adata = anndata.AnnData(counts.transpose())

adata.obs.index = cells
adata.var.index = genes```
````

---

<div class="post-metadata">

**Author:** ![adamgayoso](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/adamgayoso/32/100_2.png) [@adamgayoso](https://discourse.scverse.org/u/adamgayoso)\
**Post date:** [March 25, 2022, 12:48am UTC](https://discourse.scverse.org/t/reading-in-a-1-1-million-cell-hdf5-dataset/375/2 "2022-03-25T00:48:28Z")

</div>

> [@salwanbutrus](#):
>
> `hf.get('data').get('counts')`

What data format is this?

---

<div class="post-metadata">

**Author:** ![adamgayoso](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/adamgayoso/32/100_2.png) [@adamgayoso](https://discourse.scverse.org/u/adamgayoso)\
**Post date:** [March 25, 2022, 1:52am UTC](https://discourse.scverse.org/t/reading-in-a-1-1-million-cell-hdf5-dataset/375/3 "2022-03-25T01:52:16Z")

</div>

@salwanbutrus I checked and it’s a hdf5dataset. I’m afraid there’s not much you can do (given what I understand). The authors should have saved an hdf5 with the components to build a sparse array instead of the dense array.

Though it may work to write a for loop and make the anndata for a chunk of cells (e.g., 100 cells). Then concatenate the anndatas together (as they will all have sparse format data). That whole dataset as dense in memory is probably causing issues.

---

<div class="post-metadata">

**Author:** ![ivirshup](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/ivirshup/32/160_2.png) [@ivirshup](https://discourse.scverse.org/u/ivirshup)\
**Post date:** [March 25, 2022, 12:33pm UTC](https://discourse.scverse.org/t/reading-in-a-1-1-million-cell-hdf5-dataset/375/4 "2022-03-25T12:33:27Z")

</div>

> Though it may work to write a for loop and make the anndata for a chunk of cells (e.g., 100 cells).

I think I would probably just work with the arrays directly here, then construct the `anndata` from that. `anndata` already has some code for doing this internally, which you could just copy:

> **That code**
>
> ```python
> from scipy import sparse
> 
> def idx_chunks_along_axis(shape: tuple, axis: int, chunk_size: int):
> """\
> Gives indexer tuples chunked along an axis.
> 
> Params
> ------
> shape
> Shape of array to be chunked
> axis
> Axis to chunk along
> chunk_size
> Size of chunk along axis
> 
> Returns
> -------
> An iterator of tuples for indexing into an array of passed shape.
> """
> total = shape[axis]
> cur = 0
> mutable_idx = [slice(None) for i in range(len(shape))]
> while cur + chunk_size < total:
> mutable_idx[axis] = slice(cur, cur + chunk_size)
> yield tuple(mutable_idx)
> cur += chunk_size
> mutable_idx[axis] = slice(cur, None)
> yield tuple(mutable_idx)
> 
> def read_dense_as_csr(dataset, axis_chunk=6000):
> sub_matrices = []
> for idx in idx_chunks_along_axis(dataset.shape, 0, axis_chunk):
> dense_chunk = dataset[idx]
> sub_matrix = sparse.csr_matrix(dense_chunk)
> sub_matrices.append(sub_matrix)
> return sparse.vstack(sub_matrices, format="csr")
> 
> counts = read_dense_as_csr(hf["data"])
> 
> ```

---

<div class="post-metadata">

**Author:** ![salwanbutrus](https://yyz1.discourse-cdn.com/flex035/user_avatar/discourse.scverse.org/salwanbutrus/32/184_2.png) [@salwanbutrus](https://discourse.scverse.org/u/salwanbutrus)\
**Post date:** [March 26, 2022, 11:54pm UTC](https://discourse.scverse.org/t/reading-in-a-1-1-million-cell-hdf5-dataset/375/5 "2022-03-26T23:54:20Z")

</div>

Great, thank you both!