Assessing scRNA-seq integration with LISI and silhouette

Use local label diversity and cluster-distance summaries to examine two different questions: whether cells mix across technical batches, and whether annotated cell types remain locally distinct. These scores describe a chosen embedding and chosen labels; they do not prove that integration removed all technical effects or preserved every biological signal.

What you need before calculating scores

This workflow assumes that preprocessing and integration are complete and that you have a saved Seurat object or AnnData file. The object needs a cell-by-dimension embedding and cell-level metadata for the batch or sample label and, when available, a biological annotation. Choose the embedding used for downstream clustering or another scientifically justified high-dimensional representation, and record its name and dimensions. Compare integrated and unintegrated representations when the aim is to assess the effect of integration.

If you do not yet have such an object, first follow the official Seurat integration introduction or your AnnData-based integration workflow; this page does not download or construct biological inputs. Work from a copy and keep the original unchanged. Before scoring, confirm that the selected embedding is numeric, finite and has one row per cell, that embedding row names match metadata cell names (R), and that the requested batch, annotation and cluster columns are categorical and non-missing. The code below stops for several common schema errors; resolve them at source rather than dropping cells silently.

Use separate R and Python environments with the packages already used to create/read your objects. For R silhouette calculations, the recommended cluster package is required; install it if absent. For Python silhouettes, use pandas, numpy and scikit-learn. For R, install the package from its official repository if needed (install.packages("devtools"); devtools::install_github("immunogenomics/lisi")); Seurat must also be available to read and inspect the object. For Python silhouettes, install pandas, numpy and scikit-learn; an AnnData/Scanpy environment is needed only if you use it to export the shared CSV inputs. Check availability before running: Rscript -e 'stopifnot(requireNamespace("Seurat", quietly=TRUE), requireNamespace("lisi", quietly=TRUE), requireNamespace("cluster", quietly=TRUE)); sessionInfo()' and python -c 'import pandas, numpy, sklearn; print(pandas.__version__, numpy.__version__, sklearn.__version__)'. Save these outputs alongside results; package versions are environment-specific, not fixed versions claimed by this tutorial.

Save the first R code block as lisi_qc.R and run Rscript lisi_qc.R; save the Python block as silhouette_qc.py and run python silhouette_qc.py. Run the LISI block from the directory containing your Seurat RDS. Run either silhouette alternative from the directory containing both shared CSV files; once exported, these calculations do not read an RDS or H5AD. A successful run creates lisi_scores.csv for the R LISI workflow and silhouette_per_cell.csv for either silhouette alternative; the CSV silhouette routes leave their inputs unchanged. Inspect the printed summaries and confirm output rows match the number of embedding cells; stop if the embedding, labels or cell counts do not match. A structurally successful score calculation is not evidence that integration is biologically valid.

• iLISI: local diversity of a categorical batch label. Higher values indicate more batch categories represented in a cell's neighborhood, up to the number of observed categories.
• cLISI: local diversity of a categorical cell-type label. Lower values indicate greater local label purity; a value near 1 corresponds to a neighborhood dominated by one label. This is meaningful only to the extent that annotations are trustworthy.
• Silhouette: for a supplied set of cluster labels, compares each cell's average distance to its own cluster with its average distance to the nearest other cluster. It summarizes separation in the chosen representation, not biological correctness.

LISI and silhouette answer different questions. A metric cannot distinguish a technical effect from a real biological difference unless your design and labels make that distinction possible. For example, if each batch contains a different cell type, batch and biology are confounded; stronger batch mixing is not automatically a better result.

Calculate LISI in R

The lisi package's compute_lisi() takes a matrix with cells in rows and embedding dimensions in columns, plus metadata with one row per cell. It returns one score column for each requested metadata label. The code below uses existing Seurat reductions and reorders metadata to match embedding row names before calculation. Replace the example metadata keys with columns present in your object.

# Input: seurat_obj.rds, containing an existing PCA or integration reduction
# Outputs: lisi_scores.csv and lisi_metadata.csv; the input RDS is read only

library(lisi)
library(Seurat)

seurat_obj <- readRDS("seurat_obj.rds")
reduction_name <- if ("harmony" %in% Reductions(seurat_obj)) "harmony" else "pca"
if (!reduction_name %in% Reductions(seurat_obj)) {
  stop("No selected reduction is available; inspect Reductions(seurat_obj)")
}
emb <- Embeddings(seurat_obj, reduction_name)
meta <- seurat_obj@meta.data
if (!all(rownames(emb) %in% rownames(meta))) {
  stop("Some embedding cell names are missing from Seurat metadata")
}
meta <- meta[rownames(emb), , drop = FALSE]
labels_to_eval <- c("patient_id", "cell_type")  # replace with your metadata columns
if (!all(labels_to_eval %in% colnames(meta))) {
  stop("Update labels_to_eval: one or more metadata columns are missing")
}
if (anyNA(meta[, labels_to_eval, drop = FALSE])) {
  stop("Resolve missing values in the selected labels before calculating LISI")
}
lisi_res <- compute_lisi(emb, meta, labels_to_eval)
write.csv(lisi_res, file = "lisi_scores.csv", row.names = TRUE)
write.csv(meta, file = "lisi_metadata.csv")
summary(lisi_res)

Rows in the resulting CSV correspond to cells in the embedding. The selected batch-label column is used as iLISI and the cell-type column as cLISI; the software does not assign those biological meanings for you. The package uses a local-neighborhood calculation whose neighborhood scale affects the score. Check the installed package's documentation for its arguments, record the setting used, and compare a small, justified range if conclusions depend on it.

Summarize by cluster without losing cell alignment

Run this summary as a separate Rscript process after the LISI calculation block completes. The first block saves the matched metadata to lisi_metadata.csv; this block reloads it instead of relying on an in-memory object from another process. Keep both CSV files in the working directory. The cell-name check stops on mismatched inputs rather than silently assigning scores by row position.

# Inputs: lisi_scores.csv and lisi_metadata.csv, written by the first block
meta <- read.csv("lisi_metadata.csv", row.names = 1, check.names = FALSE)
lisi_df <- read.csv("lisi_scores.csv", row.names = 1, check.names = FALSE)
if (!setequal(rownames(meta), rownames(lisi_df))) {
  stop("LISI output cells do not match the metadata cells")
}
lisi_df <- lisi_df[rownames(meta), , drop = FALSE]
meta$ilisi <- lisi_df$patient_id
meta$clisi <- lisi_df$cell_type
if (!"seurat_clusters" %in% colnames(meta)) {
  stop("Add or select a cluster-label column before summarizing")
}
aggregate(cbind(ilisi, clisi) ~ seurat_clusters, data = meta, median)

Calculate silhouette values: equivalent R and Python alternatives

Unlike LISI above, silhouette compares supplied cluster labels against distances in an embedding. To make the language alternatives use the same data and labels, export a numeric cell-by-dimension matrix as silhouette_embedding.csv with a cell_id column, and one row per cell in silhouette_labels.csv with cell_id and cluster columns. Export the same reviewed embedding (for example, selected PCA or integration coordinates) and cluster labels from either Seurat or AnnData; cell IDs must match exactly. Both versions compute Euclidean silhouette scores, do not overwrite inputs, and require pairwise distances, which can be expensive in time and memory for large datasets. Use a justified subset or scalable alternative only if you document how it changes the calculation.

R: silhouette from the shared embedding and labels

# Inputs: silhouette_embedding.csv (cell_id, dim1, dim2, ...) and
#         silhouette_labels.csv (cell_id, cluster)
# Output: silhouette_per_cell.csv; inputs remain unchanged
library(cluster)

embedding_df <- read.csv("silhouette_embedding.csv", check.names = FALSE,
                         stringsAsFactors = FALSE)
labels_df <- read.csv("silhouette_labels.csv", check.names = FALSE,
                      stringsAsFactors = FALSE)
if (!"cell_id" %in% names(embedding_df) ||
    !all(c("cell_id", "cluster") %in% names(labels_df))) {
  stop("Expected cell_id in embedding and cell_id plus cluster in labels")
}
if (anyNA(embedding_df$cell_id) || anyNA(labels_df$cell_id) || any(!nzchar(embedding_df$cell_id)) || any(!nzchar(labels_df$cell_id))) {
  stop("Cell IDs must be present and non-empty")
}
if (anyDuplicated(embedding_df$cell_id) || anyDuplicated(labels_df$cell_id)) {
  stop("Cell IDs must be unique in both input files")
}
if (!setequal(embedding_df$cell_id, labels_df$cell_id)) {
  stop("Embedding and label inputs must contain the same cell IDs")
}
embedding_df <- embedding_df[match(labels_df$cell_id, embedding_df$cell_id), , drop = FALSE]
coords <- as.matrix(embedding_df[, setdiff(names(embedding_df), "cell_id"), drop = FALSE])
storage.mode(coords) <- "numeric"
if (nrow(coords) != nrow(labels_df) || ncol(coords) < 1L || any(!is.finite(coords))) {
  stop("Embedding must be a finite numeric matrix with one row per labeled cell")
}
cluster_label <- factor(labels_df$cluster)
if (anyNA(cluster_label) || any(!nzchar(as.character(cluster_label))) || nlevels(cluster_label) < 2L || nlevels(cluster_label) >= nrow(coords)) {
  stop("Silhouette requires non-missing labels for 2 to n_cells - 1 clusters")
}
distances <- stats::dist(coords, method = "euclidean")
sil <- cluster::silhouette(as.integer(cluster_label), distances)
result <- data.frame(cell_id = labels_df$cell_id,
                     cluster = as.character(cluster_label),
                     silhouette = sil[, "sil_width"])
write.csv(result, "silhouette_per_cell.csv", row.names = FALSE)
cat("Mean silhouette across cells:", mean(result$silhouette), "\n")

Python alternative: same embedding, labels, and distance

# Inputs: the same silhouette_embedding.csv and silhouette_labels.csv
# Output: silhouette_per_cell.csv; inputs remain unchanged
import numpy as np
import pandas as pd
from sklearn.metrics import silhouette_samples, silhouette_score

embedding_df = pd.read_csv("silhouette_embedding.csv")
labels_df = pd.read_csv("silhouette_labels.csv")
if "cell_id" not in embedding_df or not {"cell_id", "cluster"}.issubset(labels_df.columns):
    raise ValueError("Expected cell_id in embedding and cell_id plus cluster in labels")
if embedding_df["cell_id"].isna().any() or labels_df["cell_id"].isna().any():
    raise ValueError("Cell IDs must be present")
if embedding_df["cell_id"].duplicated().any() or labels_df["cell_id"].duplicated().any():
    raise ValueError("Cell IDs must be unique in both input files")
if set(embedding_df["cell_id"]) != set(labels_df["cell_id"]):
    raise ValueError("Embedding and label inputs must contain the same cell IDs")
embedding_df = embedding_df.set_index("cell_id").loc[labels_df["cell_id"]]
embedding = embedding_df.to_numpy(dtype=float)
labels = labels_df["cluster"].astype("string")
if embedding.ndim != 2 or embedding.shape[0] != len(labels) or embedding.shape[1] < 1:
    raise ValueError("Embedding must have one or more dimensions and one row per labeled cell")
if not np.isfinite(embedding).all() or labels.isna().any() or labels.eq("").any():
    raise ValueError("Embedding must be finite and cluster labels non-missing")
n_labels = labels.nunique()
if not 2 <= n_labels <= len(labels) - 1:
    raise ValueError("Silhouette requires 2 to n_cells - 1 distinct labels")
scores = silhouette_samples(embedding, labels.to_numpy(), metric="euclidean")
result = pd.DataFrame({"cell_id": labels_df["cell_id"],
                       "cluster": labels.astype(str), "silhouette": scores})
result.to_csv("silhouette_per_cell.csv", index=False)
print("Mean silhouette across cells:", silhouette_score(embedding, labels.to_numpy(), metric="euclidean"))

Per-cell scores lie between -1 and 1. Positive values mean a cell is closer, on average, to its assigned cluster than to the nearest alternative; values near zero indicate overlap; negative values flag cells whose distances favor another cluster. Treat these as prompts for inspection, not automatic misassignment calls. Silhouette depends on the chosen distance, embedding dimensions and cluster labels, and can favor compact, well-separated clusters even when that geometry is not the biological question.

Read the results as a comparison

Compare the same cells and label definitions across the unintegrated and integrated embeddings where possible. Report medians and interquartile ranges for per-cell LISI values, and cluster-wise summaries for silhouette. Keep the embedding name and dimensions, batch and cell-type columns, cluster definition, package versions and LISI neighborhood setting with the result.

• Higher batch-label LISI alongside stable cell-type-label LISI may be consistent with improved local batch mixing while retaining annotation purity. It does not demonstrate that all technical effects are gone.
• Higher cell-type-label LISI means neighborhoods contain more annotation labels. Inspect annotation quality and marker expression before attributing this pattern to overcorrection.
• Low or negative silhouette values identify cells or clusters with weak separation under the chosen distance and representation. Review cluster sizes, labels, marker expression and the underlying embedding before changing the clustering.

Do not interpret a single score as a pass/fail threshold. LISI is conditional on supplied categorical labels and neighborhood scale; silhouette is conditional on the supplied clusters and distance. If batches are confounded with cell identity, neither score alone can resolve that design limitation. These diagnostics do not establish causality or replace biological validation.

Troubleshoot common failures

• LISI reports an error or yields missing results: check that embedding rows and metadata rows refer to the same cells in the same order, the requested metadata fields exist, and labels are not missing. A missing reduction name or absent metadata field should be fixed explicitly rather than replaced with a different representation without recording the change.
• Batch LISI remains low: confirm that the selected embedding is the integrated representation you intended, inspect batch levels and sample sizes, then compare with the unintegrated embedding. Low mixing can also be expected when batches contain genuinely different cell populations.
• Cell-type LISI increases: check annotation granularity and uncertainty, and review marker expression and the study design. A score change alone does not establish loss of biology.
• Silhouette calculation fails: confirm labels are present and that the number of distinct labels is between 2 and the number of cells minus 1. Reduce the dataset only with a documented sampling strategy if full pairwise calculation is too costly; a subset estimate is not the exact full-data result.
• LISI and silhouette appear to disagree: they assess label diversity and cluster separation, respectively. Check whether they use compatible embeddings and whether cluster labels and biological annotations answer the same question before drawing a conclusion.

References

• LISI package repository describes the cell-by-coordinate matrix, metadata input, and effective number of local categories.
• scikit-learn silhouette coefficient guide documents the distance-based definition and valid label-count condition.
• AnnData write_h5ad reference documents writing the scored object to a distinct file.
• Seurat integration introduction explains integration reductions and evaluating shared cell states in context.
• Korsunsky et al. (2019), Harmony describes the integration method cited in this tutorial.