Introduction

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/.

Introduction

import numpy as np
import pandas as pd
import scanpy as sc
import os
import logging
from matplotlib import rcParams
from sklearn.model_selection import StratifiedShuffleSplit
from sklearn import metrics
from sklearn.metrics import f1_score
import seaborn as sns
import matplotlib.pyplot as plt

path = '.../Analysis/Python/Packages/PIASO'
import sys
sys.path.append(path)
path = '.../Analysis/Python/Packages/COSG'
import sys
sys.path.append(path)
import piaso
import cosg

stderr

.../site-packages/networkx/utils/backends.py:135: RuntimeWarning: networkx backend defined more than once: nx-loopback
  backends.update(_get_backends("networkx.backends"))

sc.set_figure_params(dpi=96,dpi_save=300, color_map='viridis',facecolor='white')
rcParams['font.sans-serif'] = "Arial"
rcParams['font.family'] = "Arial"
sc.settings.verbosity = 3
sc.logging.print_header()

scanpy==1.10.3 anndata==0.10.8 umap==0.5.7 numpy==1.26.4 scipy==1.13.0 pandas==2.2.3 scikit-learn==1.5.2 statsmodels==0.14.4 igraph==0.11.5 louvain==0.8.2 pynndescent==0.5.13

!.../gdrive files download --overwrite --destination .../Data/Public/PIASO 1EdRA0ECvPlEnNaOzKqj19GrEtucB7hmE

Downloading SEA-AD_RNA_MTG_subsample_excludeReference_20k_piaso.h5ad
Successfully downloaded SEA-AD_RNA_MTG_subsample_excludeReference_20k_piaso.h5ad

adata=sc.read('.../Data/Public/PIASO/SEA-AD_RNA_MTG_subsample_excludeReference_20k_piaso.h5ad')
adata

AnnData object with n_obs × n_vars = 20000 × 36601
    obs: 'sample_id', 'Neurotypical reference', 'Donor ID', 'Organism', 'Brain Region', 'Sex', 'Gender', 'Age at Death', 'Race (choice=White)', 'Race (choice=Black/ African American)', 'Race (choice=Asian)', 'Race (choice=American Indian/ Alaska Native)', 'Race (choice=Native Hawaiian or Pacific Islander)', 'Race (choice=Unknown or unreported)', 'Race (choice=Other)', 'specify other race', 'Hispanic/Latino', 'Highest level of education', 'Years of education', 'PMI', 'Fresh Brain Weight', 'Brain pH', 'Overall AD neuropathological Change', 'Thal', 'Braak', 'CERAD score', 'Overall CAA Score', 'Highest Lewy Body Disease', 'Total Microinfarcts (not observed grossly)', 'Total microinfarcts in screening sections', 'Atherosclerosis', 'Arteriolosclerosis', 'LATE', 'Cognitive Status', 'Last CASI Score', 'Interval from last CASI in months', 'Last MMSE Score', 'Interval from last MMSE in months', 'Last MOCA Score', 'Interval from last MOCA in months', 'APOE Genotype', 'Primary Study Name', 'Secondary Study Name', 'NeuN positive fraction on FANS', 'RIN', 'cell_prep_type', 'facs_population_plan', 'rna_amplification', 'sample_name', 'sample_quantity_count', 'expc_cell_capture', 'method', 'pcr_cycles', 'percent_cdna_longer_than_400bp', 'rna_amplification_pass_fail', 'amplified_quantity_ng', 'load_name', 'library_prep', 'library_input_ng', 'r1_index', 'avg_size_bp', 'quantification_fmol', 'library_prep_pass_fail', 'exp_component_vendor_name', 'batch_vendor_name', 'experiment_component_failed', 'alignment', 'Genome', 'ar_id', 'bc', 'GEX_Estimated_number_of_cells', 'GEX_number_of_reads', 'GEX_sequencing_saturation', 'GEX_Mean_raw_reads_per_cell', 'GEX_Q30_bases_in_barcode', 'GEX_Q30_bases_in_read_2', 'GEX_Q30_bases_in_UMI', 'GEX_Percent_duplicates', 'GEX_Q30_bases_in_sample_index_i1', 'GEX_Q30_bases_in_sample_index_i2', 'GEX_Reads_with_TSO', 'GEX_Sequenced_read_pairs', 'GEX_Valid_UMIs', 'GEX_Valid_barcodes', 'GEX_Reads_mapped_to_genome', 'GEX_Reads_mapped_confidently_to_genome', 'GEX_Reads_mapped_confidently_to_intergenic_regions', 'GEX_Reads_mapped_confidently_to_intronic_regions', 'GEX_Reads_mapped_confidently_to_exonic_regions', 'GEX_Reads_mapped_confidently_to_transcriptome', 'GEX_Reads_mapped_antisense_to_gene', 'GEX_Fraction_of_transcriptomic_reads_in_cells', 'GEX_Total_genes_detected', 'GEX_Median_UMI_counts_per_cell', 'GEX_Median_genes_per_cell', 'Multiome_Feature_linkages_detected', 'Multiome_Linked_genes', 'Multiome_Linked_peaks', 'ATAC_Confidently_mapped_read_pairs', 'ATAC_Fraction_of_genome_in_peaks', 'ATAC_Fraction_of_high_quality_fragments_in_cells', 'ATAC_Fraction_of_high_quality_fragments_overlapping_TSS', 'ATAC_Fraction_of_high_quality_fragments_overlapping_peaks', 'ATAC_Fraction_of_transposition_events_in_peaks_in_cells', 'ATAC_Mean_raw_read_pairs_per_cell', 'ATAC_Median_high_quality_fragments_per_cell', 'ATAC_Non-nuclear_read_pairs', 'ATAC_Number_of_peaks', 'ATAC_Percent_duplicates', 'ATAC_Q30_bases_in_barcode', 'ATAC_Q30_bases_in_read_1', 'ATAC_Q30_bases_in_read_2', 'ATAC_Q30_bases_in_sample_index_i1', 'ATAC_Sequenced_read_pairs', 'ATAC_TSS_enrichment_score', 'ATAC_Unmapped_read_pairs', 'ATAC_Valid_barcodes', 'Number of mapped reads', 'Number of unmapped reads', 'Number of multimapped reads', 'Number of reads', 'Number of UMIs', 'Genes detected', 'Doublet score', 'Fraction mitochondrial UMIs', 'Used in analysis', 'Class confidence', 'Class', 'Subclass confidence', 'Subclass', 'Supertype confidence', 'Supertype (non-expanded)', 'Supertype', 'Continuous Pseudo-progression Score', 'Severely Affected Donor'
    var: 'gene_ids'
    uns: 'APOE4 Status_colors', 'Braak_colors', 'CERAD score_colors', 'Cognitive Status_colors', 'Great Apes Metadata', 'Highest Lewy Body Disease_colors', 'LATE_colors', 'Overall AD neuropathological Change_colors', 'Sex_colors', 'Subclass_colors', 'Supertype_colors', 'Thal_colors', 'UW Clinical Metadata', 'X_normalization', 'batch_condition', 'default_embedding', 'neighbors', 'title', 'umap'
    obsm: 'X_scVI', 'X_umap'
    layers: 'UMIs'
    obsp: 'connectivities', 'distances'

adata.X=adata.layers['UMIs'].copy()

sc.pp.filter_cells(adata, min_genes=200)

adata.var['mt'] = adata.var_names.str.startswith('MT-')  # annotate the group of mitochondrial genes as 'mt'
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)

ribo_cells = adata.var_names.str.startswith('RPS','RPL')
adata.obs['pct_counts_ribo'] = np.ravel(100*np.sum(adata[:, ribo_cells].X, axis = 1) / np.sum(adata.X, axis = 1))

piaso.pl.plot_features_violin(adata,
                              ['n_genes_by_counts', 'total_counts', 'pct_counts_mt','pct_counts_ribo'],
                              groupby='Subclass')

adata.layers['raw']=adata.X.copy()

sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

adata.layers['log1p']=adata.X

normalizing counts per cell
    finished (0:00:00)

sc.experimental.pp.highly_variable_genes(adata,
                                         layer='raw',
                                         n_top_genes=3000)

pr_offset=sc.experimental.pp.normalize_pearson_residuals(adata,
                                                         layer='raw',
                                                         inplace=False)

adata.X=pr_offset['X']

del pr_offset

extracting highly variable genes
--> added
    'highly_variable', boolean vector (adata.var)
    'highly_variable_rank', float vector (adata.var)
    'highly_variable_nbatches', int vector (adata.var)
    'highly_variable_intersection', boolean vector (adata.var)
    'means', float vector (adata.var)
    'variances', float vector (adata.var)
    'residual_variances', float vector (adata.var)
computing analytic Pearson residuals on raw

stderr

.../site-packages/scanpy/experimental/pp/_normalization.py:70: RuntimeWarning: invalid value encountered in divide
  residuals = diff / np.sqrt(mu + mu**2 / theta)

finished (0:00:33)

adata.obsm['X_umap_backup']=adata.obsm['X_umap'].copy()

piaso.tl.runSVD(adata,
                   use_highly_variable=True,
                   n_components=50,
                   random_state=10,
                   key_added='X_svd')

%%time
sc.pp.neighbors(adata,
                use_rep='X_svd',
                n_neighbors=15,
                random_state=10,
                knn=True,
                method="umap")
sc.tl.umap(adata)

computing neighbors
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:38)
computing UMAP
    finished: added
    'X_umap', UMAP coordinates (adata.obsm)
    'umap', UMAP parameters (adata.uns) (0:00:36)
CPU times: user 5min 26s, sys: 843 ms, total: 5min 27s
Wall time: 1min 14s

sc.pl.umap(adata,
           color=['Subclass'],
           palette=piaso.pl.color.d_color4,
           cmap=piaso.pl.color.c_color4,
           size=10,
           frameon=True)

adata.obsm['X_umap_raw']=adata.obsm['X_umap'].copy()

sc.pl.umap(adata,
           color=['n_genes_by_counts', 'total_counts','pct_counts_mt','pct_counts_ribo'],
           cmap='Spectral_r',
           palette=piaso.pl.color.d_color4,
           ncols=4,
           size=10,
           frameon=False,)

%%time
sc.tl.leiden(adata,resolution=0.5,key_added='Leiden',flavor="leidenalg",n_iterations=-1)

running Leiden clustering

stderr

<timed eval>:1: FutureWarning: In the future, the default backend for leiden will be igraph instead of leidenalg.

 To achieve the future defaults please pass: flavor="igraph" and n_iterations=2.  directed must also be False to work with igraph's implementation.

finished: found 24 clusters and added
    'Leiden', the cluster labels (adata.obs, categorical) (0:00:01)
CPU times: user 1.57 s, sys: 40.1 ms, total: 1.61 s
Wall time: 1.61 s

logging.getLogger('matplotlib.font_manager').disabled = True

sc.pl.umap(adata,
           color=['Leiden'],
           palette=piaso.pl.color.d_color4,
           cmap=piaso.pl.color.c_color4,
           legend_fontsize=12,
           legend_fontoutline=2,
           legend_loc='on data',
           size=10,
           frameon=False)

pip install cosg

n_gene=30
cosg.cosg(adata,
          key_added='cosg',
          use_raw=False,
          layer='log1p',
          mu=100,
          expressed_pct=0.1,
          remove_lowly_expressed=True,
          n_genes_user=100,
          groupby='Leiden')

sc.tl.dendrogram(adata,groupby='Leiden',use_rep='X_svd')
df_tmp=pd.DataFrame(adata.uns['cosg']['names'][:3,]).T
df_tmp=df_tmp.reindex(adata.uns['dendrogram_'+'Leiden']['categories_ordered'])
marker_genes_list={idx: list(row.values) for idx, row in df_tmp.iterrows()}
marker_genes_list = {k: v for k, v in marker_genes_list.items() if not any(isinstance(x, float) for x in v)}

sc.pl.dotplot(adata,
              marker_genes_list,
              groupby='Leiden',
              layer='log1p',
              dendrogram=True,
              swap_axes=True,
              standard_scale='var',
              cmap='Spectral_r',
              figsize=[10,20])

Storing dendrogram info using `.uns['dendrogram_Leiden']`

marker_gene=pd.DataFrame(adata.uns['cosg']['names'])

sc.pl.umap(adata,
           color=['Leiden'],
           palette=piaso.pl.color.d_color4,
           cmap=piaso.pl.color.c_color4,
           legend_fontsize=12,
           legend_fontoutline=2,
           legend_loc='on data',
           size=10,
           frameon=False)

adata.obsm['X_umap_svd']=adata.obsm['X_umap'].copy()

piaso.tl.runGDR(adata,
                batch_key=None,
                groupby='Leiden',
                n_gene=30,
                mu=1.0,
                use_highly_variable=True,
                n_highly_variable_genes=5000,
                layer='log1p',
                score_layer='log1p',
                n_svd_dims=50,
                resolution=1.0,
                scoring_method=None,
                key_added='X_gdr',
                verbosity=0)

GDR embeddings saved to adata.obsm['X_gdr']

%%time
sc.pp.neighbors(adata,
                use_rep='X_gdr',
                n_neighbors=15,
                random_state=10,
                knn=True,
                method="umap")
sc.tl.umap(adata)

computing neighbors
    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:04)
computing UMAP
    finished: added
    'X_umap', UMAP coordinates (adata.obsm)
    'umap', UMAP parameters (adata.uns) (0:00:30)
CPU times: user 2min 46s, sys: 565 ms, total: 2min 47s
Wall time: 34.7 s

sc.pl.umap(adata,
           color=['Subclass'],
           palette=piaso.pl.color.d_color4,
           cmap=piaso.pl.color.c_color4,
           ncols=1,
           size=10,
           frameon=True)

stratified_split = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=10)
X = adata.X.copy()
y = adata.obs["Subclass"].values

stratified_split.get_n_splits(X, y)
for i, (train_index, test_index) in enumerate(stratified_split.split(X, y)):
    adata_ref = adata[train_index].copy()
    adata_test = adata[test_index].copy()

piaso.tl.predictCellTypeByGDR(
    adata_test,
    adata_ref,
    layer = 'log1p',
    layer_reference = 'log1p',
    reference_groupby  = 'Subclass',
    query_groupby = 'Leiden',
    mu = 10.0,
    n_genes= 15,
    return_integration = False,
    use_highly_variable = True,
    n_highly_variable_genes = 5000,
    n_svd_dims = 50,
    resolution= 1.0,
    scoring_method= None,
    key_added= None,
    verbosity= 0,
)

stderr

.../site-packages/piaso/tools/_predictCellType.py:137: FutureWarning: Use anndata.concat instead of AnnData.concatenate, AnnData.concatenate is deprecated and will be removed in the future. See the tutorial for concat at: https://anndata.readthedocs.io/en/latest/concatenation.html
  adata_combine=sc.AnnData.concatenate(adata_ref, adata[:,adata_ref.var_names])

Running GDR for the query dataset and the reference dataset:
GDR embeddings saved to adata.obsm['X_gdr']

stderr

2025-03-14 10:25:49,364 - harmonypy - INFO - Computing initial centroids with sklearn.KMeans...
2025-03-14 10:25:55,199 - harmonypy - INFO - sklearn.KMeans initialization complete.
2025-03-14 10:25:55,505 - harmonypy - INFO - Iteration 1 of 10
2025-03-14 10:26:02,568 - harmonypy - INFO - Iteration 2 of 10
2025-03-14 10:26:09,421 - harmonypy - INFO - Iteration 3 of 10
2025-03-14 10:26:14,227 - harmonypy - INFO - Iteration 4 of 10
2025-03-14 10:26:16,166 - harmonypy - INFO - Iteration 5 of 10
2025-03-14 10:26:18,110 - harmonypy - INFO - Iteration 6 of 10
2025-03-14 10:26:20,048 - harmonypy - INFO - Converged after 6 iterations

Predicting cell types:
All finished. The predicted cell types are saved as `CellTypes_gdr` in adata.obs.

adata_test.obs['CellTypes_gdr']=adata_test.obs['CellTypes_gdr'].astype('category')
adata_test.obs['CellTypes_gdr']=adata_test.obs['CellTypes_gdr'].cat.reorder_categories(adata_ref.obs['Subclass'].cat.categories)

sc.pl.embedding(adata_test,
                basis='X_umap',
                color=['CellTypes_gdr'],
                palette=piaso.pl.color.d_color4,
                cmap=piaso.pl.color.c_color3,
                ncols=1,
                size=10,
                frameon=False)

sc.pl.embedding(adata_test,
                basis='X_umap',
                color=['Subclass'],
                palette=piaso.pl.color.d_color4,
                cmap=piaso.pl.color.c_color3,
                ncols=1,
                size=10,
                frameon=False)

confusion_matrix = metrics.confusion_matrix(y[test_index], adata_test.obs['CellTypes_gdr'].values)

confusion_matrix_df = pd.DataFrame(confusion_matrix, columns=adata_test.obs['Subclass'].cat.categories, index=adata_test.obs['Subclass'].cat.categories)
normalized_cf_matrix_df = (confusion_matrix_df - confusion_matrix_df.mean(axis=0))/confusion_matrix_df.std(axis=0)

sns.set_style("whitegrid", {'axes.grid' : False})
sns.set(rc={'figure.figsize':(10, 6)})
sns.heatmap(normalized_cf_matrix_df,
            cmap="Purples",
            xticklabels=True,
            yticklabels=True)
plt.show()

piaso_f1_score=np.round(f1_score(adata_test.obs['Subclass'], adata_test.obs['CellTypes_gdr'], average='micro'), decimals=3)
print(f"The Micro F1 score for PIASO prediction: {piaso_f1_score}")
piaso_f1_score=np.round(f1_score(adata_test.obs['Subclass'], adata_test.obs['CellTypes_gdr'], average='macro'), decimals=3)
print(f"The Macro F1 score for PIASO prediction: {piaso_f1_score}")

The Micro F1 score for PIASO prediction: 0.963
The Macro F1 score for PIASO prediction: 0.935