刘小泽写于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的转换有兴趣，可以期待明天的推送~
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com
Welcome to our bioinfoplanet!
网友评论