Pre-processing CellRanger outputs (multiple samples, Human PBMCs)

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

Pre-processing CellRanger outputs (multiple samples, Human PBMCs)

import piaso
import cosg
import numpy as np
import pandas as pd
import scanpy as sc
import logging
from sklearn.preprocessing import StandardScaler
import warnings
import os
import anndata as ad

from matplotlib import rcParams
sc.settings.set_figure_params(scanpy=True, dpi=80, dpi_save=300, frameon=True, vector_friendly=False, fontsize=14)
import matplotlib as mpl
mpl.rcParams['pdf.fonttype'] = 42
mpl.rcParams['ps.fonttype'] = 42

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

data_dir = ".../Data/Public/PBMCMultiomeRop2023"
save_dir = ".../Results/single-cell/Methods/DataProcessing/PBMCMultiomeRop2023_SAN1_SAN2"
prefix = "HumanPBMCs_Multiome_RNA"

# !.../gdrive files download --destination .../Data/Public/PBMCMultiomeRop2023 1g1V-VTy2qHK7PJIbW-t9hRFzHCqFCe4T

!.../gdrive files download --destination .../Data/Public/PBMCMultiomeRop2023/SAN2 1q5mXAE6YWVlzTJLSG0qUgnATitml0gzX

Downloading filtered_feature_bc_matrix.h5
Successfully downloaded filtered_feature_bc_matrix.h5

!.../gdrive files download --destination .../Data/Public/PBMCMultiomeRop2023/SAN1 1Bjr81oJWto-06ak14IXB804xcLypPbRj

Downloading filtered_feature_bc_matrix.h5
Successfully downloaded filtered_feature_bc_matrix.h5

data_path='.../Data/Public/PBMCMultiomeRop2023'
samples = [entry.name for entry in os.scandir(data_path) if entry.is_dir()]
print(samples)

['SAN1', 'SAN2']

def load_multiple_anndata(
    data_path,
    samples,
    key_added):
    anndata_list = []
    for sample in samples:
        try:
            # Load the data
            file_name = os.listdir(data_path+'/'+sample+'/RNA')[0]
            print(file_name)
            h5_path=data_path+'/'+sample+'/RNA/'+file_name
            if os.path.exists(h5_path):
                adata = sc.read_10x_h5(h5_path)
                print('Loading the h5 file')

            else:
                adata = sc.read_10x_mtx(data_path+'/'+sample, cache=True)
            adata.var_names_make_unique()
            adata.obs_names_make_unique()
            adata.obs[key_added]=sample
            anndata_list.append(adata)
            print('Sample loaded: ', sample)
        except Exception as e:
            print(f"Error loading {data_path}: {e}")

    return anndata_list

adata_list=load_multiple_anndata(data_path=data_path, samples=samples, key_added='Sample')

filtered_feature_bc_matrix.h5

stderr

.../site-packages/anndata/_core/anndata.py:1758: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
  utils.warn_names_duplicates("var")
.../site-packages/anndata/_core/anndata.py:1758: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
  utils.warn_names_duplicates("var")

Loading the h5 file
Sample loaded:  SAN1
filtered_feature_bc_matrix.h5
Loading the h5 file
Sample loaded:  SAN2

stderr

.../site-packages/anndata/_core/anndata.py:1758: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
  utils.warn_names_duplicates("var")
.../site-packages/anndata/_core/anndata.py:1758: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
  utils.warn_names_duplicates("var")

adata_list

[AnnData object with n_obs × n_vars = 3545 × 36601
     obs: 'Sample'
     var: 'gene_ids', 'feature_types', 'genome', 'interval',
 AnnData object with n_obs × n_vars = 4360 × 36601
     obs: 'Sample'
     var: 'gene_ids', 'feature_types', 'genome', 'interval']

adata=ad.concat(adata_list, join='outer',index_unique="-")
adata

AnnData object with n_obs × n_vars = 7905 × 36601
    obs: 'Sample'

adata.X.data

array([ 1.,  1.,  1., ..., 43., 18., 56.],
      shape=(12509350,), dtype=float32)

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

AnnData object with n_obs × n_vars = 7905 × 36601
    obs: 'Sample'
    layers: 'raw'

adata.var

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

adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)

adata.var['ribo'] = adata.var_names.str.startswith('RPS','RPL')
sc.pp.calculate_qc_metrics(adata, qc_vars=['ribo'], percent_top=None, log1p=False, inplace=True)

stderr

/tmp/ipykernel_2699886/3834852883.py:1: FutureWarning: Allowing a non-bool 'na' in obj.str.startswith is deprecated and will raise in a future version.
  adata.var['ribo'] = adata.var_names.str.startswith('RPS','RPL')

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

experiments=np.unique(adata.obs['Sample'])
adata.obs['scrublet_score']=np.repeat(0,adata.n_obs)
adata.obs['predicted_doublets']=np.repeat(False,adata.n_obs)

import scrublet as scr

for experiment in experiments:
    print(experiment)
    adatai=adata[adata.obs['Sample']==experiment]

    scrub = scr.Scrublet(adatai.X.todense(),random_state=10)
    doublet_scores, predicted_doublets = scrub.scrub_doublets()

    adata.obs['predicted_doublets'][adatai.obs_names]=predicted_doublets

    adata.obs['scrublet_score'][adatai.obs_names]=doublet_scores

SAN1
Preprocessing...
Simulating doublets...
Embedding transcriptomes using PCA...
Calculating doublet scores...
Automatically set threshold at doublet score = 0.15
Detected doublet rate = 15.9%
Estimated detectable doublet fraction = 77.1%
Overall doublet rate:
        Expected   = 10.0%
        Estimated  = 20.6%
Elapsed time: 3.1 seconds
SAN2

stderr

/tmp/ipykernel_2699886/3464195482.py:10: FutureWarning: ChainedAssignmentError: behaviour will change in pandas 3.0!
You are setting values through chained assignment. Currently this works in certain cases, but when using Copy-on-Write (which will become the default behaviour in pandas 3.0) this will never work to update the original DataFrame or Series, because the intermediate object on which we are setting values will behave as a copy.
A typical example is when you are setting values in a column of a DataFrame, like:

df["col"][row_indexer] = value

Use `df.loc[row_indexer, "col"] = values` instead, to perform the assignment in a single step and ensure this keeps updating the original `df`.

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy

  adata.obs['predicted_doublets'][adatai.obs_names]=predicted_doublets
/tmp/ipykernel_2699886/3464195482.py:10: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  adata.obs['predicted_doublets'][adatai.obs_names]=predicted_doublets
/tmp/ipykernel_2699886/3464195482.py:12: FutureWarning: ChainedAssignmentError: behaviour will change in pandas 3.0!
You are setting values through chained assignment. Currently this works in certain cases, but when using Copy-on-Write (which will become the default behaviour in pandas 3.0) this will never work to update the original DataFrame or Series, because the intermediate object on which we are setting values will behave as a copy.
A typical example is when you are setting values in a column of a DataFrame, like:

df["col"][row_indexer] = value

Use `df.loc[row_indexer, "col"] = values` instead, to perform the assignment in a single step and ensure this keeps updating the original `df`.

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
... [9 more lines]

Preprocessing...
Simulating doublets...
Embedding transcriptomes using PCA...
Calculating doublet scores...
Automatically set threshold at doublet score = 0.14
Detected doublet rate = 17.1%
Estimated detectable doublet fraction = 79.9%
Overall doublet rate:
        Expected   = 10.0%
        Estimated  = 21.5%
Elapsed time: 5.5 seconds

stderr

/tmp/ipykernel_2699886/3464195482.py:10: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  adata.obs['predicted_doublets'][adatai.obs_names]=predicted_doublets

piaso.pl.plot_features_violin(adata,
                              ['scrublet_score'],
                              groupby='Sample',
                              width_single=2)

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

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

CPU times: user 198 ms, sys: 26.5 ms, total: 225 ms
Wall time: 124 ms

%%time
piaso.tl.infog(adata,
               copy=False,
               inplace=False,
               n_top_genes=3000,
               key_added='infog',
               key_added_highly_variable_gene='highly_variable',
               layer='raw')

The normalized data is saved as `infog` in `adata.layers`.
The highly variable genes are saved as `highly_variable` in `adata.var`.
Finished INFOG normalization.
CPU times: user 590 ms, sys: 199 ms, total: 789 ms
Wall time: 792 ms

%%time
piaso.tl.runSVD(adata,
                use_highly_variable=True,
                n_components=30,
                random_state=10,
                key_added='X_svd',
                layer='infog')

CPU times: user 4.07 s, sys: 6.2 ms, total: 4.07 s
Wall time: 493 ms

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

sc.tl.umap(adata)

CPU times: user 27.5 s, sys: 92.7 ms, total: 27.6 s
Wall time: 24.2 s

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

save_dir+'/'+prefix+'_raw_QC.h5ad'

'.../Results/single-cell/Methods/DataProcessing/PBMCMultiomeRop2023_SAN1_SAN2/HumanPBMCs_Multiome_RNA_raw_QC.h5ad'

adata.write(save_dir+'/'+prefix+'_raw_QC.h5ad')

# adata = sc.read(save_dir+'/'+prefix+'_raw_QC.h5ad')

adata=adata[adata.obs['n_genes_by_counts']>500].copy()

adata=adata[adata.obs['n_genes_by_counts']<10000].copy()

adata = adata[adata.obs['scrublet_score']<0.3].copy()
adata

AnnData object with n_obs × n_vars = 7218 × 36601
    obs: 'Sample', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'pct_counts_ribo', 'scrublet_score', 'predicted_doublets'
    var: 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'ribo', 'infog_var', 'highly_variable'
    uns: 'Sample_colors', 'log1p', 'neighbors', 'umap'
    obsm: 'X_svd', 'X_umap'
    layers: 'infog', 'log1p', 'raw'
    obsp: 'connectivities', 'distances'

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

piaso.pl.plot_embeddings_split(adata,
                               color='n_genes_by_counts',
                               cmap=piaso.pl.color.c_color1,
                               palette=piaso.pl.color.d_color1,
                               layer=None,
                               splitby='Sample',
                               size=25,
                               frameon=False,)

marker_gene_df = pd.read_csv(".../Data/Public/10x_Human_PBMC/PIASOmarkerDB_AllenHumanImmuneHealthAtlas_L2_251219.csv")
marker_gene_df

marker_gene_dict = {}
cell_types = np.unique(marker_gene_df['Cell_Type'])
for cell_type in cell_types:
    marker_gene_dict[cell_type] = list(marker_gene_df[marker_gene_df['Cell_Type']==cell_type]['Gene'])

%%time
piaso.tl.infog(adata,
               copy=False,
               inplace=False,
               n_top_genes=3000,
               key_added='infog',
               key_added_highly_variable_gene='highly_variable',
               layer='raw')

The normalized data is saved as `infog` in `adata.layers`.
The highly variable genes are saved as `highly_variable` in `adata.var`.
Finished INFOG normalization.
CPU times: user 747 ms, sys: 28.8 s, total: 29.5 s
Wall time: 29.6 s

%%time
piaso.tl.runSVD(adata,
                use_highly_variable=True,
                n_components=30,
                random_state=10,
                key_added='X_svd',
                layer='infog')

CPU times: user 5.36 s, sys: 134 ms, total: 5.5 s
Wall time: 894 ms

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

sc.tl.umap(adata)

CPU times: user 39.3 s, sys: 308 ms, total: 39.6 s
Wall time: 36.3 s

%%time
piaso.tl.predictCellTypeByMarker(adata,
                                 marker_gene_set=marker_gene_dict,
                                 score_method='piaso',
                                 use_rep='X_svd',
                                 use_score=False,
                                 smooth_prediction=True,
                                 inplace=True)

Calculating gene set scores using piaso method...

stderr

Scoring gene sets: 100%|██████████| 29/29 [00:16<00:00,  1.75set/s]

Predicting cell types based on marker gene p-values...
Smoothing cell type predictions...
Smoothing cell type predictions from 'CellTypes_predicted_raw' using 7-nearest neighbors
Smoothed predictions stored in adata.obs['CellTypes_predicted_smoothed']
Confidence scores stored in adata.obs['CellTypes_predicted_smoothed_confidence']
Modified 1632 cell labels (22.61% of total)
Cell type prediction completed. Results saved to:
  - adata.obs['CellTypes_predicted']: predicted cell types
  - adata.obsm['CellTypes_predicted_score']: full score matrix
  - adata.obsm['CellTypes_predicted_nlog10pvals']: full -log10(p-value) matrix
  - adata.obs['CellTypes_predicted_nlog10pvals']: maximum -log10(p-values)
  - adata.obs['CellTypes_predicted_raw']: original unsmoothed predictions
  - adata.obs['CellTypes_predicted_confidence_smoothed']: smoothing confidence scores
CPU times: user 6.85 s, sys: 146 ms, total: 6.99 s
Wall time: 24.6 s

category_dict = {"Human_PBMC": ['ASDC', 'CD14 monocyte', 'CD16 monocyte', 'CD56bright NK cell',
                                'CD56dim NK cell', 'CD8aa', 'DN T cell', 'Effector B cell',
                                'Erythrocyte', 'ILC', 'Intermediate monocyte', 'MAIT', 'Memory B cell',
                                'Memory CD4 T cell', 'Memory CD8 T cell', 'Naive B cell',
                                'Naive CD4 T cell', 'Naive CD8 T cell', 'Plasma cell', 'Platelet',
                                'Progenitor cell','Proliferating NK cell', 'Proliferating T cell',
                                'Transitional B cell', 'Treg', 'cDC1','cDC2', 'gdT','pDC']}

cell_types

array(['ASDC', 'CD14 monocyte', 'CD16 monocyte', 'CD56bright NK cell',
       'CD56dim NK cell', 'CD8aa', 'DN T cell', 'Effector B cell',
       'Erythrocyte', 'ILC', 'Intermediate monocyte', 'MAIT',
       'Memory B cell', 'Memory CD4 T cell', 'Memory CD8 T cell',
       'Naive B cell', 'Naive CD4 T cell', 'Naive CD8 T cell',
       'Plasma cell', 'Platelet', 'Progenitor cell',
       'Proliferating NK cell', 'Proliferating T cell',
       'Transitional B cell', 'Treg', 'cDC1', 'cDC2', 'gdT', 'pDC'],
      dtype=object)

np.unique(adata.obs['CellTypes_predicted'])

array(['CD14 monocyte', 'CD16 monocyte', 'CD56bright NK cell',
       'CD56dim NK cell', 'CD8aa', 'DN T cell', 'Effector B cell', 'ILC',
       'Intermediate monocyte', 'MAIT', 'Memory B cell',
       'Memory CD4 T cell', 'Memory CD8 T cell', 'Naive B cell',
       'Naive CD4 T cell', 'Naive CD8 T cell', 'Plasma cell', 'Platelet',
       'Progenitor cell', 'Proliferating NK cell', 'Transitional B cell',
       'Treg', 'cDC1', 'cDC2', 'gdT', 'pDC'], dtype=object)

adata.obs['CellTypes_predicted'] = adata.obs['CellTypes_predicted'].astype('category')

adata.obs['CellTypes_predicted'] = adata.obs['CellTypes_predicted'].cat.add_categories(
    [cat for cat in category_dict["Human_PBMC"]
     if cat not in adata.obs['CellTypes_predicted'].cat.categories]
)

adata.obs['CellTypes_predicted'] = adata.obs['CellTypes_predicted'].cat.reorder_categories(
    category_dict["Human_PBMC"]
)
adata

AnnData object with n_obs × n_vars = 7218 × 36601
    obs: 'Sample', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'pct_counts_ribo', 'scrublet_score', 'predicted_doublets', 'CellTypes_predicted', 'CellTypes_predicted_nlog10pvals', 'CellTypes_predicted_raw', 'CellTypes_predicted_smoothed', 'CellTypes_predicted_smoothed_confidence', 'CellTypes_predicted_confidence_smoothed'
    var: 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'ribo', 'infog_var', 'highly_variable'
    uns: 'Sample_colors', 'log1p', 'neighbors', 'umap'
    obsm: 'X_svd', 'X_umap', 'CellTypes_predicted_score', 'CellTypes_predicted_nlog10pvals'
    layers: 'infog', 'log1p', 'raw'
    obsp: 'connectivities', 'distances'

sc.pl.umap(adata,
           color=['CellTypes_predicted'],
           palette=piaso.pl.color.d_color20,
           legend_fontsize=8,
           legend_fontoutline=1,
           # legend_loc='on data',
           ncols=1,
           size=10,
           frameon=False)

%%time
sc.tl.leiden(adata,resolution=2.5,key_added='Leiden')

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.

CPU times: user 4.01 s, sys: 16.3 ms, total: 4.02 s
Wall time: 4.05 s

sc.pl.umap(adata,
           color=['CellTypes_predicted','Leiden'],
           palette=piaso.pl.color.d_color20,
           legend_fontsize=12,
           legend_fontoutline=2,
           # legend_loc='on data',
           ncols=1,
           size=10,
           frameon=False)

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

CPU times: user 915 ms, sys: 540 ms, total: 1.46 s
Wall time: 1.47 s

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='infog',
              dendrogram=True,
              swap_axes=False,
              standard_scale='var',
              cmap='Spectral_r')

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

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

adata = adata[adata.obs['pct_counts_mt']<20].copy()

adata = adata[adata.obs['pct_counts_ribo']<5].copy()

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

sc.pl.umap(adata,
           color=['Leiden'],
           groups=['29'],
           palette=piaso.pl.color.d_color20,
           legend_fontsize=12,
           legend_fontoutline=2,
           # legend_loc='on data',
           ncols=1,
           size=10,
           frameon=False)

cluster_check = '20'
marker_gene[cluster_check].values

array(['RPL19', 'RPL32', 'RPL13A', 'RPL28', 'RPS25', 'RPS8', 'RPLP2',
       'RPL14', 'RPS12', 'RPS27', 'MKI67', 'RPL22', 'RPL23A', 'RPS15A',
       'RPL35', 'RPL13', 'RPS18', 'RPL3', 'RPL17', 'RPS3', 'EEF1A1',
       'RPS20', 'RPL12', 'RPS3A', 'RPL18A', 'RPLP0', 'RPS23', 'RPL39',
       'RPS6', 'RPL31'], dtype=object)

sc.pl.umap(adata,
           color=['IGHA1', 'SDC1', 'IGLV2-11', 'JCHAIN', 'IGHG4', 'TNFRSF17',
       'IGLV2-14', 'GLDC', 'IGHV3-15', 'TXNDC5', 'IGHG2', 'MIXL1', 'MZB1'],
           cmap=piaso.pl.color.c_color1,
           palette=piaso.pl.color.d_color1,
           ncols=3,
           size=10,
           frameon=False)

sc.pl.dotplot(adata,
              marker_gene[cluster_check].values[:30],
              groupby='Leiden',
              dendrogram=False,
              swap_axes=True,
              standard_scale='var',
              cmap='Spectral_r')

sc.pl.umap(adata,
           color=['Leiden'],
           groups=['17','20','27','28'],
           palette=piaso.pl.color.d_color20,
           legend_fontsize=12,
           legend_fontoutline=2,
           # legend_loc='on data',
           ncols=1,
           size=10,
           frameon=False)

piaso.pl.plot_embeddings_split(adata,
                               color='CellTypes_predicted',
                               layer=None,
                               splitby='Sample',
                               size=25,
                               frameon=False,)

adata=adata[~adata.obs['Leiden'].isin(['17','20','27','28'])].copy()
adata

AnnData object with n_obs × n_vars = 6591 × 36601
    obs: 'Sample', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'pct_counts_ribo', 'scrublet_score', 'predicted_doublets', 'CellTypes_predicted', 'CellTypes_predicted_nlog10pvals', 'CellTypes_predicted_raw', 'CellTypes_predicted_smoothed', 'CellTypes_predicted_smoothed_confidence', 'CellTypes_predicted_confidence_smoothed', 'Leiden'
    var: 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'ribo', 'infog_var', 'highly_variable'
    uns: 'Sample_colors', 'log1p', 'neighbors', 'umap', 'CellTypes_predicted_colors', 'Leiden', 'Leiden_colors', 'cosg', 'dendrogram_Leiden'
    obsm: 'X_svd', 'X_umap', 'CellTypes_predicted_score', 'CellTypes_predicted_nlog10pvals'
    layers: 'infog', 'log1p', 'raw'
    obsp: 'connectivities', 'distances'

%%time
piaso.tl.infog(adata,
               copy=False,
               inplace=False,
               n_top_genes=3000,
               key_added='infog',
               key_added_highly_variable_gene='highly_variable',
               layer='raw')

The normalized data is saved as `infog` in `adata.layers`.
The highly variable genes are saved as `highly_variable` in `adata.var`.
Finished INFOG normalization.
CPU times: user 685 ms, sys: 4.9 s, total: 5.59 s
Wall time: 5.61 s

%%time
piaso.tl.runSVD(adata,
                use_highly_variable=True,
                n_components=30,
                random_state=10,
                key_added='X_svd',
                layer='infog')

CPU times: user 4.72 s, sys: 0 ns, total: 4.72 s
Wall time: 656 ms

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

sc.tl.umap(adata)

CPU times: user 25.4 s, sys: 11.8 ms, total: 25.4 s
Wall time: 21 s

%%time
piaso.tl.predictCellTypeByMarker(adata,
                                 marker_gene_set=marker_gene_dict,
                                 score_method='piaso',
                                 use_rep='X_svd',
                                 use_score=False,
                                 key_added = 'CellTypes_predicted',
                                 smooth_prediction=True,
                                 inplace=True)

Calculating gene set scores using piaso method...

stderr

Scoring gene sets: 100%|██████████| 29/29 [00:17<00:00,  1.69set/s]

Predicting cell types based on marker gene p-values...
Smoothing cell type predictions...
Smoothing cell type predictions from 'CellTypes_predicted_raw' using 7-nearest neighbors
Smoothed predictions stored in adata.obs['CellTypes_predicted_smoothed']
Confidence scores stored in adata.obs['CellTypes_predicted_smoothed_confidence']
Modified 1402 cell labels (21.27% of total)
Cell type prediction completed. Results saved to:
  - adata.obs['CellTypes_predicted']: predicted cell types
  - adata.obsm['CellTypes_predicted_score']: full score matrix
  - adata.obsm['CellTypes_predicted_nlog10pvals']: full -log10(p-value) matrix
  - adata.obs['CellTypes_predicted_nlog10pvals']: maximum -log10(p-values)
  - adata.obs['CellTypes_predicted_raw']: original unsmoothed predictions
  - adata.obs['CellTypes_predicted_confidence_smoothed']: smoothing confidence scores
CPU times: user 10.3 s, sys: 103 ms, total: 10.4 s
Wall time: 24.6 s

sc.pl.umap(adata,
           color=['CellTypes_predicted'],
           groups=['CD8aa'],
           palette=piaso.pl.color.d_color20,
           legend_fontsize=8,
           legend_fontoutline=1,
           # legend_loc='on data',
           ncols=1,
           size=10,
           frameon=False)

sc.pl.umap(adata,
           color=['CellTypes_predicted'],
           # groups=['CD8aa'],
           palette=piaso.pl.color.d_color20,
           legend_fontsize=8,
           legend_fontoutline=1,
           # legend_loc='on data',
           ncols=1,
           size=10,
           frameon=False)

%%time
n_gene=30
cosg.cosg(adata,
          key_added='cosg',
          use_raw=False,
          layer='infog',
          mu=100,
          expressed_pct=0.1,
          remove_lowly_expressed=True,
          n_genes_user=n_gene,
          groupby='CellTypes_predicted')

CPU times: user 723 ms, sys: 98.9 ms, total: 822 ms
Wall time: 825 ms

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

sc.tl.dendrogram(adata,groupby='CellTypes_predicted',use_rep='X_svd')
df_tmp=pd.DataFrame(adata.uns['cosg']['names'][:3,]).T
df_tmp=df_tmp.reindex(adata.uns['dendrogram_'+'CellTypes_predicted']['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='CellTypes_predicted',
              layer='infog',
              dendrogram=True,
              swap_axes=False,
              standard_scale='var',
              cmap='Spectral_r')

cluster_check = 'CD8aa'
marker_gene[cluster_check].values

array(['AC005944.1', 'AL031590.1', 'TNIP3', 'ZSCAN26', 'PDGFB',
       'AL121839.2', 'AMIGO1', 'ZSCAN2', 'AL359706.1', 'STYK1', 'ZNF681',
       'RNFT2', 'METTL18', 'TMEM116', 'BCAT2', 'VIPR1-AS1', 'ZNF268',
       'LYSMD4', 'KLRC4', 'TNFRSF9', 'STK36', 'LINC01118', 'C21orf62-AS1',
       'KLHL25', 'SLC35C1', 'DHRS13', 'MED12', 'TRPC1', 'TULP2',
       'AC137810.1'], dtype=object)

sc.pl.umap(adata,
           color=['RNF26', 'BCKDK', 'STYK1', 'ZNF268', 'RNFT2', 'BICRA-AS1',
       'ZNF510', 'AMIGO1', 'MIER3', 'AC013437.1', 'SLC24A1', 'LAGE3'],
           cmap=piaso.pl.color.c_color1,
           palette=piaso.pl.color.d_color1,
           ncols=3,
           size=10,
           frameon=False)

sc.pl.dotplot(adata,
              marker_gene[cluster_check].values[:30],
              groupby='CellTypes_predicted',
              dendrogram=False,
              swap_axes=True,
              standard_scale='var',
              cmap='Spectral_r')

np.unique(adata.obs['CellTypes_predicted'])

array(['ASDC', 'CD14 monocyte', 'CD16 monocyte', 'CD56bright NK cell',
       'CD56dim NK cell', 'CD8aa', 'DN T cell', 'Effector B cell',
       'Erythrocyte', 'ILC', 'Intermediate monocyte', 'MAIT',
       'Memory B cell', 'Memory CD4 T cell', 'Memory CD8 T cell',
       'Naive B cell', 'Naive CD4 T cell', 'Naive CD8 T cell',
       'Plasma cell', 'Proliferating NK cell', 'Transitional B cell',
       'Treg', 'cDC1', 'cDC2', 'gdT', 'pDC'], dtype=object)

adata.write(save_dir + '/' + prefix + '_QC.h5ad')

adata = sc.read(save_dir + '/' + prefix + '_QC.h5ad')

marker_gene.to_csv(save_dir + '/' + prefix + '_CellType_markerGenes.csv')

!.../gdrive files upload --parent  .../Results/single-cell/Methods/DataProcessing/PBMCMultiomeRop2023_SAN1_SAN2/HumanPBMCs_Multiome_RNA_QC.h5ad

Uploading .../Results/single-cell/Methods/DataProcessing/PBMCMultiomeRop2023_SAN1_SAN2/HumanPBMCs_Multiome_RNA_QC.h5ad
File successfully uploaded
Id:
Name: HumanPBMCs_Multiome_RNA_QC.h5ad
Mime: application/octet-stream
Size: 318.5 MB
Created: 2026-01-20 14:02:15
Modified: 2026-01-20 14:02:15
MD5: bdd76ad91aa88845d0791327f8675972
Shared: True
Parents:
ViewUrl: https://drive.google.com/file/d/1bpWCqhadqw29Av9RsqAJqWW7kmEpQOHQ/view?usp=drivesdk

!.../gdrive files upload --parent  .../Results/single-cell/Methods/DataProcessing/PBMCMultiomeRop2023_SAN1_SAN2/HumanPBMCs_Multiome_RNA_raw_QC.h5ad

Uploading .../Results/single-cell/Methods/DataProcessing/PBMCMultiomeRop2023_SAN1_SAN2/HumanPBMCs_Multiome_RNA_raw_QC.h5ad
File successfully uploaded
Id:
Name: HumanPBMCs_Multiome_RNA_raw_QC.h5ad
Mime: application/octet-stream
Size: 407.7 MB
Created: 2026-01-20 14:02:26
Modified: 2026-01-20 14:02:26
MD5: 97852a95d1b331836930c376a5aed3e5
Shared: True
Parents:
ViewUrl: https://drive.google.com/file/d/1DaeasGmgNIBuppSkqT_jiRMDCouxEuiN/view?usp=drivesdk

!.../gdrive files upload --parent  .../Results/single-cell/Methods/DataProcessing/PBMCMultiomeRop2023_SAN1_SAN2/HumanPBMCs_Multiome_RNA_CellType_markerGenes.csv

Uploading .../Results/single-cell/Methods/DataProcessing/PBMCMultiomeRop2023_SAN1_SAN2/HumanPBMCs_Multiome_RNA_CellType_markerGenes.csv
File successfully uploaded
Id:
Name: HumanPBMCs_Multiome_RNA_CellType_markerGenes.csv
Mime: text/csv
Size: 6.1 KB
Created: 2026-01-20 14:02:28
Modified: 2026-01-20 14:02:28
MD5: 1ac8ebc9c6b336588015cb2be1ed9432
Shared: True
Parents:
ViewUrl: https://drive.google.com/file/d/1BLabNMEkVlHfbooFe55jD9d1yrcl_66K/view?usp=drivesdk