查看原文
其他

狸猫换太子-"掉包"某个包装好的绘图函数

豆豆花花 生信星球 2022-06-07

 今天是生信星球陪你的第502天


   大神一句话,菜鸟跑半年。我不是大神,但我可以缩短你走弯路的半年~

   就像歌儿唱的那样,如果你不知道该往哪儿走,就留在这学点生信好不好~

   这里有豆豆和花花的学习历程,从新手到进阶,生信路上有你有我!

豆豆写于19.12.11
写这一篇的目的是考虑到使用某些R中包装好的函数得不到想要的结果,应该怎么DIY?

前言

之前改造过Seurat包的DoHeatmap函数:https://www.jianshu.com/p/b2816a72e79e

实际上是把数据重新拿来自己作图

这次来看看,如何半路”掉包“作者定义的函数

例如

在使用Scater包的时候,其中一个作图函数叫:plotHighestExprs
关于Scater的学习:https://www.jianshu.com/p/d3bf2a0bea6e

这个函数操作的是一个single cell experiment(sce)对象,然后它会帮我们找出表达量变化最大的50个基因,并对它们作图。

image-20191211115642938

很多时候,我们想用这个图来对数据有一个初步的认识(相当于一个质控的过程),但是如果数据的基因名都是Ensembl ID,结果是非常不清楚的,因为还是不能直观判断这些Ensembl ID对应哪个symbol ID。

所以需要ID转换:ID转换失败,爬起来继续high

但是问题来了

如果我的数据有4万行(也就是4万多个基因),而且使用org.xxx.db转换的结果不理想,于是利用Ensembl的biomaRt 。难不成要将这4万多个基因都进行转换,然后再做这个只需要top50的图片吗?

虽然利用biomRt转换是可行的,并且我也想出来了一个”转换速度慢“的可行的解决办法,可以看明天的推送

但是就目前的问题来讲,显然有点”花大力气办小事情“

换一种思路

思考🤔这个函数到底做了什么事?【这个思考过程对于所有函数有效】
step1:拿到sce对象的矩阵
step2:拿到top50基因
step3:根据这些基因用它的方法去作图

既然这个函数能帮我们找到这top50高变化的基因,那么我们能不能就从step2这里,对这个函数”掉包“把它的top50拿出来,然后进行一个ID转换。把这50个Ensembl ID转成Symbol ID后,重新加到这个函数中,继续进行后续的作图

那么好,第一步:查看这个包装好的函数代码

方法就是:找到这个函数,然后按住ctrl / cmd,同时在这个函数上点一下鼠标,就跳转到这个函数的具体代码上了【不过不排除有的函数封装的密闭性太好,以至于这样找不到它的源代码】

因为函数功能很多,所以代码长是很正常的,只有自己不被吓到就好

第二步 根据刚才的思路,去找对应的代码

例如:step1:拿到sce对象的矩阵

就会找到:

exprs_mat <- assay(object, exprs_values, withDimnames = FALSE)这样一行

这时,就可以把我们自己的sce对象赋值给object,另外还有开头定义的一些默认参数,也一起定义一下(保证后续代码可以顺利运行)

就这样一步步向下走,并时刻注意哪里可能得到top50

当看到这里的时候,就要警惕了,可能chosen就是得到的50个基因名,因为这个n是默认参数50(不就是对应默认提取的top50基因吗?);另外下面的也都是对矩阵取子集的过程

但是,这仅仅是可能,还需要分辨:这个chosen是位置还是名称?因为矩阵取子集可以根据三样东西:逻辑值、位置、名称

如果不是名称的话,那么还要继续向下走

直到出现了对feature_names取子集,得到sub_names

第三步 得到了想要的东西,去”掉包“换个结果

下面就是根据之前的ID转换过程,对这个sub_names 进行ID转换,将其他ID换成Symbol ID

第四步 尊重这个函数,我们只修改我们的数据

记住:不要修改别人写好的函数
因为一旦它不高兴,它就再也不给你干活了。而且与其一行一行修改别人的代码,还不如从源头开始,把我们自己的数据修改好,反正我们知道最后函数还是会对我们修改的内容进行操作

我们上面得到了top50基因对应的Symbol ID,那么我们最好是将我们数据中这些位置的这些基因名先换成Symbol ID,然后还是提交我们的数据给函数。

# chosen是top50基因的位置信息
# new_ids是sub_names利用biomart转换的结果

# 【重点】需要将ids的Ensembl与原始的sub_names位置对应起来
new_ids <- ids[match(sub_names,ids$ensembl_gene_id),]

rownames(sce[chosen] <- new_ids$hgnc_symbol

随后依然使用:plotHighestExprs(sce, exprs_values = "counts") ,什么也不用改

但实际上,我们已经完成了”掉包“的工作,函数所操作的,就是那50个包裹在Ensembl 海洋的Symbol ID

因此,最后的结果也是会显示Symbol ID了

小结

分享就到此结束了,上面👆的几步操作可以实际操作一下,然后自己也能顺利地完成”掉包“工作啦

另外,之前说过使用biomart结果很全,但转换速度可能会比较慢,直接转换几万个基因很可能会和数据库断线。
如果对使用biomart转换大量的基因ID有兴趣,可以期待明天的推送~


点击底部的“阅读原文”,获得更好的阅读体验哦😻

初学生信,很荣幸带你迈出第一步

🤓生信星球 🌎一个不拽术语、通俗易懂的生信知识平台

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

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