跳转至

174_单细胞细胞周期评分

一句话概述

在单细胞RNA-seq分析中利用细胞周期相关基因集对每个细胞进行S期和G2M期评分，从而判断细胞所处的细胞周期阶段并在下游分析中消除或保留周期效应。

核心知识点表格

知识点说明
细胞周期三阶段G1期（生长）、S期（DNA合成）、G2/M期（分裂）
评分原理基于S期和G2M期标记基因集的平均表达与背景基因的差值
Seurat方法CellCycleScoring函数，基于Tirosh 2016基因集
Scanpy方法score_genes_cell_cycle函数，原理类似
基因集来源Tirosh et al. 2016 Science论文定义的基因列表
应用场景判断增殖状态、回归细胞周期效应、肿瘤增殖分析
回归策略可选择回归S.Score和G2M.Score差异或绝对值
注意事项某些生物学场景下细胞周期效应是有意义的信号

步骤详解

第一步：理解细胞周期评分的原理

白话解释：细胞周期就像一个时钟，细胞按照G1→S→G2→M的顺序周而复始地运转。在不同阶段，细胞会表达不同的基因。通过检测这些特征基因的表达水平，我们就能判断每个细胞处于哪个阶段。

技术细节：评分方法基于Tirosh等人2016年定义的两组基因——S期基因集（43个基因）和G2/M期基因集（54个基因）。评分算法计算目标基因集的平均表达值减去一组随机背景基因的平均表达值。如果S分数和G2M分数都很低，则认为细胞处于G1期。

# Seurat中的基因集
s.genes <- cc.genes$s.genes    # 43个S期基因
g2m.genes <- cc.genes$g2m.genes # 54个G2M期基因

# 查看基因列表
head(s.genes, 10)
# [1] "MCM5"  "PCNA"  "TYMS"  "FEN1"  "MCM2"  "MCM4"  "RRM1"  "UNG"  "GINS2" "MCM6"
head(g2m.genes, 10)
# [1] "HMGB2"  "CDK1"   "NUSAP1" "UBE2C"  "BIRC5"  "TPX2"   "TOP2A"  "NDC80"  "CKS2"  "NUF2"

第二步：使用Seurat进行细胞周期评分

白话解释：Seurat的CellCycleScoring函数会自动为每个细胞计算S期和G2M期的分数，并据此判定细胞处于G1、S还是G2M期。

技术细节：函数内部使用AddModuleScore方法，该方法将基因分成若干组（bins），从每组中随机选取对照基因计算背景值，然后用目标基因集的均值减去背景均值作为得分。

library(Seurat)

# 加载数据并预处理
pbmc <- Read10X(data.dir = "filtered_feature_bc_matrix/")
pbmc <- CreateSeuratObject(counts = pbmc, min.cells = 3, min.features = 200)
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
pbmc <- ScaleData(pbmc)

# 细胞周期评分
pbmc <- CellCycleScoring(pbmc,
                          s.features = cc.genes$s.genes,
                          g2m.features = cc.genes$g2m.genes,
                          set.ident = FALSE)

# 查看结果
head(pbmc@meta.data[, c("S.Score", "G2M.Score", "Phase")])
#                    S.Score   G2M.Score Phase
# AAACATACAACCAC-1 -0.04326   0.01254    G2M
# AAACATTGAGCTAC-1 -0.08319  -0.17803    G1
# AAACATTGATCAGC-1  0.14520   0.06780    S

# 统计各阶段细胞数
table(pbmc$Phase)
# G1   G2M   S
# 5432  1230  1338

第三步：使用Scanpy进行细胞周期评分（Python）

白话解释：Scanpy的方法和Seurat类似，也是基于同一组基因来打分。不同的是它使用Python实现，适合习惯Python的用户。

import scanpy as sc
import pandas as pd

# 加载数据
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")

# 预处理
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# 定义基因集（Tirosh 2016）
s_genes = ['MCM5', 'PCNA', 'TYMS', 'FEN1', 'MCM2', 'MCM4', 'RRM1',
           'UNG', 'GINS2', 'MCM6', 'CDCA7', 'DTL', 'PRIM1', 'UHRF1',
           'MLF1IP', 'HELLS', 'RFC2', 'RPA2', 'NASP', 'RAD51AP1',
           'GMNN', 'WDR76', 'SLBP', 'CCNE2', 'UBR7', 'POLD3', 'MSH2',
           'ATAD2', 'RAD51', 'RRM2', 'CDC45', 'CDC6', 'EXO1', 'TIPIN',
           'DSCC1', 'BLM', 'CASP8AP2', 'USP1', 'CLSPN', 'POLA1',
           'CHAF1B', 'BRIP1', 'E2F8']

g2m_genes = ['HMGB2', 'CDK1', 'NUSAP1', 'UBE2C', 'BIRC5', 'TPX2',
             'TOP2A', 'NDC80', 'CKS2', 'NUF2', 'CKS1B', 'MKI67',
             'TMPO', 'CENPF', 'TACC3', 'FAM64A', 'SMC4', 'CCNB2',
             'CKAP2L', 'CKAP2', 'AURKB', 'BUB1', 'KIF11', 'ANP32E',
             'TUBB4B', 'GTSE1', 'KIF20B', 'HJURP', 'CDCA3', 'HN1',
             'CDC20', 'TTK', 'CDC25C', 'KIF2C', 'RANGAP1', 'NCAPD2',
             'DLGAP5', 'CDCA2', 'CDCA8', 'ECT2', 'KIF23', 'HMMR',
             'AURKA', 'PSRC1', 'ANLN', 'LBR', 'CKAP5', 'CENPE',
             'CTCF', 'NEK2', 'G2E3', 'GAS2L3', 'CBX5', 'CENPA']

# 评分
sc.tl.score_genes_cell_cycle(adata, s_genes=s_genes, g2m_genes=g2m_genes)

# 查看结果
print(adata.obs[['S_score', 'G2M_score', 'phase']].head())
print(adata.obs['phase'].value_counts())

第四步：回归细胞周期效应

白话解释：有时候我们不希望细胞周期的差异影响下游分析（比如聚类和差异表达）。这时可以把周期效应"回归掉"，就像统计学中控制混杂变量一样。

技术细节：有两种回归策略。完全回归（regress out S.Score和G2M.Score）会移除所有周期相关的变异。差异回归（regress out S.Score - G2M.Score的差值）只移除增殖与非增殖的差异，保留G2M和S期之间的相对关系。

# 策略1：完全回归细胞周期效应
pbmc <- ScaleData(pbmc, vars.to.regress = c("S.Score", "G2M.Score"))

# 策略2：保留增殖细胞间差异，只去除增殖vs非增殖的效应
pbmc$CC.Difference <- pbmc$S.Score - pbmc$G2M.Score
pbmc <- ScaleData(pbmc, vars.to.regress = "CC.Difference")

# 回归后重新降维聚类
pbmc <- RunPCA(pbmc)
pbmc <- RunUMAP(pbmc, dims = 1:20)
pbmc <- FindNeighbors(pbmc, dims = 1:20)
pbmc <- FindClusters(pbmc, resolution = 0.5)

第五步：可视化细胞周期分布

library(ggplot2)

# UMAP上展示细胞周期阶段
DimPlot(pbmc, group.by = "Phase", cols = c("G1" = "grey", "S" = "blue", "G2M" = "red"))

# 各cluster中细胞周期比例
ggplot(pbmc@meta.data, aes(x = seurat_clusters, fill = Phase)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("G1" = "#E8E8E8", "S" = "#4292C6", "G2M" = "#EF3B2C")) +
  theme_minimal() +
  labs(x = "Cluster", y = "Proportion", title = "Cell Cycle Phase by Cluster")

# S分数 vs G2M分数散点图
ggplot(pbmc@meta.data, aes(x = S.Score, y = G2M.Score, color = Phase)) +
  geom_point(size = 0.5, alpha = 0.5) +
  theme_minimal() +
  scale_color_manual(values = c("G1" = "grey60", "S" = "blue", "G2M" = "red"))

实战命令速查

# 一步完成评分（Seurat）
pbmc <- CellCycleScoring(pbmc, s.features = cc.genes$s.genes, g2m.features = cc.genes$g2m.genes)

# 回归周期效应
pbmc <- ScaleData(pbmc, vars.to.regress = c("S.Score", "G2M.Score"))

# 非人物种基因转换（小鼠示例）
library(biomaRt)
convertHumanGeneList <- function(x) {
  human <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
  mouse <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")
  genes <- getLDS(attributes = c("hgnc_symbol"), filters = "hgnc_symbol",
                  values = x, mart = human, attributesL = c("mgi_symbol"), martL = mouse)
  return(unique(genes[, 2]))
}
s.genes.mm <- convertHumanGeneList(cc.genes$s.genes)
g2m.genes.mm <- convertHumanGeneList(cc.genes$g2m.genes)

面试常问点

Q1: 为什么需要进行细胞周期评分？ A: 细胞周期是单细胞数据中重要的变异来源。增殖活跃的细胞会表达大量周期相关基因，导致不同周期阶段的细胞在降维空间中分开聚类。如果研究目的不是周期本身，这种效应会干扰细胞类型的鉴定。

Q2: 什么时候应该回归细胞周期效应，什么时候不应该？ A: 当研究目的是细胞分化、细胞类型鉴定时，通常应该回归。但如果研究关注增殖（如肿瘤增殖异质性）、干细胞静息vs活化、或免疫细胞活化，则不应回归，因为周期状态本身就是重要的生物学信号。

Q3: Tirosh基因集能用于非人物种吗？ A: 不能直接使用，需要进行同源基因转换。小鼠可用biomaRt或直接手动对应；对于非模式生物，需确认直系同源基因的存在和命名。

Q4: 细胞周期评分与增殖指数（如MKI67表达）有什么区别？ A: MKI67只是一个标记基因，在G2/M期高表达，灵敏度有限。细胞周期评分使用几十个基因的综合信号，更鲁棒也更准确，能区分S期和G2M期。

Q5: 如果大部分周期基因在数据中检测不到怎么办？ A: 这可能是测序深度不够或该组织/细胞类型的表达特征决定的。可以降低nbin参数、使用更小的基因集，或考虑使用专门为低深度数据设计的方法如Cyclone（scran包）。

易错点

  1. 忘记标准化：CellCycleScoring需要在NormalizeData之后使用，未标准化的原始counts会导致评分错误
  2. 物种基因名不匹配：人类基因名用于小鼠数据会导致大部分基因无法匹配，评分失效
  3. 过度回归：对所有数据无脑回归周期效应可能丢失重要的增殖相关生物学信号
  4. 基因集不在数据中：部分基因可能被质控过滤掉，应检查实际匹配的基因数量
  5. SCTransform与CellCycleScoring的顺序：如果使用SCTransform，应在SCTransform之前先用NormalizeData做评分

补充知识

替代方法对比

方法包原理优势
CellCycleScoringSeurat基因集均值评分简单快速
cyclonescran训练对比对模型对低深度数据更鲁棒
tricycleBioconductor连续周期角度估计提供连续值而非离散分类
ccRemoverccRemoverPCA方法直接在表达矩阵上操作

tricycle提供连续的细胞周期角度

library(tricycle)
# 计算连续的细胞周期角度（theta, 0-2π）
sce <- project_cycle_space(sce)
sce <- estimate_cycle_position(sce)
# theta值对应：0=G1/S, π/2=S, π=G2M, 3π/2=M/G1

这种连续估计在分析细胞周期动态变化时比离散的G1/S/G2M分类更精细。

📚 相关文章推荐

单细胞衰老评分
25% 相关
379_单细胞代谢通路活性评分
25% 相关
376_单细胞数据标准化流程
24% 相关
183_RNA-seq批次效应校正
21% 相关
单细胞数据整合方法
20% 相关