376_单细胞数据标准化流程¶

library(Seurat)     # v5.5 单细胞分析包

# 加载10x Genomics数据（CellRanger输出目录）
counts <- Read10X(data.dir = "filtered_feature_bc_matrix/")

# 创建 Seurat 对象
seurat_obj <- CreateSeuratObject(
  counts = counts,
  project = "ScRNA_tutorial",
  min.cells = 3,       # 基因至少在3个细胞中检测到
  min.features = 200   # 细胞至少表达200个基因
)
cat("细胞数:", ncol(seurat_obj), "\n")

# 质控：计算线粒体基因比例
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(
  seurat_obj,
  pattern = "^MT-"    # 人类线粒体基因以MT-开头（小鼠用mt-）
)

# 质控过滤
seurat_obj <- subset(
  seurat_obj,
  subset = nFeature_RNA > 200 &    # 至少200个基因
           nFeature_RNA < 6000 &   # 不超过6000（可能是双细胞）
           percent.mt < 15         # 线粒体比例不超过15%
)
cat("过滤后细胞数:", ncol(seurat_obj), "\n")

# ============================================================
# 方法A：经典三步法（适合理解原理）
# ============================================================
seurat_obj <- NormalizeData(
  seurat_obj,
  normalization.method = "LogNormalize",  # log(CPM/10+1)
  scale.factor = 10000                    # 归一化到1万
)

seurat_obj <- FindVariableFeatures(
  seurat_obj,
  selection.method = "vst",    # 方差稳定变换法（推荐）
  nfeatures = 2000             # 保留2000个高变基因
)

# 查看最高变基因
top10 <- head(VariableFeatures(seurat_obj), 10)
cat("Top10 高变基因:", top10, "\n")

seurat_obj <- ScaleData(
  seurat_obj,
  vars.to.regress = c("percent.mt", "nCount_RNA")  # 回归去除线粒体和测序深度
)

# ============================================================
# 方法B：SCTransform（更推荐，一步搞定）
# ============================================================
seurat_obj <- SCTransform(
  seurat_obj,
  vars.to.regress = "percent.mt",  # 回归去除线粒体比例
  verbose = FALSE,                  # 不显示进度条
  vst.flavor = "v2"                 # 使用v2（Pearson残差，Seurat v5默认）
)

import scanpy as sc    # 单细胞分析Python包
import numpy as np

# 设置随机种子（保证结果可重现）
sc.settings.verbosity = 2   # 显示关键步骤信息
sc.settings.n_jobs = 4      # 多线程

# 读取 10x 数据
adata = sc.read_10x_mtx(
    "filtered_feature_bc_matrix/",   # CellRanger输出目录
    var_names="gene_symbols",         # 用基因symbol命名（非Ensembl ID）
    cache=True                        # 缓存加速二次读取
)
adata.var_names_make_unique()         # 确保基因名唯一

# 质控
adata.var["mt"] = adata.var_names.str.startswith("MT-")  # 标记线粒体基因
sc.pp.calculate_qc_metrics(
    adata,
    qc_vars=["mt"],      # 计算线粒体比例
    inplace=True         # 结果写回 adata.obs
)

# 质控可视化
sc.pl.violin(
    adata,
    ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
    jitter=0.4,
    multi_panel=True,
    save="_qc_violin.pdf"
)

# 过滤低质量细胞
sc.pp.filter_cells(adata, min_genes=200)    # 至少200个基因
sc.pp.filter_genes(adata, min_cells=3)     # 基因至少在3个细胞中表达
adata = adata[adata.obs.pct_counts_mt < 15]   # 线粒体比例<15%

print(f"过滤后：{adata.n_obs} 细胞, {adata.n_vars} 基因")

# 保存原始 counts（后续DEG分析需要）
adata.raw = adata   # 保存未标准化数据到 .raw 层

# 标准化：总量归一化
sc.pp.normalize_total(
    adata,
    target_sum=1e4    # 归一化到每细胞10000个reads
)

# log1p 变换（log(x+1)，消除大值影响）
sc.pp.log1p(adata)

# 筛选高变基因
sc.pp.highly_variable_genes(
    adata,
    n_top_genes=2000,    # 保留2000个高变基因
    flavor="seurat_v3",  # 使用Seurat v3方法（推荐）
    batch_key=None       # 如果有批次，填批次列名
)
print(f"高变基因数: {adata.var.highly_variable.sum()}")

# z-score 标准化
sc.pp.scale(adata, max_value=10)    # 截断到最大值10，防止极端值