Data initialization with setup_adata

In contrast to default Scanpy, where preprocessing parameters must be specified at each call and not centrally stored, ScToolkit introduces a standardized initialization step via the setup_adata function. This augments the AnnData object with: 1.

Configuration dictionary (adata.uns[“config”])

A unified registry covering quality control, normalization, clustering, plotting, and downstream analyses parameters. It encodes both Scanpy arguments with proper defaults and project-specific arguments, including curated marker-gene handling, organism-specific geneset handling, and file output policies. GPU acceleration is integrated into the config, enabling automatic switching between CPU and GPU backends for PCA, nearest-neighbor search, and clustering.

Existing configurations and statistics are preserved unless explicitly overwritten, facilitating reproducibility and iterative refinement of analyses. By consolidating the config, state-tracking, and GPU initialization into a single entry point, setup_adata transforms a raw AnnData object into a fully contextualized analysis container. This contrasts with unmodified Scanpy, where users must manually manage parameters, track preprocessing steps, and configure hardware acceleration. From now on if config is stated, it refers to adata.uns[‘config’].

Quality control preprocessing with all_qc_preprocessing

After initialization, ScToolkit provides a dedicated entry point for systematic quality control via the all_qc_preprocessing function. This standardizes QC annotation, reduces defective cell likelihood, and ensures reproducibility and interoparability across datasets.

Key operations include: •

Gene family annotation using curated organism-specific definitions for mitochondrial, ribosomal, hemoglobins, pseudogenes, and immunoglobulin/T-cell receptor loci. Gene families are projected into adata.var as boolean flags for downstream filtering and visualization. Databases bundled with the package (Ensembl-GTF-derived, ¹⁵ HGNC/MGI ^{16, 17} ) curated categories, and housekeeping gene sets ¹⁸ ) extend Scanpy’s defaults by systematically tracking additional biologically relevant gene families.

Pre-clustering normalization with Run_all_prep_steps_clustering

Following QC, Run_all_prep_steps_clustering provides centralized preprocessing and normalization. This integrates normalization, log transformation, scaling, and gene selection into a reproducible pipeline controlled entirely by the config.

Extensions beyond Scanpy include: •

Mixture of Regex-based and gene-list gene exclusion and centralized control of preprocessing steps.

By consolidating these steps, the function guarantees consistent inputs for clustering, avoiding redundant calls and user errors.

Clustering orchestration with cluster_SC_scanpy_like

cluster_SC_scanpy_like structures PCA, neighborhood graph construction, Leiden clustering, and UMAP embedding in a single highly configurable function, providing a failsafe and streamlined workflow beyond Scanpy’s separate steps. Each step can be perfomed on a CPU or GPU.

Enhancements include: •

PCA with variance-based selection of informative components, capped by data rank, to estimate the appropriate number of components, n_pcs.

This ensures reproducible and biologically meaningful clustering without manual trial-and-error.

Downstream analysis and visualization

run_downstream integrates core post-clustering functionality: •

Differential expression (DEG) with parallelized Wilcoxon tests (standard z-score approximation and JIT-compiled permutation-based variants), and group-specific DEG handling.

These downstream steps extend Scanpy’s minimal plotting and DEG support into a comprehensive, automated framework for reproducible structured interpretation.

ScToolkit is distributed as a Docker container bundling both Python and R environments to ensure reproducibility and straightforward deployment across platforms as far as possible for GPU usage. For users preferring local installation, the GitHub repository provides multiple conda environment definitions (.yml files) tailored to different use cases: •

Default environment: GPU-enabled with RAPIDS and CUDA libraries preconfigured.

Documentation is generated with Sphinx and kept synchronized with the source code, ensuring that both API references and tutorials remain current. The repository includes extensive Jupyter notebooks that demonstrate full end-to-end analyses, from raw data initialization to clustering and downstream interpretation. Additional notebooks showcase advanced features such as curated gene family handling, GPU acceleration, and automated DEG export.

This dual strategy, providing both containerized execution and conda-based environments, ensures accessibility for beginners while maintaining flexibility for advanced users and developers. Together with reproducible notebooks, this structure guarantees that ScToolkit can be deployed consistently in diverse computational settings, from laptops to GPU clusters.

For use case analysis, we accessed publicly available single-cell RNA sequencing (scRNA-seq) resources. Specifically, we retrieved datasets from the DISCO database ( https://disco.bii.a-star.edu.sg/) ¹⁹ including scRNA-seq profiles from human tonsil (n = 92,114 cells) and peripheral blood (n = 169,686 cells). In addition, we obtained a tonsil scRNA-seq dataset from the Tonsil Atlas database ¹¹ encopassing 462,352 cells. The full use-case analysis described below is available in the Supplementary Data as a PDF export of the corresponding Jupyter Notebooks.

Next, we applied a rigorous quality control workflow to isolate high-quality single-cell transcriptomes and used violin plots to assess data quality across conditions and tissues ( Figure 1A). Cells were excluded if they expressed fewer than 500 genes, fewer than 700 total counts, or more than 18,000 counts. To minimize the inclusion of low-complexity libraries, cells with fewer than 10 unique genes were removed. Similarly cells with more than 5,000 genes were removed to prevent doublets. Cells with >10% mitochondrial content were excluded to reduce bias from dying or stressed cells. In addition, cells with <1.5% ribosomal transcripts or >1% HGNC-annotated pseudogene expression were filtered out ( Figure 1A). At the gene level, only features expressed in at least 2,000 cells, with a minimum of seven unique counts, a minimum of 2,000 total counts, and at least 10 maximum counts across cells, were retained. Beyond the standard parameters (n_genes and n_counts), we included n_unique as a measure of distributional spread, providing a more general indicator of transcript diversity than mean- or variance-based metrics that assume normality. The percentages of housekeeping (pct_HK_gtex), intergenic (pct_GTF_IG_genes), and pseudogene (pct_HGNC_pseudogene) expression offer complementary insights into sample composition and annotation fidelity. The displayed metrics represent a selected subset of all available quality control measures, shown here to demonstrate the extended functionality and interpretability of the workflow. In Figure 1A, only the top 10% quantile range of each metric is displayed, emphasizing the high end of the distributions where duplicate or outlier cells typically accumulate. This focused view enables refined inspection of potential overrepresented or low-complexity transcriptomes that define the upper QC thresholds. The shown metrics represent a subset of all available QC measures, selected to demonstrate the workflow’s extended functionality and interpretability.

The results reveal that parameters such as blood of solid tumor patient (BOSTP), non-tumor disease tissue (NTDT), and Tonsil removed during surgery for oropharyngeal papillomatosis (TRDSFOP) show generally lower levels across n_genes, n_counts, and n_unique, lowering the expectation of doublets within these samples. In contrast, the consistent percent of house-keeping genes (pct_HK_gtex) values across all conditions indicate viable and transcriptionally active cells, suggesting comparable baseline activity among cell types. The pct_GTF_IG_genes distributions show a distinct multimodal pattern, where the upper mode indicates the presence of plasma cells. Finally, elevated pct_HGNC_pseudogene values in BONTDP, BOSTP, NTDT, and Disco_blood_normal flag the presence of subsets of low-quality or degraded cells. The remaining groups display lower pseudogene percentages consistent with higher-quality transcriptomic profiles.

To explore the initial annotation of cell clusters, we utilized a curated library of cell type -specific gene signatures from four complementary resources: PanglaoDB, ²⁰ CellMarker, ²¹ SCType, ²² and CellTypist (v2). ²³ To increase tissue relevance, each database was then subset signatures associated with tonsil and blood, ensuring that the resulting markers reflected the diversity of resident immune cell populations within these compartments. This process demonstrates the database integration and refinement functionality of the workflow, we used three core functions: get_ref_db(), search_database(), and get_element_and_counts_with_positional_sum(). First, get_ref_db() collects all available cell type-specific genesets from the supported databases (PanglaoDB, CellTypist, CellMarker, and SCType). Next, search_database() extracts the subset of genesets matching a specified cell type, defined by the user using regular expression.

Notably, substantial heterogeneity was observed across gene-signatures for the same annotated cell type. To address this, we generated consensus signatures by integrating overlapping features. We calculated the frequency of occurrence of each gene across signatures for a given cell type, ranked genes by their frequency, and defined the top 20 most recurrent genes as the consensus signature ( Figure 1B). Specifically, we used get_element_and_counts_with_positional_sum() to rank all genes by their recurrence across the retrieved genesets, generating consensus signature by taking the top 20 genes. As an exampler, we retrieved cross-database overlap of B-cell gene sets and subsequently generated their consensus signature. For each B cell signature from the databases and the B cell consensus signature, the summed Jaccard index was calculated against all nonB cell gene sets in the reference databases. Results showed that our B cell consensus signature exhibits a minimum similarity (relatively low Jaccard index) to nonB cell signatures compared to the database redived B cell signatures ( Figure 1B). This prevents the selection of unspecific gene-set signature from many existing signature ressources. The refined scToolkit reference represents a consensus-based, high-confidence collection of cell-type-specific gene signatures, illustrating the workflow’s ability to generalize this procedure to any cell-type-specific refinement task.

Following quality control, we performed unsupervised clustering on the combined dataset using the Leiden algorithm. Cells were projected into UMAP space, where coloring by tissue origin revealed intermixed clustering of cells from Blood and Tonsil ( Figure 2A). The UMAP density distributions show a continuous embedding without visible batch effects between blood- and tonsil-derived cells. Cells from both tissues occupy overlapping yet structured regions, indicating comparable transcriptomic profiles and proper integration of datasets ( Figure 2A). To assign cell identities, we applied single-cell level gene set enrichment analysis using over-representation analysis algorithm (ORA) with our consensus gene signatures. Each cell was assigned the identity for which its transcriptome showed the maximum ORA enrichment score (Idxmax). The major cell types identified included B cells, dendritic cells (DCs), granulocytes, macrophages, NK cells, T cells, and plasma cells ( Figure 2B). To validate these annotations, we compared our results with reference annotations (referred to as ground-truth annotation from hereafter) provided by the DISCO and Tonsil Atlas databases ( Figure 2C). The idxmax ORA assignments align with the broad ground-truth annotation ( Figure 2C), demonstrating that clusters capture major immune and myeloid populations while maintaining distinct lineage boundaries. Notably, cell types of similar identity cluster in close proximity within the joint embedding. This confirms that the refinement and preprocessing preserve biological relationships rather than dataset origin. For instance, cells belonging to the same category (e.g., B cells or T cells) cluster together regardless of tissue origin. This highlights the robustness and reliability of our workflow, demonstrating its ability to successfully integrate datasets from different tissues without the need for specialized batch-correction algorithms ( Figure 2A-B).

To quantitatively validate/support our annotations, we computed the Jaccard similarity index between our ORA-based assignments and the reference (ground-truth) labels. The resulting heatmap ( Figure 2D) quantitatively supports this correspondence, with high overlap values between identical or closely related cell types and minor off-diagonal similarities reflecting lineage proximity. Concordance as reflected by the Jaccard index ( Figure 2D), demonstrated the strongest agreement (Jaccard index = 1) for the most abundant populations, including B cells, T/NK cells, granulocytes, plasma cells, DCs, and macrophages. Of note, in ground-truth annotation T and NK cells were combined into a single category, and kept separate in T and NK cells in Idxmax ORA annotations. Interestingly, cells classified as “unknown” within Idxmax ORA annotations, showed a distribution of Jaccard Scores across many ground-truth cell types with the highest score for T/NK cell categroy. To annotate the clusters we selected the predominant cell types of the ground truth ( Figure 2E). The Leiden clustering labeled with ground-truth cell types ( Figure 2E) further shows that clusters are composed of cells from both tissues and databases, confirming that integration was successful and not driven by dataset-specific biases.

To investigate the functional roles of the annotated cell types, we performed a pathway enrichment analysis based on differentially expressed genes (DEGs). DEGs were first identified by comparing each cell type specific cluster against all remaining clusters and stringently filtered to ensure robustness. DEGs were retained only if expressed in at least 15% of cells in one of the compared cell type resolved cluster, with an adjusted p-value < 0.0001 and an absolute log ₂ fold change ≥ 0.5. Following filtering, DEGs are intersected with hallmark gene sets from MSigDB. For each pathway, the z-scores of Wilcoxon rank-test u-statistics from overlapping genes are summed and scaled to –1 to 1. This analysis revealed distinct pathway overlap patterns across cell types. Myeloid cell cluster showed strong overlap of inflammatory signaling pathways, including inflammatory response, TNF-α signaling, and complement activation ( Figure 3A). Notably, differences also emerged between T-cell subclusters where the C3: T-cell cluster showed selective overlap of Allograft Rejection and TNF-α signaling genes, which were absent in C1: T-cell cluster, suggesting functional specialization among T-cell subtypes ( Figure 3A).

Another visualization option is to show the top DEGs for each cluster as split dot plots, where the log2 fold change is color coded ( Figure 3B). The percent expression of DEGs is displayed as split circles, with the left half representing the expression of the cluster of interest and the right half representing expression of all other cells combined. Corroborating with earlier findings, DEG analysis revealed two transcriptionally distinct T cell populations. C3: T-cell cluster was characterized by high expression of canonical T cell receptor (TCR) components ( CD3E, CD3D, CD3G) together with IL32, consistent with a population of conventional, activated T-cells ( Figure 3B). In contrast, the C1: T-cell cluster lacked strong expression of CD3 genes and was instead enriched for signaling-related transcripts, including INPP4B (PI3K/AKT signaling), RASGRF2 (Ras/MAPK activation), and PLCL1 (calcium signaling regulation), as well as CATSPERB and OPRM1, suggesting a specialized, non-classical functional state ( Figure 3B). The two Bcell clusters (C0: B-cell and C4: B-cell) also showed differential marker experssion. C4: B-cell cluster express genes characteristic of mature B-cells, including B-cell receptor (BCR) components ( HLA-DRA , CD79B, CD79A), a B-cell surface antigen ( MS4A1), and a molecule involved in antigen presentation that’s also expressed in B cells ( CD74) ( Figure 3B). C0: B-cell cluster consists of genes related to B-cell activation, survival, and differentiation, with KHDRBS2 and TRIO involved in signaling, CDK14 and PRKCE playing roles in cell cycle and signaling respectively, and KCNQ5 being less specific to B-cells but also involved in cell processes ( Figure 3B).

Next, we investigated compartment-specific transcriptional differences within the same immune cell types. For this analysis, we selected clusters that contained a comparable number of cells from both blood and tonsil, including C2: Myeloid-cell, C3: T-cell, and C4: B-cell clusters. Differential expression analysis (DEG) was then performed between the two compartments (blood and tonsil) to identify genes driving tissue-specific adaptations within each cell type. DEG analysis revealed a distinct enrichment of genes associated with immune regulation and tissue adaptation in the tonsil ( Figure 3C). A comparative analysis of myeloid cell clusters from blood and tonsil revealed distinct transcriptional signatures. Blood myeloid cells were characterized by high expression of S100A4, S100A6, S100A8, and S100A9, genes, which are associated with inflammation, cytoskeletal dynamics, and acute immune activation ( Figure 3C). In contrast, tonsil myeloid cells showed elevated levels of HLA-DRA , HLA-DPB1, and HLA-DQA1, reflecting enhanced antigen-presenting capacity, together with increased expression of SERPINF1, and DNASE1L3, indicative of tissue adaptation and immune regulatory functions ( Figure 3C). Similarly, DEG analysis revealed distinct transcriptional features between blood and tonsillar T-cells. Compared with tonsillar T cells, blood-derived T cells showed higher expression of the interferon-stimulated genes IFITM1 and RARRES3, together with increased AES, PTPRCAP, and SEPT7, indicative of an interferon-responsive and activated circulating T-cell state ( Figure 3C). In contrast, tonsillar T cells exhibited elevated levels of TRAC and IGKC, reflecting active T-cell receptor signaling and local immune activation within the lymphoid tissue ( Figure 3C). Analysis of distinct transcriptional profiles between blood and tonsillar B-cells revealed that blood B-cells exhibited higher expression of JUN, JUNB, DUSP1, and TSC22D3, indicating activation of immediate-early response and stress-regulatory pathways characteristic of circulating immune cells ( Figure 3C). In contrast, tonsillar B-cells showed increased expression of ACTG1, consistent with enhanced cytoskeletal remodeling and local activation within the lymphoid microenvironment ( Figure 3C). Together, these results indicate that tonsil immune cells adopt a tissue-adapted, immunoregulatory phenotype that differs markedly from blood-derived immune populations.

Analysis of HLA gene expression revealed compartment-specific patterns across immune cell types ( Figure 3D). HLA-C , a class I molecule, showed elevated expression in blood myeloid, B-cells, and T-cells, reflecting stronger class I–mediated antigen presentation in circulating immune populations. In contrast, HLA-DRA and HLA-DQA1, both class II molecules, were highly expressed in tonsillar myeloid cells ( Figure 3D), consistent with enhanced antigen-presenting activity within the lymphoid tissue. These findings highlight differential regulation of antigen presentation pathways between blood and tonsil compartments.