单细胞分析流程：Seurat v5 与 Scanpy 从质控到细胞注释

# Seurat v5 / SeuratObject 5.4: replacements for v4 code
GetAssayData(obj, slot = "counts")   # Error: The `slot` argument of `GetAssayData()` was deprecated in SeuratObject 5.0.0 and is now defunct.
obj@assays$RNA@counts                # Error: no slot of name "counts" for this object of class "Assay5"

obj[["RNA"]]$counts                  # v5 accessors
LayerData(obj, assay = "RNA", layer = "data")
GetAssayData(obj, layer = "data")    # works only while the assay has one "data" layer

# Old v4 object read into v5: convert the assay class when a tool needs Assay5 (e.g. BPCells, IntegrateLayers)
obj[["RNA"]] <- as(obj[["RNA"]], "Assay5")
# Tool that only accepts v4 assays: convert back (merge layers first)
obj <- JoinLayers(obj); obj[["RNA"]] <- as(obj[["RNA"]], "Assay")

# Scanpy 1.12 (Python >= 3.12): PBMC 3k from Cell Ranger output to annotated clusters
import numpy as np, scanpy as sc
from scipy.stats import median_abs_deviation

adata = sc.read_10x_mtx("filtered_gene_bc_matrices/hg19", var_names="gene_symbols")
# Cell Ranger >= 3: sc.read_10x_h5("filtered_feature_bc_matrix.h5")
adata.var_names_make_unique()

# QC metrics: mitochondrial / ribosomal / hemoglobin genes (mouse: "mt-", ("Rps","Rpl"), "^Hb[^(p)]")
adata.var["mt"] = adata.var_names.str.startswith("MT-")
adata.var["ribo"] = adata.var_names.str.startswith(("RPS", "RPL"))
adata.var["hb"] = adata.var_names.str.contains("^HB[^(P)]")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt", "ribo", "hb"], percent_top=[20], log1p=True, inplace=True)

# Doublets on raw counts, per capture (batch_key="sample" when several samples are concatenated)
sc.pp.scrublet(adata, random_state=0)

# MAD outliers: 5 MAD on counts / genes / top-20 share, mito 3 MAD AND above a fixed floor
def outlier(x, nmads, upper_only=False):
    med, mad = np.median(x), median_abs_deviation(x)
    return (x > med + nmads * mad) if upper_only else ((x < med - nmads * mad) | (x > med + nmads * mad))
o = adata.obs
bad = (outlier(o.log1p_total_counts, 5) | outlier(o.log1p_n_genes_by_counts, 5)
       | outlier(o.pct_counts_in_top_20_genes, 5)
       | (outlier(o.pct_counts_mt, 3, upper_only=True) & (o.pct_counts_mt > 8)))
adata = adata[~bad & ~adata.obs.predicted_doublet].copy()
sc.pp.filter_genes(adata, min_cells=3)

# Normalize, HVG, PCA
adata.layers["counts"] = adata.X.copy()
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat_v3", layer="counts")  # add batch_key="sample" for multi-sample data
adata.raw = adata                                     # keep all genes (log-normalized) for markers and plots
hv = adata[:, adata.var.highly_variable].copy()       # scale only HVGs; scaling all genes makes X dense
sc.pp.scale(hv, max_value=10)
sc.tl.pca(hv, n_comps=50)

# Graph, clustering, UMAP
sc.pp.neighbors(hv, n_neighbors=15, n_pcs=30)
for r in (0.2, 0.4, 0.6, 0.8, 1.0):
    sc.tl.leiden(hv, resolution=r, flavor="igraph", n_iterations=2, key_added=f"leiden_{r}")
sc.tl.umap(hv)

# Markers on log-normalized expression of all genes
adata.obs["leiden"] = hv.obs["leiden_0.6"].values
sc.tl.rank_genes_groups(adata, "leiden", method="wilcoxon", use_raw=False)
sc.pl.rank_genes_groups_dotplot(adata, n_genes=5)

# Seurat v5: same workflow in R
library(Seurat)
counts <- Read10X("filtered_gene_bc_matrices/hg19")      # Cell Ranger >= 3: Read10X_h5("filtered_feature_bc_matrix.h5")
obj <- CreateSeuratObject(counts, project = "pbmc3k", min.cells = 3, min.features = 200)
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")   # mouse: "^mt-"

# Doublets: scDblFinder on the raw counts of each capture, before strict QC
library(scDblFinder); library(SingleCellExperiment)
set.seed(1)
sce <- scDblFinder(SingleCellExperiment(list(counts = LayerData(obj, layer = "counts"))))
obj$dbl_class <- sce$scDblFinder.class               # several captures: scDblFinder(sce, samples = "sample")

# MAD thresholds instead of the 200 / 2500 / 5% values of the PBMC3k tutorial
mad_out <- function(x, n, upper = FALSE) {
  m <- median(x); d <- mad(x)
  if (upper) x > m + n * d else x < m - n * d | x > m + n * d
}
bad <- mad_out(log1p(obj$nCount_RNA), 5) | mad_out(log1p(obj$nFeature_RNA), 5) |
       (mad_out(obj$percent.mt, 3, upper = TRUE) & obj$percent.mt > 8)
obj <- subset(obj, cells = colnames(obj)[!bad & obj$dbl_class == "singlet"])

obj <- NormalizeData(obj)                            # LogNormalize, scale.factor = 1e4
obj <- FindVariableFeatures(obj, nfeatures = 2000)
obj <- ScaleData(obj)                                # HVGs only by default
obj <- RunPCA(obj, npcs = 50)
ElbowPlot(obj, ndims = 50)
obj <- FindNeighbors(obj, dims = 1:20)
obj <- FindClusters(obj, resolution = c(0.2, 0.4, 0.6, 0.8, 1.0))
obj <- RunUMAP(obj, dims = 1:20)

library(clustree)                                    # attach it; clustree::clustree() alone fails at ggsave
clustree(obj, prefix = "RNA_snn_res.")
Idents(obj) <- "RNA_snn_res.0.6"
markers <- FindAllMarkers(obj, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

# Seurat v5 multi-sample integration: split -> per-layer preprocessing -> IntegrateLayers -> JoinLayers
obj <- merge(s1, y = list(s2, s3), add.cell.ids = c("s1", "s2", "s3"))   # already one layer per object: counts.1, counts.2, counts.3
# Object built from one matrix, or after JoinLayers: split by sample first
# obj[["RNA"]] <- split(obj[["RNA"]], f = obj$sample)
# Splitting a merged object again fails: The following layers are already split: 'counts.1', 'counts.2' Please join before splitting
obj <- NormalizeData(obj); obj <- FindVariableFeatures(obj); obj <- ScaleData(obj); obj <- RunPCA(obj)

# Unintegrated baseline first
obj <- FindNeighbors(obj, dims = 1:30, reduction = "pca")
obj <- FindClusters(obj, resolution = 0.5, cluster.name = "unintegrated_clusters")
obj <- RunUMAP(obj, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")

obj <- IntegrateLayers(obj, method = HarmonyIntegration, orig.reduction = "pca", new.reduction = "harmony")
# alternatives: method = RPCAIntegration (more conservative), CCAIntegration, FastMNNIntegration, scVIIntegration
obj <- FindNeighbors(obj, reduction = "harmony", dims = 1:30)
obj <- FindClusters(obj, resolution = 0.5, cluster.name = "harmony_clusters")
obj <- RunUMAP(obj, reduction = "harmony", dims = 1:30, reduction.name = "umap.harmony")

obj <- JoinLayers(obj)          # before FindMarkers / FindAllMarkers / GetAssayData

# Re-clustering a subset: layers are joined now, so split again before IntegrateLayers
sub <- subset(obj, idents = c("0", "3"))
sub[["RNA"]] <- split(sub[["RNA"]], f = sub$sample)
sub <- NormalizeData(sub); sub <- FindVariableFeatures(sub); sub <- ScaleData(sub); sub <- RunPCA(sub)
sub <- IntegrateLayers(sub, method = HarmonyIntegration, orig.reduction = "pca", new.reduction = "harmony_sub")
sub <- FindNeighbors(sub, reduction = "harmony_sub", dims = 1:20); sub <- FindClusters(sub, resolution = 0.3)
sub <- JoinLayers(sub)

# Scanpy: Harmony (harmonypy) and scVI on the same AnnData
import scanpy as sc
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat_v3", layer="counts", batch_key="sample")
hv = adata[:, adata.var.highly_variable].copy()
sc.pp.scale(hv, max_value=10); sc.tl.pca(hv, n_comps=50)
# sc.external.pp.harmony_integrate(hv, key="sample") fails with harmonypy 2.x:
#   ValueError: Value passed for key 'X_pca_harmony' is of incorrect shape ... (50,) while it should have had (n_cells,)
# harmonypy 2.x already returns Z_corr as cells x PCs, the Scanpy 1.12.4 wrapper transposes it again
import harmonypy
ho = harmonypy.run_harmony(hv.obsm["X_pca"], hv.obs, "sample")
Z = ho.Z_corr
hv.obsm["X_pca_harmony"] = Z if Z.shape[0] == hv.n_obs else Z.T
sc.pp.neighbors(hv, use_rep="X_pca_harmony", n_pcs=30)
sc.tl.leiden(hv, resolution=0.5, flavor="igraph", n_iterations=2)

# scVI: needs raw counts; a GPU shortens training, CPU works for tens of thousands of cells
import scvi
scvi.model.SCVI.setup_anndata(hv, layer="counts", batch_key="sample")
m = scvi.model.SCVI(hv, n_latent=30)
m.train()                                                      # max_epochs scales down automatically for large data
hv.obsm["X_scVI"] = m.get_latent_representation()
sc.pp.neighbors(hv, use_rep="X_scVI")

# CellTypist on log1p data normalized to 10,000 counts per cell (adata after normalize_total + log1p)
import celltypist
from celltypist import models
models.download_models(model=["Immune_All_Low.pkl"])
pred = celltypist.annotate(adata, model="Immune_All_Low.pkl",
                           majority_voting=True, over_clustering=adata.obs["leiden"])
adata.obs["celltypist"] = pred.predicted_labels["majority_voting"]
adata.obs["celltypist_conf"] = pred.probability_matrix.max(axis=1)   # low values = label not in the model

# SingleR + celldex (Bioconductor). On macOS arm64 bioconda has no SingleR build:
# install with BiocManager::install(c("SingleR", "celldex")) inside the R environment
library(SingleR); library(celldex)
ref <- celldex::MonacoImmuneData()                    # immune; HumanPrimaryCellAtlasData() for broad tissue types
pred <- SingleR(test = LayerData(obj, layer = "data"), ref = ref,
                labels = ref$label.fine, clusters = obj$seurat_clusters)   # per-cluster mode
obj$singler <- pred$pruned.labels[match(obj$seurat_clusters, rownames(pred))]
plotScoreHeatmap(pred)                                # check that the best score clearly beats the second

# R: anndataR 1.2 (Bioconductor). rhdf5 is a separate install, otherwise:
#   Error: HDF5AnnData requires the rhdf5 package
# BiocManager::install(c("anndataR", "rhdf5"))
library(Seurat); library(anndataR)

obj <- read_h5ad("pbmc3k_processed.h5ad", as = "Seurat")   # h5ad -> Seurat: X and layers become layers, obsm -> reductions
obj <- JoinLayers(obj)                                       # split layers are skipped with a warning only
write_h5ad(obj, "from_seurat.h5ad")                          # Seurat -> h5ad: layers "counts" and "data", X left empty

# Python: X is empty after the Seurat -> h5ad direction, point it at a layer
import anndata as ad
adata = ad.read_h5ad("from_seurat.h5ad")
adata.X = adata.layers["counts"].copy()                      # or layers["data"] for log-normalized values

# Seurat v5 + BPCells: counts stay on disk, analyze a 50,000-cell sketch, project back
library(Seurat); library(BPCells)
options(future.globals.maxSize = 4e9)
mat <- open_matrix_10x_hdf5("filtered_feature_bc_matrix.h5")
write_matrix_dir(mat, dir = "bpcells_counts")              # one-time conversion to the on-disk format
obj <- CreateSeuratObject(open_matrix_dir("bpcells_counts"))
obj <- NormalizeData(obj); obj <- FindVariableFeatures(obj)
obj <- SketchData(obj, ncells = 50000, method = "LeverageScore", sketched.assay = "sketch")
DefaultAssay(obj) <- "sketch"
obj <- FindVariableFeatures(obj); obj <- ScaleData(obj); obj <- RunPCA(obj)
obj <- FindNeighbors(obj, dims = 1:50); obj <- FindClusters(obj, resolution = 1)
obj <- RunUMAP(obj, dims = 1:50, return.model = TRUE)
obj <- ProjectData(obj, assay = "RNA", full.reduction = "pca.full", sketched.assay = "sketch",
                   sketched.reduction = "pca", umap.model = "umap", dims = 1:50,
                   refdata = list(cluster_full = "seurat_clusters"))

# Python side: open .h5ad without loading X (anndata >= 0.12)
import anndata as ad
a = ad.experimental.read_lazy("big.h5ad")   # lazy obs/var/X; subset first, then .to_memory()
a = ad.read_h5ad("big.h5ad", backed="r")    # older backed mode: read-only X on disk