查看原文
其他

拟时序分析的热图提取基因问题

生信技能树 单细胞天地 2022-06-07


昨天我在单细胞天地讲解了使用monocle2进行拟时序分析的方法,基本上跟着我的代码走一波就可以学会了,当然具体参数理解需要自行发力哦,见:使用monocle做拟时序分析(单细胞谱系发育) 用法只是最基础的知识而已,更多的时候,我们需要活学活用,比如课程学员提到的问题,就是因为做不到活学活用,他想知道下面的拟时序分析的热图提取基因,学员把基因按照发育顺序绘制了热图,而这些基因被他分成了3组,想拿基因去做GO/KEGG等数据库进行功能注释,不知道如何获取基因名字。

我这里不能拿学员真实项目数据来演示,所以还是用我们的老朋友,拿scRNAseq包的表达矩阵测试,见:使用monocle做拟时序分析(单细胞谱系发育)

首先根据细胞发育谱系来绘制热图

因为前面的教程 使用monocle做拟时序分析(单细胞谱系发育) 我们已经把细胞发育情况做出来了,就是NPC细胞跟另外3种细胞从生理上就不一样,所以是单独的发育轨迹,而 “GW16” and “GW21” ,“GW21+3” 这种孕期细胞,就可以很清晰的看到时间被反映在我们的拟时序分析结果了。

大家可以重温教程  使用monocle做拟时序分析(单细胞谱系发育)  里面的4个绘图代码:

plot_cell_trajectory(cds, color_by = "Biological_Condition")  
# 可以很明显看到细胞的发育轨迹 
plot_cell_trajectory(cds, color_by = "State")
plot_cell_trajectory(cds, color_by = "Pseudotime")
plot_cell_trajectory(cds, color_by = "State") +
  facet_wrap(~State, nrow = 1)

我们前面构建细胞发育谱系,使用的是不同Biological_Condition的细胞类型之间使用monocle找到的两千多个基因。

简单的一个函数就可以绘制热图:

plot_pseudotime_heatmap(cds[ordering_genes,],
                        num_clusters = 3,
                        cores = 1,
                        show_rownames = T)

可以看到, 这个图就和学员求助的图一模一样啦,因为基因数量两千多个,所以画出来肯定是看不清晰的啦。


既然基因被分成了3组,想拿基因去做GO/KEGG等数据库进行功能注释,就需要获取基因名字。

这个做不出来,不怪学员,因为正常人很难想到,这个绘图函数其实是可以调整返回数据对象的,而且官网例子也没有提到。

需要看函数的帮助文档,如下:

hclust_method    
The method used by pheatmap to perform hirearchical clustering of the rows.
hclust_method = "ward.D2"
return_heatmap    
Whether to return the pheatmap object to the user.

很明显,这个函数其实就是pheatmap的一个包装罢了,本质上也是调用  hclust 而已,使用的是ward.D2距离。

然后解析热图函数返回对象

根据帮助文档,我们修改参数,这样monocle的plot_pseudotime_heatmap函数就有返回值了,是一个对象。

p=plot_pseudotime_heatmap(cds[ordering_genes,],
                        num_clusters = 3,
                        cores = 1,return_heatmap=T,
                        show_rownames = T)

从pheatmap的对象里面提取基因名字就很简单了,就是在p$tree_row里面

> p$tree_row

Call:
hclust(d = d, method = method)

Cluster method   : ward.D2 
Number of objects: 2200 

就可以拿到基因名对应的cluster啦,代码如下:

clusters <- cutree(p$tree_row, k = 3)
clustering <- data.frame(clusters)
clustering[,1] <- as.character(clustering[,1])
colnames(clustering) <- "Gene_Clusters"
table(clustering)

实际上学员提问是有问题的

因为学员之间丢出这个热图,然后咨询如何在图片里面提取基因名字,所以大家只能是问是pdf还是png的图片呢?是不是可以AI或者PS解析它,拿到基因名字呢?

如果学员是直接问:使用monocle的plot_pseudotime_heatmap函数绘制的热图里面的基因聚集成为3类,该如何提取基因名字,其实就很简单了。

这个答疑虽然是给基础课程的学员,但实际上这个内容在单细胞转录组的进阶课程,如果学了就不需要提问了。

今天是教师节,我放一个收款二维码吧,如果你曾经受益于我的教程,也算是我的学生,意思一下即可!

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

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