Gene-set activity scoring primer

# Import tools for writing a temporary GMT file.
from pathlib import Path
from tempfile import TemporaryDirectory

# Import plotting and table tools for inspecting the scores.
import matplotlib.pyplot as plt
import pandas as pd

# Import CyteArc's datastore and gene-set scoring methods.
import cytearc

# Keep warnings visible while hiding routine progress messages.
cytearc.configure_output(level='WARNING', progress=False)

# Download the prepared PBMC store and its saved analysis.
dataset = cytearc.cytebase.connect("cytearc_docs").download_dataset(
    'tenx_5K_pbmc_rnaseq',
    destination='cytearc_datasets',
    zarr=True,
)

# Open the store so scoring can save its results.
ds = cytearc.DataStore(f'{dataset}/data.zarr', nthreads=4)

# Reuse the saved analysis and its selected cells.
run = ds.pipeline.open(label='docs_default')

# Inspect the opened assays and their dimensions.
ds

# Keep the demonstration GMT file in a temporary folder.
input_directory = TemporaryDirectory()

# Choose the filename for the three PBMC signatures.
gmt_path = Path(input_directory.name) / 'pbmc_signatures.gmt'

# Write five target genes for each demonstration signature.
gmt_path.write_text(
    'T_cell\tna\tCD3D\tCD3E\tTRAC\tLTB\tIL7R\n'
    'B_cell\tna\tMS4A1\tCD79A\tCD37\tCD74\tHLA-DRA\n'
    'Myeloid\tna\tLST1\tS100A8\tS100A9\tCTSS\tFCER1G\n',
    encoding='utf-8',
)

# Read the GMT as a table of signature names and target genes.
gene_sets = cytearc.read_gmt(gmt_path)

# Inspect the genes assigned to each signature.
gene_sets

# Normalize gene-name case before checking signature overlap.
available = {str(name).upper() for name in ds.RNA.feats.fetch_all('names')}

# Mark each target whose name occurs in this assay.
matched_targets = gene_sets.assign(matched=gene_sets['target'].str.upper().isin(available))

# Compare matched and total target counts for each signature.
matched_targets.groupby('source')['matched'].agg(['sum', 'count'])

# Reuse the cells selected by the saved analysis.
cell_selection = run['analysis_cell_selection']

# Include every RNA feature in the scoring universe.
all_features = ds.features.universe(from_assay='RNA')

# Show the number of cells and genes used by both methods.
{'cells': int(run.cells.fetch_all('I').sum()), 'genes': ds.RNA.feats.N}

# Aggregate normalized expression for each demonstration signature.
waggr = ds.scores.waggr(
    gene_sets,
    cell_selection,
    features=all_features,
    tmin=3,
)

# Keep the same signature order throughout the comparison.
score_sources = ['T_cell', 'B_cell', 'Myeloid']

# Load only the three signatures compared in this tutorial.
waggr_result = ds.scores.load(waggr, sources=score_sources)

# Materialize the three score columns as a table.
waggr_scores = pd.DataFrame(
    waggr_result.data.compute(),
    columns=list(waggr_result.source_names),
)

# Inspect the range and median of each signature's WAGGR scores.
waggr_scores.describe().loc[['min', '50%', 'max']]

# Score recovery of each signature among the top-ranked genes.
aucell = ds.scores.aucell(
    gene_sets,
    cell_selection,
    features=all_features,
    tmin=3,
)

# Load AUCell scores in the same signature order as WAGGR.
aucell_result = ds.scores.load(aucell, sources=score_sources)

# Materialize the three AUCell score columns as a table.
aucell_scores = pd.DataFrame(
    aucell_result.data.compute(),
    columns=list(aucell_result.source_names),
)

# Inspect the range and median of each signature's AUCell scores.
aucell_scores.describe().loc[['min', '50%', 'max']]

# Read the saved UMAP coordinates in the score table's cell order.
umap = run.cells.to_pandas_dataframe(['umap_1', 'umap_2'])

# Create one panel per signature.
figure, axes = plt.subplots(1, 3, figsize=(12, 4))

# Map each signature's AUCell score onto the same UMAP.
for axis, source in zip(axes, score_sources, strict=True):

    # Use the common zero-to-one scale for every signature.
    points = axis.scatter(
        umap['umap_1'],
        umap['umap_2'],
        c=aucell_scores[source],
        s=3,
        vmin=0,
        vmax=1,
    )

    # Identify the signature shown in the panel.
    axis.set_title(f'{source} AUCell')

    # Label the score scale used to color the cells.
    figure.colorbar(points, ax=axis, label='AUCell score')

# Make room for panel titles and color scales.
figure.tight_layout()

# Display the three signature maps.
plt.show()

# Align the Myeloid signature's scores from both methods by cell.
myeloid_compare = pd.DataFrame(
    {
        'Myeloid_WAGGR': waggr_scores['Myeloid'],
        'Myeloid_AUCell': aucell_scores['Myeloid'],
    }
)

# Compare the distributions without assuming a shared numerical scale.
myeloid_compare.describe()

# Create an axis for the comparison between scoring methods.
figure, axis = plt.subplots(figsize=(4, 4))

# Plot one point per cell using its two Myeloid signature scores.
axis.scatter(
    myeloid_compare['Myeloid_WAGGR'],
    myeloid_compare['Myeloid_AUCell'],
    s=4,
    alpha=0.35,
)

# Label the expression-based score on the horizontal axis.
axis.set_xlabel('Myeloid WAGGR')

# Label the rank-recovery score on the vertical axis.
axis.set_ylabel('Myeloid AUCell')

# Display the comparison between the two scoring methods.
plt.show()