如果讀過一些單細胞的文獻，應該會經常看到一群名為"cycling cell"的亞群，T細胞、B細胞、上皮細胞等等在分亞群的時候，都可能碰到。實際上，這群細胞只是因為高表達一些與細胞周期相關的基因，才會被單獨聚成一群，里面可能包含多種細胞亞群的混合，只要它們處於細胞周期/增殖狀態。
當這群細胞占比很少時，可以不做處理；占比較大，或者處於不同細胞周期時間點(G1, S, G2, M)的細胞相互分開，這時就需要消除細胞周期的影響，留下有意義的生物學差異。下圖是比較極端的例子，所有細胞都處於細胞周期，且不同時間點的細胞分開：
這張圖中的細胞是什么類型，有幾類亞群暫且不論，僅僅是細胞周期相關基因主導聚類，干擾聚類結果，就挺影響后續分析的。這個例子來自Seurat官網，我也借此例子進行具體演示，我改了其中一些流程代碼，希望能說明得清楚一些，后台回復20210316可獲取全部代碼和測試數據集。
總的來說，Seurat中的函數可以根據經典的細胞周期基因計算周期得分，推斷細胞所處的cell-cycle時間點，然后在預處理過程中將其回歸掉。
先把標准流程跑下來，
library(Seurat)
library(tidyverse)
test.mat=read.table("nestorawa_forcellcycle_expressionMatrix.remove_some_lines.txt",header = T,row.names = 1,sep = "\t",stringsAsFactors = F)
test.seu=CreateSeuratObject(counts = test.mat)
test.seu <- NormalizeData(test.seu)
test.seu <- FindVariableFeatures(test.seu, selection.method = "vst")
test.seu <- ScaleData(test.seu, features = rownames(test.seu))
test.seu <- RunPCA(test.seu, npcs = 50, verbose = FALSE)
test.seu <- FindNeighbors(test.seu, dims = 1:30)
test.seu <- FindClusters(test.seu, resolution = 0.4)
test.seu <- RunUMAP(test.seu, dims = 1:30)
test.seu <- RunTSNE(test.seu, dims = 1:30)
DimPlot(test.seu,reduction = "umap")
FeaturePlot(test.seu,features = c("PCNA","MKI67"),reduction = "umap")
s期基因PCNA，G2/M期基因MKI67均高表達，說明細胞均處於細胞周期
接下來判斷這些細胞分別處於什么期，用到Seurat的CellCycleScoring()函數，依據的基因集已內置在Seurat中
s.genes=Seurat::cc.genes.updated.2019$s.genes
g2m.genes=Seurat::cc.genes.updated.2019$g2m.genes
test.seu <- CellCycleScoring(test.seu, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
set.ident參數是將時間點的推斷結果賦值為每個細胞的身份，之前聚類之后每個細胞的身份是seurat_clusters
這一步之后test.seu@meta.data數據框會多4列：S.Score、G2M.Score、Phase、old.ident
> head(test.seu@meta.data,2)
orig.ident nCount_RNA nFeature_RNA RNA_snn_res.0.4 seurat_clusters S.Score
Prog_013 Prog 2562453 10206 0 0 -0.1230880
Prog_019 Prog 3030553 9987 3 3 -0.1407003
G2M.Score Phase old.ident
Prog_013 -0.4613545 G1 0
Prog_019 0.5763395 G2M 3
我們來看一下兩個評分和細胞周期phase的關系
test.seu@meta.data %>% ggplot(aes(S.Score,G2M.Score))+geom_point(aes(color=Phase))+
theme_minimal()
從這個圖中，我們可以看出S.Score較高的為S期，G2M.Score較高的為G2M期，都比較低的為G1期
再回過頭來看看前面得到的降維圖
DimPlot(test.seu,reduction = "umap")
右上角混得還行，左下角這些基本就是細胞周期基因表達譜引起的聚類了，這就是我們要消除的影響
Seurat主要是在數據縮放過程中回歸細胞周期分數，用到的仍然是ScaleData()函數，這一步中，對於每個基因，Seurat建模了基因表達值與S、G2M細胞周期分數之間的關系。模型的殘差表示一個矯正后的表達矩陣，用於下游降維。
test.seu <- ScaleData(test.seu, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(test.seu))
#耗時較長
之后重新跑標准流程
test.seu <- RunPCA(test.seu, npcs = 50, verbose = FALSE)
test.seu <- FindNeighbors(test.seu, dims = 1:30)
test.seu <- FindClusters(test.seu, resolution = 0.4)
test.seu <- RunUMAP(test.seu, dims = 1:30)
test.seu <- RunTSNE(test.seu, dims = 1:30)
DimPlot(test.seu,reduction = "umap",group.by = "Phase")
將新圖和之前的圖對比可以發現，此時細胞周期對聚類結果的影響已經減少了，不同時間點的細胞混合得更均勻，但也不是100%，沒有工具可以保證這個效果
因水平有限，有錯誤的地方，歡迎批評指正！