查看原文
其他

把基因表达量画在拟时序结果图上

生信技能树 生信技能树 2022-08-10

学员在练习单细胞时遇到了一个拟时序难题,如下所示的右图,因为它在monocle这样的拟时序R包的官方文档并没有介绍,而我以前没有这样的需求,所以就没有分享这方面笔记。

 

其实实现起来非常的简单,我觉得本质上仍然是对每个环节的R包的对象的熟悉程度。其实我一直强调,大家的基础课程学完后需要完成作业:单细胞基础视频课程结业考核20题 , 也就是熟练掌握5个R包,分别是: scater,monocle,Seurat,scran,M3Drop 需要熟练掌握它们的对象,而且分析流程也大同小异:

  • step1: 创建对象
  • step2: 质量控制
  • step3: 表达量的标准化和归一化
  • step4: 去除干扰因素(多个样本整合)
  • step5: 判断重要的基因
  • step6: 多种降维算法
  • step7: 可视化降维结果
  • step8: 多种聚类算法
  • step9: 聚类后找每个细胞亚群的标志基因
  • step10: 继续分类

前面的教程:拟时序分析就是差异分析的细节剖析,我们产生了  output_of_phe2_monocle.Rdata 文件,就是拟时序分析的结果, 就可以进行各种各样的可视化啦!

而且我们演示了3种可视化选择,包括:Cluster,Pseudotime,State ,请务必自行看前面的教程:拟时序分析就是差异分析的细节剖析,可视化的代码如下所示:

library(ggsci)
p1=plot_cell_trajectory(cds, color_by = "Cluster")  + scale_color_nejm() 
p1
ggsave('trajectory_by_cluster.pdf')
plot_cell_trajectory(cds, color_by = "celltype")  
p2=plot_cell_trajectory(cds, color_by = "Pseudotime")  
p2
ggsave('trajectory_by_Pseudotime.pdf')
p3=plot_cell_trajectory(cds, color_by = "State")  + scale_color_npg()
p3
ggsave('trajectory_by_State.pdf')
library(patchwork)
p1+p2/p3

两个单核细胞的亚群拟时序

其中,Cluster,以及,State 是离散变量,而Pseudotime是连续性变量,既然Pseudotime可以被可视化,那么基因表达量当然可以同样的可视化啊。

因为每个细胞都有一个推断好的Pseudotime值,同样的,每个细胞都有一个基因的表达量值。所以很简单的代码即可:

library(ggsci)
colnames(pData(cds))
pData(cds)$FCGR3A = log2( exprs(cds)['FCGR3A',]+1)
p1=plot_cell_trajectory(cds, color_by = "FCGR3A")  + scale_color_gsea()
pData(cds)$CD14 = log2(exprs(cds)['CD14',]+1)
p2=plot_cell_trajectory(cds, color_by = "CD14")    + scale_color_gsea()
library(patchwork)
p1+p2

如下所示:

表达量就跟拟时序的结果结合

表达量就跟拟时序的结果结合起来啦!

这样的基础认知,也可以看基础10讲:

最基础的往往是降维聚类分群,参考前面的例子:人人都能学会的单细胞聚类分群注释

文末友情推荐

做教学我们是认真的,如果你对我们的马拉松授课(直播一个月互动教学)有疑问,可以看完我们从2000多个提问互动交流里面精选的200个问答! 2021第二期_生信入门班_微信群答疑整理,以及 2021第二期_数据挖掘班_微信群答疑笔记

与十万人一起学生信,你值得拥有下面的学习班:


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存