Scanpy (python) tutorial¶

# Read in test dataset - subsetted to existing features
with resources.path("ouroboros.data", "neurogen_example.h5ad") as path:
    adata = sc.read_h5ad(path)
    
adata

# See if any genes are missing from the annadata 
missing_genes = obo.check_genes(adata)

# Continue with nomal scanpy filtering

# mitochondrial genes, "MT-" for human, "Mt-" for mouse
adata.var["mt"] = adata.var_names.str.startswith("MT-")
# ribosomal genes
adata.var["ribo"] = adata.var_names.str.startswith(("RPS", "RPL"))
# hemoglobin genes
adata.var["hb"] = adata.var_names.str.contains("^HB[^(P)]")

# Calculate QC metrics with scanpy 
sc.pp.calculate_qc_metrics(
    adata, qc_vars=["mt", "ribo", "hb"], inplace=True, percent_top=[20], log1p=True
)

# Plot QC metrics
p1 = sns.displot(adata.obs["total_counts"], bins=100, kde=False)
p2 = sc.pl.violin(adata, "pct_counts_mt")
p3 = sc.pl.scatter(adata, "total_counts", "n_genes_by_counts", color="pct_counts_mt")

sc.pp.filter_cells(adata, min_genes=20)
# Usually you would use Scanpy to filter genes with sc.pp.filter_genes(), but if any of our feature genes have low counts they will be removed and then the model will retrain
# If you want Ouroboros to be as accurate as possible, it's best to leave all feature genes in the matrix, even if they have low counts 

# Alternative way to filter genes to keep feature genes present:
# Read in feature genes
feature_genes = obo.read_in_features()
# Count how many cells each gene is expressed in
gene_counts = np.array((adata.X > 0).sum(axis=0)).flatten()
# Apply min_cells threshold
gene_mask = gene_counts >= 3
# Get gene names
gene_names = adata.var_names.to_numpy()
# Force-keep genes in keep_genes
gene_mask |= np.isin(gene_names, list(feature_genes))
# Subset the AnnData
adata = adata[:, gene_mask].copy()
adata

# Then you can run Ouroboros on the saved h5ad 
z_df = obo.run_ouroboros('output_example/example.h5ad', 'h5ad', species = 'human', outdir = '.')

# This is just for rendering in the wiki 
IFrame(src="../output_example/ouroboros_KNN_sphere.html", width="100%", height=600)

# This is just for rendering in the wiki 
IFrame(src="../output_example/ouroboros_KNN_sphere_no_refembed.html", width="100%", height=600)

IFrame(src="../output_example/example_sphere_ct.html", width="100%", height=600)

#obo.plot_gene_sphere(z_df, adata, gene_name = 'TOP2A', ref = 'default', show = True, outpath = 'output_example/example_sphere_top2a.html')

IFrame(src="../output_example/example_sphere_top2a.html", width="100%", height=600)

# Make your snapshot
obo.sphere_snapshot(0, 0, z_df = z_df, colour_by='dormancy_pseudotime', radius = 1.4, ref_embed = 'default', vel_df = None, save_as_png='output_example/dormancy_depth.png')

# Make plot with discrete variable 
# set palette
phase_pal_transition = {
    'G1':'#1f77b4',
       'S': '#ff7f0e',
       'G2M': '#2ca02c', 
       'G0': 'black', 
       'G1-G0 transition': '#d62728'}

# plot
obo.plot_robinson_projection(
    z_df=z_df,
    colour_by='KNN_phase',
    palette=phase_pal_transition,
    ref='default',
    title="Robinson",
    scale=15,
    save_fig="output_example/robinson.png"
)

# Make plot with continous variable 
obo.plot_robinson_projection(
    z_df=z_df,
    colour_by='cell_cycle_pseudotime',
    palette="rocket",
    ref=None,
    title="Robinson",
    scale=15,
    save_fig="output_example/robinson.png")

# Code from Scanpy's clustering tutorial [https://scanpy.readthedocs.io/en/stable/tutorials/basics/clustering.html]

# Normalizing to median total counts
sc.pp.normalize_total(adata)
# Logarithmize the data
sc.pp.log1p(adata)

# Feature selection
sc.pp.highly_variable_genes(adata)

#PCA
sc.tl.pca(adata)

sc.pl.pca(
    adata,
    color=["cell_type", "KNN_phase", 'cell_cycle_pseudotime', 'dormancy_pseudotime'],
    ncols=2,
    wspace=0.4
)

# Nearest neighbour graph construction and visualization 
sc.pp.neighbors(adata)

sc.tl.umap(adata)

sc.pl.umap(
    adata,
    color=["cell_type", "KNN_phase", 'cell_cycle_pseudotime', 'dormancy_pseudotime'],
    ncols=2,
    wspace=0.4
)

from matplotlib.colors import Normalize
import seaborn as sns
import matplotlib.pyplot as plt

# Colormaps and ranges (same as the Robinson plot)
cmap_dorm = sns.color_palette('crest_r', as_cmap=True)
cmap_cyc  = sns.color_palette('flare',   as_cmap=True)
norm_dorm = Normalize(vmin=-1, vmax=0)
norm_cyc  = Normalize(vmin=0,  vmax=1)

# Split cells: whichever pseudotime each cell has
non_cyc = adata[adata.obs['dormancy_pseudotime'].notna()]
cyc     = adata[adata.obs['cell_cycle_pseudotime'].notna()]

fig, ax = plt.subplots(figsize=(6, 6))

# Background: all cells in light grey so the full UMAP shape is visible
sc.pl.umap(adata, ax=ax, show=False, size=10, alpha=0.3, frameon=False)
# (grey is the default when no colour is given)

sc.pl.umap(non_cyc, color='dormancy_pseudotime', ax=ax, show=False,
           cmap=cmap_dorm, vmin=-1, vmax=0, alpha=0.9,
           colorbar_loc=None, frameon=False, title='')
sc.pl.umap(cyc, color='cell_cycle_pseudotime', ax=ax, show=False,
           cmap=cmap_cyc, vmin=0, vmax=1, alpha=0.9,
           colorbar_loc=None, frameon=False, title='Ouroboros pseudotimes')

# Colourbars underneath
fig.subplots_adjust(bottom=0.15)

sm2 = plt.cm.ScalarMappable(cmap=cmap_dorm, norm=norm_dorm)
cbar2_ax = fig.add_axes([0.15, 0.06, 0.33, 0.02])  # [left, bottom, width, height]
cbar2 = fig.colorbar(sm2, cax=cbar2_ax, orientation='horizontal')
cbar2.ax.set_title('Φ', fontsize=10, pad=6)
cbar2.set_ticks([-0.2, -0.4, -0.6, -0.8])

sm1 = plt.cm.ScalarMappable(cmap=cmap_cyc, norm=norm_cyc)
cbar1_ax = fig.add_axes([0.55, 0.06, 0.33, 0.02])
cbar1 = fig.colorbar(sm1, cax=cbar1_ax, orientation='horizontal')
cbar1.ax.set_title('θ', fontsize=10, pad=6)

plt.show()