单细胞学习---细胞周期

看人和动物单细胞相关的文章时候，都会有个关于细胞周期的分析，但是植物相关文章的话，这块还很少涉及。我们今天就来学习一下怎么从单细胞数据里面获得细胞周期的故事。

我们先看下维基百科上对于cell cycle的定义：

A cell cycle is a series of events that takes place in a cell as it grows and divides.

即描述细胞生长、分裂整个过程中细胞变化过程。最重要的两个特点就是DNA复制、分裂成两个一样的子细胞。

如上图所示，一般分成4个阶段：

G1(gap1)：Cell increases in size(Cellular contents duplicated)

S(synthesis) ：DNA replication, each of the 46 chromosomes (23 pairs) is replicated by the cell

G2(gap2)：Cell grows more，organelles and proteins develop in preparation for cell division，为分裂做准备

M(mitosis)：'Old' cell partitions the two copies of the genetic material into the two daughter cells.

And the cell cycle can begin again.

因此，在分析单细胞数据时，同一类型的细胞往往来自于不同的细胞周期阶段，这可能对下游聚类分析，细胞类型注释产生混淆；

由于细胞周期也是通过cell cycle related protein 调控，即每个阶段有显著的marker基因；

通过分析细胞周期有关基因的表达情况，就有望对细胞所处周期阶段进行注释。

注：但是，在单细胞周期分析时，通常只考虑三个阶段：G1、S、G2M。(即把G2和M当做一个phase)

下面，我们来测试如何通过单细胞数据来推断细胞周期。

=======scran包cyclone注释cell cycle=======

1、方法简介

http://bioconductor.org/packages/release/bioc/manuals/scran/man/scran.pdf

scran包cyclone函数是利用‘marker基因对’表达来对细胞所在周期阶段进行预测的方法Scialdone (2015)

“maker基因对”由作者根据训练集细胞（已注释了cell cycle）的基因表达特征产生，我们可以直接使用。对于每一细胞周期阶段（人/鼠）都有一组“maker基因对”集合。

library(scran)

scran包里面包含人和小鼠的marker基因对。

# hs.pairs <- readRDS(system.file("exdata", "human_cycle_markers.rds", package="scran"))

mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", package="scran"))

str(mm.pairs)

head(mm.pairs$G1)

如上图所示，就是小鼠G1的marker基因对。具体来说，即比较某个细胞的ENSMUSG00000000001基因表达值是否大于ENSMUSG00000001785基因表达值。

如果对于所有G1期marker基因对，某个细胞的“first”列基因表达量大于对应的“second”基因的情况越多，则越有把握认为该细胞就是处于G1期（因为越符合训练集特征）。

而cyclone则是通过计算score，即对于某个细胞，符合上述比较关系的marker基因对数占全部marker基因对数的比值。

=======示例数据一=======

http://bioconductor.org/books/release/OSCA/cell-cycle-assignment.html#motivation-12

第一个测试数据：单细胞数据来自scRNAseq包的LunSpikeInData。选这个数据集是因为It is known to contain actively cycling cells after oncogene induction.

该数据包含192 cells，46,604 genes；实验及部分注释信息已经提供：是研究embryonic stem cell(ESC)胚胎肝细胞经癌基因诱导前后的变化。

数据前处理主要(主要基于sce对象的操作)

library(scRNAseq)

sce.416b <- LunSpikeInData(which="416b")

sce.416b$block <- factor(sce.416b$block)

下面将ID进行转化：（因为scran默认提供marker基因对是ensemble格式，因此其它类基因ID需要转换一下）

library(AnnotationHub)

ens.mm.v97 <- AnnotationHub()[["AH73905"]]

rowData(sce.416b)$ENSEMBL <- rownames(sce.416b)

rowData(sce.416b)$SYMBOL <- mapIds(ens.mm.v97, keys=rownames(sce.416b),

                                  keytype="GENEID", column="SYMBOL")

下面是单细胞数据的常规操作，和seurat处理一样：

library(scater)

sce.416b <- logNormCounts(sce.416b)

#--- variance-modelling 挑选高变基因---#

dec.416b <- modelGeneVarWithSpikes(sce.416b, "ERCC", block=sce.416b$block)

chosen.hvgs <- getTopHVGs(dec.416b, prop=0.1)

#--- dimensionality-reduction PCA降维---#

sce.416b <- runPCA(sce.416b, ncomponents=10, subset_row=chosen.hvgs)

head(head(colData(sce.416b)))

relabel <- c("onco", "WT")[factor(sce.416b$phenotype)]

plotPCA(sce.416b, colour_by=I(relabel), shape_by=I(relabel))

如下图所示，是各个cell的二维PCA图，其中WT代表wild type phenotype组；onco代表induced CBFB-MYH11 oncogene expression组。

下面，我们利用scran包来注释不同细胞的细胞周期。

mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", package="scran"))

assignments <- cyclone(sce.416b, pairs=mm.pairs)

结果包含三个部分，phases（最终分配的时期）；scores（每个时期的得分）；已经标准化的scores。

head(assignments$scores)

所以，我们根据每一个细胞对于三个周期阶段的scores，可进行判断；

1. 如果G1或者G2M的值大于0.5，则为G1或者G2M时期；

2. 如果G1和G2M均小余0.5，则为S时期；

3. 如果G1和G2M都大于0.5，则那个高是那个时期。

其实逻辑就是，看下G1/G2M最大值是否大于0.5，如果是，则比较G1和G2M大小，大的为分配时期；反之，为S时期。

画个图就好理解了。

plot(assignments$score$G1, assignments$score$G2M,

    xlab="G1 score", ylab="G2/M score", pch=16)

abline(h = 0.5, col='red')

abline(v = 0.5, col='red')

如下图所示：除了左下角是S时期，其余就看G1和G2M谁大，谁为主。

table(assignments$phases)

plotPCA(sce.416b, colour_by=I(relabel), shape_by=I(assignments$phases))

我们也能查看不同组细胞中细胞周期的情况。从图中来看，对于WT组，S和G2M几乎很近，分不开。但是在onco组，G2M则可以完全分开，S和G1时期比较近。

========示例数据二========

还是用以前的pbmc的数据。还是因为scran记录的ID都是ENSEMBL ID，所以要转换，有两种方式。

第一种方法：改变库里的ID，也就是转换hs.pairs里面的ID为symbol。

htrs <- lapply(names(hs.pairs), function(x){

df <- hs.pairs[[x]]

first <- bitr(hs.pairs[[x]][,1],fromType= "ENSEMBL", toType="SYMBOL",OrgDb="org.Hs.eg.db")

second <- bitr(hs.pairs[[x]][,2],fromType= "ENSEMBL", toType="SYMBOL",OrgDb="org.Hs.eg.db")

df <- df[first$ENSEMBL %in% df$first,]

df <- df[second$ENSEMBL %in% df$second,]

df$first <-  lapply(df$first, function(x){

    first[first$ENSEMBL == x,2][1]

  })  %>% unlist()

df$second <-  lapply(df$second, function(x){

    second[second$ENSEMBL == x,2][1]

  })  %>% unlist()

  return(df)

})

names(htrs) <- names(hs.pairs)

把hs.pairs里面的ENSEMBL ID转化成常见的symbol。

然后，我们就可以常规的利用

pbmc <- readRDS("pbmc.rds")

library(scran)

assigndata <- cyclone(pbmc@assays[["RNA"]]@data, pairs=htrs)

table(assigndata$phases)

plot(assigndata$score$G1, assigndata$score$G2M, xlab="G1 score", ylab="G2/M score", pch=16)

abline(h = 0.5, col='red')

abline(v = 0.5, col='red')

也可以像前面一样，看下每个cell 的G1和G2M score的分布情况。

然后把phase加入到pbmc的meta.data中：

pbmc@meta.data$scan_phase<-assigndata[["phases"]]

第二种方法：转化Seurat数据里面的symbol，也就是转化symbol为ID。

ensembl <- mapIds(org.Hs.eg.db, keys=rownames(pbmc), keytype="SYMBOL", column="ENSEMBL")

assigndata <- cyclone(pbmc@assays[["RNA"]]@data, pairs=hs.pairs, gene.names=ensembl)

注意：这两种转换方式最终计算的结果是不完全一致的！！！至于选择哪种，我个人倾向于第一种，因为第二种好像很像很多ID都不能识别。当然，也可以自己转化一下。

=======使用Seurat的CellCycleScoring来注释cell cycle=======

Seurat包CellCycleScoring较scran包cyclone函数最主要的区别是直接根据每个cycle，一组marker基因表达值判断。

library(Seurat)

str(cc.genes)

还有个2019年版本的数据集：

s.genes <- cc.genes.updated.2019$s.genes

g2m.genes <- cc.genes.updated.2019$g2m.genes

上面的cc.genes就是Seurat包提供的人的细胞中分别与S期、G2M期直接相关的marker基因。

CellCycleScoring即根据此，对每个细胞的S期、G2M期可能性进行打分；具体如何计算的，暂时在Seurat官方文档中提及的较少。

查资料后，大致是这样子的：

根据每个细胞的S期(或者G2/M期)基因集是否显著高表达，对应的score就是表示在该细胞中，S期(或者G2/M期)基因集高表达的程度(如果是负数，就认为不属于该phase)。

区别于scran包的另外重要的一点就是Seurat包仅提供了人类细胞有关的cell cycle related gene，没有小鼠的。小鼠的需要通过同源基因去获取。

下面，我们来测试一下Seurat执行细胞周期鉴定的结果。

srat <- CellCycleScoring(pbmc, s.features = s.genes, g2m.features = g2m.genes)

plot(srat@meta.data$S.Score,srat@meta.data$G2M.Score,col=factor(srat@meta.data$Phase),main="CellCycleScoring")

legend("topleft",inset=.05,title = "cell cycle",  c("G1","S","G2M"), pch = c(1),col=c("black","green","red"))

如下图，分类条件如上面所说：细胞的S.Score与G2M.Score均小于0时，则为G1期；否则那个值大，就是属于哪个phase。

但是，从结果来看，Scran和Seurat的结果更是差别巨大!!!最终选择哪种，我个人的原则是去看一下这些原文的数据集是来自什么细胞。哪种和自己细胞相近用哪个，如果都不相近，选多的那个。