Standardize metadata on-the-fly¶

This use cases runs on a LaminDB instance with populated CellType and Pathway registries. Make sure you run the GO Ontology notebook before executing this use case.

Here, we demonstrate how to standardize the metadata on-the-fly during cell type annotation and pathway enrichment analysis using these two registries.

For more information, see:

!lamin load use-cases-registries

Entity has to be a laminhub URL or 'artifact' or 'transform'

import lamindb as ln
import bionty as bt
from lamin_usecases import datasets as ds
import scanpy as sc
import matplotlib.pyplot as plt
import celltypist
import gseapy as gp

→ connected lamindb: testuser1/use-cases-registries

sc.settings.set_figure_params(dpi=50, facecolor="white")

ln.track("hsPU1OENv0LS0000")

→ created Transform('hsPU1OEN'), started new Run('bg8vT3Xj') at 2024-12-20 15:05:36 UTC

→ notebook imports: bionty==0.53.2 celltypist==1.6.3 gseapy==1.1.4 lamin_usecases==0.0.1 lamindb==0.77.3 matplotlib==3.10.0 scanpy==1.10.4

An interferon-beta treated dataset¶

A small peripheral blood mononuclear cell dataset that is split into control and stimulated groups. The stimulated group was treated with interferon beta.

Let’s load the dataset and perform some preprocessing:

adata = ds.anndata_seurat_ifnb(preprocess=False, populate_registries=True)
adata

AnnData object with n_obs × n_vars = 13999 × 9943
    obs: 'stim'
    var: 'symbol'

sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
sc.pp.pca(adata, n_comps=20)
sc.pp.neighbors(adata, n_pcs=10)
sc.tl.umap(adata)

Analysis: cell type annotation using CellTypist¶

model = celltypist.models.Model.load(model="Immune_All_Low.pkl")

Show code cell output Hide code cell output

🔎 No available models. Downloading...

📜 Retrieving model list from server https://celltypist.cog.sanger.ac.uk/models/models.json

📚 Total models in list: 53

📂 Storing models in /home/runner/.celltypist/data/models

💾 Downloading model [1/53]: Immune_All_Low.pkl

💾 Downloading model [2/53]: Immune_All_High.pkl

💾 Downloading model [3/53]: Adult_COVID19_PBMC.pkl

💾 Downloading model [4/53]: Adult_CynomolgusMacaque_Hippocampus.pkl

💾 Downloading model [5/53]: Adult_Human_MTG.pkl

💾 Downloading model [6/53]: Adult_Human_PancreaticIslet.pkl

💾 Downloading model [7/53]: Adult_Human_PrefrontalCortex.pkl

💾 Downloading model [8/53]: Adult_Human_Skin.pkl

💾 Downloading model [9/53]: Adult_Human_Vascular.pkl

💾 Downloading model [10/53]: Adult_Mouse_Gut.pkl

💾 Downloading model [11/53]: Adult_Mouse_OlfactoryBulb.pkl

💾 Downloading model [12/53]: Adult_Pig_Hippocampus.pkl

💾 Downloading model [13/53]: Adult_RhesusMacaque_Hippocampus.pkl

💾 Downloading model [14/53]: Autopsy_COVID19_Lung.pkl

💾 Downloading model [15/53]: COVID19_HumanChallenge_Blood.pkl

💾 Downloading model [16/53]: COVID19_Immune_Landscape.pkl

💾 Downloading model [17/53]: Cells_Adult_Breast.pkl

💾 Downloading model [18/53]: Cells_Fetal_Lung.pkl

💾 Downloading model [19/53]: Cells_Human_Tonsil.pkl

💾 Downloading model [20/53]: Cells_Intestinal_Tract.pkl

💾 Downloading model [21/53]: Cells_Lung_Airway.pkl

💾 Downloading model [22/53]: Developing_Human_Brain.pkl

💾 Downloading model [23/53]: Developing_Human_Gonads.pkl

💾 Downloading model [24/53]: Developing_Human_Hippocampus.pkl

💾 Downloading model [25/53]: Developing_Human_Organs.pkl

💾 Downloading model [26/53]: Developing_Human_Thymus.pkl

💾 Downloading model [27/53]: Developing_Mouse_Brain.pkl

💾 Downloading model [28/53]: Developing_Mouse_Hippocampus.pkl

💾 Downloading model [29/53]: Fetal_Human_AdrenalGlands.pkl

💾 Downloading model [30/53]: Fetal_Human_Pancreas.pkl

💾 Downloading model [31/53]: Fetal_Human_Pituitary.pkl

💾 Downloading model [32/53]: Fetal_Human_Retina.pkl

💾 Downloading model [33/53]: Fetal_Human_Skin.pkl

💾 Downloading model [34/53]: Healthy_Adult_Heart.pkl

💾 Downloading model [35/53]: Healthy_COVID19_PBMC.pkl

💾 Downloading model [36/53]: Healthy_Human_Liver.pkl

💾 Downloading model [37/53]: Healthy_Mouse_Liver.pkl

💾 Downloading model [38/53]: Human_AdultAged_Hippocampus.pkl

💾 Downloading model [39/53]: Human_Colorectal_Cancer.pkl

💾 Downloading model [40/53]: Human_Developmental_Retina.pkl

💾 Downloading model [41/53]: Human_Embryonic_YolkSac.pkl

💾 Downloading model [42/53]: Human_IPF_Lung.pkl

💾 Downloading model [43/53]: Human_Longitudinal_Hippocampus.pkl

💾 Downloading model [44/53]: Human_Lung_Atlas.pkl

💾 Downloading model [45/53]: Human_PF_Lung.pkl

💾 Downloading model [46/53]: Human_Placenta_Decidua.pkl

💾 Downloading model [47/53]: Lethal_COVID19_Lung.pkl

💾 Downloading model [48/53]: Mouse_Dentate_Gyrus.pkl

💾 Downloading model [49/53]: Mouse_Isocortex_Hippocampus.pkl

💾 Downloading model [50/53]: Mouse_Postnatal_DentateGyrus.pkl

💾 Downloading model [51/53]: Mouse_Whole_Brain.pkl

💾 Downloading model [52/53]: Nuclei_Lung_Airway.pkl

💾 Downloading model [53/53]: Pan_Fetal_Human.pkl

predictions = celltypist.annotate(
    adata, model="Immune_All_Low.pkl", majority_voting=True
)
adata.obs["cell_type_celltypist"] = predictions.predicted_labels.majority_voting

🔬 Input data has 13999 cells and 9943 genes

🔗 Matching reference genes in the model

🧬 3699 features used for prediction

⚖️ Scaling input data

🖋️ Predicting labels

✅ Prediction done!

👀 Detected a neighborhood graph in the input object, will run over-clustering on the basis of it

⛓️ Over-clustering input data with resolution set to 10

🗳️ Majority voting the predictions

✅ Majority voting done!

adata.obs["cell_type_celltypist"] = bt.CellType.standardize(
    adata.obs["cell_type_celltypist"]
)

sc.pl.umap(
    adata,
    color=["cell_type_celltypist", "stim"],
    frameon=False,
    legend_fontsize=10,
    wspace=0.4,
)

... storing 'cell_type_celltypist' as categorical

_images/54bf9e5faa046818a0891d3030549bdbb49344c5ab6c935b0f5daa321ed79df4.png

Analysis: Pathway enrichment analysis using Enrichr¶

This analysis is based on the GSEApy scRNA-seq Example.

First, we compute differentially expressed genes using a Wilcoxon test between stimulated and control cells.

# compute differentially expressed genes
sc.tl.rank_genes_groups(
    adata,
    groupby="stim",
    use_raw=False,
    method="wilcoxon",
    groups=["STIM"],
    reference="CTRL",
)

rank_genes_groups_df = sc.get.rank_genes_groups_df(adata, "STIM")
rank_genes_groups_df.head()

	names	scores	logfoldchanges
0	ISG15	99.456627	7.132626
1	ISG20	96.736649	5.074122
2	IFI6	94.972954	5.828700
3	IFIT3	92.482460	7.432328
4	IFIT1	90.699036	8.053429

Next, we filter out up/down-regulated differentially expressed gene sets:

degs_up = rank_genes_groups_df[
    (rank_genes_groups_df["logfoldchanges"] > 0)
    & (rank_genes_groups_df["pvals_adj"] < 0.05)
]
degs_dw = rank_genes_groups_df[
    (rank_genes_groups_df["logfoldchanges"] < 0)
    & (rank_genes_groups_df["pvals_adj"] < 0.05)
]

degs_up.shape, degs_dw.shape

((541, 5), (937, 5))

Run pathway enrichment analysis on DEGs and plot top 10 pathways:

enr_up = gp.enrichr(degs_up.names, gene_sets="GO_Biological_Process_2023").res2d
gp.dotplot(enr_up, figsize=(2, 3), title="Up", cmap=plt.cm.autumn_r);

_images/79a3cb68f8132a838f5128fc51c9832d140b7c55fd460e257ef557bb50f96678.png

enr_dw = gp.enrichr(degs_dw.names, gene_sets="GO_Biological_Process_2023").res2d
gp.dotplot(enr_dw, figsize=(2, 3), title="Down", cmap=plt.cm.winter_r);

_images/3e5f72b4e7027683f6cdbae3eaf3bb9bcc5326d79b4c4c69215412d2c7433099.png

Annotate & save dataset¶

gRegister new features and labels (check out more details here):

new_features = ln.Feature.from_df(adata.obs)
ln.save(new_features)
new_labels = [ln.ULabel(name=i) for i in adata.obs["stim"].unique()]
ln.save(new_labels)

! You have few permissible values for feature stim, consider dtype 'cat' instead of 'str'

! You have few permissible values for feature cell_type_celltypist, consider dtype 'cat' instead of 'str'

features = ln.Feature.lookup()

Register dataset using a Artifact object:

artifact = ln.Artifact.from_anndata(
    adata,
    description="seurat_ifnb_activated_Bcells",
).save()

# TODO: rewrite based on ln.Curator.from_anndata()
# artifact.features._add_set_from_anndata(
#     var_field=bt.Gene.symbol,
#     organism="human",  # optionally, globally set organism via bt.settings.organism = "human"
# )
# cell_type_records = bt.CellType.from_values(adata.obs["cell_type_celltypist"])
# artifact.labels.add(cell_type_records, features.cell_type_celltypist)
# stim_records = ln.ULabel.from_values(adata.obs["stim"])
# artifact.labels.add(stim_records, features.stim)

Link pathway labels¶

Let’s enable tracking of the current notebook as the transform of this artifact:

We further create two feature sets for degs_up and degs_dw which we can later associate with the associated pathways:

degs_up_featureset = ln.FeatureSet.from_values(
    degs_up.names,
    bt.Gene.symbol,
    name="Up-regulated DEGs STIM vs CTRL",
    type="category",
    organism=(  # optionally, globally set organism via bt.settings.organism = "human"
        "human"
    ),
)
degs_dw_featureset = ln.FeatureSet.from_values(
    degs_dw.names,
    bt.Gene.symbol,
    name="Down-regulated DEGs STIM vs CTRL",
    type="category",
    organism=(  # optionally, globally set organism via bt.settings.organism = "human"
        "human"
    ),
)

# Link feature sets to artifact
artifact.features.add_feature_set(degs_up_featureset, slot="STIM-up-DEGs")
artifact.features.add_feature_set(degs_dw_featureset, slot="STIM-down-DEGs")

Link the top 10 pathways to the corresponding differentially expressed genes:

def parse_ontology_id_from_enrichr_results(key):
    """Parse out the ontology id.

    "ATF6-mediated Unfolded Protein Response (GO:0036500)" -> ("GO:0036500", "ATF6-mediated Unfolded Protein Response")
    """
    id = key.split(" ")[-1].replace("(", "").replace(")", "")
    name = key.replace(f" ({id})", "")
    return (id, name)


# get ontology ids for the top 10 pathways
enr_up_top10 = [
    pw_id[0]
    for pw_id in enr_up.head(10).Term.apply(parse_ontology_id_from_enrichr_results)
]
enr_dw_top10 = [
    pw_id[0]
    for pw_id in enr_dw.head(10).Term.apply(parse_ontology_id_from_enrichr_results)
]

# get pathway records
enr_up_top10_pathways = bt.Pathway.from_values(enr_up_top10, bt.Pathway.ontology_id)
enr_dw_top10_pathways = bt.Pathway.from_values(enr_dw_top10, bt.Pathway.ontology_id)

Associate the pathways to the differentially expressed genes:

degs_up_featureset.pathways.set(enr_up_top10_pathways)
degs_dw_featureset.pathways.set(enr_dw_top10_pathways)

degs_up_featureset.pathways.list("name")

['cellular response to cytokine stimulus',
 'defense response to symbiont',
 'defense response to virus',
 'negative regulation of viral genome replication',
 'negative regulation of viral process',
 'positive regulation of cytokine production',
 'regulation of viral genome replication',
 'response to cytokine',
 'response to interferon-beta',
 'response to type II interferon']

Querying pathways¶

Querying for pathways contains “interferon-beta” in the name:

bt.Pathway.filter(name__contains="interferon-beta").df()

	uid	name	ontology_id	abbr	synonyms	description	source_id	run_id	created_at	created_by_id
id
684	1l4z0v8W	cellular response to interferon-beta	GO:0035458	None	cellular response to fiblaferon\|cellular respo...	Any Process That Results In A Change In State ...	86	None	2024-12-20 15:04:53.293956+00:00	1
2130	1NzHDJDi	negative regulation of interferon-beta production	GO:0032688	None	down regulation of interferon-beta production\|...	Any Process That Stops, Prevents, Or Reduces T...	86	None	2024-12-20 15:04:53.352050+00:00	1
3127	3x0xmK1y	positive regulation of interferon-beta production	GO:0032728	None	up-regulation of interferon-beta production\|up...	Any Process That Activates Or Increases The Fr...	86	None	2024-12-20 15:04:53.403539+00:00	1
4334	54R2a0el	regulation of interferon-beta production	GO:0032648	None	regulation of IFN-beta production	Any Process That Modulates The Frequency, Rate...	86	None	2024-12-20 15:04:53.456787+00:00	1
4953	3VZq4dMe	response to interferon-beta	GO:0035456	None	response to fiblaferon\|response to fibroblast ...	Any Process That Results In A Change In State ...	86	None	2024-12-20 15:04:53.484359+00:00	1

Query pathways from a gene:

bt.Pathway.filter(genes__symbol="KIR2DL1").df()

	uid	name	ontology_id	abbr	synonyms	description	source_id	run_id	created_at	created_by_id
id
1346	7S7qlEkG	immune response-inhibiting cell surface recept...	GO:0002767	None	immune response-inhibiting cell surface recept...	The Series Of Molecular Signals Initiated By A...	86	None	2024-12-20 15:04:53.319836+00:00	1

Query artifacts from a pathway:

ln.Artifact.filter(feature_sets__pathways__name__icontains="interferon-beta").first()

Artifact(uid='umqzYDd43eJeS3zf0000', is_latest=True, description='seurat_ifnb_activated_Bcells', suffix='.h5ad', type='dataset', size=215085769, hash='J9oa_JXX4wZ-frPnlk2slL', _hash_type='sha1-fl', _accessor='AnnData', visibility=1, _key_is_virtual=True, storage_id=1, transform_id=1, run_id=1, created_by_id=1, created_at=2024-12-20 15:09:05 UTC)

Query featuresets from a pathway to learn from which geneset this pathway was computed:

pathway = bt.Pathway.get(ontology_id="GO:0035456")
pathway

Pathway(uid='3VZq4dMe', name='response to interferon-beta', ontology_id='GO:0035456', synonyms='response to fiblaferon|response to fibroblast interferon|response to interferon beta', description='Any Process That Results In A Change In State Or Activity Of A Cell Or An Organism (In Terms Of Movement, Secretion, Enzyme Production, Gene Expression, Etc.) As A Result Of An Interferon-Beta Stimulus. Interferon-Beta Is A Type I Interferon.', created_by_id=1, source_id=86, created_at=2024-12-20 15:04:53 UTC)

degs = ln.FeatureSet.get(pathways__ontology_id=pathway.ontology_id)

Now we can get the list of genes that are differentially expressed and belong to this pathway:

contributing_genes = pathway.genes.all() & degs.genes.all()
contributing_genes.list("symbol")

['PLSCR1',
 'IFITM2',
 'IFITM1',
 'XAF1',
 'SHFL',
 'MNDA',
 'OAS1',
 'IFI16',
 'STAT1',
 'AIM2',
 'PNPT1',
 'IFITM3',
 'CALM1',
 'BST2',
 'IRF1']

# clean up test instance
!lamin delete --force use-cases-registries
!rm -r ./use-cases-registries

Show code cell output Hide code cell output

Traceback (most recent call last):
  File "/opt/hostedtoolcache/Python/3.11.11/x64/bin/lamin", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages/rich_click/rich_command.py", line 367, in __call__
    return super().__call__(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages/click/core.py", line 1157, in __call__
    return self.main(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages/rich_click/rich_command.py", line 152, in main
    rv = self.invoke(ctx)
         ^^^^^^^^^^^^^^^^
  File "/opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages/click/core.py", line 1688, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages/click/core.py", line 1434, in invoke
    return ctx.invoke(self.callback, **ctx.params)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages/click/core.py", line 783, in invoke
    return __callback(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages/lamin_cli/__main__.py", line 209, in delete
    return delete(instance, force=force)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages/lamindb_setup/_delete.py", line 102, in delete
    n_objects = check_storage_is_empty(
                ^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages/lamindb_setup/core/upath.py", line 836, in check_storage_is_empty
    raise InstanceNotEmpty(message)
lamindb_setup.core.upath.InstanceNotEmpty: Storage '/home/runner/work/lamin-usecases/lamin-usecases/docs/use-cases-registries/.lamindb' contains 1 objects - delete them prior to deleting the instance