Seurat亮点之细胞周期评分和回归
作者：张虎
校对：生信宝典斯坦福大学Satija lab的 Seurat v3.1 guidline于近日更新啦！其中包括许多个性化的模块，其中我个人比较感兴趣的是Cell-Cycle Scoring and Regression模块,因为在条件干预的情况下，部分细胞处于非稳定状态下，如增殖类细胞出现由于细胞周期相关基因的不同导致细胞聚类发生一定的偏移。
其实在许多已经发表的文献中，我们也可以看到由于细胞周期不同所导致的分类偏移。如在Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing 中 vCAFs和cCAFs中只有cell cycle genes的表达差异。
当作者使用SC3重新进行聚类时，两种亚型又重新聚在了一起 (重点关注图中的cluster1和cluster3，不要被颜色误导了。同一个cluster，不同的配色方案，也是服了，还是看看史上最全的图表色彩运用原理吧)。
系统介绍
通过计算各个细胞可能所处的细胞周期阶段的得分 (G1,S,G2期），并在预处理过程 (主要是ScaleData步骤)中将细胞周期得分作为混杂因素移除，从而排除单细胞所处细胞周期不同对基因表达和细胞分型的影响。作者在小鼠造血祖细胞的数据集上证明了该观点 （Nestorowa et al. Blood 2016.），其实只要是Seurat v3对象，自己的数据都是可以跑得通的。
细胞周期相关基因集使用的是人的基因，用小鼠进行测试，说明该细胞周期相关基因数据集适合人和小鼠；如果是其它物种，准备方法见 https://github.com/satijalab/seurat/issues/462。
下面是操作代码：
library(Seurat)
# 读取表达矩阵， The first row is a header row, the first column is rownames
exp.mat <- read.table(file = "../data/nestorawa_forcellcycle_expressionMatrix.txt",
header = TRUE, as.is = TRUE, row.names = 1)
# 一系列的细胞周期相关的markers,其中包括处于S期的43个细胞周期相关基因，54个G2M期的细胞周期相关基因，
# from Tirosh et al, 2015, is loaded with Seurat. We can
# segregate this list into markers of G2/M phase and markers of S phase
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes
# 创建Seurat对象并进行标准化；marrow <- CreateSeuratObject(counts = exp.mat)
marrow <- NormalizeData(marrow)
marrow <- FindVariableFeatures(marrow, selection.method = "vst")
marrow <- ScaleData(marrow, features = rownames(marrow))如果我们在Seurat对象上进行PCA分析，使用FindVariableFeatures中找到高可变基因，进行PCA分析，并展示对各个主成分贡献最大的基因。
PCA分析有不少需要注意的，具体见用了这么多年的PCA可视化竟然是错的！！！ 和 PCA主成分分析实战和可视化 | 附R代码和测试数据。
marrow <- RunPCA(marrow, features = VariableFeatures(marrow), ndims.print = 6:10, nfeatures.print = 10)我们看到对PC8和PC10贡献最大的基因中一部分是细胞周期相关基因，包括TOP2A和MKI67。
从下面的热图也可以看出来。
DimHeatmap(marrow, dims = c(8, 10))#热图表示
我们将尝试从数据中去除该组分，从而确保细胞周期异质性对PCA或下游分析没有贡献。
分配细胞周期分数
首先，根据其G2/M和S期标记基因的表达为每个细胞分配一个所处周期的分数。这些标记基因的表达水平应该是反相关的，而不表达这些标记基因的细胞可能处于G1期。
我们用CellCycleScoring函数计算细胞周期分数，并在metadata中存储S和G2/M分数，以及G2M，S或G1阶段中每个细胞的预测分类。如果指定set.ident=T，则CellCycleScoring将Seurat对象中每个细胞的分组信息设置为其所处的细胞周期阶段。
marrow <- CellCycleScoring(marrow, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
# view cell cycle scores and phase assignments
head(marrow[[]])
观察细胞周期基因的表达情况 （评估计算的准确性）
注：R语言可视化学习笔记之ggridges包可绘制类似图形。
# 观察细胞周期基因的表达情况
RidgePlot(marrow, features = c("PCNA", "TOP2A", "MCM6", "MKI67"), ncol = 2)
我们用CellCycleScoring函数计算细胞周期分数，并在metadata中存储S和G2/M分数，以及G2M，S或G1阶段中每个细胞的预测分类。如果指定set.ident=T，则CellCycleScoring将Seurat对象中每个细胞的分组信息设置为其所处的细胞周期阶段。利用细胞周期基因进行PCA分析 （这个例子可以拓展，使用任意指定的基因集进行细胞周期分析）
PCA分析有不少需要注意的，具体见用了这么多年的PCA可视化竟然是错的！！！ 和 PCA主成分分析实战和可视化 | 附R代码和测试数据。
# Running a PCA on cell cycle genes reveals, unsurprisingly, that cells separate entirely by
# phase
marrow <- RunPCA(marrow, features = c(s.genes, g2m.genes))
DimPlot(marrow)
在数据标归一化时去除细胞周期影响
marrow <- ScaleData(marrow, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(marrow))再次做PCA时，就看不到细胞周期相关基因对主成分的贡献了。
# 我们可以看到组分中不再出现细胞周期相关的基因
marrow <- RunPCA(marrow, features = VariableFeatures(marrow), nfeatures.print = 10)
再次根据细胞周期相关基因进行PCA分析时，也不分不出群了，说明移除细胞周期影响的效果还是比较好的。
# When running a PCA on only cell cycle genes, cells no longer separate by cell-cycle phase
marrow <- RunPCA(marrow, features = c(s.genes, g2m.genes))
DimPlot(marrow)
如果细胞周期不合适时怎么办？
上述过程去除了与细胞周期相关的所有信息。在某些情况下，我们发现这会对下游分析产生负面影响，特别是在分化过程（如小鼠血细胞分化生成过程中 hematopoiesis）中，干细胞处于静止状态，分化细胞正在增殖（反之亦然）。在这种情况下，消除所有细胞周期效应也会模糊干细胞和前体细胞之间的区别。
作为替代方案，我们建议消除G2M和S期分数之间的差异。这意味着将保持非周期细胞和周期细胞的组分差异，但是增殖细胞之间的细胞周期阶段的差异将从数据中去除。
# 计算并移除分数差异
marrow$CC.Difference <- marrow$S.Score - marrow$G2M.Score
marrow <- ScaleData(marrow, vars.to.regress = "CC.Difference", features = rownames(marrow))细胞周期相关的基因对各个主成分贡献减小
# cell cycle effects strongly mitigated in PCA
# 在PCA中不再出现大量的细胞周期相关的基因
marrow <- RunPCA(marrow, features = VariableFeatures(marrow), nfeatures.print = 10)
G1期区分开，G2/M和S期聚在一起。
# when running a PCA on cell cycle genes, actively proliferating cells remain distinct from G1
# cells.
# however, within actively proliferating cells, G2M and S phase cells group together
marrow <- RunPCA(marrow, features = c(s.genes, g2m.genes))
DimPlot(marrow)
单细胞是目前很火的领域，分析工具很多，而且也还在不断发展中。Seurat是其中一个，虽然好用，却不一定是最好的。而且运用好工具，需要对原理有很好的理解。