Methods

ARCHS4 expression data and gene annotation

Human gene-level expression data was obtained from ARCHS4. The expression matrix contained genes as rows and RNA sequencing samples as columns. Sample metadata distributed within ARCHS4 was used to distinguish bulk RNA sequencing samples from samples predicted to represent single-cell experiments. Samples were classified as bulk RNA sequencing when their single-cell probability was below 0.05.

Gene annotations were used to split the matrix according to gene biotype. Rows annotated as snoRNAs were assigned to the snoRNA expression matrix. Rows annotated as protein-coding genes were assigned to the protein-coding expression matrix. Both matrices retained the same sample columns but contained different gene rows.

Gene symbols, Ensembl gene identifiers, biotype assignments, and selected sample indices were preserved. These identifiers were retained in a consistent order across the raw-expression, normalized-expression, correlation, pathway-ranking, and overlapping-gene files.

Bulk-sample selection and expression filtering

Sample selection was performed after excluding samples that did not satisfy the bulk RNA sequencing criterion (Single-cell probability > 0.05). For every retained bulk sample, the mean snoRNA expression was calculated across all annotated snoRNA rows. The median of the resulting sample-level mean snoRNA expression values was used as the sample-selection threshold. Samples with mean snoRNA expression above this median were included in the eligible high-snoRNA-expression sample pool.

A random subset of 10,000 samples was selected without replacement from the eligible sample pool. A fixed random seed was used so that the same sample indices could be reproduced in repeated executions of the workflow. The selected sample indices were used for both the snoRNA and protein-coding matrices.

Feature filtering was performed within the selected sample subset. A snoRNA or protein-coding gene was retained only when its raw expression value was greater than zero in at least 50 selected samples. This filter was applied separately to the snoRNA and protein-coding matrices.

The original nonzero-expression status of each retained feature in every selected sample was preserved as a Boolean expression mask. These masks are later used to define the samples included in each pairwise snoRNA-gene correlation calculation.

Expression transformation and normalization

Retained expression values were transformed using a base-2 logarithm after adding one to each value. Adding one allowed zero-valued observations to remain defined after the logarithmic transformation. The same transformation was applied to both the snoRNA and protein-coding expression matrices.

Quantile normalization was applied after the logarithmic transformation. Within each expression matrix, values in every selected sample were sorted by rank. The mean value for each rank was calculated across samples. Each observation was then replaced by the mean value corresponding to its rank before values were returned to their original gene ordering.

The resulting sample columns had a common empirical expression distribution. The normalized snoRNA matrix, normalized protein-coding matrix, expression masks, gene identifiers, gene symbols, and selected sample indices were stored in an HDF5 file for downstream analysis.

Pairwise snoRNA to protein-coding gene coexpression

Pearson correlation coefficients were calculated between every retained snoRNA and every retained protein-coding gene. Each calculation was restricted to samples in which both features had nonzero raw expression.

Samples in which either member of the pair had zero raw expression were excluded from the calculation for that specific pair. The overlapping sample set therefore differed between snoRNA-gene pairs according to their expression patterns.

For each valid pair, Pearson correlation was calculated using the normalized expression values from the shared nonzero samples.

The Pearson correlation coefficient, corresponding p-value, and number of shared nonzero samples were recorded for every valid pair. The resulting HDF5 correlation file contained protein-coding genes as rows and snoRNAs as columns.

Gene-set libraries and term matching

Functional terms were obtained from 6 gene set libraries from Enrichr. The analysis included the following libraries:

• ChEA 2022
• GWAS Catalog 2025
• MGI Mammalian Phenotype Level 4 2024
• KEGG 2026
• WikiPathways 2024 Human
• Gene Ontology Biological Process 2025

Each library consisted of named terms and the protein-coding gene symbols assigned to those terms. Gene symbols were standardized between the Enrichr libraries and the correlation matrix before matching. For every term, the original gene list was intersected with the protein-coding genes represented in the snoRNA-gene correlation matrix.

The original gene-set size and the number of matched protein-coding genes were recorded separately. The original gene-set size represents the complete number of genes supplied for the term by the source library. The overlap size represents the number of genes from that term that were present in the processed protein-coding correlation matrix.

Terms with fewer than 5 overlapping protein-coding genes were excluded from functional scoring. This minimum-overlap criterion was applied independently to every term in every library.

Functional association scoring and ranking

For each retained snoRNA and each eligible gene-set term, correlations for the matched protein-coding genes were extracted from the correlation matrix.

The snoRNA-term association score was calculated as the unweighted arithmetic mean of the available snoRNA-gene Pearson correlations for those matched genes.

Every matched gene-level correlation contributed equally to the term score. No gene-specific weighting, enrichment p-value, pathway-size correction, or variance-based adjustment was applied to the mean.

The pathway score table contained the snoRNA symbol, snoRNA Ensembl identifier, library name, term name, overlap size, original gene-set size, matched matrix-row count, overlapping gene symbols, and mean Pearson correlation.

Terms were ranked separately for each snoRNA and gene-set library. Descending sorting produced the highest-correlated ranking. Ascending sorting produced the lowest-correlated ranking. Direction and rank were stored with every selected term in the final ranking file.

The current web report displays the 10 highest-correlated terms across all included libraries and the 10 highest-correlated terms within each individual library.

Overall term selection and overlapping-gene reports

Overall highest-correlated terms were selected by collecting the highest-direction term rankings across all included libraries for a queried snoRNA. The collected terms were then sorted by mean Pearson correlation in descending order. Duplicate library-term combinations were removed, and the first 10 terms were retained.

Library-specific displays were generated independently. For each library, terms in the highest direction were sorted by mean Pearson correlation in descending order. Duplicate term names were removed, and up to 10 terms were retained.

Gene-level results were generated for each selected overall term. Protein-coding genes assigned to the term were matched to the correlation matrix, and their pairwise correlations with the selected snoRNA were retrieved.

When multiple matrix rows mapped to the same gene symbol, the available correlations were averaged to produce one reported value for that gene symbol. Correlation distance was calculated by subtracting the Pearson correlation from one.

Genes were ordered by ascending correlation distance. This ordering is equivalent to sorting by descending Pearson correlation. Duplicate gene symbols were removed after gene-level aggregation, and up to 25 genes were retained for each term.

The resulting overlapping-gene tables contain gene rank, gene symbol, Pearson correlation, and correlation distance. The overall and library-specific term tables contain term rank, library or term name, mean Pearson correlation, overlap size, and original gene-set size.

Tables and figures are are built in the browser. Table downloads are written as UTF-8 CSV files. Figure downloads are rendered as PNG files with a white background.

Perturb-Seqr analysis

Positively correlated protein-coding genes for each snoRNA were submitted to the Columbia Perturb-Seqr API. Up to 500 positively correlated genes were used when at least that many were available; otherwise, all available positively correlated genes were submitted.

Each supported perturbation library was queried separately so that a snoRNA could retain results from libraries that returned matches even when other libraries returned no results. Unrecognized submitted identifiers were retained in the query gene list rather than removed before submission.

The report groups results by perturbation type and library, then separates up-regulated and down-regulated perturbations. Each library panel reports query status, query-gene count, total enrichment-result count, perturbation rank, nominal and adjusted p-values, odds ratio, overlap, and significance status.

The website reads the combined all_snornas_perturbseqr_results.csv file. The current analysis retains up to 10 up-regulated and 10 down-regulated results per snoRNA-library query, for a maximum of 20 displayed results per successful library.

Software and data storage

Data extraction and numerical processing were implemented in Python. HDF5 files were read and written using h5py. Matrix operations were performed using NumPy. Tabular processing, identifier matching, grouping, ranking, and CSV generation were performed using pandas.

Pearson statistics were calculated using SciPy.

Large expression and correlation matrices were stored in HDF5 format to permit incremental access without loading the complete source dataset into memory. Final snoRNA-term rankings and ranked overlapping-gene results were stored as CSV files.