scrna2/6 Jupyter Notebook lamindata

Append a new batch of data#

Here, we’ll learn

  • how to standardize a less well curated dataset

  • how to append it as a new batch of data to the growing versioned dataset

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

ln.track()

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

💡 notebook imports: lamindb==0.57.2 lnschema_bionty==0.33.0 pandas==1.5.3

💡 Transform(uid='ManDYgmftZ8Cz8', name='Append a new batch of data', short_name='scrna2', version='0', type=notebook, updated_at=2023-10-23 17:47:41, created_by_id=1)

💡 Run(uid='ZDcOJl5E15UoqnTxZJ7t', run_at=2023-10-23 17:47:41, transform_id=2, created_by_id=1)

Standardize a less-well curated dataset#

Access #

Let’s now consider a dataset with less-well curated features:

adata = ln.dev.datasets.anndata_pbmc68k_reduced()
adata

AnnData object with n_obs × n_vars = 70 × 765
    obs: 'cell_type', 'n_genes', 'percent_mito', 'louvain'
    var: 'n_counts', 'highly_variable'
    uns: 'louvain', 'louvain_colors', 'neighbors', 'pca'
    obsm: 'X_pca', 'X_umap'
    varm: 'PCs'
    obsp: 'connectivities', 'distances'

We see that this dataset is indexed by gene symbols. Because we assume that in-house, we index all datasets by Ensembl IDs, we’ll need to re-curate:

adata.var.head()
n_counts highly_variable
index
HES4 1153.387451 True
TNFRSF4 304.358154 True
SSU72 2530.272705 False
PARK7 7451.664062 False
RBP7 272.811035 True

We are still working with human data, and can globally instruct bionty to assume human:

lb.settings.organism = "human"

Validate #

Curate & validate genes#

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

❗ 70 terms (9.20%) are not validated for symbol: ATPIF1, C1orf228, CCBL2, RP11-782C8.1, RP11-277L2.3, RP11-156E8.1, AC079767.4, GPX1, H1FX, SELT, ATP5I, IGJ, CCDC109B, FYB, H2AFY, FAM65B, HIST1H4C, HIST1H1E, ZNRD1, C6orf48, ...

lb.Gene.inspect(adata.var.index, lb.Gene.symbol);

❗ 70 terms (9.20%) are not validated for symbol: ATPIF1, C1orf228, CCBL2, RP11-782C8.1, RP11-277L2.3, RP11-156E8.1, AC079767.4, GPX1, H1FX, SELT, ATP5I, IGJ, CCDC109B, FYB, H2AFY, FAM65B, HIST1H4C, HIST1H1E, ZNRD1, C6orf48, ...

   detected 54 terms with synonyms: ATPIF1, C1orf228, CCBL2, AC079767.4, H1FX, SELT, ATP5I, IGJ, CCDC109B, FYB, H2AFY, FAM65B, HIST1H4C, HIST1H1E, ZNRD1, C6orf48, SEPT7, WBSCR22, RSBN1L-AS1, CCDC132, ...

→  standardize terms via .standardize()

   detected 5 Gene terms in Bionty for symbol: 'RN7SL1', 'GPX1', 'SNORD3B-2', 'SOD2', 'IGLL5'

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

   couldn't validate 11 terms: 'RP11-291B21.2', 'RP11-277L2.3', 'RP3-467N11.1', 'TMBIM4-1', 'RP11-156E8.1', 'RP11-390E23.6', 'AC084018.1', 'RP11-782C8.1', 'CTD-3138B18.5', 'RP11-620J15.3', 'RP11-489E7.4'

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

Standardize symbols and register additional symbols from Bionty:

adata.var.index = lb.Gene.standardize(adata.var.index, lb.Gene.symbol)
gene_records = lb.Gene.from_values(adata.var.index, lb.Gene.symbol)
ln.save(gene_records)

❗ did not create Gene records for 11 non-validated symbols: 'AC084018.1', 'CTD-3138B18.5', 'RP11-156E8.1', 'RP11-277L2.3', 'RP11-291B21.2', 'RP11-390E23.6', 'RP11-489E7.4', 'RP11-620J15.3', 'RP11-782C8.1', 'RP3-467N11.1', 'TMBIM4-1'

We only want to register data with validated genes: data related to other features wouldn’t be useful to us, anyway.

Hence, we subset the AnnData object to the validated genes:

validated = lb.Gene.validate(adata.var.index, lb.Gene.symbol)
adata_validated = adata[:, validated].copy()

❗ 11 terms (1.40%) are not validated for symbol: RP11-782C8.1, RP11-277L2.3, RP11-156E8.1, RP3-467N11.1, RP11-390E23.6, RP11-489E7.4, RP11-291B21.2, RP11-620J15.3, TMBIM4-1, AC084018.1, CTD-3138B18.5

We also subset raw of the anndata object to the validated genes

adata_validated.raw = adata.raw[:, validated].to_adata()

Now, we need to convert gene symbols into ensembl gene ids:

records = lb.Gene.filter(id__in=[record.id for record in gene_records])
mapper = pd.DataFrame(records.values_list("symbol", "ensembl_gene_id")).set_index(0)[1]
adata_validated.var.insert(0, "gene_symbol", adata_validated.var.index)
adata_validated.var.rename(index=mapper, inplace=True)

adata_validated.var.head()
gene_symbol n_counts highly_variable
ENSG00000188290 HES4 1153.387451 True
ENSG00000186827 TNFRSF4 304.358154 True
ENSG00000160075 SSU72 2530.272705 False
ENSG00000116288 PARK7 7451.664062 False
ENSG00000162444 RBP7 272.811035 True

Raw has the same genes, so set them also.

adata_validated.raw.var.index = adata_validated.var.index

Curate & validate cell types#

Inspection shows none of the terms are validated:

inspector = lb.CellType.inspect(adata_validated.obs.cell_type)

❗ received 9 unique terms, 61 empty/duplicated terms are ignored

❗ 9 terms (100.00%) are not validated for name: Dendritic cells, CD19+ B, CD4+/CD45RO+ Memory, CD8+ Cytotoxic T, CD4+/CD25 T Reg, CD14+ Monocytes, CD56+ NK, CD8+/CD45RA+ Naive Cytotoxic, CD34+

   couldn't validate 9 terms: 'CD4+/CD45RO+ Memory', 'CD4+/CD25 T Reg', 'CD56+ NK', 'CD34+', 'CD19+ B', 'CD8+ Cytotoxic T', 'Dendritic cells', 'CD8+/CD45RA+ Naive Cytotoxic', 'CD14+ Monocytes'

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

Let us search the cell type names from the public ontology, and add the name value found in the AnnData object as a synonym to the top match found in the public ontology.

bionty = lb.CellType.bionty()  # access the public ontology through bionty
name_mapper = {}
for name in adata_validated.obs.cell_type.unique():
    ontology_id = (
        bionty.search(name).iloc[0].ontology_id
    )  # search the public ontology and use the ontology id of the top match
    record = lb.CellType.from_bionty(
        ontology_id=ontology_id
    )  # create a record by loading the top match from bionty
    name_mapper[name] = record.name  # map the original name to standardized name
    record.save()  # save the record
    record.add_synonym(
        name
    )  # add the original name as a synonym, so that next time, we can just run .standardize()

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

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

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

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

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

We can now standardize cell type names using the search-based mapper:

adata_validated.obs.cell_type = adata_validated.obs.cell_type.map(name_mapper)

Now, all cell types are validated:

validated = lb.CellType.validate(adata_validated.obs.cell_type)
assert all(validated)

We don’t want to store any of the other metadata columns:

for column in ["n_genes", "percent_mito", "louvain"]:
    adata.obs.drop(column, axis=1)

Register #

modalities = ln.Modality.lookup()
experimental_factors = lb.ExperimentalFactor.lookup()
organism = lb.Organism.lookup()
features = ln.Feature.lookup()

file = ln.File.from_anndata(
    adata_validated,
    description="10x reference adata",
    field=lb.Gene.ensembl_gene_id,
    modality=modalities.rna,
)

❗    3 terms (75.00%) are not validated for name: n_genes, percent_mito, louvain

As we do not want to manage the remaining unvalidated terms in registries, we can save the file.

file.save()

file.labels.add(adata_validated.obs.cell_type, features.cell_type)
file.labels.add(organism.human, feature=features.organism)
file.labels.add(experimental_factors.single_cell_rna_sequencing, feature=features.assay)

file.describe()

File(uid='LnH0qcPKHapOHPwc6LrM', suffix='.h5ad', accessor='AnnData', description='10x reference adata', size=857752, hash='U4UmKxr_rN8KlwVK4WYKjw', hash_type='md5', updated_at=2023-10-23 17:48:08)

Provenance:
  🗃️ storage: Storage(uid='ZIqsTwSZ', root='/home/runner/work/lamin-usecases/lamin-usecases/docs/test-scrna', type='local', updated_at=2023-10-23 17:46:37, created_by_id=1)
  💫 transform: Transform(uid='ManDYgmftZ8Cz8', name='Append a new batch of data', short_name='scrna2', version='0', type=notebook, updated_at=2023-10-23 17:47:41, created_by_id=1)
  👣 run: Run(uid='ZDcOJl5E15UoqnTxZJ7t', run_at=2023-10-23 17:47:41, transform_id=2, created_by_id=1)
  👤 created_by: User(uid='DzTjkKse', handle='testuser1', name='Test User1', updated_at=2023-10-23 17:46:37)
Features:
  var: FeatureSet(uid='JP1OvnEx2ITWEjurSjHk', n=754, type='number', registry='bionty.Gene', hash='WMDxN7253SdzGwmznV5d', updated_at=2023-10-23 17:48:08, modality_id=1, created_by_id=1)
    'IL18', 'NPM3', 'S100A9', 'S100A8', 'CNN2', 'ARHGAP45', 'RNF34', 'GPX4', 'S100A6', 'ADISSP', 'S100A4', 'FAM174C', 'SIT1', 'CCDC107', 'RSL1D1', 'TLN1', 'HES4', 'TNFRSF17', 'PCNA', 'RAB13', ...
  obs: FeatureSet(uid='kq35k5DEtE4DkfG02hgo', n=1, registry='core.Feature', hash='iuyPNK3smIuvcnCIlFSV', updated_at=2023-10-23 17:48:08, modality_id=2, created_by_id=1)
    🔗 cell_type (9, bionty.CellType): 'B cell, CD19-positive', 'dendritic cell', 'CD24-positive, CD4 single-positive thymocyte', 'CD8-positive, CD25-positive, alpha-beta regulatory T cell', 'CD16-positive, CD56-dim natural killer cell, human', 'gamma-delta T cell', 'monocyte', 'cytotoxic T cell', 'CD4-positive, alpha-beta T cell'
  external: FeatureSet(uid='4tILjB6W30HNtBGOI5uS', n=2, registry='core.Feature', hash='LNTvx6NLtiHp0Nw8h_T5', updated_at=2023-10-23 17:48:08, modality_id=2, created_by_id=1)
    🔗 assay (1, bionty.ExperimentalFactor): 'single-cell RNA sequencing'
    🔗 organism (1, bionty.Organism): 'human'
Labels:
  🏷️ organism (1, bionty.Organism): 'human'
  🏷️ cell_types (9, bionty.CellType): 'B cell, CD19-positive', 'dendritic cell', 'CD24-positive, CD4 single-positive thymocyte', 'CD8-positive, CD25-positive, alpha-beta regulatory T cell', 'CD16-positive, CD56-dim natural killer cell, human', 'gamma-delta T cell', 'monocyte', 'cytotoxic T cell', 'CD4-positive, alpha-beta T cell'
  🏷️ experimental_factors (1, bionty.ExperimentalFactor): 'single-cell RNA sequencing'

file.view_flow()
_images/8d253affb2d099da2f03f128551e3b13244e0057db5f5bd2b453fa7206890193.svg

Append a file to the growing dataset#

Query the previous dataset:

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

Create a new version of the dataset by sharding it across the new file file and the file underlying version 1 of the dataset:

dataset_v2 = ln.Dataset(
    [file, dataset_v1.file],
    is_new_version_of=dataset_v1,
)
dataset_v2.save()

# annotate the dataset
dataset_v2.labels.add_from(file)
dataset_v2.labels.add_from(dataset_v1)

Version 2 of the dataset covers significantly more conditions.

dataset_v2.describe()

Dataset(uid='V8nTN6y8N38rnSsSMk70', name='My versioned scRNA-seq dataset', version='2', hash='xTgBfXxl15bqEVzTaMxs', updated_at=2023-10-23 17:48:11)

Provenance:
  💫 transform: Transform(uid='ManDYgmftZ8Cz8', name='Append a new batch of data', short_name='scrna2', version='0', type=notebook, updated_at=2023-10-23 17:47:41, created_by_id=1)
  👣 run: Run(uid='ZDcOJl5E15UoqnTxZJ7t', run_at=2023-10-23 17:47:41, transform_id=2, created_by_id=1)
  🔖 initial_version: Dataset(uid='KHDuUmCEv6SwM8UNQulB', name='My versioned scRNA-seq dataset', version='1', hash='6Hu1BywwK6bfIU2Dpku2xZ', updated_at=2023-10-23 17:47:33, transform_id=1, run_id=1, file_id=1, created_by_id=1)
  👤 created_by: User(uid='DzTjkKse', handle='testuser1', name='Test User1', updated_at=2023-10-23 17:46:37)
Features:
  var: FeatureSet(uid='aB5o5Aj22CBiUsXkn7j0', n=37257, type='number', registry='bionty.Gene', hash='bWfNZ3hy-yGnPj65T3kc', updated_at=2023-10-23 17:48:10, created_by_id=1)
    'MIR1302-2HG', 'FAM138A', 'OR4F5', 'None', 'None', 'None', 'None', 'None', 'None', 'None', 'OR4F29', 'None', 'OR4F16', 'None', 'LINC01409', 'FAM87B', 'LINC01128', 'LINC00115', 'FAM41C', 'None', ...
  obs: FeatureSet(uid='YMf0jITlhFBIzo97KFNP', n=4, registry='core.Feature', hash='sTKN4H5ua2p0aUB7Qf41', updated_at=2023-10-23 17:47:28, modality_id=2, created_by_id=1)
    🔗 cell_type (39, bionty.CellType): 'B cell, CD19-positive', 'dendritic cell', 'CD24-positive, CD4 single-positive thymocyte', 'CD8-positive, CD25-positive, alpha-beta regulatory T cell', 'CD16-positive, CD56-dim natural killer cell, human', 'gamma-delta T cell', 'monocyte', 'cytotoxic T cell', 'CD4-positive, alpha-beta T cell', 'classical monocyte', ...
    🔗 assay (4, bionty.ExperimentalFactor): 'single-cell RNA sequencing', '10x 3' v3', '10x 5' v2', '10x 5' v1'
    🔗 tissue (17, bionty.Tissue): 'blood', 'thoracic lymph node', 'spleen', 'lung', 'mesenteric lymph node', 'lamina propria', 'liver', 'jejunal epithelium', 'omentum', 'bone marrow', ...
    🔗 donor (12, core.ULabel): 'D496', '621B', 'A29', 'A36', 'A35', '637C', 'A52', 'A37', 'D503', '640C', ...
  external: FeatureSet(uid='4tILjB6W30HNtBGOI5uS', n=2, registry='core.Feature', hash='LNTvx6NLtiHp0Nw8h_T5', updated_at=2023-10-23 17:48:08, modality_id=2, created_by_id=1)
    🔗 assay (4, bionty.ExperimentalFactor): 'single-cell RNA sequencing', '10x 3' v3', '10x 5' v2', '10x 5' v1'
    🔗 organism (1, bionty.Organism): 'human'
Labels:
  🏷️ organism (1, bionty.Organism): 'human'
  🏷️ tissues (17, bionty.Tissue): 'blood', 'thoracic lymph node', 'spleen', 'lung', 'mesenteric lymph node', 'lamina propria', 'liver', 'jejunal epithelium', 'omentum', 'bone marrow', ...
  🏷️ cell_types (39, bionty.CellType): 'B cell, CD19-positive', 'dendritic cell', 'CD24-positive, CD4 single-positive thymocyte', 'CD8-positive, CD25-positive, alpha-beta regulatory T cell', 'CD16-positive, CD56-dim natural killer cell, human', 'gamma-delta T cell', 'monocyte', 'cytotoxic T cell', 'CD4-positive, alpha-beta T cell', 'classical monocyte', ...
  🏷️ experimental_factors (4, bionty.ExperimentalFactor): 'single-cell RNA sequencing', '10x 3' v3', '10x 5' v2', '10x 5' v1'
  🏷️ ulabels (12, core.ULabel): 'D496', '621B', 'A29', 'A36', 'A35', '637C', 'A52', 'A37', 'D503', '640C', ...

View the flow:

dataset_v2.view_flow()
_images/6b23a4543212ebe2b7e1512d2eef64a4d33f75a82f94d18cb4d276eb95cc5102.svg