scanpy in Single-Cell-Projekten: QC, Clustering und robuste Markerlogikscanpy in Single-cell Projects: QC, Clustering, and Robust Marker Logic

Abstract

scanpy ist das Python-Referenzpaket für die Analyse von Single-Cell-RNA-Seq-Daten (scRNA-Seq). Es deckt die gesamte Pipeline ab – von der Qualitätskontrolle und Normalisierung über Dimensionsreduktion (PCA, UMAP, t-SNE) und Clustering (Leiden, Louvain) bis zur Identifikation von Markergenen und Zelltyp-Annotation. Auf Basis des AnnData-Objekts als zentraler Datenstruktur integriert scanpy sich nahtlos mit scvi-tools, CellRank und dem gesamten scverse-Ökosystem. Für Bioinformatik-Teams, die große Single-Cell-Datensätze (10&sup5;–10&sup7; Zellen) effizient verarbeiten müssen, ist scanpy die erste Wahl im Python-Stack.

Typisches Projektszenario

Ein Immunologie-Labor untersucht die Heterogenität von CD8⁺ T-Zellen in Tumorgewebe (TILs) aus 8 Melanom-Patienten. Pro Patient werden ~5.000 Zellen per 10x Genomics Chromium erfasst (40.000 Zellen total). Die CellRanger-Ausgabe (gefilterte Feature-Barcode-Matrix) wird in scanpy geladen. Ziel: Subpopulationen identifizieren (naïv, effector, exhausted, progenitor exhausted), deren Marker definieren und die klonale Expansion über TCR-Daten mit dem Erschöpfungsphänotyp verbinden.

Welches Problem löst scanpy?

• Sparse, hochdimensionale Daten: Eine typische scRNA-Seq-Matrix hat 20.000 Gene × 40.000 Zellen, wobei >90% der Einträge Null sind (Dropout). scanpy arbeitet nativ mit Sparse-Matrizen (scipy.sparse) und skaliert effizient auf Millionen von Zellen.
• Batch-Effekte zwischen Patienten: Technische und biologische Variabilität überlagern sich. scanpy integriert mit Harmony, scVI, bbknn und anderen Batch-Korrektur-Methoden, um patientenspezifische Effekte zu entfernen, während biologische Variation erhalten bleibt.
• Zelltyp-Identifikation: Die unbekannte Zellzusammensetzung muss aus den Daten selbst erschlossen werden – über Clustering, Markergene und Vergleich mit Referenzatlanten.

Warum Teams scanpy einsetzen

1. AnnData-Objekt: Das zentrale Datenobjekt speichert Expressionsmatrix, Zell-Metadaten, Gen-Annotationen, Dimensionsreduktionen und Graphen in einer einzigen Struktur – keine separaten Dateien oder Konvertierungen nötig.
2. GPU-Beschleunigung: Über rapids-singlecell können rechenintensive Schritte (PCA, Neighbor-Graphen, UMAP) auf NVIDIA-GPUs ausgeführt werden – 10–100× Speedup.
3. Modularer Aufbau: Jeder Analyseschritt ist eine eigenständige Funktion (sc.pp.*, sc.tl.*, sc.pl.*), die das AnnData-Objekt in-place modifiziert – das macht Pipelines reproduzierbar und leicht anpassbar.
4. scverse-Integration: squidpy (spatial), scvi-tools (probabilistische Modelle), CellRank (Trajektorien), muon (Multiomics) – alle bauen auf AnnData/scanpy auf.
5. Community: >1.500 GitHub-Stars, aktive Entwicklung, umfangreiche Tutorials und ein wachsendes Plugin-Ökosystem.

scanpy-Pipeline im Detail

Die Qualitätskontrolle bildet das Fundament jeder scRNA-Seq-Analyse. scanpy berechnet per-Zell-Metriken (Gene pro Zelle, Counts pro Zelle, Anteil mitochondrialer Reads) über sc.pp.calculate_qc_metrics(). Typische Filter: Zellen mit <200 Genen oder >20% mitochondrialen Reads werden entfernt – erstere sind leere Droplets oder Debris, letztere sterbende Zellen mit Zytoplasma-Verlust. Gene, die in weniger als 3 Zellen exprimiert sind, werden ebenfalls gefiltert.

Die Normalisierung mit sc.pp.normalize_total() (Library-Size-Normalisierung zu 10.000 Counts pro Zelle) gefolgt von sc.pp.log1p() (Log-Transformation) erzeugt die Grundlage für alle nachfolgenden Schritte. Alternativ bietet scvi-tools eine probabilistische Normalisierung, die biologische und technische Varianz trennt.

Die Feature-Selektion mit sc.pp.highly_variable_genes() (typisch: 2.000–4.000 HVGs) reduziert das Rauschen und fokussiert die Analyse auf biologisch informative Gene. Der Seurat-v3-Flavor (flavor='seurat_v3') ist der robusteste Ansatz für UMI-basierte Daten.

Python-Code: scanpy-Pipeline

import scanpy as sc
import numpy as np

sc.settings.verbosity = 2
sc.settings.figdir = './figures/'

# Daten laden (10x CellRanger Output)
adata = sc.read_10x_mtx('./filtered_feature_bc_matrix/',
                         var_names='gene_symbols', cache=True)
adata.var_names_make_unique()

# QC-Metriken + Filtern
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'],
                            percent_top=None, inplace=True)
adata = adata[adata.obs.n_genes_by_counts.between(200, 5000) &
              (adata.obs.pct_counts_mt < 20)].copy()

# Normalisierung + Log-Transformation
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# Highly Variable Genes
sc.pp.highly_variable_genes(adata, n_top_genes=3000,
                             flavor='seurat_v3')
adata.raw = adata
adata = adata[:, adata.var.highly_variable]

# PCA + Nachbar-Graph + Clustering
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, n_comps=50, svd_solver='arpack')
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)
sc.tl.leiden(adata, resolution=0.8)

# UMAP + Markergene
sc.tl.umap(adata, min_dist=0.3)
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')

# Visualisierung
sc.pl.umap(adata, color=['leiden', 'CD8A', 'PDCD1', 'TCF7'],
           save='_til_markers.png')
sc.pl.rank_genes_groups_dotplot(adata, n_genes=5,
                                 save='_markers_dotplot.png')

Beispielausgabe

# Loaded 42,318 cells x 33,694 genes
# Filtered to 38,456 cells (QC)
# 3,000 highly variable genes selected
# PCA: 50 components, 30 used for neighbors
# Leiden clustering: 14 clusters at resolution 0.8

# Top markers per cluster (Wilcoxon):
# Cluster 0: TCF7, LEF1, CCR7, SELL     (naive)
# Cluster 1: GZMB, PRF1, GNLY, NKG7     (effector)
# Cluster 3: PDCD1, HAVCR2, LAG3, TOX   (exhausted)
# Cluster 5: TCF7, PDCD1, SLAMF6        (prog. exhausted)

Diagnostische Plots

Vergleich mit Alternativen

Merkmal	scanpy	Seurat (R)	scvi-tools
Sprache	Python	R	Python (PyTorch)
Datenobjekt	AnnData	SeuratObject	AnnData
Max. Zellzahl	>107	~106	>107
GPU-Support	Ja (rapids)	Nein (nativ)	Ja (PyTorch)
Batch-Korrektur	bbknn, Harmony	CCA, RPCA, Harmony	scVI (probabilistisch)
Normalisierung	Shift-Log, scran	LogNormalize, SCTransform	scVI (latent)
Clustering	Leiden, Louvain	Louvain, Leiden	Über scanpy

Statistische Vertiefung

Die Leiden-Algorithmus-Implementierung in scanpy nutzt den modularity-basierten Community-Detection-Ansatz auf dem k-Nearest-Neighbor-Graphen. Der Resolution-Parameter γ steuert die Granularität: γ=0.5 liefert wenige große Cluster, γ=2.0 viele kleine. Die Wahl von γ ist kritisch und sollte über Silhouette-Scores oder biologische Marker validiert werden – nicht blind optimiert. Der Leiden-Algorithmus garantiert im Gegensatz zu Louvain, dass alle resultierenden Communities zusammenhängend sind.

Die UMAP-Embeddings sind Visualisierungen, keine quantitativen Analysetools. Absolute Distanzen zwischen Clustern in der UMAP sind nicht interpretierbar – nur die relative Nachbarschaft ist informativ. Zwei Cluster, die im UMAP weit auseinander liegen, können im hochdimensionalen Raum durchaus ähnlich sein. Alle quantitativen Analysen (DE-Tests, Trajektorie-Inferenz) sollten auf den PCA-Daten oder dem Nachbar-Graphen basieren, nicht auf UMAP-Koordinaten.

Der Wilcoxon-Rangsummentest in rank_genes_groups() testet für jedes Gen, ob es in einem Cluster höher exprimiert ist als in allen anderen Zellen. Der Test ist nicht-parametrisch und robust gegen die Verletzung von Normalverteilungsannahmen in scRNA-Seq-Daten. Die Benjamini-Hochberg-Korrektur kontrolliert die False Discovery Rate. Alternative: der t-Test mit Overestimation-Correction (method='t-test_overestim_var') ist konservativer und vermeidet falsch-positive Marker bei großen Zellzahlen.

Zitationen

1. Wolf FA, Angerer P, Theis FJ (2018). “SCANPY: large-scale single-cell gene expression data analysis.” Genome Biology, 19, 15.
2. Virshup I et al. (2023). “The scverse project provides a computational framework for single-cell omics data analysis.” Nature Biotechnology, 41, 604–606.
3. Traag VA, Waltman L, van Eck NJ (2019). “From Louvain to Leiden: guaranteeing well-connected communities.” Scientific Reports, 9, 5233.

Fazit

scanpy ist das Arbeitspferd der Single-Cell-Analyse im Python-Ökosystem. Es vereint Geschwindigkeit (Sparse-Matrizen, GPU-Support), Modularität (scverse-Plugins) und eine aktive Community. Für Teams, die bereits in Python arbeiten, ist scanpy die natürliche Wahl gegenüber Seurat. Limitierungen: (1) Die Qualität der Ergebnisse hängt stark von der Parameterwahl ab (n_neighbors, resolution, n_pcs). (2) Die Log-Normalisierung ist statistisch suboptimal – probabilistische Methoden (scVI) liefern bessere latente Repräsentationen. (3) Zelltyp-Annotation erfordert nach wie vor Expertenwissen und manuelle Kuration – automatisierte Tools (CellTypist, scArches) können unterstützen, aber nicht ersetzen.

Dokumentation

• scanpy Dokumentation
• scanpy Tutorials

import scanpy as sc
import numpy as np

sc.settings.verbosity = 2
sc.settings.figdir = './figures/'

# Load data (10x CellRanger output)
adata = sc.read_10x_mtx('./filtered_feature_bc_matrix/',
                         var_names='gene_symbols', cache=True)
adata.var_names_make_unique()

# QC metrics + filtering
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'],
                            percent_top=None, inplace=True)
adata = adata[adata.obs.n_genes_by_counts.between(200, 5000) &
              (adata.obs.pct_counts_mt < 20)].copy()

# Normalization + log transformation
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# Highly variable genes
sc.pp.highly_variable_genes(adata, n_top_genes=3000,
                             flavor='seurat_v3')
adata.raw = adata
adata = adata[:, adata.var.highly_variable]

# PCA + neighbor graph + clustering
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, n_comps=50, svd_solver='arpack')
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)
sc.tl.leiden(adata, resolution=0.8)

# UMAP + marker genes
sc.tl.umap(adata, min_dist=0.3)
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')

# Visualization
sc.pl.umap(adata, color=['leiden', 'CD8A', 'PDCD1', 'TCF7'],
           save='_til_markers.png')
sc.pl.rank_genes_groups_dotplot(adata, n_genes=5,
                                 save='_markers_dotplot.png')

# Loaded 42,318 cells x 33,694 genes
# Filtered to 38,456 cells (QC)
# 3,000 highly variable genes selected
# PCA: 50 components, 30 used for neighbors
# Leiden clustering: 14 clusters at resolution 0.8

# Top markers per cluster (Wilcoxon):
# Cluster 0: TCF7, LEF1, CCR7, SELL     (naive)
# Cluster 1: GZMB, PRF1, GNLY, NKG7     (effector)
# Cluster 3: PDCD1, HAVCR2, LAG3, TOX   (exhausted)
# Cluster 5: TCF7, PDCD1, SLAMF6        (prog. exhausted)