1 前言
有时我们希望从scRNA数据中获得细胞周期的信息，一般来讲，细胞在各个周期的分布可能与真实情况还是存在差距，但是可以用这种方法看看不同亚群之间或者各个处理之间的扩增差异。许多细胞周期相关的关键节点（比如通过检查点）属于转录后调控，在转录组数据中不可见，但通过表达量的变化也可以进行初步的推断。
另外还是先要放一张细胞周期的图片
下面就用416B数据来看看
load('../0-data/clusted.sce.416b.RData')
sce.416b
## class: SingleCellExperiment
## dim: 46604 185
## metadata(0):
## assays(3): counts logcounts corrected
## rownames(46604): 4933401J01Rik Gm26206 ... CAAA01147332.1
## CBFB-MYH11-mcherry
## rowData names(4): Length ENSEMBL SYMBOL SEQNAME
## colnames(185): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
## colData names(11): Source Name cell line ... sizeFactor label
## reducedDimNames(2): PCA TSNE
## altExpNames(2): ERCC SIRV
2 使用细胞周期蛋白
2.1 先找细胞周期蛋白基因，再对照数据集
细胞周期蛋白（cyclins）控制着整个细胞周期的进展，并在细胞周期各阶段的表达模式具有一定特点。例如：
cyclin B在G2后期和M（mitosis有丝分裂期）表达量最高
cyclin D在整个过程都表达，但在G1期表达最高；
根据相关基因在各个cluster的表达可以帮助判断：
cyclin.genes <- grep("^Ccn[abde][0-9]$", rowData(sce.416b)$SYMBOL)
cyclin.genes <- rownames(sce.416b)[cyclin.genes]
cyclin.genes
## [1] "Ccnb3" "Ccna2" "Ccna1" "Ccne2" "Ccnd2" "Ccne1" "Ccnd1" "Ccnb2" "Ccnb1"
## [10] "Ccnd3"
画个热图
library(scater)
plotHeatmap(sce.416b, order_columns_by="label",
cluster_rows=FALSE, features=sort(cyclin.genes))
图中可以看到：cluster1很明显ccnd系列的基因表达量高，说明它可能是G1期，而其他几个clusters分散在后面的时期
2.2 先找marker基因，再对应细胞周期蛋白
不管怎样，416b数据的对应关系还是很不错的，界限清晰，但并非任何数据都像它一样，尤其是一些异质性较高的数据中，细胞周期可能不是驱动基因表达变化的主力。
不过好在，即使不非常清楚地知道每个cluster的每个细胞所处的周期，也还是能回答一下细胞周期相关的问题。例如：我们可以利用常规的寻找差异基因的方法，找找有没有基因在细胞周期蛋白中上调，来帮助判断某个cluster中是否有更多细胞属于某个细胞周期。
首先寻找marker基因
library(scran)
markers <- findMarkers(sce.416b, subset.row=cyclin.genes,
test.type="wilcox", direction="up")
# 4个cluster都各自找了细胞周期蛋白基因的排名
> markers
List of length 4
names(4): 1 2 3 4
然后对应细胞周期蛋白
比如看cluster4，它的Ccnb系列基因AUC指标高于其余的基因，说明在cyclin B中表达量较高，存在更多的细胞处于G2/M期
这种直接通过表达量判断的方法很直观，也很好理解。但需要保证得到每个细胞周期蛋白基因表达量，而且需要一个假设：细胞周期蛋白基因的表达不受其他生物因素的影响，尤其对于那些非常容易混淆正常细胞周期的事件，如恶性癌细胞的扩增。
3 辅以参考数据
就是使用别人已经做好的细胞周期注释数据，来辅助我们自己的数据进行注释【有没有很像之前细胞类型注释的环节？】加入参考注释的目的是：更精确地寻找细胞周期相关的基因
这里使用的参考数据是来自 Buettner et al. 2015 的小鼠胚胎干细胞结果
library(scRNAseq)
sce.ref <- BuettnerESCData()
sce.ref
## class: SingleCellExperiment
## dim: 38293 288
## metadata(0):
## assays(1): counts
## rownames(38293): ENSMUSG00000000001 ENSMUSG00000000003 ...
## ENSMUSG00000097934 ENSMUSG00000097935
## rowData names(3): EnsemblTranscriptID AssociatedGeneName GeneLength
## colnames(288): G1_cell1_count G1_cell2_count ... G2M_cell95_count
## G2M_cell96_count
## colData names(1): phase
## reducedDimNames(0):
## altExpNames(1): ERCC
> table(sce.ref$phase)
G1 G2M S
96 96 96
为了提高判断的精准度，可以除了参考数据集的注释结果外，再加上GO中与细胞周期相关的基因，再与我们自己的基因取个交集
# 细胞周期通路（https://www.ebi.ac.uk/QuickGO/term/GO:0007049）基因
library(org.Mm.eg.db)
cycle.anno <- select(org.Mm.eg.db, keytype="GOALL", keys="GO:0007049",
columns="ENSEMBL")[,"ENSEMBL"]
# 再取交集
candidates <- Reduce(intersect,
list(rownames(sce.ref), rowData(sce.416b)$ENSEMBL, cycle.anno))
str(candidates)
## chr [1:1605] "ENSMUSG00000000001" "ENSMUSG00000000028" ...
用三者交集的基因，代入到注释好的参考数据集中，继续找marker基因（相当于再过滤一遍）
sce.ref <- logNormCounts(sce.ref)
phase.stats <- pairwiseWilcox(logcounts(sce.ref), sce.ref$phase,
direction="up", subset.row=candidates)
cycle.markers <- getTopMarkers(phase.stats[[1]], phase.stats[[2]])
将得到的marker基因给到416b数据
# 由于参考数据中的基因ID是Ensembl，所以这里也使用Ensembl
test.data <- logcounts(sce.416b)
rownames(test.data) <- rowData(sce.416b)$ENSEMBL
library(SingleR)
assignments <- SingleR(test.data, ref=sce.ref,
label=sce.ref$phase, genes=cycle.markers)
tab <- table(assignments$labels, colLabels(sce.416b))
tab
##
## 1 2 3 4
## G1 71 7 19 1
## G2M 2 60 1 13
## S 5 2 4 0
# 可以再添加一步检验，看看两个cluster之间细胞分布的差异大小
chisq.test(tab[,1:2])
##
## Pearson's Chi-squared test
##
## data: tab[, 1:2]
## X-squared = 107.91, df = 2, p-value < 2.2e-16
可以看到，cluster1大部分的细胞都在G1期，和之前利用细胞周期蛋白得到的结果差不多。但和之前的方法不同的是，这次我们是清楚看到有多少细胞处于哪个时期，而不是通过图去观察
4 用已发表的分类器
例如Scialdone et al. (2015)发表了一个cyclone的分类器，属于scran包，其中也是包括了小鼠和人的训练数据集
# 加载内置数据
set.seed(100)
library(scran)
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds",
package="scran"))
> class(mm.pairs)
[1] "list"
> names(mm.pairs)
[1] "G1" "S" "G2M"
> head(mm.pairs$G1)
first second
1 ENSMUSG00000000001 ENSMUSG00000001785
2 ENSMUSG00000000001 ENSMUSG00000005470
3 ENSMUSG00000000001 ENSMUSG00000012443
4 ENSMUSG00000000001 ENSMUSG00000015120
5 ENSMUSG00000000001 ENSMUSG00000022033
6 ENSMUSG00000000001 ENSMUSG00000023015
# 进行比较
assignments <- cyclone(sce.416b, mm.pairs, gene.names=rowData(sce.416b)$ENSEMBL)
> class(assignments)
[1] "list"
> names(assignments)
[1] "phases" "scores"
[3] "normalized.scores"
> dim(assignments$scores)
[1] 185 3
> head(assignments$scores)
G1 S G2M
1 0.663 0.995 0.000
2 0.960 0.456 0.003
3 0.055 0.517 0.144
4 0.000 0.768 1.000
5 0.927 0.473 0.002
6 0.914 0.944 0.000
看这个score结果
其中对416b数据中的每个细胞都划定了一个分数，分数越高就越倾向于那个周期，比如看看G1和G2/M期
plot(assignments$score$G1, assignments$score$G2M,
xlab="G1 score", ylab="G2/M score", pch=16)
判断标准
如果G1 分数高于0.5并且大于G2/M的分数，那么就判断为G1期；同样的方法可以判断G2/M期；如果二者分数都不大于0.5，那就是S期
table(assignments$phases, colLabels(sce.416b))
##
## 1 2 3 4
## G1 74 8 20 0
## G2M 1 48 0 13
## S 3 13 4 1
这种方法使用简单，就是实际使用过程可能会耗时
5 移除细胞周期导致的差异
有时担心细胞周期带来的数据差异会影响下游分析，如果要去掉这部分影响，可以先按照上面的方法推断出细胞周期，然后将每个周期当成一个独立的批次，再进行批次矫正。最简单的方法是使用线性模型（如regressBatches()）来移除这个影响
library(batchelor)
sce.nocycle <- regressBatches(sce.416b, batch=assignments$phases)
# 另外其他支持block的函数也是可以移除这个影响的
dec.nocycle <- modelGeneVarWithSpikes(sce.416b, "ERCC",
block=assignments$phases)
marker.nocycle <- findMarkers(sce.416b, block=assignments$phases)
不过，这个操作并非常规的分析流程。因为大部分数据中，细胞周期导致的差异是次要的，不会超过细胞真实生物差异的影响。
