Seurat之细胞周期评分

1. 细胞周期背景介绍

细胞周期(cell cycle)是指细胞从一次分裂完成开始到下一次分裂结束所经历的全过程，分为间期与分裂期两个阶段。在细胞周期过程中，细胞的遗传物质复制并均等地分配给两个子细胞。
（一） 间期
间期又分为三期、即DNA合成前期（G1期）、DNA合成期（S期）与DNA合成后期（G2期）。

1. G1期（first gap）从有丝分裂到DNA复制前的一段时期，又称合成前期，此期主要合成RNA和核糖体。该期特点是物质代谢活跃，迅速合成RNA和蛋白质，细胞体积显著增大。这一期的主要意义在于为下阶段S期的DNA复制作好物质和能量的准备。细胞进入G1期后，并不是毫无例外地都进入下一期继续增殖，在此时可能会出现三种不同前景的细胞：①增殖细胞：这种细胞能及时从G1期进入S期，并保持旺盛的分裂能力。例如消化道上皮细胞及骨髓细胞等；②暂不增殖细胞或休止细胞：这类细胞进入G1期后不立即转入S期，在需要时，如损伤、手术等，才进入S期继续增殖。例如肝细胞及肾小管上皮细胞等；③不增殖细胞：此种细胞进入G1期后，失去分裂能力，终身处于G1期，最后通过分化、衰老直至死亡。例如高度分化的神经细胞、肌细胞及成熟的红细胞等。

2. S期（synthesis）即DNA合成期。在此期，除了合成DNA外，同时还要合成组蛋白。DNA复制所需要的酶都在这一时期合成。

3. G2期（second gap）为DNA合成后期，是有丝分裂的准备期。在这一时期，DNA合成终止，大量合成RNA及蛋白质，包括微管蛋白和促成熟因子等。

（二）细胞分裂期，也就是M期

细胞的有丝分裂（mitosis）期又分为前期，中期，后期，末期四个阶段，是一个连续变化过程，由一个母细胞分裂成为两个子细胞。一般需1～2小时。

1. 前期（prophase）染色质丝高度螺旋化，逐渐形成染色体（chromosome）。染色体短而粗，强嗜碱性。两个中心体向相反方向移动，在细胞中形成两极；而后以中心粒随体为起始点开始合成微管，形成纺锤体。随着核仁相随染色质的螺旋化，核仁逐渐消失。核被膜开始瓦解为离散的囊泡状内质网。

2. 中期（metaphase）细胞变为球形，核仁与核被膜已完全消失。染色体均移到细胞的赤道平面，从纺锤体两极发出的微管附着于每一个染色体的着丝点上。从中期细胞可分离得到完整的染色体群，共46个，其中44个为常染色体，2个为性染色体。男性的染色体组型为44+XY，女性为44+XX。分离的染色体呈短粗棒状或发夹状，均由两个染色单体借狭窄的着丝点连接构成。

3．后期（anaphase）由于纺锤体微管的活动，着丝点纵裂，每一染色体的两个染色单体分开，并向相反方向移动，接近各自的中心体，染色单体遂分为两组。与此同时，细胞被拉长，并由于赤道部细胞膜下方环行微丝束的活动，该部缩窄，细胞遂呈哑 铃形。

4．末期（telophase）染色单体逐渐解螺旋，重新出现染色质丝与核仁；内质网囊泡组合为核被膜；细胞赤道部缩窄加深，最后完全分裂为两个2倍体的子细胞。

（三）G0期
暂时离开细胞周期，停止细胞分裂，去执行一定生物学功能的细胞所处的时期。

参考：https://baike.so.com/doc/2392736-2529977.html

2. Seurat对细胞周期进行评分

参考：https://satijalab.org/seurat/archive/v3.1/cell_cycle_vignette.html

原理：首先，我们基于G2/M和S期的经典marker基因的表达，计算每个细胞可能所处的细胞周期的分数。这些marker基因集应该与它们的表达水平是成反相关的关系，而那些都不表达这些marker基因的细胞可能处于G1期。
我们使用CellCycleScoring函数计算每个细胞的细胞周期得分，并将计算出的S期和G2/M期的评分保存在metadata中，以及细胞处于G2M，S或G1期的预测分类。通过设置set.ident = TRUE，则CellCycleScoring将Seurat对象中每个细胞的分组信息设置为其所处的细胞周期阶段。

示例数据下载：pbmc3k（注 ：pbmc多为分化成熟的细胞，受细胞周期影响较小。该数据仅用于演示，以熟悉细胞周期分析流程。）

• 读入数据，创建seurat对象，并进行质控和标准化

pbmc <- Read10X('./filtered_gene_bc_matrices/hg19/')
pbmc <- CreateSeuratObject(pbmc,project = 'pbmc3k',min.cells = 3,min.features = 200)
#质控
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
# 标准化
pbmc <- NormalizeData(pbmc)

• 细胞周期评分
  cc.genes是Seurat自带的一个包含两个vectors的list。
  这里包含的s.genes和g2m.genes也可按需求换成自己想要的基因，然后进行CellCycleSorting

View(cc.genes)

pbmc <- NormalizeData(pbmc)
g2m_genes <- cc.genes$g2m.genes ## 获取G2M期marker基因
g2m_genes <- CaseMatch(search=g2m_genes, match=rownames(pbmc)) #提取pbmc矩阵中的G2M期marker基因
s_genes <- cc.genes$s.genes   #获取S期marker基因 
s_genes <- CaseMatch(search=s_genes, match=rownames(pbmc)) #提取pbmc矩阵中的S期marker基因 
#通过提取到的g2m期基因和s期基因，使用CellCycleScoring函数，对pbmc进行细胞周期评分
pbmc <- CellCycleScoring(pbmc, g2m.features=g2m_genes, s.features=s_genes)

⚠️细胞周期评分用的是Normalize Data（经log转换后的矩阵）

查看细胞周期评分结果
在metadata中出现了"S.Score" "G2M.Score" 和"Phase" （CellCycleScoring的结果）
预测结果处于G1期细胞有1174个，G2M期有474个，S期有990个。

colnames(pbmc@meta.data)
[1] "orig.ident"   "nCount_RNA"   "nFeature_RNA" "percent.mt"  "S.Score"  "G2M.Score"  
[7]   "Phase"  
table(pbmc$Phase)
  G1  G2M    S 
1174  474  990 

• 可视化

# Visualize the distribution of cell cycle markers across 
RidgePlot(pbmc, features = c("MCM4", "TYMS", "MCM5", "MCM2"), ncol = 2) 
#查看的这几个基因表达量较低，可以更换其他基因查看

# Running a PCA on cell cycle genes reveals, unsurprisingly, that cells separate entirely by phase 
pbmc <- RunPCA(pbmc, features = c(s_genes, g2m_genes)) 
# 数据可视化，可以看到细胞按不同的细胞周期进行了分群 
DimPlot(pbmc,group.by = 'Phase')

可以看到，该组数据受细胞周期影响较小（因为是pbmc的演示数据）

若数据受细胞周期影响较大，pca图三个时期的重合度会较差。
如：

这时就要考虑消除细胞周期对后续分析的影响。

3. 消除细胞周期对数据的影响

方法1：使用ScaleData函数进行数据标准化，并设置vars.to.regress参数指定对细胞周期评分进行回归处理，消除细胞周期异质性的影响

pbmc <- ScaleData(pbmc, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(pbmc))

pbmc <- RunPCA(pbmc, features = VariableFeatures(pbmc))
pbmc <- RunPCA(pbmc, features = c(s_genes, g2m_genes))
DimPlot(pbmc,group.by = 'Phase')

方法2：可选择的替代方法
在上述的分析过程中，我们消除了与细胞周期相关的所有信号。但是，在某些情况下，我们发现这会对下游的分析产生一定的负面影响，尤其是在细胞分化过程中（如鼠类造血过程）。在此过程中干细胞处于静止状态，而分化的细胞正在增殖。在这种情况下，清除所有细胞周期效应也会使干细胞和祖细胞之间的区别模糊。作为替代方案，建议可以逐步消除G2M和S期评分之间的差异。这意味着将保持非周期细胞和周期细胞的组分差异，但是增殖细胞之间的细胞周期阶段的差异将从数据中去除。

# 计算分数差异 
pbmc$CC.Difference <- pbmc$S.Score - pbmc$G2M.Score 
# 数据标准化消除分数差异 
pbmc <- ScaleData(pbmc, vars.to.regress = "CC.Difference", features = rownames(pbmc)) 
# cell cycle effects strongly mitigated in PCA 
pbmc <- RunPCA(pbmc, features = VariableFeatures(pbmc), nfeatures.print = 10) 

pbmc <- RunPCA(pbmc, features = c(s_genes, g2m_genes)) 
DimPlot(pbmc,group.by = 'Phase') 

方法3：SCTransform

pbmc <- SCTransform(pbmc, vars.to.regress = c("S.Score","G2M.Score"))
pbmc <- RunPCA(pbmc, features = c(s_genes, g2m_genes)) 
DimPlot(pbmc,group.by = 'Phase') 

⚠️vars.to.regress可以用于消除细胞周期、线粒体含量等的影响，但一般不用于消除批次影响。（效果有限，且有专门的，更好的去批次方法）