ixxmu
commented
1 month ago 单细胞细胞周期矫正分析流程学习(Seurat)以及关于是否应该矫正的思考 by 生信方舟
单细胞细胞周期矫正的目的从宽泛的角度来说是因为它可以消除由于细胞周期阶段不同而导致的基因表达异质性。在单细胞数据中，不同细胞可能处于不同的细胞周期阶段(如G1、S、G2或M期)，这些阶段会影响细胞的基因表达模式，尤其是那些与细胞周期密切相关的基因。如果不进行细胞周期矫正，这些周期性变化的基因可能会误导后续的数据分析，如聚类和差异表达分析，导致分析结果不能准确反映细胞的生物学状态和异质性。
但这个时候问题就来了，如果可以消除由于细胞周期阶段不同而导致的基因表达异质性，那么换句话来说是不是也就可能消除了本身的生物学差异，而这种差异也许是研究者正是所想要研究的呢？
笔者认为正确也不一定正确，这需要依情况而定。
我们在做单细胞分析的时候是需要根据样本的"特征"(这里可以泛指是各种重要的且表达有差异的基因)去对数据做分析。那么如何更好的去找出这些特征基因且必须把一些表达有差异但其实是没有意义的基因给筛去？在刚开始的质控环节就已经开始了，比如限制线粒体基因的百分比，nCount和nFeature的值等等。
在单细胞时代之下数据的分辨率已经达到了细胞级，不同细胞应看做不同的样本，那么它们之间是存在极大异质性的。这种异质性从表象上都可归类于生物学差异，但更深层次的，研究者真正想要知道是生物学差异背后的样本之间的功能差异！
而细胞周期既可以看做是生物学上的混杂因素(不同细胞的细胞周期情况在取样的那一瞬间就已经被定格了，而他们的周期情况绝对是不可能完全一样的)，也可以在控制大部分混杂因素(控制变量这件事要完美是不可能的)后的后续分析中认为是功能上的差异(不同处理下，细胞周期通路在实验组和对照组之间存在的差异)。如果说是生物学的差异，我们就需要在正式分析之前就要进行矫正，尽量减少这种生物学差异导致的影响！
但这里还有一个问题，细胞周期的影响真的那么大吗？是的，在既往的文献中(见参考资料)已经明确提到了细胞周期这个生物学差异会导致很大的异质性，从笔者的理解来看，这种异质性可能会掩盖功能学上的差异！(这个掩盖可以解释为由于细胞周期导致的生物学误差过大之后，其他非细胞周期相关的功能学差异反而显得没那么有差异了..)
当然，我们还是要依情况而定，如果研究方向就是跟细胞周期有直接且密切关系的，那就不可以去矫正这个因素了，比如现在热门—衰老/干细胞(涉及静息态等情况？这块不太懂)等相关的研究。如果是非直接关系的，应当要去矫正，换个角度去想，如果一个矫正能把"所谓的差异"给矫正没了，那请问这个"差异"到底真有那么大的研究价值吗？
步骤流程
1、导入
sce.all是已经做过质量控制之后的数据
rm(list=ls())
sce <- qread("./Data.qs")
dim(sce.all)
2、细胞周期分析
# 标准化数据
sce <- NormalizeData(sce,
normalization.method = "LogNormalize",
scale.factor = 10000)
# 找到高变基因
sce <- FindVariableFeatures(sce, selection.method = "vst")
# ScaleData数据
all.genes <- rownames(sce)
sce <- ScaleData(sce,features = all.genes)
# s和g2m期基因和数据集的基因取交集
s.genes <- intersect(cc.genes$s.genes,rownames(sce))
g2m.genes <- intersect(cc.genes$g2m.genes,rownames(sce))
# 细胞周期分析
sce <- CellCycleScoring(sce,
s.features = s.genes,
g2m.features = g2m.genes,
set.ident = T)
# PCA可视化
sce <- RunPCA(sce,features = c(s.genes,g2m.genes))
PCAPlot(sce,group.by = "Phase")
# S.Score评分
p1 = VlnPlot(ch1,"S.Score",pt.size = 0,group.by = "Phase");p1
# G2M.Score评分
p2 = VlnPlot(ch1,"G2M.Score",pt.size = 0,group.by = "Phase");p2
wrap_plots(p1,p2,nrow = 1)
ggsave("Cell_cycle_analysis.png",width = 18,height = 6,dpi = 300)
这张图主要看数据是不是集中在0的附近，从图片看还是比较集中的。下图是官网上的示例数据，左边是正常的，右边是不正常的！如果出现右边这样的就一定需要矫正哦~
3、矫正之后检查数据
# 根据上面的结果决定是否要去矫正
# 重新ScaleData这次用S.Score/G2M.Score数据去矫正
sce <- ScaleData(sce,
vars.to.regress = c("S.Score", "G2M.Score"),
features = all.genes)
# check一下
# PCA可视化
sce <- RunPCA(sce,features = c(s.genes,g2m.genes))
PCAPlot(sce,group.by = "Phase")
# S.Score评分
p1 = VlnPlot(ch1,"S.Score",pt.size = 0,group.by = "Phase");p1
# G2M.Score评分
p2 = VlnPlot(ch1,"G2M.Score",pt.size = 0,group.by = "Phase");p2
wrap_plots(p1,p2,nrow = 1)
ggsave("Cell_cycle_analysis_after.png",width = 18,height = 6,dpi = 300)
关键代码是 sce <- ScaleData(sce,vars.to.regress = c("S.Score", "G2M.Score"), features = all.genes)。
其他接下来就是常规流程了，harmony整合，注释等等。
参考资料：
1、Seurat ：https://satijalab.org/seurat/articles/cell_cycle_vignette
2、A single-cell resolution map of mouse hematopoietic stem and progenitor cell differentiation. Blood. 2016 Aug 25;128(8):e20-31.
3、Complex Analysis of Single-Cell RNA Sequencing Data. Biochemistry (Mosc). 2023 Feb;88(2):231-252.
4、Single-cell RNA-seq reveals changes in cell cycle and differentiation programs upon aging of hematopoietic stem cells. Genome Res. 2015 Dec;25(12):1860-72. doi: 10.1101/gr.192237.115
5、生信星球：https://mp.weixin.qq.com/s/bJ0tNj4w5p4R1MX5RHG9Ng
6、观科研：https://mp.weixin.qq.com/s/4KtwcV_eJnkdLfO47Yo_XA
7、Top生物信息：https://mp.weixin.qq.com/s/tjP4sOXxGswxU8aOY00rjg
8、老俊俊的生信笔记：https://mp.weixin.qq.com/s/i5ByW8L_X759hVoAYCWEsg
9、生信百宝箱：https://mp.weixin.qq.com/s/kgNUodM18cCUls7ZQXYgEQ
致谢：感谢曾老师、小洁老师以及生信技能树团队全体成员。
注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟
- END -
https://mp.weixin.qq.com/s/mlAhviD-7cyi891Xob1x3A