标准单细胞分析流程：scRNA-seq分步指南（10x实战）｜生信云服务

本部分是本指南的实践核心，我们将以最经典的10x Genomics 3k PBMC（外周血单个核细胞）数据集为例，一步步拆解标准的下游分析流程。这个数据集也是Seurat和Scanpy官方教程中使用的范例，非常适合初学者上手实践 。

3.1 上游分析：从原始序列到表达矩阵

在开始下游分析之前，测序仪产出的原始FASTQ文件需要经过上游流程的处理。对于10x Genomics平台的数据，最常用的工具是其官方软件Cell Ranger。Cell Ranger会自动完成序列比对、细胞条形码的识别与过滤、UMI的定量等一系列复杂步骤，最终生成我们下游分析所需的输入文件 。

通常，Cell Ranger的输出是一个文件夹，其中包含三个核心文件：matrix.mtx（记录了非零表达值的稀疏矩阵）、barcodes.tsv（细胞条形码列表）和features.tsv或genes.tsv（基因列表） 。这三者共同构成了我们的起点——基因-细胞表达矩阵。

3.2 下游分析：两大工具包实战演练 (Seurat & Scanpy)

对于下游分析的每一步，我们将阐明其目标 (Goal)、方法 (Method)，并提供在Seurat和Scanpy中的核心代码 (Code)。

1. 加载数据与创建对象

• 目标：将上游分析生成的表达矩阵读入内存，并创建一个专门用于存储和管理单细胞数据的对象。

• 方法：两大工具包都使用一个核心对象来统一管理数据。这个对象不仅存储了原始表达矩阵，还会陆续存入标准化后的数据、细胞和基因的元数据（metadata）、PCA降维结果、聚类信息等。

Seurat (R)：使用Seurat对象。通过Read10X()函数读取Cell Ranger的输出，然后用CreateSeuratObject()创建核心对象 。

• R

# library(Seurat)
# library(dplyr)
pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)

Scanpy (Python)：使用AnnData对象。通过sc.read_10x_mtx()函数直接读取数据并创建对象 。

• Python

# import scanpy as sc
adata = sc.read_10x_mtx(
    'filtered_gene_bc_matrices/hg19/',
    var_names='gene_symbols',
    cache=True)

2. 质量控制 (QC) 与细胞过滤

目标：识别并移除低质量的细胞（如死细胞、空油滴）和潜在的“双包”，确保后续分析的准确性。

方法：计算每个细胞的QC指标，并根据设定的阈值进行过滤。三个核心QC指标是：

1. 检测到的基因数 (nFeature_RNA / n_genes_by_counts)：过低可能为空油滴，过高可能为双包。

2. 总UMI数 (nCount_RNA / total_counts)：与基因数高度相关，反映测序深度。

3. 线粒体基因表达比例：比例过高通常意味着细胞状态不佳或已死亡，因为细胞质mRNA降解后，线粒体内的mRNA会相对富集 。

Seurat (R)：使用PercentageFeatureSet()计算线粒体基因比例。使用subset()函数进行过滤。通过VlnPlot可视化QC指标分布 。

• R

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
# VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

Scanpy (Python)：使用sc.pp.calculate_qc_metrics()计算QC指标。使用pandas的索引方式进行过滤。通过sc.pl.violin可视化 。

• Python

adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
# sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'], jitter=0.4, multi_panel=True)
adata = adata[adata.obs.n_genes_by_counts < 2500, :]
adata = adata[adata.obs.pct_counts_mt < 5, :]

3. 数据标准化

• 目标：校正细胞间因测序深度不同而引入的技术偏差。

• 方法：最常用的方法是“LogNormalize”。具体操作为：将每个细胞中每个基因的UMI计数除以该细胞的总UMI数，再乘以一个固定的缩放因子（通常是10,000），最后对结果进行自然对数转换（log(1+x)） 。

Seurat (R)：一步完成。

• R

pbmc <- NormalizeData(pbmc)

Scanpy (Python)：分两步完成。

• Python

sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

4. 特征选择 (鉴定高变基因)

• 目标：从数万个基因中，筛选出那些在不同细胞间表达量差异最大的基因（Highly Variable Genes, HVGs）。这些基因被认为携带了最主要的生物学信号，后续的降维和聚类将主要基于这些基因进行，以减少噪音并提高计算效率。

• 方法：计算每个基因的表达均值和离散度（方差），并识别出那些相对于其表达水平均值而言，离散度异常高的基因。

Seurat (R)：

• R

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

Scanpy (Python)：

• Python

sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata = adata[:, adata.var.highly_variable] # 过滤数据，只保留高变基因

5. 数据缩放

• 目标：对高变基因的表达值进行线性变换，使得每个基因在所有细胞中的表达均值为0，方差为1。这一步的目的是为了在后续的PCA分析中，所有基因的权重都相同，避免那些高表达基因（即使其变异不大）不成比例地主导分析结果。

• 方法：对每个高变基因的表达值应用标准的Z-score变换。

Seurat (R)：ScaleData()函数还可以同时“回归掉”不想要的变异来源，如线粒体比例或细胞周期评分 。

• R

all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

Scanpy (Python)：

• Python

sc.pp.scale(adata, max_value=10)

6. 线性降维 (PCA)

• 目标：将高维的高变基因空间（约2000维）压缩到低维的主成分空间（通常10-50维），捕获数据中的主要变异模式。

• 方法：在缩放后的数据上运行主成分分析。一个关键的决策是选择使用多少个主成分（PCs）进行下游分析。通常使用“肘部图”（Elbow Plot）来辅助判断，寻找p值显著下降趋势变缓的“拐点” 。

Seurat (R)：

• R

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
# ElbowPlot(pbmc)

Scanpy (Python)：

• Python

sc.tl.pca(adata, svd_solver='arpack')
# sc.pl.pca_variance_ratio(adata, log=True)

7. 细胞聚类

• 目标：基于细胞在PCA空间中的相似性，将它们划分成不同的群组（clusters）。

• 方法：这是一个两步过程。首先，在选定的PC空间中构建一个k-近邻图（k-Nearest Neighbor graph），其中每个细胞是一个节点，与它最相似的k个细胞之间有边相连。然后，在这个图上运行社区发现算法（如Louvain或Leiden算法），将图分割成密集连接的社区，每个社区即为一个细胞簇。

Seurat (R)：

• R

pbmc <- FindNeighbors(pbmc, dims = 1:10) # 使用前10个PCs
pbmc <- FindClusters(pbmc, resolution = 0.5) # resolution参数控制聚类粒度

Scanpy (Python)：Leiden算法是Louvain算法的改进版，通常效果更好 。

• Python

sc.pp.neighbors(adata, n_neighbors=10, n_pcs=10) # 使用前10个PCs
sc.tl.leiden(adata)

8. 非线性降维 (可视化)

• 目标：将高维的PCA空间进一步降维到二维，生成直观的可视化图像，以便观察细胞的分布和聚类结果。

• 方法：在选定的PC空间上运行UMAP（Uniform Manifold Approximation and Projection）或t-SNE算法。UMAP目前更受青睐，因为它通常能更好地保留数据的全局结构。

Seurat (R)：

• R

pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap", label = TRUE)

Scanpy (Python)：

• Python

sc.tl.umap(adata)
sc.pl.umap(adata, color=['leiden'])

9. 寻找差异表达基因 (聚类标记物)

• 目标：为每个细胞簇找到特异性高表达的基因（marker genes）。这些基因是后续进行细胞类型注释的关键生物学证据。

• 方法：对每个细胞簇，将其与所有其他细胞簇进行比较，执行差异表达分析，找出在该簇中表达水平显著上调的基因。

Seurat (R)：

• R

# 寻找所有簇的marker genes
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
# pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)

Scanpy (Python)：

• Python

sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')
# sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

表2：核心分析流程：Seurat vs. Scanpy 命令参考

这张“罗塞塔石碑”式的表格，将帮助你在两个生态系统之间自由切换，理解它们在实现相同科学逻辑时所使用的不同语法。