Manage a pathway registry

Background

Pathways in single-cell analysis represent the interconnected networks of molecular signaling cascades that govern critical cellular processes. They are of utmost importance as they offer a comprehensive understanding of the intricate regulatory mechanisms underlying cellular behavior, providing insights into disease pathogenesis, therapeutic responses, and the identification of potential targets for precision medicine and intervention strategies.

Managing pathways across different datasets is crucial in a biotech company to gain a comprehensive understanding of complex biological processes and facilitate efficient research and development.

In this notebook we are registering the 2023 GO Biological process pathway ontology with Lamin. Afterwards, we are linking the pathways to genes and conducting a pathway enrichment analysis on an interferon-beta treated dataset. Finally, we will demonstrate how to fetch datasets with pathway queries using Lamin.

Setup

Warning

Please ensure that you have created or loaded a LaminDB instance before running the remaining part of this notebook!

Show code cell content Hide code cell content

# A lamindb instance containing bionty schema (skip if you already loaded your instance)
!lamin init --storage ./enrichr --schema bionty

import lamindb as ln
import lnschema_bionty as lb

lb.settings.species = "human"  # globally set species

import gseapy as gp
import scanpy as sc
import matplotlib.pyplot as plt

from lamin_usecases import datasets as ds


lb.settings.species = "human"  # globally set species

Fetch GO_Biological_Process_2023 pathways annotated with human genes using Enrichr

First we fetch the “GO_Biological_Process_2023” pathways for humans using GSEApy which wraps GSEA and Enrichr.

go_bp = gp.get_library(name="GO_Biological_Process_2023", organism="Human")
print(f"Number of pathways {len(go_bp)}")

go_bp["ATF6-mediated Unfolded Protein Response (GO:0036500)"]

Parse out the ontology_id from keys, convert into the format of {ontology_id: (name, genes)}

def parse_ontology_id_from_keys(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)

go_bp_parsed = {}

for key, genes in go_bp.items():
    id, name = parse_ontology_id_from_keys(key)
    go_bp_parsed[id] = (name, genes)

go_bp_parsed["GO:0036500"]

Register pathway ontology in LaminDB

pathway_bionty = lb.Pathway.bionty()  # equals to bionty.Pathway()
pathway_bionty

Next, we register all the pathways and genes in LaminDB to finally link pathways to genes.

Register pathway terms

To register the pathways we make use of .from_values to directly parse the annotated GO pathway ontology IDs into LaminDB.

pathway_records = lb.Pathway.from_values(go_bp_parsed.keys(), lb.Pathway.ontology_id)

lb.Pathway.from_bionty(ontology_id="GO:0015868")

ln.save(pathway_records, parents=False)  # not recursing through parents

Register gene symbols

Similarly, we use .from_values for all Pathway associated genes to register them with LaminDB.

all_genes = {g for genes in go_bp.values() for g in genes}

gene_records = lb.Gene.from_values(all_genes, lb.Gene.symbol)

gene_records[:3]

ln.save(gene_records);

Link pathway to genes

Now that we are tracking all pathways and genes records, we can link both of them to make the pathways even more queryable.

gene_records_ids = {record.symbol: record for record in gene_records}

for pathway_record in pathway_records:
    pathway_genes = go_bp_parsed.get(pathway_record.ontology_id)[1]
    pathway_genes_records = [gene_records_ids.get(gene) for gene in pathway_genes]
    pathway_record.genes.set(pathway_genes_records)

Now genes are linked to pathways:

pathway_record.genes.list("symbol")

A interferon-beta treated dataset

We will now conduct a pathway enrichment analysis on a small peripheral blood mononuclear cell dataset that is split into control and stimulated groups. The stimulated group was treated with interferon beta.

The dataset was initially obtained using From "SeuratData::ifnb".

Let’s load the dataset and look at the cell type annotations.

adata = ds.anndata_seurat_ifnb()
adata

adata.obs["seurat_annotations"].value_counts()

For simplicity, we subset to “B Activated” cells:

adata_ba = adata[adata.obs.seurat_annotations == "B Activated"].copy()
adata_ba

Pathway enrichment analysis using Enrichr

This analysis is based on: https://gseapy.readthedocs.io/en/master/singlecell_usecase.html

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_ba,
    groupby="stim",
    use_raw=False,
    method="wilcoxon",
    groups=["STIM"],
    reference="CTRL",
)

rank_genes_groups_df = sc.get.rank_genes_groups_df(adata_ba, "STIM")

rank_genes_groups_df.head()

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

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);

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, size=10);

Track datasets containing annotated pathways in LaminDB

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

ln.track()

We further create a File object to track the dataset.

file = ln.File.from_anndata(
    adata_ba, description="seurat_ifnb_activated_Bcells", var_ref=lb.Gene.symbol
)

ln.save(file)

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, lb.Gene.symbol)

degs_dw_featureset = ln.FeatureSet.from_values(degs_dw.names, lb.Gene.symbol)

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

# 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_keys)
]
enr_dw_top10 = [
    pw_id[0] for pw_id in enr_dw.head(10).Term.apply(parse_ontology_id_from_keys)
]

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

Link feature sets to file:

file.features.add_feature_set(degs_up_featureset, slot="up-DEGs")
file.features.add_feature_set(degs_dw_featureset, slot="down-DEGs")

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")

Querying for pathways

Querying for pathways is now simple with .filter:

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

Query pathways from a gene:

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

Query files from a pathway:

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

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

pathway = lb.Pathway.filter(ontology_id="GO:0035456").one()
pathway

degs = ln.FeatureSet.filter(pathways__ontology_id=pathway.ontology_id).one()

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

pathway_genes = set(pathway.genes.list("symbol"))
degs_genes = set(degs.genes.list("symbol"))

pathway_genes.intersection(degs_genes)

Conclusion

Registering pathways and associated gene sets is made simple with .from_values that ensures that all parsed objects are linked to ontology IDs.Linking both sets is possible with FeatureSet to facilitate simple querying for datasets that contain specific pathways. Since the pathways are linked to genes, Lamin also enables fetching the associated genes of a registered pathway to, for usecase, retrieve sets of differentially expressed genes that are a part of a specific pathway.

Try it yourself

This notebook is available at laminlabs/lamin-usecases.

Show code cell content Hide code cell content

!lamin delete enrichr
!rm -r ./enrichr