153-”掉包“某个包装好的绘图函数

Dec 11, 2019 4 min read cnposts

刘小泽写于19.12.11 写这一篇的目的是考虑到使用某些R中包装好的函数得不到想要的结果，应该怎么办？

前言

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

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

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

例如

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

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

很多时候，我们想用这个图来对数据有一个初步的认识（相当于一个质控的过程），但是如果数据的基因名都是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结果很全，但转换速度可能会比较慢，直接转换几万个基因很可能会和数据库断线。如果对大量的基因ID的转换有兴趣，可以期待明天的推送~

R tips

153-”掉包“某个包装好的绘图函数