Normalize

import scanpy as sc
import cupy as cp # enable GPU-accelerated computing

import time
import rapids_singlecell as rsc

import warnings

warnings.filterwarnings("ignore")

import rmm # The Rapids Memory Manager (RMM) Python package that allows users to customize the allocation of memory for devices and hosts.
from rmm.allocators.cupy import rmm_cupy_allocator

rmm.reinitialize(
    managed_memory=False,
    pool_allocator=False,
    devices=0, # GPU device IDs to register, By default registers only GPU 0.
)

cp.cuda.set_allocator(rmm_cupy_allocator)

data_load_start=time.time()

%%time
#adata = sc.read_h5ad("/stomics_data/liminData/LiLi_projects/D02672A2/seurat/DAPI/addimage/tissue_sc.h5ad")
adata = sc.read_h5ad("/stomics_data/liminData/stereopy_demo/formatConvertion/scanpy/S135TLD1_addSCimg.h5ad")
adata

CPU times: user 110 ms, sys: 403 ms, total: 513 ms
Wall time: 511 ms

AnnData object with n_obs × n_vars = 2322 × 24289
    obs: 'X_index', 'nCount_Spatial', 'nFeature_Spatial', 'orig.ident', 'in_tissue', 'array_row', 'array_col'
    var: 'n_cells', 'n_counts', 'mean_umi'
    uns: 'spatial'
    obsm: 'X_spatial', 'spatial'

%%time
rsc.get.anndata_to_GPU(adata)

CPU times: user 59.5 ms, sys: 14.3 ms, total: 73.8 ms
Wall time: 71.5 ms

%%time
rsc.pp.flag_gene_family(adata, gene_family_name="MT", gene_family_prefix="mt-")

CPU times: user 11.6 ms, sys: 1.35 ms, total: 13 ms
Wall time: 12 ms

%%time
rsc.pp.calculate_qc_metrics(adata, qc_vars=["MT"])

CPU times: user 17.1 ms, sys: 7.03 ms, total: 24.1 ms
Wall time: 22.5 ms

adata

AnnData object with n_obs × n_vars = 2322 × 24289
    obs: 'X_index', 'nCount_Spatial', 'nFeature_Spatial', 'orig.ident', 'in_tissue', 'array_row', 'array_col', 'n_genes_by_counts', 'total_counts', 'log1p_n_genes_by_counts', 'log1p_total_counts', 'total_counts_MT', 'pct_counts_MT', 'log1p_total_counts_MT'
    var: 'n_cells', 'n_counts', 'mean_umi', 'MT', 'n_cells_by_counts', 'total_counts', 'mean_counts', 'pct_dropout_by_counts', 'log1p_total_counts', 'log1p_mean_counts'
    uns: 'spatial'
    obsm: 'X_spatial', 'spatial'

%%time
adata = adata[adata.obs["n_genes_by_counts"] >= 50] # filter out low quality cells
adata = adata[adata.obs["n_genes_by_counts"] <= 80000]
adata

CPU times: user 82.2 ms, sys: 134 ms, total: 216 ms
Wall time: 214 ms

View of AnnData object with n_obs × n_vars = 2316 × 24289
    obs: 'X_index', 'nCount_Spatial', 'nFeature_Spatial', 'orig.ident', 'in_tissue', 'array_row', 'array_col', 'n_genes_by_counts', 'total_counts', 'log1p_n_genes_by_counts', 'log1p_total_counts', 'total_counts_MT', 'pct_counts_MT', 'log1p_total_counts_MT'
    var: 'n_cells', 'n_counts', 'mean_umi', 'MT', 'n_cells_by_counts', 'total_counts', 'mean_counts', 'pct_dropout_by_counts', 'log1p_total_counts', 'log1p_mean_counts'
    uns: 'spatial'
    obsm: 'X_spatial', 'spatial'

%%time
adata = adata[adata.obs["pct_counts_MT"] < 20]

CPU times: user 53.5 ms, sys: 71 ms, total: 124 ms
Wall time: 122 ms

%%time
rsc.pp.filter_genes(adata, min_count=3) # filter out genes that are expressed in less than 3 cells

filtered out 2616 genes based on n_cells_by_counts
CPU times: user 232 ms, sys: 538 ms, total: 770 ms
Wall time: 768 ms

adata

AnnData object with n_obs × n_vars = 2316 × 21673
    obs: 'X_index', 'nCount_Spatial', 'nFeature_Spatial', 'orig.ident', 'in_tissue', 'array_row', 'array_col', 'n_genes_by_counts', 'total_counts', 'log1p_n_genes_by_counts', 'log1p_total_counts', 'total_counts_MT', 'pct_counts_MT', 'log1p_total_counts_MT'
    var: 'n_cells', 'n_counts', 'mean_umi', 'MT', 'n_cells_by_counts', 'total_counts', 'mean_counts', 'pct_dropout_by_counts', 'log1p_total_counts', 'log1p_mean_counts'
    uns: 'spatial'
    obsm: 'X_spatial', 'spatial'

# store the raw expression counts (after filtering) in the `.layer["counts"]`
adata.layers["counts"] = adata.X.copy()
adata.shape

(2316, 21673)

adata

AnnData object with n_obs × n_vars = 2316 × 21673
    obs: 'X_index', 'nCount_Spatial', 'nFeature_Spatial', 'orig.ident', 'in_tissue', 'array_row', 'array_col', 'n_genes_by_counts', 'total_counts', 'log1p_n_genes_by_counts', 'log1p_total_counts', 'total_counts_MT', 'pct_counts_MT', 'log1p_total_counts_MT'
    var: 'n_cells', 'n_counts', 'mean_umi', 'MT', 'n_cells_by_counts', 'total_counts', 'mean_counts', 'pct_dropout_by_counts', 'log1p_total_counts', 'log1p_mean_counts'
    uns: 'spatial'
    obsm: 'X_spatial', 'spatial'
    layers: 'counts'

%%time
rsc.pp.normalize_total(adata, target_sum=1e4)
rsc.pp.log1p(adata)
rsc.pp.highly_variable_genes(
    adata,
    n_top_genes=5000,
    flavor="seurat_v3",
    #batch_key="PatientNumber",
    layer="counts",
)

CPU times: user 78 ms, sys: 11.8 ms, total: 89.8 ms
Wall time: 87 ms

adata.raw=adata
adata = adata[:, adata.var["highly_variable"]]

%%time
rsc.pp.scale(adata, max_value=10)

CPU times: user 7.15 ms, sys: 1.38 ms, total: 8.53 ms
Wall time: 6.65 ms

%%time
rsc.pp.pca(adata, n_comps=50)

CPU times: user 803 ms, sys: 235 ms, total: 1.04 s
Wall time: 1.05 s

%%time
rsc.pp.neighbors(adata, n_neighbors=15, n_pcs=50)

CPU times: user 97.7 ms, sys: 94.6 ms, total: 192 ms
Wall time: 190 ms

%%time
rsc.tl.umap(adata)

CPU times: user 123 ms, sys: 86.8 ms, total: 210 ms
Wall time: 206 ms

%%time
rsc.tl.leiden(adata, resolution=0.6,key_added="clusters")

CPU times: user 45.9 ms, sys: 15.6 ms, total: 61.4 ms
Wall time: 59.2 ms

adata

AnnData object with n_obs × n_vars = 2316 × 5000
    obs: 'X_index', 'nCount_Spatial', 'nFeature_Spatial', 'orig.ident', 'in_tissue', 'array_row', 'array_col', 'n_genes_by_counts', 'total_counts', 'log1p_n_genes_by_counts', 'log1p_total_counts', 'total_counts_MT', 'pct_counts_MT', 'log1p_total_counts_MT', 'leiden', 'louvain', 'clusters'
    var: 'n_cells', 'n_counts', 'mean_umi', 'MT', 'n_cells_by_counts', 'total_counts', 'mean_counts', 'pct_dropout_by_counts', 'log1p_total_counts', 'log1p_mean_counts', 'highly_variable', 'highly_variable_rank', 'means', 'variances', 'variances_norm', 'mean', 'std'
    uns: 'spatial', 'log1p', 'hvg', 'pca', 'neighbors', 'umap', 'leiden', 'louvain'
    obsm: 'X_spatial', 'spatial', 'X_pca', 'X_umap'
    varm: 'PCs'
    layers: 'counts'
    obsp: 'distances', 'connectivities'

%%time
rsc.tl.louvain(adata, resolution=0.6)

CPU times: user 58.9 ms, sys: 17.8 ms, total: 76.7 ms
Wall time: 74.5 ms

help(rsc.tl.louvain)

Help on function louvain in module rapids_singlecell.tools._clustering:

louvain(adata: 'AnnData', resolution: 'float' = 1.0, *, restrict_to: 'tuple[str, Sequence[str]] | None' = None, key_added: 'str' = 'louvain', adjacency: 'sparse.spmatrix | None' = None, n_iterations: 'int' = 100, threshold: 'float' = 1e-07, use_weights: 'bool' = True, neighbors_key: 'int | None' = None, obsp: 'str | None' = None, copy: 'bool' = False) -> 'AnnData | None'
    Performs Louvain clustering using cuGraph, which implements the method
    described in:
    
    Blondel, V.D., Guillaume, J.-L., Lambiotte, R., & Lefebvre, E. (2008).
    Fast unfolding of community hierarchies in large networks, J. Stat.
    Mech., P10008. DOI: 10.1088/1742-5468/2008/10/P10008
    
    Parameters
    ----------
        adata :
            annData object
    
        resolution
            A parameter value controlling the coarseness of the clustering
            (called gamma in the modularity formula). Higher values lead to
            more clusters.
    
        restrict_to
            Restrict the clustering to the categories within the key for
            sample annotation, tuple needs to contain
            `(obs_key, list_of_categories)`.
    
        key_added
            `adata.obs` key under which to add the cluster labels.
    
        adjacency
            Sparse adjacency matrix of the graph, defaults to neighbors
            connectivities.
    
        n_iterations
            This controls the maximum number of levels/iterations of the
            Louvain algorithm. When specified the algorithm will terminate
            after no more than the specified number of iterations. No error
            occurs when the algorithm terminates early in this manner.
            Capped at 500 to prevent excessive runtime.
    
        threshold
            Modularity gain threshold for each level/iteration. If the gain
            of modularity between two levels of the algorithm is less than
            the given threshold then the algorithm stops and returns the
            resulting communities. Defaults to 1e-7.
    
        use_weights
            If `True`, edge weights from the graph are used in the
            computation (placing more emphasis on stronger edges).
    
        neighbors_key
            If not specified, `louvain` looks at `.obsp['connectivities']`
            for neighbors connectivities. If specified, `louvain` looks at
            `.obsp[.uns[neighbors_key]['connectivities_key']]` for neighbors
            connectivities.
    
        obsp
            Use `.obsp[obsp]` as adjacency. You can't specify both `obsp`
            and `neighbors_key` at the same time.
    
        copy
            Whether to copy `adata` or modify it in place.