使用scran包的cyclone函数进行单细胞转录组的细胞周期状态推断

多次在CNS文章里面看到：Gene sets reflecting five phases of the HeLa cell cycle (G1/S, S, G2/M, M and M/G1) were referred from Whitfield et al. (2002). 其中一篇提到了scran包的cyclone函数进行单细胞转录组的细胞周期状态推断，学习一下并分享给大家。

模拟单细胞转录组数据

首先模拟一个单细胞转录组表达矩阵，代码如下：

library(scran)
set.seed(100)
## ------------------------------------------------------------------------
ngenes <- 10000
ncells <- 200
mu <- 2^runif(ngenes, -1, 5)
gene.counts <- matrix(rnbinom(ngenes*ncells, mu=mu, size=10), nrow=ngenes)

## ------------------------------------------------------------------------
library(org.Mm.eg.db)
all.ensembl <- unique(toTable(org.Mm.egENSEMBL)$ensembl_id)
rownames(gene.counts) <- sample(all.ensembl, ngenes)

## ------------------------------------------------------------------------
nspikes <- 100
ncells <- 200
mu <- 2^runif(nspikes, -1, 5)
spike.counts <- matrix(rnbinom(nspikes*ncells, mu=mu, size=10), nrow=nspikes)
rownames(spike.counts) <- paste0("ERCC-", seq_len(nspikes))
all.counts <- rbind(gene.counts, spike.counts)


一万个基因在200个细胞的表达矩阵，还模拟了100个ERCC的表达量，用来做QC。

使用cyclone进行细胞周期推断

代码如下：

library(scran)
sce <- SingleCellExperiment(list(counts=all.counts))
isSpike(sce, "MySpike") <- grep("^ERCC", rownames(sce))
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", package="scran"))
assigned <- cyclone(sce, pairs=mm.pairs)
head(assigned$scores)
table(assigned$phases)

或者对真实测序数据进行推断

set.seed(6473)
library(scRNAseq)
library(RColorBrewer)
library(scran)
all.counts=assays(fluidigm)$rsem_counts
dim(all.counts)
head(rownames(all.counts))
sce <- SingleCellExperiment(list(counts=all.counts))
library(org.Hs.eg.db)
ensembl <- mapIds(org.Hs.eg.db, keys=rownames(sce), keytype="SYMBOL", column="ENSEMBL")
mm.pairs <- readRDS(system.file("exdata", "human_cycle_markers.rds", package="scran"))
assigned <- cyclone(sce, pairs=mm.pairs, gene.names=ensembl)
head(assigned$scores)
table(assigned$phases)
plot(assigned$score$G1, assigned$score$G2M,
    xlab="G1 score", ylab="G2/M score", pch=16)

最后得到了每个细胞的细胞周期推断情况：

得到了每个细胞的各种周期得分，那么就可以制作热图啦，我随意演示如下：

> head(assigned$scores)
    G1     S   G2M
1 0.889 0.336 0.247
2 0.880 0.271 0.366
3 0.844 0.296 0.473
4 0.861 0.037 0.436
5 0.512 0.187 0.613
6 0.696 0.327 0.093
> pheatmap::pheatmap(t(assigned$scores))

因为是模拟数据，而且我是随意写的绘图代码，所以必然不如SCI文章配图，但是意思到了。

而且，有趣的是这个文章用的是5个阶段，而scran包的cyclone函数只提供了3个阶段，所以，要完全重现这个流程，还需要再去搜索看看。