1、什么是细胞周期
大家在生物课上应该都学过细胞周期的相关知识，这里带大家简单回顾一下。细胞周期（cell cycle）是指细胞从一次分裂完成开始到下一次分裂结束所经历的全过程，分为间期与分裂期两个阶段。
间期又分为G1，S和G2期。G1期是DNA合成准备期，主要合成RNA和核糖体，S期是DNA合成期，合成DNA，组蛋白和DNA复制所需要的酶，G2期是分裂准备期，DNA合成终止，大量合成RNA及蛋白质，包括微管蛋白和促成熟因子等。
分裂期（M期）是细胞形态结构发生急速变化的时期，包括一系列核的变化、染色质的浓缩、纺锤体的出现，以及染色体精确均等地分配到两个子细胞中的过程。蛋白质合成明显下降，RNA合成及其他代谢周转停止。
2、细胞周期对 scRNA 数据的影响
很多研究都发现[1]，相同类型的细胞在细胞周期中不是同步的，可能处于不同的细胞周期阶段，这会导致相同类型细胞的异质性。不同细胞类型的细胞可能会因为处于相同的细胞周期阶段，而具有相同的基因表达特征，在进行 scRNA 细胞聚类分析时聚到同一个类中，但实际上它们是属于不同的细胞类型，如果进行亚群分析，可以将它们再分为不同的细胞类型。
下图为 ccRemover[1]使用的一个模拟数据集的基因表达水平的 PCA 降维图。这个数据集有 50 个细胞和2000个基因，不同的形状表示不同细胞类型，不同颜色表示不同的细胞周期阶段，可以看到细胞按细胞类型和细胞周期阶段分成了 6 个类。
3 细胞周期 marker 基因
细胞繁殖的 marker 基因包括Cdk1，Foxm1和Pbk等，G1/S 期的典型marker基因包括Cdk2,Ccne1, Mcm 和 E2f 家族成员等，G2/M 期的典型 marker基因包括Cdk1, *Ccnb *和纺锤体形成相关的基因等[2]。当在关键分析结果中发现这些基因时，就需要注意检查自己的数据是否受细胞周期影响啦。
下图为Seurat [3]内置的小鼠细胞周期相关基因，可以使用 cc.genes 命令进行查看。
4 推断细胞周期
4.1 seurat 包的CellCycleScoring函数
使用 seurat 函数CellCycleScoring可以预测每个细胞所处的细胞周期阶段。其计算原理为：
- 1、 对每个细胞计算S期基因和 G2/M期基因的平均表达值；
- 2、把所有基因基于平均表达水平划分到不同 bin 里去，然后从每个 bin 中随机抽取作为背景的 control 基因，计算这些 control 基因的平均表达水平；
- 3、从 S 期基因和 G2/M 期基因的平均表达水平中减去control 基因的平均表达水平，得到 S.Score和 G2M.Score；
- 4、 S.Score <0 且 G2M.Score <0 判断为 G1 期，否则，哪个分数大就判断为哪个期。
分析代码如下：
# 创建一个 Seurat 对象，counts为单细胞基因表达矩阵，行名为基因，列名为细胞barcode，每一个值为该细胞中该基因的UMI count数
seuratObj <- CreateSeuratObject(counts = counts, project = "seurat")
# 对表达量进行标准化
seuratObj <- NormalizeData(object = seuratObj, verbose = FALSE)
# 计算细胞周期分数
seuratObj <- CellCycleScoring(seuratObj, s.features = s.genes, g2m.features = g2m.genes)
s.genes 为一个S期基因名的向量，g2m.genes 为G2/M期基因的向量。如果物种是小鼠，可以用seurat内置的cc.genesg2m.genes。也可以使用从文献中找到的细胞周期相关基因进行计算。
计算后得到的S.Score，G2M.Score和推断出的细胞周期结果如下图所示。Phase就是由细胞周期分数推断得到的该细胞所处的细胞周期阶段。
4.2 scran 包的 cyclone 函数
使用 scran [4]函数cyclone也可以预测每个细胞所处的细胞周期阶段。
分析代码如下：
# 创建一个 SingleCellExperiment 对象，输入的counts为单细胞基因表达矩阵，行名为基因，列名为细胞barcode，每一个值为该细胞中该基因的UMI count数
sce <- SingleCellExperiment(assays = list(counts = counts))
# 计算细胞周期分数
result <- cyclone(sce, pairs= ccGenes, gene.names= rownames(sce))
ccGenes 是细胞周期 marker 基因，在 scran 包的 exdata 文件夹中内置了人和小鼠的细胞周期 marker 基因，可以通过以下方法读取小鼠的细胞周期 marker 基因。如果想读取人的细胞周期marker基因，只需要将mouse替换成human就行。
ccGenes <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", package="scran"))
scran 内置的G1期基因：
5 消除细胞周期影响的方法
5.1 消除细胞周期效应
- seurat regress 细胞周期分数
如果采用seurat进行分析，在 ScaleData 或 SCTransform 中，可以对通过 seurat CellCycleScoring 计算出的细胞周期分数进行 regress，从而消除细胞周期对细胞分类的影响。
分析代码如下：
seuratObj <- ScaleData(seuratObj, vars.to.regress = c("S.Score", "G2M.Score"), verbose = FALSE)
# 或
seuratObj <- SCTransform(seuratObj, vars.to.regress =c("S.Score", "G2M.Score"), verbose = FALSE)
5.1.2 ccRemover
ccRemover是一个用于去除单细胞数据中的细胞周期效应的R包。其原理为，计算细胞周期相关基因的在 PCA 降维的结果中对各 PC 的贡献分数，如果细胞周期基因的贡献分数平均值大于背景基因，则判断该 PC 为细胞周期的影响，从总的 PC 中减去这个 PC，最后返回去掉细胞周期效应的表达矩阵。
分析代码如下：
# 获取每个基因的平均表达值
mean_gene_exp <- rowMeans(as.matrix(exp))
# 对每个基因的表达值进行中心化
exp_center <- exp - mean_gene_exp
# 获取基因名
gene_names <- rownames(exp_center)
# 利用包内置的细胞周期相关基因，从我们的基因中提取出细胞周期相关基因
cell_cycle_gene_indices <- gene_indexer(gene_names, species = "mouse", name_type = "symbols" )
# 查看我们的基因中细胞周期相关基因的个数
length(cell_cycle_gene_indices)
# 生成一个向量，里面存放每个基因是否是细胞周期相关基因的信息，如果是，则为 TRUE，否则为 FALSE
is_cc_gene <- rep(FALSE, nrow(exp_center))
is_cc_gene[cell_cycle_gene_indices] <- TRUE
summary(is_cc_gene)
# 把表达矩阵和 is_cc_gene 放到一起
dat <- list(x=exp_center, if_cc=is_cc_gene)
# 执行 ccRemover
xhat <- ccRemover(dat, bar=FALSE)
# 把平均值加回来
xhat <- xhat + mean_gene_exp
5.2 去掉受细胞周期影响严重的细胞
在有些文献[2]的分析过程中，会去掉不感兴趣的，或者是对细胞分类产生混淆的细胞群。如果对细胞周期的研究不感兴趣，并且样本中只有少部分细胞受到细胞周期影响而聚成一类，可以考虑去除这个细胞群体。一般，大部分细胞会处于G1期，也可以采用去掉 G1 期以外的细胞的方法。
5.3 去掉细胞周期基因
在进行细胞分类前，去掉在各细胞间表达量差异较大的细胞周期基因，在细胞分类时，不同类型的细胞就不会因为，在这些高差异的细胞周期基因中表达模式相似，而聚在一起了。
例如有些文献[5]中采用层次聚类的方法，从所有细胞的 top 2000 差异基因中找到细胞周期相关基因，去掉这些基因，然后再进行PCA降维，并且在后续分析（如差异分析）结果中也去掉这些基因。
6 总结
细胞周期对不同样本的影响程度不一定是一样的，针对不同数据的特点应该考虑采用不同的方法。
比如，有的数据中细胞周期效应不那么明显，对细胞分类的影响不大，可以把细胞周期当成一种批次效应，在分析时消除就行。
有的数据受细胞周期影响很大，从数据中去掉所有细胞周期基因，也不一定能消除细胞周期效应，因为还有很多潜在的受细胞周期影响的基因可能会影响细胞分类，此时可能需要综合考虑细胞的整体状态，是否需要将细胞周期效应作为研究的一部分，而不是进行消除。
参考文献：
[1] Barron, M. and Li, J. Identifying and removing the cell-cycle effect from single-cell RNA-Sequencing data. Sci. Rep. 6, 33892; doi: 10.1038/srep33892 (2016).
[2] Wei-Lin Qiu, Yu-Wei Zhang, Ye Feng et al. Deciphering Pancreatic Islet β Cell and α Cell Maturation Pathways and Characteristic Features at the Single-Cell Level. Cell Metabolism 25, 1194–1205.e1–e4, May 2, 2017.
[3] Butler, A., Hoffman, P., Smibert, P. et al. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol 36, 411–420 (2018).
[4] Lun ATL, McCarthy DJ, Marioni JC (2016). “A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor.” F1000Res., 5, 2122.
[5] Ye Feng, Wei-Lin Qiu, Xin-Xin Yu et al. Characterizing pancreatic b-cell heterogeneity in the streptozotocin model by single-cell transcriptomic analysis. Molecular Metabolism 37 (2020) 100982.
文章转自微信公众号：嘉因生物
原创作者： 冲啊蜗牛
网友评论