Seurat v5 完整教學：單細胞 RNA-seq 分析從零到進階

1. 專案定位與生態系統

1.1 什麼是 Seurat？

Seurat 是一套以 R 語言開發的 single-cell genomics (單細胞基因體學) 分析工具套件，由 New York Genome Center 的 Satija Lab 維護。它提供了從原始 count matrix (計數矩陣) 到生物學詮釋的完整分析流程。

為什麼叫 Seurat？ 名稱取自法國後印象派畫家 Georges Seurat，他以 pointillism (點描法) 聞名——由無數小點組成完整圖像，正如單細胞分析從個別細胞拼湊出組織全貌。

1.2 Seurat 在單細胞分析生態系中的位置

flowchart TB
    subgraph INPUT["資料來源"]
        10X["10X Genomics\nCell Ranger"]
        STAR["STARsolo"]
        KB["kallisto|bustools"]
        CXG["CellxGene\nDiscover"]
    end

    subgraph SEURAT["Seurat v5 核心"]
        QC["Quality Control\n品質控制"]
        NORM["Normalization\n正規化"]
        HVG["Feature Selection\n特徵選擇"]
        DIM["Dimensionality Reduction\nPCA / UMAP / t-SNE"]
        CLUST["Clustering\n聚類分析"]
        DE["Differential Expression\n差異表現分析"]
        INT["Integration\n多批次整合"]
    end

    subgraph DOWNSTREAM["下游分析"]
        TRAJ["Monocle 3\nTrajectory Analysis"]
        COMM["CellChat\nCell Communication"]
        REG["SCENIC\nGene Regulatory Network"]
        SPATIAL["Spatial Analysis\n空間轉錄體"]
    end

    subgraph OUTPUT["輸出"]
        VIZ["互動視覺化"]
        MARKERS["Marker Genes\n標記基因"]
        ANNOT["Cell Type\nAnnotation"]
    end

    INPUT --> SEURAT
    SEURAT --> DOWNSTREAM
    SEURAT --> OUTPUT

1.3 版本歷史與重要里程碑

目前最新穩定版：v5.5.0（2026-04-22 發佈）

2. 安裝指南

2.1 系統需求

2.2 安裝步驟

從 CRAN 安裝（穩定版）：

1install.packages("Seurat")

從 GitHub 安裝（開發版）：

1if (!requireNamespace("remotes", quietly = TRUE)) {
2  install.packages("remotes")
3}
4remotes::install_github("satijalab/seurat", "seurat5", quiet = TRUE)

額外推薦套件：

 1# SCTransform v2 加速
 2install.packages("glmGamPoi")
 3
 4# 快速差異表現分析
 5install.packages("presto")
 6
 7# BPCells 大型資料集磁碟儲存
 8install.packages("BPCells", repos = "https://bnprks.r-universe.dev")
 9
10# Leiden 聚類演算法
11install.packages("leidenbase")

2.3 驗證安裝

1library(Seurat)
2packageVersion("Seurat")
3# [1] '5.5.0'

3. Seurat Object 核心架構

3.1 SeuratObject 資料結構

Seurat 的核心是 SeuratObject——一個多層級的 S4 class (S4 類別) 物件，包裝了所有分析資料：

classDiagram
    class SeuratObject {
        +assays : list~Assay5~
        +meta.data : data.frame
        +reductions : list~DimReduc~
        +graphs : list~Graph~
        +images : list~SpatialImage~
        +active.assay : character
        +active.ident : factor
        +project.name : character
    }

    class Assay5 {
        +layers : list
        +counts : dgCMatrix
        +data : dgCMatrix
        +scale.data : matrix
        +features : character
        +cells : character
        +var.features : character
        +meta.data : data.frame
    }

    class DimReduc {
        +cell.embeddings : matrix
        +feature.loadings : matrix
        +stdev : numeric
        +key : character
    }

    class Graph {
        +adjacency : dgCMatrix
    }

    SeuratObject "1" *-- "1..*" Assay5 : assays
    SeuratObject "1" *-- "0..*" DimReduc : reductions
    SeuratObject "1" *-- "0..*" Graph : graphs

3.2 Layers 系統（v5 新功能）

Seurat v5 引入了 layers (層) 概念，取代 v4 的 slots：

v5 重要改變：多批次資料可 split() 成獨立 layers，各自正規化後再 JoinLayers() 合併，大幅提升整合分析的彈性。

3.3 資料存取方式速查

 1# 細胞名稱
 2Cells(obj)           # 或 colnames(obj)
 3
 4# 基因名稱
 5Features(obj)        # 或 rownames(obj)
 6
 7# 細胞數 / 基因數
 8ncol(obj)            # 細胞數
 9nrow(obj)            # 基因數
10
11# 取得表現資料
12obj[["RNA"]]$counts          # 原始 counts
13obj[["RNA"]]$data            # 正規化後
14obj[["RNA"]]$scale.data      # 縮放後
15
16# 取得 metadata
17obj[[]]                      # 完整 metadata data.frame
18obj$nCount_RNA               # 單一欄位
19
20# 取得降維結果
21Embeddings(obj, "pca")       # PCA cell embeddings
22Loadings(obj, "pca")         # PCA feature loadings
23
24# 彈性取值
25FetchData(obj, vars = c("PC_1", "nFeature_RNA", "MS4A1"))

4. 單細胞 RNA-seq 分析完整流程

4.1 標準分析 Pipeline 總覽

flowchart TD
    A["原始 Count Matrix\n(genes x cells)"] --> B["Step 1: CreateSeuratObject\n建立 Seurat 物件"]
    B --> C["Step 2: QC Filtering\n品質控制過濾"]
    C --> D{"正規化方法選擇"}
    D -->|"標準流程"| E1["Step 3a: NormalizeData\nLogNormalize"]
    D -->|"進階流程"| E2["Step 3b: SCTransform\n(取代 3a+4+5)"]
    E1 --> F["Step 4: FindVariableFeatures\n高變異基因選擇"]
    F --> G["Step 5: ScaleData\n資料縮放"]
    G --> H["Step 6: RunPCA\n主成分分析"]
    E2 --> H
    H --> I["Step 7: ElbowPlot\n決定 PC 數量"]
    I --> J["Step 8: FindNeighbors\nKNN 圖建構"]
    J --> K["Step 9: FindClusters\n社群偵測聚類"]
    K --> L["Step 10: RunUMAP\n非線性降維視覺化"]
    L --> M["Step 11: FindAllMarkers\n差異表現分析"]
    M --> N["Step 12: Cell Type Annotation\n細胞類型註釋"]

    style D fill:#f9f,stroke:#333,stroke-width:2px
    style E2 fill:#bbf,stroke:#333,stroke-width:2px

4.2 兩種正規化路線比較

建議：新手先學標準流程理解每步原理，實戰時改用 SCTransform 獲得更好結果。

5. Quality Control (QC; 品質控制)

5.1 三大 QC 指標

單細胞資料的品質控制是分析成功的關鍵第一步：

5.2 為什麼粒線體比例重要？

當細胞膜 (cell membrane) 破裂時，細胞質 mRNA (cytoplasmic mRNA) 會洩漏，但 mitochondrial mRNA (粒線體 mRNA) 因被雙層膜包裹而保留。因此 percent.mt 高 = 細胞受損。

5.3 QC 實作

 1# 計算粒線體 RNA 比例
 2pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
 3
 4# 視覺化 QC 指標
 5VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
 6
 7# 散點圖檢查相關性
 8FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
 9FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
10
11# 過濾
12pbmc <- subset(pbmc,
13  subset = nFeature_RNA > 200 &    # 排除空液滴
14           nFeature_RNA < 2500 &    # 排除 doublets
15           percent.mt < 5           # 排除瀕死細胞
16)

flowchart LR
    RAW["原始細胞\n(所有 barcodes)"]
    RAW --> F1{"nFeature_RNA > 200?"}
    F1 -->|"No"| DROP1["丟棄\n(空液滴 / 死細胞)"]
    F1 -->|"Yes"| F2{"nFeature_RNA < 2500?"}
    F2 -->|"No"| DROP2["丟棄\n(Doublet 雙細胞)"]
    F2 -->|"Yes"| F3{"percent.mt < 5%?"}
    F3 -->|"No"| DROP3["丟棄\n(瀕死 / 破裂細胞)"]
    F3 -->|"Yes"| PASS["通過 QC\n(高品質細胞)"]

    style PASS fill:#9f9,stroke:#333
    style DROP1 fill:#f99,stroke:#333
    style DROP2 fill:#f99,stroke:#333
    style DROP3 fill:#f99,stroke:#333

5.4 參數選擇指引

6. Normalization (正規化)

6.1 為什麼需要正規化？

不同細胞的 sequencing depth (定序深度) 不同——有些細胞被讀了 10,000 次，有些只有 1,000 次。不正規化的話，高 sequencing depth 的細胞所有基因看起來都「高表現」，但這是技術假象而非生物差異。

6.2 方法一：LogNormalize

1pbmc <- NormalizeData(pbmc,
2  normalization.method = "LogNormalize",
3  scale.factor = 10000
4)

數學公式：

$$normalized_{i,j} = \log\left(\frac{count_{i,j}}{\sum_i count_{i,j}} \times 10000 + 1\right)$$

其中 $i$ = gene, $j$ = cell。步驟：

1. 每個細胞除以總 UMI 數（消除 sequencing depth 差異）
2. 乘以 scale factor (10,000)
3. 取 log1p（壓縮極端值分布）

結果存於 obj[["RNA"]]$data

6.3 方法二：SCTransform（推薦）

1pbmc <- SCTransform(pbmc,
2  vars.to.regress = "percent.mt",  # 迴歸移除粒線體效應
3  verbose = FALSE
4)

SCTransform 演算法原理：

flowchart TD
    A["UMI Count Matrix"] --> B["Step 1: Regularized NB Regression\n對每個基因建立\nNegative Binomial 模型"]
    B --> C["Step 2: 從所有基因學習\n全域 mean-variance 趨勢"]
    C --> D["Step 3: 正則化參數\n(用全域趨勢修正個別基因)"]
    D --> E["Step 4: 計算 Pearson Residuals\n(觀察值 - 預期值) / 標準差"]
    E --> F["輸出：Variance-stabilized\nresiduals 作為 PCA 輸入"]

    G["傳統 LogNormalize"] -.->|"分三步\n各自獨立"| H["NormalizeData\n+ FindVariableFeatures\n+ ScaleData"]
    
    style F fill:#bbf,stroke:#333

SCTransform 優勢：

• 不需要 pseudocount 或 log transformation (偽計數或對數轉換) 等 heuristic (啟發式) 步驟
• 更好的 variance stabilization (變異穩定化)
• 參數更穩健（可安全用 30+ PCs）
• 更細緻的 cluster resolution (聚類解析度)

7. Feature Selection (特徵選擇; 高變異基因)

7.1 為什麼需要選擇高變異基因？

人類基因組有 ~20,000 protein-coding genes (蛋白質編碼基因)，但大多數基因在所有細胞間表現穩定（housekeeping genes; 管家基因）。只有一小部分基因在不同細胞間有顯著差異——這些 highly variable features (HVF; 高變異特徵) 才能揭示細胞間的真實生物差異。

7.2 VST 方法

 1pbmc <- FindVariableFeatures(pbmc,
 2  selection.method = "vst",    # Variance Stabilizing Transformation
 3  nfeatures = 2000             # 選擇前 2,000 個高變異基因
 4)
 5
 6# 查看 top 10 HVF
 7top10 <- head(VariableFeatures(pbmc), 10)
 8print(top10)
 9
10# 視覺化
11VariableFeaturePlot(pbmc)
12LabelPoints(plot = VariableFeaturePlot(pbmc),
13            points = top10, repel = TRUE)

VST 演算法：

1. 對每個基因計算 mean expression (平均表現量) 和 variance (變異量)
2. 用 LOESS 回歸建立 mean-variance relationship (均值-變異關係)
3. 計算每個基因的 standardized variance (標準化變異)
4. 排序後取前 N 個（預設 2,000）

8. Scaling (資料縮放)

8.1 為什麼需要縮放？

PCA 是 variance-based (基於變異) 的方法。如果不做 scaling (縮放)，高表現基因會主導 PCA 結果——不是因為它們生物學上更重要，而是因為它們的數值更大。

1pbmc <- ScaleData(pbmc, features = rownames(pbmc))

ScaleData 做了什麼：

1. Centering (置中)：每個基因減去平均值 → mean = 0
2. Scaling (縮放)：每個基因除以標準差 → variance = 1

結果存於 obj[["RNA"]]$scale.data

選擇性迴歸（移除混淆變數）：

1pbmc <- ScaleData(pbmc, vars.to.regress = c("percent.mt", "nCount_RNA"))

9. Dimensionality Reduction (降維分析)

9.1 PCA — Principal Component Analysis (主成分分析)

為什麼需要 PCA？ 2,000 個高變異基因 = 2,000 維空間。人無法直觀理解高維空間，且 noise (雜訊) 在高維中被放大。PCA 找出資料中變異最大的方向（principal components; 主成分），將 2,000 維壓縮到 10-50 維，保留主要訊號並移除雜訊。

 1pbmc <- RunPCA(pbmc, features = VariableFeatures(pbmc))
 2
 3# 檢視各 PC 的 top genes
 4print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
 5
 6# 視覺化 PC loadings
 7VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
 8
 9# PCA 散點圖
10DimPlot(pbmc, reduction = "pca")
11
12# Heatmap 檢視 PC 與基因的關聯
13DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)

PCA 演算法概念：

flowchart LR
    subgraph INPUT["輸入"]
        M["Scale.data Matrix\n2000 genes x 2638 cells"]
    end

    subgraph PCA["PCA 過程"]
        COV["計算 Covariance Matrix\n(共變異矩陣)"]
        EIG["Eigendecomposition\n(特徵值分解)"]
        SORT["按 Eigenvalue 排序\n(變異量大 → 小)"]
    end

    subgraph OUTPUT["輸出"]
        EMB["Cell Embeddings\n(每個細胞在 PC 空間的座標)"]
        LOAD["Feature Loadings\n(每個基因對 PC 的貢獻)"]
        SDEV["Standard Deviations\n(各 PC 解釋的變異量)"]
    end

    INPUT --> PCA --> OUTPUT

9.2 決定使用多少 PCs

1ElbowPlot(pbmc)

ElbowPlot 顯示每個 PC 解釋的 variance percentage (變異百分比)。找到 elbow point (肘部點)——曲線開始趨於平坦的位置：

• Elbow 之前的 PCs = signal (訊號)
• Elbow 之後的 PCs = noise (雜訊)

經驗法則：寧可多選幾個（10-15），不要少選。選多了只是加一點雜訊，選少了可能遺失生物訊號。

9.3 UMAP vs t-SNE 視覺化

在 PCA 降維後，再用非線性方法將 10-50 維降到 2 維以便視覺化：

1# UMAP（推薦，速度快、保留全域結構）
2pbmc <- RunUMAP(pbmc, dims = 1:10)
3
4# t-SNE（替代選項，保留局部結構）
5pbmc <- RunTSNE(pbmc, dims = 1:10)

⚠️ 重要提醒：UMAP/t-SNE 僅供視覺化，不應基於這些圖做生物學結論。所有定量分析（聚類、差異表現）都在 PCA 空間完成。

10. Clustering (聚類分析)

10.1 圖論聚類演算法

Seurat 使用 graph-based clustering (圖論聚類)，這是目前單細胞領域最主流的方法：

1# Step 1: 建構 KNN 圖
2pbmc <- FindNeighbors(pbmc, dims = 1:10)
3
4# Step 2: 社群偵測
5pbmc <- FindClusters(pbmc, resolution = 0.5)
6
7# 查看聚類結果
8head(Idents(pbmc))
9table(Idents(pbmc))

flowchart TD
    A["PCA Embeddings\n(cells in PC space)"] --> B["Step 1: KNN Graph\n找每個細胞的 K 個最近鄰居\n(Euclidean distance in PC space)"]
    B --> C["Step 2: SNN Graph\nShared Nearest Neighbor\n計算任兩細胞共享的鄰居數\n(Jaccard similarity)"]
    C --> D{"聚類演算法"}
    D -->|"預設"| E1["Louvain Algorithm\n最大化 modularity\n(社群內部連結密度)"]
    D -->|"替代"| E2["Leiden Algorithm\n改良版 Louvain\n(保證連通性)"]
    D -->|"替代"| E3["SLM Algorithm\nSmart Local Moving"]
    E1 --> F["輸出：每個細胞的\nCluster Identity"]
    E2 --> F
    E3 --> F

    style C fill:#ffd,stroke:#333
    style E1 fill:#dfd,stroke:#333

10.2 Resolution 參數

resolution 控制聚類的粒度——值越大，cluster 越多：

1# 可以測試多個 resolution
2pbmc <- FindClusters(pbmc, resolution = c(0.1, 0.3, 0.5, 0.8, 1.0))
3
4# 結果都存在 metadata 中
5head(pbmc@meta.data[, grep("RNA_snn_res", colnames(pbmc@meta.data))])

10.3 Louvain 演算法簡介

Louvain algorithm (Louvain 演算法) 是一種 modularity optimization (模組化最佳化) 方法：

1. 初始化：每個 node (節點/細胞) 自成一個 community (社群)
2. Phase 1 — 局部移動：
   • 對每個 node，計算將它移到鄰居社群後 modularity 的變化 (ΔQ)
   • 將 node 移到 ΔQ 最大的社群（若 ΔQ > 0）
   • 重複直到沒有 node 可以移動
3. Phase 2 — 聚合：
   • 將同一社群的 nodes 合併為一個 super-node
   • 重建新的 graph
4. 迭代：回到 Phase 1，直到 modularity 不再增加

Modularity Q 衡量的是：社群內部的連結密度是否顯著高於隨機 graph 的預期？

11. Differential Expression (DE; 差異表現分析)

11.1 尋找 Marker Genes (標記基因)

 1# 找特定 cluster 的 markers（vs. 所有其他 cluster）
 2cluster2.markers <- FindMarkers(pbmc, ident.1 = 2)
 3head(cluster2.markers, n = 5)
 4
 5# 比較兩組 clusters
 6cluster5vs03.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3))
 7
 8# 自動找所有 cluster 的 markers
 9pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE)
10
11# 使用 presto 加速（建議）
12pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

11.2 輸出欄位說明

11.3 可用的統計檢定

11.4 Marker 視覺化

 1# Violin plot — 顯示基因表現分布
 2VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
 3
 4# Feature plot — 在 UMAP 上顯示基因表現
 5FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14",
 6                                "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
 7
 8# Heatmap — 各 cluster top markers
 9top10 <- pbmc.markers %>%
10  group_by(cluster) %>%
11  dplyr::filter(avg_log2FC > 1) %>%
12  slice_head(n = 10) %>%
13  ungroup()
14DoHeatmap(pbmc, features = top10$gene)
15
16# Dot plot — 同時顯示表現量與表現比例
17DotPlot(pbmc, features = unique(top10$gene)) + RotatedAxis()

12. Cell Type Annotation (細胞類型註釋)

12.1 手動註釋策略

根據已知的 canonical marker genes (經典標記基因)，對照聚類結果進行註釋：

1new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B",
2                      "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
3names(new.cluster.ids) <- levels(pbmc)
4pbmc <- RenameIdents(pbmc, new.cluster.ids)
5
6DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

12.2 自動註釋工具

13. PBMC 3K 完整實戰範例

13.1 範例概述

這是 Seurat 最經典的入門範例——分析 10X Genomics 提供的 2,700 個 Peripheral Blood Mononuclear Cells (PBMC; 周邊血液單核細胞)。

資料集：

• 來源：10X Genomics（免費下載）
• 細胞數：2,700
• 技術：Illumina NextSeq 500
• 初始基因數：13,714

13.2 完整程式碼（可直接執行）

 1# ── 載入套件 ──
 2library(dplyr)
 3library(Seurat)
 4library(patchwork)
 5
 6# ── Step 1: 讀取 10X 資料 ──
 7pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")
 8pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k",
 9                           min.cells = 3, min.features = 200)
10# min.cells = 3  → 只保留至少在 3 個細胞中表現的基因
11# min.features = 200 → 只保留至少偵測到 200 個基因的細胞
12
13# ── Step 2: QC ──
14pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
15VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
16pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
17
18# ── Step 3: 正規化 ──
19pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
20
21# ── Step 4: 高變異基因 ──
22pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
23
24# ── Step 5: 縮放 ──
25pbmc <- ScaleData(pbmc, features = rownames(pbmc))
26
27# ── Step 6: PCA ──
28pbmc <- RunPCA(pbmc, features = VariableFeatures(pbmc))
29
30# ── Step 7: 決定 PC 數 ──
31ElbowPlot(pbmc)
32
33# ── Step 8-9: 聚類 ──
34pbmc <- FindNeighbors(pbmc, dims = 1:10)
35pbmc <- FindClusters(pbmc, resolution = 0.5)
36
37# ── Step 10: UMAP ──
38pbmc <- RunUMAP(pbmc, dims = 1:10)
39DimPlot(pbmc, reduction = "umap")
40
41# ── Step 11: 差異表現 ──
42pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE,
43                                min.pct = 0.25, logfc.threshold = 0.25)
44
45# ── Step 12: 細胞類型註釋 ──
46new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B",
47                      "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
48names(new.cluster.ids) <- levels(pbmc)
49pbmc <- RenameIdents(pbmc, new.cluster.ids)
50
51# ── 最終視覺化 ──
52DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
53
54# ── 儲存結果 ──
55saveRDS(pbmc, file = "pbmc3k_final.rds")

13.3 SCTransform 簡化版（推薦）

 1library(Seurat)
 2library(sctransform)
 3
 4pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")
 5pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k",
 6                           min.cells = 3, min.features = 200)
 7pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
 8pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
 9
10# SCTransform 一步完成正規化 + 特徵選擇 + 縮放
11pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE) %>%
12  RunPCA(verbose = FALSE) %>%
13  FindNeighbors(dims = 1:30, verbose = FALSE) %>%
14  FindClusters(verbose = FALSE) %>%
15  RunUMAP(dims = 1:30, verbose = FALSE)
16
17DimPlot(pbmc, label = TRUE)

14. Data Integration (資料整合)

14.1 為什麼需要整合？

當合併來自不同 batch (批次)、donor (捐贈者) 或 condition (實驗條件) 的單細胞資料時，技術性的 batch effects (批次效應) 會導致：

• 同類型細胞因來源不同而被分到不同 cluster
• 聚類結果反映的是「來自哪個批次」而非「什麼細胞類型」

14.2 可用整合方法

flowchart TD
    A["多批次 Seurat Objects"] --> B["merge() 合併"]
    B --> C["NormalizeData + FindVariableFeatures + ScaleData + RunPCA"]
    C --> D{"選擇整合方法"}
    D --> E1["CCAIntegration\n(Canonical Correlation Analysis)"]
    D --> E2["RPCAIntegration\n(Reciprocal PCA)\n★ 推薦大型資料集"]
    D --> E3["HarmonyIntegration\n(Harmony)"]
    D --> E4["scVIIntegration\n(scVI, 需 Python)"]
    E1 --> F["IntegrateLayers()"]
    E2 --> F
    E3 --> F
    E4 --> F
    F --> G["JoinLayers() → 後續分析"]

    style E2 fill:#bfb,stroke:#333,stroke-width:2px

14.3 整合程式碼範例

 1# 假設已有合併後的物件 merged_obj，含 split layers
 2merged_obj <- NormalizeData(merged_obj)
 3merged_obj <- FindVariableFeatures(merged_obj)
 4merged_obj <- ScaleData(merged_obj)
 5merged_obj <- RunPCA(merged_obj)
 6
 7# 整合（以 RPCA 為例，適合大型資料集）
 8merged_obj <- IntegrateLayers(
 9  object = merged_obj,
10  method = RPCAIntegration,
11  orig.reduction = "pca",
12  new.reduction = "integrated.rpca",
13  verbose = FALSE
14)
15
16# 合併 layers 後繼續分析
17merged_obj[["RNA"]] <- JoinLayers(merged_obj[["RNA"]])
18merged_obj <- FindNeighbors(merged_obj, reduction = "integrated.rpca", dims = 1:30)
19merged_obj <- FindClusters(merged_obj, resolution = 0.5)
20merged_obj <- RunUMAP(merged_obj, reduction = "integrated.rpca", dims = 1:30)

14.4 整合方法比較

15. CellxGene 資料庫介紹

15.1 什麼是 CellxGene？

CZ CELLxGENE Discover 是由 Chan Zuckerberg Initiative (CZI) 建立的單細胞資料平台，提供：

• 數千個 公開單細胞資料集的瀏覽與下載
• 互動式視覺化 探索工具（CellxGene Annotate）
• Census API 程式化存取全資料庫
• Gene Expression 跨資料集的基因表現查詢
• Differential Expression 線上差異表現分析

GitHub：https://github.com/chanzuckerberg/cellxgene（778 ⭐、MIT License）

15.2 資料格式

15.3 從 CellxGene 下載資料到 Seurat

方法一：網頁下載

1. 前往 https://cellxgene.cziscience.com/datasets
2. 用 filter 篩選（Organism / Tissue / Disease / Cell Count）
3. 選擇資料集 → Download → 選 .rds 格式
4. R 中載入：

1obj <- readRDS("downloaded_dataset.rds")

方法二：Census API（程式化）

 1# 安裝 cellxgene.census
 2if (!requireNamespace("cellxgene.census", quietly = TRUE)) {
 3  install.packages("cellxgene.census", repos = "https://chanzuckerberg.r-universe.dev")
 4}
 5
 6library(cellxgene.census)
 7
 8# 開啟 census
 9census <- open_soma()
10
11# 查詢特定組織的小型資料集
12adata <- get_seurat(
13  census,
14  organism = "Homo sapiens",
15  obs_value_filter = "tissue_general == 'blood' && disease == 'normal'",
16  obs_column_names = c("cell_type", "tissue", "donor_id")
17)
18
19census$close()

15.4 CellxGene → Seurat 資料流

flowchart LR
    CXG["CellxGene\nDiscover"] -->|".rds 下載"| RDS["readRDS()"]
    CXG -->|".h5ad 下載"| H5AD["SeuratDisk::LoadH5Seurat()\n或 anndata + reticulate"]
    CXG -->|"Census API"| API["cellxgene.census\nget_seurat()"]
    RDS --> SEURAT["Seurat Object\n(可直接分析)"]
    H5AD --> SEURAT
    API --> SEURAT
    SEURAT --> ANALYSIS["標準 Seurat\n分析流程"]

15.5 CellxGene Annotate 本地視覺化

1# 安裝
2pip install cellxgene
3
4# 啟動（用 h5ad 檔案）
5cellxgene launch my_dataset.h5ad
6# 瀏覽器自動開啟互動式探索介面

16. 常用指令速查表

16.1 基本分析 Pipeline

 1# 標準流程（7 行）
 2obj <- NormalizeData(obj)
 3obj <- FindVariableFeatures(obj)
 4obj <- ScaleData(obj)
 5obj <- RunPCA(obj)
 6obj <- FindNeighbors(obj, dims = 1:30)
 7obj <- FindClusters(obj)
 8obj <- RunUMAP(obj, dims = 1:30)
 9
10# SCTransform 流程（5 行）
11obj <- SCTransform(obj) %>%
12  RunPCA() %>%
13  FindNeighbors(dims = 1:30) %>%
14  FindClusters() %>%
15  RunUMAP(dims = 1:30)

16.2 Subsetting (子集操作)

 1# 依 cluster identity 取子集
 2t_cells <- subset(obj, idents = c("CD4 T", "CD8 T"))
 3
 4# 依基因表現過濾
 5high_cd3 <- subset(obj, subset = CD3E > 2)
 6
 7# 依 metadata 過濾
 8treated <- subset(obj, subset = condition == "treated")
 9
10# Downsample（每個 identity 最多 100 cells）
11downsampled <- subset(obj, downsample = 100)

16.3 視覺化指令

16.4 資料存取

 1# 物件資訊
 2Cells(obj)                           # 細胞名
 3Features(obj)                        # 基因名
 4ncol(obj)                            # 細胞數
 5nrow(obj)                            # 基因數
 6Layers(obj)                          # 列出所有 layers
 7Assays(obj)                          # 列出所有 assays
 8VariableFeatures(obj)                # 高變異基因
 9Idents(obj)                          # 當前 identity
10obj[[]]                              # 完整 metadata
11Embeddings(obj, "pca")               # PCA embeddings
12Loadings(obj, "pca")                 # PCA loadings
13FetchData(obj, vars = c("PC_1", "gene1", "nCount_RNA"))  # 彈性取值

17. 資安掃描報告

17.1 掃描結果

對 Seurat repo（v5.5.0.9006）的 R/、src/、vignettes/ 目錄進行模式掃描：

17.2 安全評估

整體評級：🟢 低風險

• 所有 eval() 呼叫皆為內部函式名稱解析，不接受使用者輸入
• rlang::exec() 是安全的程式化函式呼叫機制
• HTTP URLs 為文件參考連結，不涉及資料傳輸
• 無密碼、API key、token 等敏感資訊
• 作為 CRAN 套件，通過 CRAN 審查標準
• MIT License，無授權風險

18. FAQ (常見問題)

Q1: Seurat v4 和 v5 有什麼差異？

最大差異是 layers 系統：v5 的 assay 內部用 layers 儲存 counts/data/scale.data，支援 split 與 join 操作。原本的 GetAssayData(obj, slot = "counts") 改為 obj[["RNA"]]$counts。

Q2: 該用 LogNormalize 還是 SCTransform？

• 學習用：先用 LogNormalize 理解每步原理
• 實戰用：SCTransform v2 更穩健，推薦用於正式分析
• 整合多批次：兩者皆可，但 SCTransform 需特殊處理（PrepSCTFindMarkers()）

Q3: Resolution 該設多少？

沒有標準答案。建議用 clustree 套件或測試多個 resolution 後，選擇能分出已知生物學差異的值。

Q4: 多大的資料集 Seurat 能處理？

• < 100K cells：標準記憶體即可
• 100K–1M cells：建議使用 sketch-based analysis
• > 1M cells：使用 BPCells on-disk 模式

Q5: 如何從 h5ad 轉到 Seurat？

 1# 方法 1: 使用 SeuratDisk（deprecated 但仍可用）
 2library(SeuratDisk)
 3Convert("data.h5ad", dest = "h5seurat")
 4obj <- LoadH5Seurat("data.h5seurat")
 5
 6# 方法 2: 使用 anndata + reticulate
 7library(reticulate)
 8ad <- import("anndata")
 9adata <- ad$read_h5ad("data.h5ad")
10# 手動轉換...

Q6: 如何保存和載入 Seurat 物件？

1saveRDS(obj, file = "my_analysis.rds")
2obj <- readRDS("my_analysis.rds")

19. 進階技巧與效能優化

19.1 並行運算

1library(future)
2plan("multisession", workers = 4)  # 使用 4 個 workers
3# 之後的 Seurat 操作自動並行

19.2 大規模資料集策略

 1# Sketch-based analysis（抽樣代表子集）
 2obj <- SketchData(
 3  object = obj,
 4  ncells = 5000,
 5  method = "LeverageScore",
 6  sketched.assay = "sketch"
 7)
 8
 9# BPCells on-disk（磁碟儲存）
10library(BPCells)
11# 將 counts 存到磁碟
12write_matrix_dir(mat = obj[["RNA"]]$counts, dir = "bpcells_counts")
13bpcells_mat <- open_matrix_dir("bpcells_counts")
14obj[["RNA"]]$counts <- bpcells_mat

19.3 Cell Cycle Regression (細胞週期迴歸)

1# Seurat 提供 cell cycle genes list
2s.genes <- cc.genes.updated.2019$s.genes
3g2m.genes <- cc.genes.updated.2019$g2m.genes
4
5# 評分
6obj <- CellCycleScoring(obj, s.features = s.genes, g2m.features = g2m.genes)
7
8# 迴歸移除（若 cell cycle 不是研究目標）
9obj <- ScaleData(obj, vars.to.regress = c("S.Score", "G2M.Score"))

19.4 Pseudobulk Analysis (偽批量分析)

1# 將單細胞資料聚合為偽批量（用於嚴格的統計檢定）
2bulk <- AggregateExpression(obj,
3  group.by = c("celltype", "donor_id"),
4  return.seurat = TRUE
5)

20. Downsample 策略概覽

Downsampling (降取樣) 在單細胞分析中極為常見，用於：

• 加速分析：大型資料集的快速原型
• 平衡組別：避免大 cluster 主導分析結果
• 跨資料集比較：統一各批次的細胞數

20.1 Seurat 內建 downsample

1# 每個 identity 最多取 100 個 cells
2obj_down <- subset(obj, downsample = 100)
3
4# 隨機抽樣（不依 identity）
5set.seed(42)
6cells_to_keep <- sample(Cells(obj), size = 5000)
7obj_down <- subset(obj, cells = cells_to_keep)

20.2 Sketch-based Analysis（進階 downsample）

1# 使用 leverage score sampling（保留稀有細胞類型的代表性）
2obj <- SketchData(obj,
3  ncells = 5000,
4  method = "LeverageScore",
5  sketched.assay = "sketch"
6)

21. 重點摘要 Checklist

• 安裝：R >= 4.0.0 + install.packages("Seurat") + glmGamPoi + presto
• QC 三指標：nFeature_RNA / nCount_RNA / percent.mt
• 正規化選擇：LogNormalize（教學用）vs SCTransform（實戰推薦）
• PCA：用 ElbowPlot 決定 PC 數，寧多勿少
• 聚類：FindNeighbors → FindClusters，resolution 依需調整
• UMAP：僅供視覺化，不做定量結論
• DE 分析：FindAllMarkers + Wilcoxon（預設）
• 註釋：用 canonical markers 或 Azimuth 自動標記
• 整合：多批次用 IntegrateLayers（RPCA 推薦大資料集）
• CellxGene：下載 .rds 直接 readRDS，或用 Census API
• 大規模：> 100K cells 用 sketch analysis 或 BPCells

22. 進一步閱讀與引用

22.1 官方資源

22.2 核心論文

1. Seurat v5: Hao Y, et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nature Biotechnology (2023). DOI: 10.1038/s41587-023-01767-y
2. Seurat v4: Hao Y, et al. Integrated analysis of multimodal single-cell data. Cell (2021). DOI: 10.1016/j.cell.2021.04.048
3. Seurat v3: Stuart T, et al. Comprehensive Integration of Single-Cell Data. Cell (2019). DOI: 10.1016/j.cell.2019.05.031
4. SCTransform: Hafemeister C & Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biology (2019).
5. SCTransform v2: Choudhary S & Satija R. Comparison and evaluation of statistical error models for scRNA-seq. Genome Biology (2022).

22.3 相關教學

• Orchestrating Single-Cell Analysis with Bioconductor (OSCA): https://bioconductor.org/books/release/OSCA/
• Single-cell best practices: https://www.sc-best-practices.org/