单细胞对比：Seurat/Scanpy/Monocle¶

# === Seurat v5 标准流程 ===
library(Seurat)

# 读取10x数据
data <- Read10X(data.dir="filtered_feature_bc_matrix/")  # 读取10x目录
obj <- CreateSeuratObject(
  counts = data,                           # 计数矩阵
  project = "my_project",                  # 项目名
  min.cells = 3,                           # 基因至少在3个细胞中表达
  min.features = 200                       # 细胞至少表达200个基因
)

# 质控
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern="^MT-")  # 线粒体比例
obj <- subset(obj,
  nFeature_RNA > 200 &                    # 基因数>200
  nFeature_RNA < 5000 &                   # 基因数<5000（去双细胞）
  percent.mt < 20                          # 线粒体<20%
)

# 标准化+高变基因+PCA
obj <- NormalizeData(obj) %>%              # log标准化
  FindVariableFeatures(nfeatures=2000) %>% # 找2000个高变基因
  ScaleData() %>%                          # 缩放
  RunPCA(npcs=50)                          # PCA降维

# 聚类+UMAP
obj <- FindNeighbors(obj, dims=1:30) %>%   # KNN图
  FindClusters(resolution=0.5) %>%          # Leiden聚类
  RunUMAP(dims=1:30)                       # UMAP可视化

# 可视化
DimPlot(obj, reduction="umap", label=TRUE)  # UMAP图
FeaturePlot(obj, features=c("CD3D","CD14"))  # 标记基因

# 差异基因
markers <- FindAllMarkers(obj,
  only.pos = TRUE,                         # 只看上调
  min.pct = 0.25,                          # 至少25%细胞表达
  logfc.threshold = 0.25                   # logFC阈值
)

# === Scanpy 标准流程 ===
import scanpy as sc

# 读取数据
adata = sc.read_10x_mtx("filtered_feature_bc_matrix/")  # 读取10x数据

# 质控
sc.pp.filter_cells(adata, min_genes=200)   # 细胞最少200个基因
sc.pp.filter_genes(adata, min_cells=3)     # 基因至少在3个细胞
adata.var['mt'] = adata.var_names.str.startswith('MT-')  # 标记线粒体基因
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)
adata = adata[adata.obs.pct_counts_mt < 20, :]  # 过滤线粒体>20%

# 标准化
sc.pp.normalize_total(adata, target_sum=1e4)  # 总计数标准化
sc.pp.log1p(adata)                          # log转换
sc.pp.highly_variable_genes(adata, n_top_genes=2000)  # 高变基因
adata = adata[:, adata.var.highly_variable]  # 子集

# PCA+聚类+UMAP
sc.pp.scale(adata)                          # 缩放
sc.tl.pca(adata, n_comps=50)               # PCA
sc.pp.neighbors(adata, n_pcs=30)           # KNN图
sc.tl.leiden(adata, resolution=0.5)        # Leiden聚类
sc.tl.umap(adata)                          # UMAP

# 可视化
sc.pl.umap(adata, color='leiden')          # UMAP图
sc.pl.umap(adata, color=['CD3D','CD14'])   # 标记基因

# 差异基因
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=20)

# === Monocle3 轨迹分析 ===
library(monocle3)

# 从Seurat转换
cds <- as.cell_data_set(obj)                # Seurat→Monocle3
cds <- cluster_cells(cds)                   # 聚类
cds <- learn_graph(cds)                     # ★学习轨迹图
cds <- order_cells(cds)                     # ★排序细胞（需指定根节点）

# 可视化轨迹
plot_cells(cds,
  color_cells_by = "pseudotime",           # 用拟时间着色
  label_groups_by_cluster = FALSE,
  label_leaves = FALSE,
  label_branch_points = FALSE
)

# 找沿轨迹变化的基因
track_genes <- graph_test(cds, neighbor_graph="principal_graph")
sig_genes <- track_genes[track_genes$q_value < 0.05, ]

# === Seurat ↔ Scanpy 互转 ===
# Seurat → h5ad（给Scanpy用）
library(SeuratDisk)
SaveH5Seurat(obj, filename="data.h5seurat")
Convert("data.h5seurat", dest="h5ad")       # 生成data.h5ad

# h5ad → Seurat
Convert("data.h5ad", dest="h5seurat")
obj <- LoadH5Seurat("data.h5seurat")

# Python侧
import scanpy as sc
# 读取h5ad
adata = sc.read_h5ad("data.h5ad")
# 保存
adata.write("data.h5ad")

# === 单细胞速查 ===
# Seurat: Read10X → CreateSeuratObject → NormalizeData → FindVariableFeatures
#   → ScaleData → RunPCA → FindNeighbors → FindClusters → RunUMAP → FindAllMarkers

# Scanpy: read_10x_mtx → filter → normalize_total → log1p → highly_variable_genes
#   → scale → pca → neighbors → leiden → umap → rank_genes_groups

# Monocle3: as.cell_data_set → cluster_cells → learn_graph → order_cells

# 互转: SaveH5Seurat → Convert("h5seurat","h5ad")
# 选择: 全流程→Seurat | Python→Scanpy | 轨迹→Monocle3 | 大数据→Scanpy