Analyzing KEGG and chembl (PBMC)

From the previous release

This page was recovered from the PIASO v1.1.0 documentation and has not been re-run against the current release. The code is preserved as published; some calls may have changed. The pages under All tutorials are the ones executed against this release.

The original v1.1.0 site is archived at https://genecell.github.io/PIASO/.

Analyzing KEGG and chembl (PBMC)

import piaso
import cosg

import os
import numpy as np
import pandas as pd
import scanpy as sc
import anndata as ad
import logging
from matplotlib import rcParams
import warnings

# modifying the default figure size, use rcParams.
rcParams['figure.figsize'] = 4, 4
rcParams['font.sans-serif'] = "Arial"
rcParams['font.family'] = "Arial"
sc.settings.verbosity = 3
sc.logging.print_header()
sc.set_figure_params(dpi=80,dpi_save=300, color_map='viridis',facecolor='white')

warnings.simplefilter(action='ignore', category=FutureWarning)

stderr

.../site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

adata = sc.read_h5ad('.../SEA-AD_RNA_MTG_subsample_excludeReference_20k_piaso.h5ad')

import gseapy as gp

names = gp.get_library_name()
kegg_libs = [n for n in names if "KEGG" in n]
print(kegg_libs)

['KEGG_2013', 'KEGG_2015', 'KEGG_2016', 'KEGG_2019_Human', 'KEGG_2019_Mouse', 'KEGG_2021_Human', 'KEGG_2026']

# choosing 'KEGG_2021_Human'

kegg = gp.parser.get_library('KEGG_2021_Human')

%%time
score_matrix, gene_set_names = piaso.tl.calculateScoreParallel(adata,
                                                               kegg)

CPU times: user 2.14 s, sys: 991 ms, total: 3.13 s
Wall time: 2min 21s

adata_score = ad.AnnData(
    X=score_matrix,
    obs=adata.obs.copy(),
    var=pd.DataFrame(index=gene_set_names)
)

adata_score.obsm = adata.obsm.copy()

cosg.cosg(
    adata_score,
    key_added='cosg',
    groupby='Subclass',
    n_genes_user=100,
)

# save current adata
adata_score.write_h5ad('.../kegg2021_human_sea20k.h5ad')

adata_score = sc.read_h5ad('.../kegg2021_human_sea20k.h5ad')

cosg_names = pd.DataFrame(adata_score.uns['cosg']['names'])
cosg_scores = pd.DataFrame(adata_score.uns['cosg']['scores'])

import matplotlib.pyplot as plt

names_df = cosg_names.iloc[:5]
scores_df = cosg_scores.iloc[:5]

top_drugs = names_df.values.flatten().tolist()

# Use scanpy dotplot directly
sc.pl.dotplot(
    adata_score,
    var_names=names_df.to_dict(orient='list'),  # grouped by cell type
    groupby='Subclass',
    cmap='Spectral_r',
    standard_scale='var',
    swap_axes=True,
    dendrogram=False,
    show=False
)

## save the plot as pdf
#plt.savefig('.../kegg_sea20k_dotplot.pdf', bbox_inches='tight')

plt.show()

%%time
score_matrix, gene_set_names, p_score_matrix = piaso.tl.calculateScoreParallel(adata,
                                                                               kegg,
                                                                               return_pvals = True)

CPU times: user 2.52 s, sys: 973 ms, total: 3.49 s
Wall time: 2min 22s

p_adata = ad.AnnData(
    X=p_score_matrix,
    obs=adata.obs.copy(),
    var=pd.DataFrame(index=gene_set_names)
)

p_adata.obsm = adata.obsm.copy()

cosg.cosg(
    p_adata,
    key_added='cosg',
    groupby='Subclass',
    n_genes_user=100,
)

p_adata.write_h5ad('.../p_kegg2021_human_sea20k')

p_adata = sc.read_h5ad('.../p_kegg2021_human_sea20k')

cosg_names = pd.DataFrame(p_adata.uns['cosg']['names'])
cosg_scores = pd.DataFrame(p_adata.uns['cosg']['scores'])

names_df = cosg_names.iloc[:5]
scores_df = cosg_scores.iloc[:5]

top_drugs = names_df.values.flatten().tolist()

# Use scanpy dotplot directly
sc.pl.dotplot(
    p_adata,
    var_names=names_df.to_dict(orient='list'),  # grouped by cell type
    groupby='Subclass',
    cmap='Spectral_r',
    standard_scale='var',
    swap_axes=True,
    dendrogram=False,
    show=False
)
#plt.savefig('.../p_kegg2021_human_sea20k_dotplot.pdf', bbox_inches='tight')
plt.show()

# !wget ftp://ftp.sanger.ac.uk/pub/users/kp9/chembl_30_merged_genesymbols_humans.pkl

import warnings
import numpy as np

warnings.filterwarnings('ignore', category=RuntimeWarning)
warnings.simplefilter("ignore")

import pandas as pd
import drug2cell as d2c

# original = pd.read_pickle("chembl_30_merged_genesymbols_humans.pkl
original = pd.read_pickle(".../chembl_30_merged_genesymbols_humans.pkl")

#pChEMBL is -log10() as per https://chembl.gitbook.io/chembl-interface-documentation/frequently-asked-questions/chembl-data-questions#what-is-pchembl
#the threshold values were updated on 21Dec2024
thresholds_dict={
    'none':6, #1uM
    'NHR':7, #100nM
    'GPCR':7, #100nM
    'Ion Channel':5, #10uM
    'Kinase':7.53, #30nM
}

filtered_df = d2c.chembl.filter_activities(
    dataframe=original,
    drug_max_phase=4,
    assay_type='F',
    add_drug_mechanism=True,
    remove_inactive=True,
    include_active=True,
    pchembl_target_column="target_class",
    pchembl_threshold=thresholds_dict
)
print(filtered_df.shape)

(39660, 55)

chembldict = d2c.chembl.create_drug_dictionary(
    filtered_df,
    drug_grouping='ATC_level'
)
new_chembldict = {}
for outer_key, inner_dict in chembldict.items():
    for inner_key, gene_list in inner_dict.items():
        if inner_key in new_chembldict:
            new_chembldict[inner_key] = list(set(new_chembldict[inner_key]) | set(gene_list))
        else:
            new_chembldict[inner_key] = gene_list

%%time

score_matrix, gene_set_names = piaso.tl.calculateScoreParallel(adata, new_chembldict)

CPU times: user 16.3 s, sys: 2.89 s, total: 19.2 s
Wall time: 17min 46s

adata_score = ad.AnnData(
    X=score_matrix,
    obs=adata.obs.copy(),
    var=pd.DataFrame(index=gene_set_names)
)

adata_score.obsm = adata.obsm.copy()

cosg.cosg(
    adata_score,
    key_added='cosg',
    groupby='Subclass',
    n_genes_user=100,
    expressed_pct=0.1,
    remove_lowly_expressed=True,
)

adata_score.write_h5ad('.../chembl_sea20k.h5ad')

adata_score = sc.read_h5ad('.../chembl_sea20k.h5ad')

cosg_names = pd.DataFrame(adata_score.uns['cosg']['names'])
cosg_scores = pd.DataFrame(adata_score.uns['cosg']['scores'])

names_df = cosg_names.iloc[:5]
scores_df = cosg_scores.iloc[:5]

top_drugs = names_df.values.flatten().tolist()

# Use scanpy dotplot directly
sc.pl.dotplot(
    adata_score,
    var_names=names_df.to_dict(orient='list'),  # grouped by cell type
    groupby='Subclass',
    cmap='Spectral_r',
    standard_scale='var',
    swap_axes=True,
    dendrogram=False,
    show=False
)

plt.savefig('.../chembl_sea20k_dotplot.pdf', bbox_inches='tight')
plt.show()

%%time
score_matrix, gene_set_names, p_score_matrix = piaso.tl.calculateScoreParallel(adata,
                                                                               new_chembldict,
                                                                               return_pvals = True)

CPU times: user 19.6 s, sys: 2.42 s, total: 22 s
Wall time: 18min

p_adata = ad.AnnData(
    X=p_score_matrix,
    obs=adata.obs.copy(),
    var=pd.DataFrame(index=gene_set_names)
)

p_adata.obsm = adata.obsm.copy()

cosg.cosg(
    p_adata,
    key_added='cosg',
    groupby='Subclass',
    n_genes_user=100,
    expressed_pct=0.1,
    remove_lowly_expressed=True,
)

# Save current adata

p_adata.write_h5ad('.../p_chembl_sea20k.h5ad')

p_adata = sc.read_h5ad('.../p_chembl_sea20k.h5ad')

cosg_names = pd.DataFrame(p_adata.uns['cosg']['names'])
cosg_scores = pd.DataFrame(p_adata.uns['cosg']['scores'])

names_df = cosg_names.iloc[:5]
scores_df = cosg_scores.iloc[:5]

top_drugs = names_df.values.flatten().tolist()

# Use scanpy dotplot directly
sc.pl.dotplot(
    p_adata,
    var_names=names_df.to_dict(orient='list'),  # grouped by cell type
    groupby='Subclass',
    cmap='Spectral_r',
    standard_scale='var',
    swap_axes=True,
    dendrogram=False,
    show=False
)

# plt.savefig('.../p_chembl_sea20k_dotplot.pdf', bbox_inches='tight')
plt.show()