Scanorama

Scanorama is designed for the integration and batch correction of heterogeneous scRNA-seq datasets. It is intended to be used in pipelines after noise-reduction steps (like imputation or highly-variable gene filtering).

There are two main modes of operation:

1. Integration: Focuses on finding low-dimensional embeddings (X_scanorama) that align different datasets, useful for clustering and visualization.
2. Batch Correction: Transforms the actual cell-by-gene expression matrices (adata.X) to remove batch effects while preserving biological signals.

For large-scale datasets, Scanorama can be accelerated using geometric sketching (via the geosketch tool) to compactly summarize the transcriptomic landscape.

You can install Scanorama using pip. If you are using an Anaconda environment, you must also install python-annoy from the conda-forge channel to ensure compatibility.

# Standard installation
pip install scanorama

# If using Anaconda, install annoy first
conda install -c conda-forge python-annoy

The param_sensitivity.py script provides a framework for evaluating how different Scanorama hyperparameters affect the quality of data integration. It uses the Silhouette Coefficient as a metric to quantify how well-clustered the integrated data is.

Users can test the sensitivity of the following assemble parameters:

• knn: Number of nearest neighbors used in the integration process.
• sigma: The bandwidth parameter for the kernel.
• alpha: A parameter controlling the integration strength.
• approx: Whether to use approximate nearest neighbors (set to False for exact NN testing).

For each parameter tested, the script performs a t-test (ttest_ind) comparing the current parameter's performance against previous results and generates a boxplot visualization saved as an .svg file (e.g., param_sensitivity_knn.svg).

After integration, the quality of downstream visualization (like t-SNE) can be assessed by testing the sensitivity of t-SNE hyperparameters. The script evaluates how these parameters affect the Silhouette Coefficient of the resulting embedding.

Parameters tested include:

• perplexity: Controls the balance between local and global aspects of the data.
• learn_rate: The learning rate for the t-SNE optimization.

Results are visualized using boxplots saved as .svg files (e.g., param_sensitivity_perplexity.svg).

Common issues and solutions:

• MemoryError during large integration: Lower the batch_size parameter or use the sketch parameter in integrate() to enable sketch-based acceleration.
• Illegal instruction or Segfault: This is often caused by recent versions of the annoy package. Fix this by passing approx=False to the function to use scikit-learn's nearest neighbors matching instead.
• Incorrect Data Shape: Ensure your input matrices are cells-by-genes, not the transpose.
• t-SNE performance: The t-SNE visualization step can be very slow and memory-intensive for >100k cells. Consider using faster implementations like bhtsne, Multicore-TSNE, or net-SNE in your pipeline.

Scanorama provides wrappers for seamless integration with scanpy.AnnData objects.

• integrate_scanpy(adatas): Adds an entry to adata.obsm called 'X_scanorama' for each object in the list. This contains the low-dimensional embeddings for downstream analysis like KNN graph construction or visualization.
• correct_scanpy(adatas): Creates new AnnData objects where adata.X is replaced with the Scanorama-transformed cell-by-gene matrix, while preserving other metadata.
• correct_scanpy(adatas, return_dimred=True): Performs both integration and batch correction.

import scanpy
import scanorama

# adatas: list of scanpy.AnnData

# Integration
scanorama.integrate_scanpy(adatas)

# Batch correction
corrected = scanorama.correct_scanpy(adatas)

# Integration and batch correction
corrected = scanorama.correct_scanpy(adatas, return_dimred=True)

You can call Scanorama from R using the reticulate package.

Important: Because reticulate has difficulty returning sparse matrices, you should set return_dense=TRUE when using the correct() method. This returns the corrected data as R matrix objects, though it will increase memory usage for large datasets.

# datasets: list of list of matrix
# genes_list: list of list of list of string

library(reticulate)
scanorama <- import('scanorama')

# Integration
integrated.data <- scanorama$integrate(datasets, genes_list)

# Batch correction (use return_dense=TRUE for R compatibility)
corrected.data <- scanorama$correct(datasets, genes_list, return_dense=TRUE)

# Integration and batch correction
integrated.corrected.data <- scanorama$correct(datasets, genes_list,
                                               return_dimred=TRUE, return_dense=TRUE)

The load_alignments_from_log function parses log files containing alignment results to reconstruct the alignment matrix and dataset names.

Expected Log Format:

• Lines starting with Loaded indicate the name of a dataset (e.g., Loaded dataset_name).
• Data rows are enclosed in brackets [...] and contain space-separated floating-point values.
• The function handles multi-line rows and ignores lines containing t-SNE within brackets.

Returns:

• A: A numpy array representing the alignment matrix.
• names: A list of strings containing the names of the datasets.

Parameters:

• fname: Path to the log file.

from scanorama import plt

# Example usage
A, names = load_alignments_from_log('path/to/alignment_log.txt')

The time_align_correlate function calculates the correlation between temporal distance and alignment scores. It computes the Euclidean distance between time points and compares it against the provided alignment matrix to determine if datasets that are closer in time are also more similar in their alignment scores.

It prints both the Spearman rho and Pearson rho values to the console.

Parameters:

• alignments: A square matrix of alignment scores.
• time: An array of time values corresponding to the datasets in the alignment matrix.

Requirement: The shape of the time_dist (derived from time) must match the shape of alignments.

from scanorama import plt
import numpy as np

# Example usage
alignments = np.random.rand(5, 5)
time = np.array([0, 1, 2, 3, 4])

time_align_correlate(alignments, time)

Use silhouette_samples to obtain the Silhouette Coefficient for each individual sample rather than the mean. This is useful for identifying specific samples that are poorly clustered.

Parameters:

• X: Array of shape [n_samples, n_features] or [n_samples, n_samples] (if metric='precomputed').
• labels: Array of label values for each sample.
• metric: String or callable for distance calculation.
• **kwds: Additional keyword arguments passed to the distance function.

Returns:

• silhouette: An array of shape [n_samples] containing the coefficient for each sample.

from bin.unsupervised import silhouette_samples
import numpy as np

X = np.array([[1, 1], [-1, -1], [1, -1], [-1, 1]])
labels = np.array([0, 0, 1, 1])

samples_scores = silhouette_samples(X, labels)
print(f'Scores per sample: {samples_scores}')