scrna1/6 Jupyter Notebook lamindata Binder

scRNA-seq#

You’ll learn how to manage a growing number of scRNA-seq data batches as a single queryable dataset.

Along the way, you’ll see how to create reports, leverage data lineage, and query statistics of individual data batches stored as files.

Specifically, you will:

  1. read a single .h5ad file as an AnnData and seed a growing dataset with it (scrna1/6, currently page)

  2. append a new data batch (a new .h5ad file) and create a new version of this dataset (scrna2/6)

  3. query & inspect files by metadata individually (scrna3/6)

  4. load the dataset into memory and save analytical results as plots (scrna4/6)

  5. iterate over the dataset, train a model, store a derived representation (scrna5/6)

  6. discuss converting a number of files to a single TileDB SOMA store of the same data (scrna6/6)

Setup#

!lamin init --storage ./test-scrna --schema bionty
Hide code cell output 
✅ saved: User(id='DzTjkKse', handle='testuser1', email='testuser1@lamin.ai', name='Test User1', updated_at=2023-10-10 15:41:26)
✅ saved: Storage(id='YpFBBtr4', root='/home/runner/work/lamin-usecases/lamin-usecases/docs/test-scrna', type='local', updated_at=2023-10-10 15:41:26, created_by_id='DzTjkKse')

💡 loaded instance: testuser1/test-scrna
💡 did not register local instance on hub (if you want, call `lamin register`)


import lamindb as ln
import lnschema_bionty as lb
import pandas as pd

ln.track()

💡 loaded instance: testuser1/test-scrna (lamindb 0.55.2)

💡 notebook imports: lamindb==0.55.2 lnschema_bionty==0.31.2 pandas==1.5.3

💡 Transform(id='Nv48yAceNSh8z8', name='scRNA-seq', short_name='scrna', version='0', type=notebook, updated_at=2023-10-10 15:41:34, created_by_id='DzTjkKse')

💡 Run(id='RNgU2xx7TeUrL3d83b4F', run_at=2023-10-10 15:41:34, transform_id='Nv48yAceNSh8z8', created_by_id='DzTjkKse')

Access #

Let us look at the data of Conde et al., Science (2022).

These data are available in standardized form from the CellxGene data portal.

Here, we’ll use it to seed a growing in-house store of scRNA-seq data managed with the corresponding metadata in LaminDB registries.

Note

If you’re not interested in managing large collections of in-house data and you’d just like to query public data, please take a look at CellxGene census, which exposes all datasets hosted in the data portal as a concatenated TileDB SOMA store.

lb.settings.species = "human"

By calling ln.dev.datasets.anndata_human_immune_cells below, we download the dataset from the CellxGene portal here and pre-populate some LaminDB registries.

adata = ln.dev.datasets.anndata_human_immune_cells(
    populate_registries=True  # this pre-populates registries
)

adata

AnnData object with n_obs × n_vars = 1648 × 36503
    obs: 'donor', 'tissue', 'cell_type', 'assay'
    var: 'feature_is_filtered', 'feature_reference', 'feature_biotype'
    uns: 'cell_type_ontology_term_id_colors', 'default_embedding', 'schema_version', 'title'
    obsm: 'X_umap'

This AnnData is already standardized using the same public ontologies underlying lnschema-bionty, hence, we expect validation to be simple.

Nonetheless, LaminDB focuses on building clean in-house registries

Note

In the next notebook, we’ll look at the more difficult case of a non-standardized dataset that requires curation.

Validate #

Validate genes in .var#

lb.Gene.validate(adata.var.index, lb.Gene.ensembl_gene_id);

❗ 148 terms (0.40%) are not validated for ensembl_gene_id: ENSG00000269933, ENSG00000261737, ENSG00000259834, ENSG00000256374, ENSG00000263464, ENSG00000203812, ENSG00000272196, ENSG00000272880, ENSG00000270188, ENSG00000287116, ENSG00000237133, ENSG00000224739, ENSG00000227902, ENSG00000239467, ENSG00000272551, ENSG00000280374, ENSG00000236886, ENSG00000229352, ENSG00000286601, ENSG00000227021, ...

148 gene identifiers can’t be validated (not currently in the Gene registry). Let’s inspect them to see what to do:

inspector = lb.Gene.inspect(adata.var.index, lb.Gene.ensembl_gene_id)

❗ 148 terms (0.40%) are not validated for ensembl_gene_id: ENSG00000269933, ENSG00000261737, ENSG00000259834, ENSG00000256374, ENSG00000263464, ENSG00000203812, ENSG00000272196, ENSG00000272880, ENSG00000270188, ENSG00000287116, ENSG00000237133, ENSG00000224739, ENSG00000227902, ENSG00000239467, ENSG00000272551, ENSG00000280374, ENSG00000236886, ENSG00000229352, ENSG00000286601, ENSG00000227021, ...

   detected 35 Gene terms in Bionty for ensembl_gene_id: 'ENSG00000275249', 'ENSG00000275063', 'ENSG00000198899', 'ENSG00000198886', 'ENSG00000273554', 'ENSG00000198840', 'ENSG00000276017', 'ENSG00000273748', 'ENSG00000212907', 'ENSG00000278633', 'ENSG00000275869', 'ENSG00000277836', 'ENSG00000278817', 'ENSG00000268674', 'ENSG00000228253', 'ENSG00000277856', 'ENSG00000274175', 'ENSG00000198695', 'ENSG00000277630', 'ENSG00000276256', ...

→  add records from Bionty to your Gene registry via .from_values()

   couldn't validate 113 terms: 'ENSG00000272567', 'ENSG00000270394', 'ENSG00000255823', 'ENSG00000280095', 'ENSG00000225932', 'ENSG00000182230', 'ENSG00000280374', 'ENSG00000254561', 'ENSG00000268955', 'ENSG00000272267', 'ENSG00000224745', 'ENSG00000256222', 'ENSG00000221995', 'ENSG00000261737', 'ENSG00000237133', 'ENSG00000254740', 'ENSG00000270672', 'ENSG00000237838', 'ENSG00000256892', 'ENSG00000227902', ...

→  if you are sure, create new records via ln.Gene() and save to your registry

Logging says 35 of the non-validated ids can be found in the Bionty reference. Let’s register them:

records = lb.Gene.from_values(inspector.non_validated, lb.Gene.ensembl_gene_id)
ln.save(records)

❗ did not create Gene records for 113 non-validated ensembl_gene_ids: 'ENSG00000112096', 'ENSG00000182230', 'ENSG00000203812', 'ENSG00000204092', 'ENSG00000215271', 'ENSG00000221995', 'ENSG00000224739', 'ENSG00000224745', 'ENSG00000225932', 'ENSG00000226377', 'ENSG00000226380', 'ENSG00000226403', 'ENSG00000227021', 'ENSG00000227220', 'ENSG00000227902', 'ENSG00000228139', 'ENSG00000228906', 'ENSG00000229352', 'ENSG00000231575', 'ENSG00000232196', ...

The remaining 113 are legacy IDs, not present in the current Ensembl assembly (e.g. ENSG00000112096).

We’d still like to register them, but won’t dive into the details of converting them from an old Ensembl version to the current one.

validated = lb.Gene.validate(adata.var.index, lb.Gene.ensembl_gene_id)
records = [lb.Gene(ensembl_gene_id=id) for id in adata.var.index[~validated]]
ln.save(records)

❗ 113 terms (0.30%) are not validated for ensembl_gene_id: ENSG00000269933, ENSG00000261737, ENSG00000259834, ENSG00000256374, ENSG00000263464, ENSG00000203812, ENSG00000272196, ENSG00000272880, ENSG00000270188, ENSG00000287116, ENSG00000237133, ENSG00000224739, ENSG00000227902, ENSG00000239467, ENSG00000272551, ENSG00000280374, ENSG00000236886, ENSG00000229352, ENSG00000286601, ENSG00000227021, ...

Now all genes pass validation:

lb.Gene.validate(adata.var.index, lb.Gene.ensembl_gene_id);

Our in-house Gene registry provides rich metadata for each gene measured in the AnnData:

lb.Gene.filter().df().head(10)
symbol stable_id ensembl_gene_id ncbi_gene_ids biotype description synonyms species_id bionty_source_id updated_at created_by_id
id
8W6kKBmqo0wm TAF11L2 None ENSG00000284373 391742 protein_coding TATA-box binding protein associated factor 11 ... uHJU pVGg 2023-10-10 15:41:49 DzTjkKse
PjNo1Yh905M8 PGAP6 None ENSG00000129925 58986 protein_coding post-GPI attachment to proteins 6 [Source:HGNC... TMEM6|TMEM8A|TMEM8|M83 uHJU pVGg 2023-10-10 15:41:49 DzTjkKse
J0pAIXXk0IV2 None None ENSG00000236838 lncRNA novel transcript, antisense to SMG6 uHJU pVGg 2023-10-10 15:41:49 DzTjkKse
Vc0vNMkflIkG PTBP2 None ENSG00000117569 58155 protein_coding polypyrimidine tract binding protein 2 [Source... PTB|BRPTB|PTBLP|NPTB uHJU pVGg 2023-10-10 15:41:49 DzTjkKse
PTQ7E0RWszBX C5orf34-AS1 None ENSG00000248554 lncRNA C5orf34 antisense RNA 1 [Source:HGNC Symbol;Ac... uHJU pVGg 2023-10-10 15:41:49 DzTjkKse
Hh0Gk3vjs5HC B4GALNT4 None ENSG00000182272 338707 protein_coding beta-1,4-N-acetyl-galactosaminyltransferase 4 ... FLJ25045|NGALNAC-T1 uHJU pVGg 2023-10-10 15:41:49 DzTjkKse
ww29UrWYeh9C LINC02958 None ENSG00000265752 lncRNA long intergenic non-protein coding RNA 2958 [S... uHJU pVGg 2023-10-10 15:41:49 DzTjkKse
6bm0vySsNcMU DMD None ENSG00000198947 1756 protein_coding dystrophin [Source:HGNC Symbol;Acc:HGNC:2928] DXS239|DXS269|DXS142|DXS164|DXS270|DXS268|DXS2... uHJU pVGg 2023-10-10 15:41:49 DzTjkKse
qSw8XnsDNynf LINC00706 None ENSG00000281186 100652997 lncRNA long intergenic non-protein coding RNA 706 [So... uHJU pVGg 2023-10-10 15:41:49 DzTjkKse
0XDePDNFMzLw EEF1AKMT1 None ENSG00000150456 221143 protein_coding EEF1A lysine methyltransferase 1 [Source:HGNC ... N6AMT2 uHJU pVGg 2023-10-10 15:41:49 DzTjkKse

There are about 36k genes in the registry, all for species “human”.

lb.Gene.filter().df().shape

(36503, 11)

Validate metadata in .obs#

adata.obs.columns

Index(['donor', 'tissue', 'cell_type', 'assay'], dtype='object')

ln.Feature.validate(adata.obs.columns)

❗ 1 term (25.00%) is not validated for name: donor

array([False,  True,  True,  True])

1 feature is not validated: "donor". Let’s register it:

feature = ln.Feature(name="donor", type="category", registries=[ln.ULabel])
ln.save(feature)

Tip

You can also use features = ln.Feature.from_df(df) to bulk create features with types.

All metadata columns are now validated:

ln.Feature.validate(adata.obs.columns)

array([ True,  True,  True,  True])

Next, let’s validate the corresponding labels of each feature.

Some of the metadata labels can be typed using dedicated registries like CellType:

validated = lb.CellType.validate(adata.obs.cell_type)

❗ received 32 unique terms, 1616 empty/duplicated terms are ignored

❗ 2 terms (6.20%) are not validated for name: germinal center B cell, megakaryocyte

Register non-validated cell types - they can all be loaded from a public ontology through Bionty:

records = lb.CellType.from_values(adata.obs.cell_type[~validated], "name")
ln.save(records)

❗ now recursing through parents: this only happens once, but is much slower than bulk saving

lb.ExperimentalFactor.validate(adata.obs.assay)
lb.Tissue.validate(adata.obs.tissue);

Because we didn’t mount a custom schema that contains a Donor registry, we use the ULabel registry to track donor ids:

ln.ULabel.validate(adata.obs.donor);

❗ received 12 unique terms, 1636 empty/duplicated terms are ignored

❗ 12 terms (100.00%) are not validated for name: D496, 621B, A29, A36, A35, 637C, A52, A37, D503, 640C, A31, 582C

Donor labels are not validated, so let’s register them:

donors = [ln.ULabel(name=name) for name in adata.obs.donor.unique()]
ln.save(donors)

ln.ULabel.validate(adata.obs.donor);

Register #

modalities = ln.Modality.lookup()
experimental_factors = lb.ExperimentalFactor.lookup()
species = lb.Species.lookup()
features = ln.Feature.lookup()

Register data#

When we create a File object from an AnnData, we’ll automatically link its feature sets and get information about unmapped categories:

file = ln.File.from_anndata(
    adata, description="Conde22", field=lb.Gene.ensembl_gene_id, modality=modalities.rna
)

file.save()

The file has the following 2 linked feature sets:

file.features

Features:
  var: FeatureSet(id='XGD5bALjzlxe3dYm2tcM', n=36503, type='number', registry='bionty.Gene', hash='dnRexHCtxtmOU81_EpoJ', updated_at=2023-10-10 15:42:29, modality_id='nG6MZ3aj', created_by_id='DzTjkKse')
    'TAF11L2', 'PGAP6', 'None', 'PTBP2', 'C5orf34-AS1', 'B4GALNT4', 'LINC02958', 'DMD', 'LINC00706', 'EEF1AKMT1', 'None', 'None', 'METRNL', 'MPND', 'NOBOX', 'LINC02706', 'TRIM50', 'IGKV6D-21', 'ZFHX4', 'AHCYL1', ...
  obs: FeatureSet(id='KrYPEOnuTBTRi4WqoelO', n=4, registry='core.Feature', hash='Z0BvFRBSIr9xpTLjV1nb', updated_at=2023-10-10 15:42:35, modality_id='NIjDnou1', created_by_id='DzTjkKse')
    🔗 donor (0, core.ULabel): 
    🔗 tissue (0, bionty.Tissue): 
    🔗 assay (0, bionty.ExperimentalFactor): 
    🔗 cell_type (0, bionty.CellType):

Create a dataset from the file#

dataset = ln.Dataset(file, name="My versioned scRNA-seq dataset", version="1")

dataset

Dataset(id='oYS8NR45oBcEPgCISfCm', name='My versioned scRNA-seq dataset', version='1', hash='6Hu1BywwK6bfIU2Dpku2xZ', transform_id='Nv48yAceNSh8z8', run_id='RNgU2xx7TeUrL3d83b4F', file_id='oYS8NR45oBcEPgCISfCm', created_by_id='DzTjkKse')

Let’s inspect the features measured in this dataset which were inherited from the file:

dataset.features

Features:
  var: FeatureSet(id='XGD5bALjzlxe3dYm2tcM', n=36503, type='number', registry='bionty.Gene', hash='dnRexHCtxtmOU81_EpoJ', updated_at=2023-10-10 15:42:29, modality_id='nG6MZ3aj', created_by_id='DzTjkKse')
    'TAF11L2', 'PGAP6', 'None', 'PTBP2', 'C5orf34-AS1', 'B4GALNT4', 'LINC02958', 'DMD', 'LINC00706', 'EEF1AKMT1', 'None', 'None', 'METRNL', 'MPND', 'NOBOX', 'LINC02706', 'TRIM50', 'IGKV6D-21', 'ZFHX4', 'AHCYL1', ...
  obs: FeatureSet(id='KrYPEOnuTBTRi4WqoelO', n=4, registry='core.Feature', hash='Z0BvFRBSIr9xpTLjV1nb', updated_at=2023-10-10 15:42:35, modality_id='NIjDnou1', created_by_id='DzTjkKse')
    🔗 donor (0, core.ULabel): 
    🔗 tissue (0, bionty.Tissue): 
    🔗 assay (0, bionty.ExperimentalFactor): 
    🔗 cell_type (0, bionty.CellType): 
  external: FeatureSet(id='BHYzxR0GLmF7wf6wZ9Od', n=1, registry='core.Feature', hash='6xnLXE9IYVLlzfQ9kHw6', updated_at=2023-10-10 15:42:36, modality_id='NIjDnou1', created_by_id='DzTjkKse')
    🔗 species (0, bionty.Species):

This looks all good, hence, let’s save it:

dataset.save()

Annotate by linking labels:

dataset.labels.add(experimental_factors.single_cell_rna_sequencing, features.assay)
dataset.labels.add(species.human, features.species)
dataset.labels.add(adata.obs.cell_type, feature=features.cell_type)
dataset.labels.add(adata.obs.assay, feature=features.assay)
dataset.labels.add(adata.obs.tissue, feature=features.tissue)
dataset.labels.add(adata.obs.donor, feature=features.donor)

For this version 1 of the dataset, dataset and file match each other. But they’re independently tracked and queryable through their registries.

dataset.describe()

Dataset(id='oYS8NR45oBcEPgCISfCm', name='My versioned scRNA-seq dataset', version='1', hash='6Hu1BywwK6bfIU2Dpku2xZ', updated_at=2023-10-10 15:42:41)

Provenance:
  💫 transform: Transform(id='Nv48yAceNSh8z8', name='scRNA-seq', short_name='scrna', version='0', type=notebook, updated_at=2023-10-10 15:41:34, created_by_id='DzTjkKse')
  👣 run: Run(id='RNgU2xx7TeUrL3d83b4F', run_at=2023-10-10 15:41:34, transform_id='Nv48yAceNSh8z8', created_by_id='DzTjkKse')
  📄 file: File(id='oYS8NR45oBcEPgCISfCm', suffix='.h5ad', accessor='AnnData', description='Conde22', size=57615999, hash='6Hu1BywwK6bfIU2Dpku2xZ', hash_type='sha1-fl', updated_at=2023-10-10 15:42:41, storage_id='YpFBBtr4', transform_id='Nv48yAceNSh8z8', run_id='RNgU2xx7TeUrL3d83b4F', created_by_id='DzTjkKse')
  👤 created_by: User(id='DzTjkKse', handle='testuser1', email='testuser1@lamin.ai', name='Test User1', updated_at=2023-10-10 15:41:26)
Features:
  var: FeatureSet(id='XGD5bALjzlxe3dYm2tcM', n=36503, type='number', registry='bionty.Gene', hash='dnRexHCtxtmOU81_EpoJ', updated_at=2023-10-10 15:42:29, modality_id='nG6MZ3aj', created_by_id='DzTjkKse')
    'TAF11L2', 'PGAP6', 'None', 'PTBP2', 'C5orf34-AS1', 'B4GALNT4', 'LINC02958', 'DMD', 'LINC00706', 'EEF1AKMT1', 'None', 'None', 'METRNL', 'MPND', 'NOBOX', 'LINC02706', 'TRIM50', 'IGKV6D-21', 'ZFHX4', 'AHCYL1', ...
  obs: FeatureSet(id='KrYPEOnuTBTRi4WqoelO', n=4, registry='core.Feature', hash='Z0BvFRBSIr9xpTLjV1nb', updated_at=2023-10-10 15:42:35, modality_id='NIjDnou1', created_by_id='DzTjkKse')
    🔗 donor (12, core.ULabel): '582C', 'A36', 'D503', 'A37', 'A29', '640C', 'D496', 'A52', 'A35', 'A31', ...
    🔗 tissue (17, bionty.Tissue): 'jejunal epithelium', 'caecum', 'ileum', 'lamina propria', 'thymus', 'duodenum', 'thoracic lymph node', 'skeletal muscle tissue', 'omentum', 'mesenteric lymph node', ...
    🔗 assay (4, bionty.ExperimentalFactor): 'single-cell RNA sequencing', '10x 3' v3', '10x 5' v2', '10x 5' v1'
    🔗 cell_type (32, bionty.CellType): 'naive thymus-derived CD8-positive, alpha-beta T cell', 'dendritic cell, human', 'non-classical monocyte', 'effector memory CD4-positive, alpha-beta T cell', 'megakaryocyte', 'naive thymus-derived CD4-positive, alpha-beta T cell', 'germinal center B cell', 'mast cell', 'alveolar macrophage', 'T follicular helper cell', ...
  external: FeatureSet(id='BHYzxR0GLmF7wf6wZ9Od', n=1, registry='core.Feature', hash='6xnLXE9IYVLlzfQ9kHw6', updated_at=2023-10-10 15:42:36, modality_id='NIjDnou1', created_by_id='DzTjkKse')
    🔗 species (1, bionty.Species): 'human'
Labels:
  🏷️ species (1, bionty.Species): 'human'
  🏷️ tissues (17, bionty.Tissue): 'jejunal epithelium', 'caecum', 'ileum', 'lamina propria', 'thymus', 'duodenum', 'thoracic lymph node', 'skeletal muscle tissue', 'omentum', 'mesenteric lymph node', ...
  🏷️ cell_types (32, bionty.CellType): 'naive thymus-derived CD8-positive, alpha-beta T cell', 'dendritic cell, human', 'non-classical monocyte', 'effector memory CD4-positive, alpha-beta T cell', 'megakaryocyte', 'naive thymus-derived CD4-positive, alpha-beta T cell', 'germinal center B cell', 'mast cell', 'alveolar macrophage', 'T follicular helper cell', ...
  🏷️ experimental_factors (4, bionty.ExperimentalFactor): 'single-cell RNA sequencing', '10x 3' v3', '10x 5' v2', '10x 5' v1'
  🏷️ ulabels (12, core.ULabel): '582C', 'A36', 'D503', 'A37', 'A29', '640C', 'D496', 'A52', 'A35', 'A31', ...

And we can access the file like so:

dataset.file

File(id='oYS8NR45oBcEPgCISfCm', suffix='.h5ad', accessor='AnnData', description='Conde22', size=57615999, hash='6Hu1BywwK6bfIU2Dpku2xZ', hash_type='sha1-fl', updated_at=2023-10-10 15:42:41, storage_id='YpFBBtr4', transform_id='Nv48yAceNSh8z8', run_id='RNgU2xx7TeUrL3d83b4F', created_by_id='DzTjkKse')

dataset.view_flow()
https://d33wubrfki0l68.cloudfront.net/202e5509c355ecafae134932e44ffbbd80469f85/bc1f1/_images/4c7cb6c4bd0ebee78dc8a8b29e4a5998de35908d93f4f6e912980ff1b52f92aa.svg