学到了！别人发了CNS的数据，我还能发吗？教你一招，做出新意！（附代码）

单细胞数据整合细胞周期分析

大家好，我是风风。欢迎来到风风的从零开始单细胞系列。在bulk-seq的分析中，“挑、圈、联、靠”中的联，我们常常会结合一些特定表型的数据进行联合分析，而从生信文章切入角度来说，特定表型作为创新进行生物信息学分析的文章也越来越多，像免疫基因、铁死亡、缺氧、自噬等。在这么多表型中，细胞周期是最经典那几种表型之一。在二十四型的课程中，细胞周期也是最先进行讲解的表型之一。那么单细胞数据分析可以做细胞周期方面的分析吗？答案当然是肯定的，而且是非常适合！

我们可以从scRNA seq数据中确定细胞周期活动。就其本身而言，细胞在细胞周期各阶段的分布通常不具有信息量，但我们可以通过它来确定亚群之间或不同治疗条件下的增殖是否存在差异。细胞周期中的许多关键进程，例如通过细胞周期检查点，便是由翻译后机制驱动的，所以不能在转录组数据中直接看到。但是我们也可以利用相关基因的表达变化来确定细胞周期的各个阶段。

Cyclin蛋白表达变化确定细胞周期阶段

细胞周期蛋白（Cyclin）控制细胞周期的进展，并在细胞周期各阶段有明确的表达模式。

其中，Cyclin D的表达贯穿整个细胞周期，但在G1期达到最高的峰值；Cyclin E在G1/S转换期表达最高；Cyclin A在S期和G2期表达；Cyclin B在G2后期和有丝分裂期表达最高（关于细胞周期的Cyclin蛋白更多背景知识可以查看这篇文章：The Cell Cycle: Principles of Control. New Science Press）。因此，我们可以利用Cyclin蛋白的表达可以帮助确定每个簇中的相对细胞周期活性。接下来我们看一下相应的分析内容：

rm(list = ls())# 安装并加载R包if(!require("scater")) install.packages("scater", update = F, ask = F)if(!require("scran")) install.packages("scran", update = F, ask = F)if(!require("batchelor")) install.packages("batchelor", update = F, ask = F)if(!require("SingleR")) install.packages("SingleR", update = F, ask = F)if(!require("org.Mm.eg.db")) BiocManager::install("org.Mm.eg.db")# 加载分析数据，这里使用我们前面推文分析处理后的sce.416b数据load(file = "inputdata.Rda")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## mainExpName: endogenous## altExpNames(2): ERCC SIRV# 获取Cyclin蛋白对应的基因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"# cyclin蛋白基因表达值的热图plotHeatmap(sce.416b, order_columns_by = "label", cluster_rows = FALSE, features = sort(cyclin.genes))

# 这里我们用标准的scran差异分析方法对这些观察结果进行量化，以测试cluster之间每个细胞周期蛋白的上调，这将意味着一个亚群包含更多处于相应细胞周期阶段的细胞，比如，我们可以根据cyclin蛋白A2和B1的高表达，推断出cluster 4在S期和G2期的细胞比例最高

markers <- findMarkers(sce.416b, subset.row = cyclin.genes, test.type = "wilcox", direction = "up")markers## List of length 4## names(4): 1 2 3 4# 查看结果markers@listData[["1"]]; markers@listData[["2"]]; markers@listData[["3"]]; markers@listData[["4"]]## DataFrame with 10 rows and 7 columns## Top p.value FDR summary.AUC AUC.2 AUC.3## <integer> <numeric> <numeric> <numeric> <numeric> <numeric>## Ccnd2 1 3.27150e-10 1.63575e-09 0.803976 0.8039762 0.723825## Ccnd1 1 1.02655e-12 1.02655e-11 0.843924 0.8439242 0.795406## Ccnd3 3 3.77192e-01 1.00000e+00 0.597070 0.4329246 0.439637## Ccna1 4 1.00000e+00 1.00000e+00 0.464286 0.4565217 0.500000## Ccne1 4 9.99999e-01 1.00000e+00 0.435897 0.2788926 0.390491## Ccne2 5 9.99917e-01 1.00000e+00 0.358059 0.3381643 0.341880## Ccnb1 5 1.00000e+00 1.00000e+00 0.353098 0.0575994 0.353098## Ccnb2 7 1.00000e+00 1.00000e+00 0.237447 0.1044221 0.237447## Ccna2 8 1.00000e+00 1.00000e+00 0.229701 0.0123560 0.229701## Ccnb3 9 1.00000e+00 1.00000e+00 0.500000 0.5000000 0.500000## AUC.4## <numeric>## Ccnd2 0.8653846## Ccnd1 0.6318681## Ccnd3 0.5970696## Ccna1 0.4642857## Ccne1 0.4358974## Ccne2 0.3580586## Ccnb1 0.0503663## Ccnb2 0.0659341## Ccna2 0.0036630## Ccnb3 0.5000000## DataFrame with 10 rows and 7 columns## Top p.value FDR summary.AUC AUC.1 AUC.3## <integer> <numeric> <numeric> <numeric> <numeric> <numeric>## Ccna2 1 2.76781e-24 2.76781e-23 0.987644 0.987644 0.918478## Ccnd3 1 1.68417e-02 2.40595e-02 0.716356 0.567075 0.515700## Ccnd2 2 6.43358e-02 8.04197e-02 0.672878 0.196024 0.449275## Ccnb1 2 3.20876e-20 1.60438e-19 0.942401 0.942401 0.917874## Ccne1 3 1.72663e-06 4.31657e-06 0.721107 0.721107 0.631643## Ccnb2 3 2.12607e-16 7.08688e-16 0.895578 0.895578 0.679952## Ccne2 4 2.52964e-04 5.05929e-04 0.661836 0.661836 0.507850## Ccna1 5 1.23004e-02 2.05006e-02 0.543478 0.543478 0.543478## Ccnd1 8 1.00000e+00 1.00000e+00 0.474638 0.156076 0.474638## Ccnb3 10 1.00000e+00 1.00000e+00 0.500000 0.500000 0.500000## AUC.4## <numeric>## Ccna2 0.358178## Ccnd3 0.716356## Ccnd2 0.672878## Ccnb1 0.480331## Ccne1 0.633540## Ccnb2 0.218427## Ccne2 0.552795## Ccna1 0.504658## Ccnd1 0.177019## Ccnb3 0.500000## DataFrame with 10 rows and 7 columns## Top p.value FDR summary.AUC AUC.1 AUC.2## <integer> <numeric> <numeric> <numeric> <numeric> <numeric>## Ccna2 1 8.82266e-05 0.000794417 0.770299 0.770299 0.0815217## Ccnd2 1 7.64325e-02 0.127387429 0.693452 0.276175 0.5507246## Ccnd3 1 3.34927e-02 0.083731626 0.726190 0.560363 0.4842995## Ccnd1 2 9.99994e-01 1.000000000 0.525362 0.204594 0.5253623## Ccnb2 2 1.58883e-04 0.000794417 0.762553 0.762553 0.3200483## Ccne2 3 9.49766e-03 0.031658861 0.658120 0.658120 0.4921498## Ccne1 4 1.05985e-01 0.151407476 0.609509 0.609509 0.3683575## Ccnb1 4 4.39155e-02 0.087830964 0.646902 0.646902 0.0821256## Ccna1 5 1.00000e+00 1.000000000 0.464286 0.500000 0.4565217## Ccnb3 9 1.00000e+00 1.000000000 0.500000 0.500000 0.5000000## AUC.4## <numeric>## Ccna2 0.0744048## Ccnd2 0.6934524## Ccnd3 0.7261905## Ccnd1 0.2232143## Ccnb2 0.1011905## Ccne2 0.5446429## Ccne1 0.5267857## Ccnb1 0.0654762## Ccna1 0.4642857## Ccnb3 0.5000000## DataFrame with 10 rows and 7 columns## Top p.value FDR summary.AUC AUC.1 AUC.2## <integer> <numeric> <numeric> <numeric> <numeric> <numeric>## Ccna2 1 4.47082e-09 4.47082e-08 0.996337 0.996337 0.641822## Ccnd1 1 2.27713e-04 5.69283e-04 0.822981 0.368132 0.822981## Ccnb1 1 1.19027e-07 5.95137e-07 0.949634 0.949634 0.519669## Ccnb2 2 3.87799e-07 1.29266e-06 0.934066 0.934066 0.781573## Ccna1 4 2.96992e-02 5.93985e-02 0.535714 0.535714 0.495342## Ccne2 5 6.56983e-02 1.09497e-01 0.641941 0.641941 0.447205## Ccne1 6 5.85979e-01 8.37113e-01 0.564103 0.564103 0.366460## Ccnd3 7 9.94578e-01 1.00000e+00 0.402930 0.402930 0.283644## Ccnd2 8 9.99993e-01 1.00000e+00 0.306548 0.134615 0.327122## Ccnb3 10 1.00000e+00 1.00000e+00 0.500000 0.500000 0.500000## AUC.3## <numeric>## Ccna2 0.925595## Ccnd1 0.776786## Ccnb1 0.934524## Ccnb2 0.898810## Ccna1 0.535714## Ccne2 0.455357## Ccne1 0.473214## Ccnd3 0.273810## Ccnd2 0.306548## Ccnb3 0.500000# 使用grun.hsc数据进行细胞周期分析，发现分类的HSCs中cyclin蛋白D2上调load(file = "sce.grun.hsc.Rda")# 改变行名rownames(sce.grun.hsc) <- uniquifyFeatureNames(rownames(sce.grun.hsc), rowData(sce.grun.hsc)$SYMBOL)cyclin.genes <- grep("^Ccn[abde][0-9]$", rowData(sce.grun.hsc)$SYMBOL)cyclin.genes <- rownames(sce.grun.hsc)[cyclin.genes]# 热图可视化plotHeatmap(sce.grun.hsc, order_columns_by = "label", cluster_rows = FALSE, features = sort(cyclin.genes), colour_columns_by = "protocol")

# 需要注意的是，这种方法是假定了细胞周期蛋白的表达不受细胞周期以外的生物过程的影响这个前提，但是细胞周期并不独立于细胞中发生的其他过程，所以，细胞周期推断最好用于密切相关的细胞类型之间的比较，因为其他细胞过程的变化较少，导致结果也相对更加可靠。

使用背景数据集确定细胞周期阶段

细胞周期分配可以被认为是细胞注释的一个特殊例子，因此，我们也可以将前面用过的细胞注释的方法用在这里，给定一个包含已知细胞周期阶段的背景数据，我们可以使用SingleR等方法来确定测试数据集中每个细胞的细胞周期阶段，具体操作如下：

load(file = "sce.ref.Rda")# 确定与细胞周期相关的基因library(org.Mm.eg.db)cycle.anno <- select(org.Mm.eg.db, keytype = "GOALL", keys = "GO:0007049", columns = "ENSEMBL")[,"ENSEMBL"]## 'select()' returned 1:many mapping between keys and columnsstr(cycle.anno)## chr [1:2812] "ENSMUSG00000026842" "ENSMUSG00000026842" ...# 转换ID以匹配引用test.data <- logcounts(sce.416b)rownames(test.data) <- rowData(sce.416b)$ENSEMBL# 我们将使用SingleR()函数，根据ESC参考中的细胞周期阶段为数据分配标签library(SingleR)assignments <- SingleR(test.data, ref = sce.ref, label = sce.ref$phase, de.method = "wilcox", restrict = cycle.anno)tab <- table(assignments$labels, colLabels(sce.416b))tab## ## 1 2 3 4## G1 71 5 18 1## G2M 2 60 2 13## S 5 4 4 0

# 结果发现，cluster 1主要由G1细胞组成，而其他cluster则有更多的细胞处于其他阶段，这与我们基于细胞周期蛋白的分析结论大致一致。与基于细胞周期蛋白的分析不同，这种方法产生了细胞周期阶段的 "绝对 "分配，不需要相对于同一数据集中的其他细胞进行解释

gridExtra::grid.arrange( plotExpression(sce.ref, features = "ENSMUSG00000080800", x = "phase"), plotExpression(sce.416b, features = "Lef1", # 区分G1和G2/M的标志基因 x="label"), ncol=2)

# 卡方检验分析cluster 1 和 cluster 2的分布差异chisq.test(tab[,1:2])## Warning in chisq.test(tab[, 1:2]): Chi-squared approximation may be incorrect## ## Pearson's Chi-squared test## ## data: tab[, 1:2]## X-squared = 111.55, df = 2, p-value < 2.2e-16# 结果发现两者具有显著统计学意义

使用scran的cyclone()分类确定细胞周期阶段

Scialdone等人在2015年提出了一种新的细胞分为细胞周期阶段的又一方法：首先使用参考数据集，计算每对基因之间的表达差异的上下调。调控方向在不同细胞周期阶段有变化的一对基因被选为标志物。然后，根据观察到的每对标记物的方向是否与一个阶段或另一个阶段一致，可以将测试数据集中的细胞分类到合适的细胞周期阶段。这种方法的详细内容可以参考下面这篇文献：Computational assignment of cell-cycle stage from single-cell transcriptome data。我们可以使用scran包的cyclone()函数实现该方法：

library(scran)mm.pairs <- readRDS(system.file("exdata","mouse_cycle_markers.rds", package = "scran"))mm.pairs# 使用 EnsemblID匹配mm.pairs进行注释assignments <- cyclone(sce.416b, mm.pairs, gene.names = rowData(sce.416b)$ENSEMBL)

# 对于每个细胞，一个细胞周期阶段的得分越高，对应于该细胞处于该阶段的概率就越高，这里我们将重点放在G1和G2/M分数上

plot(assignments$score$G1, assignments$score$G2M, xlab = "G1 score", ylab = "G2/M score", pch = 21, cex = 2.5, col = "steelblue")

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

去除细胞周期的影响

如果我们不是专门研究细胞周期的话，我们可能更想做的是去除细胞周期对于数据的影响（前面的推文也提到过相应去除办法）。在大多数情况下，细胞周期跟细胞类型特征不一样，更像是一个次要的变异因素。此外，当细胞周期活动在不同的细胞类型或条件下发生变化时，大多数去除的策略都会遇到问题，例如研究T细胞在激活后增殖，我们可能就想要去除细胞周期带来的影响。去除细胞周期影响的方法有很多，我们以前使用过直接去除细胞周期相关的基因，这里使用线性回归的方法进行去除：

# 我们假设每一种细胞类型包含相同的细胞在不同阶段的分布，并且细胞周期对表达的影响恒定，接着进行回归调整具有差异的过程，例如如果两个亚群的细胞周期阶段分布不同，回归将总是对这些亚群之间的所有差异基因进行调整

library(batchelor)dec.nocycle<- modelGeneVarWithSpikes(sce.416b, "ERCC",block = assignments$phases)reg.nocycle<- regressBatches(sce.416b, batch = assignments$phases)reg.nocycle<- runPCA(reg.nocycle,exprs_values = "corrected",subset_row = getTopHVGs(dec.nocycle, prop = 0.1))relabel<- c("onco", "WT")[factor(sce.416b$phenotype)]scaled<- scale_shape_manual(values = c(onco = 4, WT = 16))# 可视化gridExtra::grid.arrange(plotPCA(sce.416b,colour_by = I(assignments$phases),shape_by = I(relabel)) + ggtitle("Before")+ scaled,plotPCA(reg.nocycle,colour_by = I(assignments$phases),shape_by = I(relabel)) + ggtitle("After")+ scaled,ncol = 2)

# 构建线性模型design<- model.matrix( ~ as.matrix(assignments$scores))dec.nocycle2<- modelGeneVarWithSpikes(sce.416b, "ERCC",design = design)reg.nocycle2<- regressBatches(sce.416b, design = design)reg.nocycle2<- runPCA(reg.nocycle2, exprs_values = "corrected",subset_row = getTopHVGs(dec.nocycle2, prop = 0.1))plotPCA(reg.nocycle2,colour_by = I(assignments$phases), point_size = 5, shape_by = I(relabel)) + scaled

reg.nocycle3 <-fastMNN(sce.416b,batch = assignments$phases)plotReducedDim(reg.nocycle3, dimred = "corrected",point_size = 5,colour_by = I(assignments$phases),shape_by = I(relabel)) + scaled

# PCA结果展示fastMNN()后的校正PCA与来自cyclone()的细胞周期分配的关系图，每个点是一个细胞，根据其推断的细胞周期阶段添加颜色，并根据肿瘤基因的诱导状态调整形状

好了，本期的内容到此结束！在本期中，我们利用几种不同的方法来对单细胞数据中的细胞判断细胞周期的阶段，这可能有助于我们研究一些和细胞周期活动有关的表型，比如侵袭等。当然我们以往在看到单细胞教程甚至包括我们前面的推文，也说过单细胞数据的预处理中，要去除掉细胞周期带来的影响，除了以前说过的去除细胞周期相关基因，这里也提供了利用线性回归去除细胞周期带来的影响的做法。具体的情况还要根据大家研究的方向来进行使用。当然，既然我们可以通过标志物markers或者相应背景数据集来确定细胞周期的阶段，我们是不是也可以用类似的方法来确定T细胞的功能阶段，或者确定免疫循环的各个阶段呢？

后台回复“feng 55”获取代码和数据，我们下期见！

单细胞分析专栏传送门

大年初一，单细胞分析从nature communication开始

8+单细胞文章从解读到套路总结全部送给你

30分多的单细胞分析文章是怎么做的？

不用R，点点点就能get的单细胞公共数据下载，你“学废”了吗

单细胞都有两个“对象”，你一个都没有

8成生信审稿人必问！单细胞第一张可视化图，你漏了吗？

我从Science中，“偷”到了这个实用的聚类分析技能！真舍不得分享！（附代码）

不骗你！学了这个R包，搞定九成高大上的单细胞分析！包治百病！一包封神！小白赶紧收起来！

Cell的生信分析能有多高级？单细胞研究必备分析！不会做么？简单，我来教你！（附代码）

放大招！单细胞筛分子标志物，让你的课题逼格up! up! up! 你一定要试试！（附详细代码）