分裂期的细胞分析时需要过滤吗？

分享是一种态度

作者 | 周运来

男，

一个长大了才会遇到的帅哥，

稳健，潇洒，大方，靠谱。

一段生信缘，一棵技能树。

生信技能树核心成员，单细胞天地特约撰稿人，简书创作者，单细胞数据科学家。

前段时间朋友问我细胞分裂时期的细胞在分析之前是不是应该去掉，因为它们不是正常的细胞：DNA加倍了。

我抓紧时间看看RNA是不是也加倍了。

请出我们的Seurat和pbmc3k数据集：

library(Seurat)
library(SeuratData)
pbmc3k.final 

An object of class Seurat 
13714 features across 2638 samples within 1 assay 
Active assay: RNA (13714 features, 2000 variable features)
 2 dimensional reductions calculated: pca, umap

用Seurat自带的功推断一下细胞的状态：

pbmc3k.final <- CellCycleScoring(pbmc3k.final,
                                   g2m.features = cc.genes.updated.2019$g2m.genes,
                                   s.features = cc.genes.updated.2019$s.genes)

VlnPlot(pbmc3k.final,features = c("nCount_RNA","nFeature_RNA"),group.by = "Phase",pt.size = 0)

在分布上并没有什么差别，但是G2M还是有几个细胞较高：

pbmc3k.final@meta.data %>% filter(nCount_RNA > 7500,nFeature_RNA> 2000)
               orig.ident nCount_RNA nFeature_RNA seurat_annotations percent.mt RNA_snn_res.0.5 seurat_clusters     S.Score  G2M.Score Phase
ACGAACTGGCTATG     pbmc3k       8875         2413           Platelet   1.577465               8               8 -0.02272569 0.02334952   G2M
GGGCCAACCTTGGA     pbmc3k       8415         2013                 DC   1.057635               7               7 -0.01562669 0.01747374   G2M

就俩细胞啊。

我想起一种带凹口的箱型图可以比较差异：

pbmc3k.final@meta.data -> metta
metta[,c("nCount_RNA","nFeature_RNA","Phase")] -> metta
library(reshape2)
head(melt(metta) -> metta)

g <- ggplot(metta, aes(Phase, value))
g + geom_boxplot(varwidth=T,notch = T)  + facet_grid(.~variable) + theme_bw() +theme(
  strip.background = element_rect(
    color="black", fill="#FC4E07", size=1.5, linetype="solid"
  )
)

然而并没有什么差异。

我以Phase分割开来呢？会不会有细胞类型特异的细胞周期状态呢?

DimPlot(pbmc3k.final,
        label = TRUE, 
        split.by = "Phase")  + NoLegend()

DimPlot(pbmc3k.final,group.by = "Phase",
        label = TRUE,repel = T)  

于是，我很想知道在图谱上细胞周期的score值是怎的。

Idents(pbmc3k.final) <- 'seurat_annotations'
metrics <-  c("nCount_RNA", "nFeature_RNA", "S.Score", "G2M.Score",'percent.mt')

FeaturePlot(pbmc3k.final, 
            reduction = "umap", 
            features = metrics,
            pt.size = 0.4, 
            order = TRUE,
            min.cutoff = 'q10',
            label = TRUE)

这个Score值真的很低啊，我开始怀疑这个打分体系了：它默认细胞都处在分裂期，里面没有细胞的静息状态G0期。

于是，我想做个差异基因看看预测得到的每个时期的差异基因：

Idents(pbmc3k.final) <- 'Phase'
pbmc3k.final@misc$ccmk<- FindAllMarkers(pbmc3k.final,only.pos = T)
pbmc3k.final@misc$ccmk %>% group_by(cluster) %>%top_n(20,avg_logFC) -> ccmk 
ccseo <- ScaleData(pbmc3k.final,features = ccmk$gene)

intersect(ccmk$gene,c(cc.genes.updated.2019$s.genes,cc.genes.updated.2019$g2m.genes))
 [1] "NASP"   "SLBP"   "MRPL36" "PCNA"   "USP1"   "HMGB2"  "TUBB4B" "TACC3"  "CBX5"   "ANP32E" "LBR"    "TMPO"   "CKS1B"  "SMC4"   "CTCF"  

DoHeatmap(ccseo,features = ccmk$gene)

还是没看出来分裂期的RNA怎么变化的，问题是细胞分裂的时候它是如何保持其身份的。在分裂的时候，比如B细胞，在分裂的时候还是B细胞：RNA转录模式没有变。最后，我不得不请教浏览器：

细胞分裂前，细胞内基因数量翻倍，细胞如何维持mRNA水平不变?

想象一下在你开手动挡的时候，想要保持恒定的速度，这当然很简单，但如果你只被允许间歇地踩油门踏板呢？这种驱动技术看起来很奇怪，但细胞在转录基因以产生信使RNA (mRNA)分子时，通常表现出类似的行为方式。当被称为启动子的基因区域打开时，信使rna分子就产生了：由于这些启动子在开放和封闭状态之间随意切换，信使rna分子就会大量产生。细胞也会降解信使rna，就像道路摩擦会降低汽车的速度一样。因此，mRNA水平会随着时间上升或下降

现在让我们回到我们的汽车类比，并添加另一个复杂的东西。如果另一个司机奇迹般地出现在你旁边（副驾），她独立控制你们的油门，那该怎么办?现在你怎么可能不超速呢？加速增加一倍后，你必须通过减少踩油门踏板的频率或力度来进行补偿。

细胞在复制时也面临着类似的挑战(图1)。在细胞周期的某些阶段，细胞中每个基因的副本是正常情况的两倍，这可能会使产生的mRNA数量增加一倍。复制细胞如何补偿这一点，使信使rna水平保持不变?现在，在eLife中，Ido Golding和他的同事——包括第一作者Samuel Skinner——提出了一个可以用来解决这个问题的优雅框架。

在整个细胞周期中维持恒定的mRNA水平。在细胞周期的G1期(蓝色)，特定基因的启动子打开和关闭，以产生突发的mRNA分子(黑色波)。然而，在细胞周期的S期和G2期(绿色)，细胞中每个基因的副本是复制的两倍。因此，细胞必须补偿这些额外的拷贝，以便在整个细胞周期中保持恒定的mRNA水平。Skinner等人发现，这种补偿是通过减少每小时启动子打开的次数来实现的。维持一个恒定水平的信使rna的挑战类似于在一辆汽车中维持一个恒定速度的挑战，其中有另一个司机定期出现。

使用一种称为单分子信使rna荧光原位杂交的技术，斯金纳et al .,贝勒医学院,莱斯大学,伊利诺伊大学香槟分校和伊坎在西奈山医学院——在单独的老鼠胚胎干细胞中量化的成熟水平和新转录信使rna分子。他们还测量了每个细胞含有多少DNA，并以此来计算出细胞处于细胞周期的哪个阶段。这些测量允许在基因复制之前和之后分析每个特定基因拷贝的活性。

Skinner等人发现在胚胎干细胞中重要的两个基因Nanog和Oct4以“突发性”的方式转录。此外，这种活性在DNA复制后减少，以补偿细胞中额外的DNA，从而在细胞周期中平衡mRNA的水平。类似的“剂量补偿”机制改变转录速率后，DNA复制已在其他生物系统中被记录。

可以通过降低三个参数来实现剂量补偿：

• 每小时启动子打开的频率(脉冲频率);
• 开放多久
• 以及信使rna从开放启动子产生的速率(裂解大小)。为了确定这些参数中哪些与小鼠干细胞相关，Skinner等人开发了一个模型，该模型可以解释基因拷贝数在细胞周期中发生变化的事实。

在一个简化的模型中，脉冲频率是细胞周期中唯一改变的参数，结果与实验数据非常吻合。这意味着，细胞似乎是通过减少每个拷贝切换到开放状态的频率来补偿基因复制，而不是通过减少开放状态时产生的mrna数量。Skinner等人还发现，在他们研究的两个基因Oct4和Nanog中，脉冲频率是差异最大的参数。这些结果加强了最近在哺乳动物细胞中的其他研究发现，即通过突变频率的改变来调节转录率的重要性。这种调节可以在整个细胞周期中发生(Padovan-Merhar et al.， 2015)，响应转录因子水平(Senecal et al.， 2014)，或在肝脏基因中响应代谢刺激。

细胞周期给细胞带来了额外的挑战。首先，细胞分裂会立即使mRNA的产量减半，因此mRNA的水平必须提高来为这一事件做准备。这是通过一种机制来实现的，该机制监测细胞的体积，并相应地增加爆发的大小(Padovan-Merhar等，2015)。此外，不同基因之间的mRNA降解率差异很大，这对burst参数及其补偿机制造成了不同的约束。这与开车有相似之处：在土路上要保持匀速行驶要比在高速公路上困难得多。

最后，有些情况下，细胞放大mRNA水平的变异性，从而产生细胞间的差异。细菌应激反应就是一个例子，其中关键基因的不同mRNA水平有助于改变反应(Eldar和Elowitz, 2010)。干细胞分化策略也可能会根据细胞周期所达到的阶段而改变(Pauklin和Vallier, 2013)。斯金纳等人现在提供了一种工具，用于识别在整个细胞周期中对任何感兴趣的基因和生物系统都不同的脉冲参数。

最后，我们注意到在数据分析过程中，如何判断细胞周期的影响呢？一般的方法是在PCA空间或者umap空间中看看是不是有因周期而聚集的群：

DimPlot(pbmc3k.final,
        reduction = "pca",
        group.by= "Phase")

这个是没有啊，如果某一群都是某周期的细胞就要考虑了。如果真的不想把细胞周期考虑在分析框架内，可以在选择高变基因的时候去除周期相关基因，以及在标准化的时候回归掉，如：

pbmc3k.final <-  ScaleData(pbmc3k.final,vars.to.regress = c("S.Score", "G2M.Score"))

http://shalevlab.weizmann.ac.il/wp-content/uploads/2016/03/shani_elife.pdf