为R包写一本书（像Y叔致敬）

Posted on by

最近收到粉丝求助，本来是询问富集分析的时候，我给的参考代码，为什么pvaluecutoff值和qvaluecutoff值设到0.9，其实是怕大家富集不到结果。然后Y叔在自己的微信公众号中提到“富集不到结果才是正确的结果”，采用了更加稳妥和可靠的方法来判断富集结果，而粉丝的数据在DAVID中能有结果，可在Y叔的包里，结果就少了一些，如何决定采取哪个？最后又讨论到DAVID结果可视化，网上资源少，他只能做成条图，需要我给指条方向！

这个时候我还是推荐了Y叔的clusterProfiler ，就去找了找其官网，的确可视化方法又多了几个：

  • barplot
  • cnetplot
  • dotplot
  • emapplot
  • gseaplot
  • goplot
  • upsetplot

好几个都是以前没有介绍过的，有趣的是我准备浏览这些可视化函数的帮助文档的时候，看到了这样的话：Please go to https://yulab-smu.github.io/clusterProfiler-book/ for the full vignette.

重点来了，Y叔特意为其包写了一本书来介绍其用法。

镇书之图

打开Y叔的书，映入眼帘的就是流程图：

image-20191210114608881

可以说是，跟着这些函数和数据库，都可以入门生物信息学数据分析了。

做完差异分析

GEO数据挖掘代码，很容易得到上下调基因，而且转为ENTREZID，后续分析都以这个为主线。

load(file = 'deg.Rdata')
head(deg)
## 不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。
logFC_t=1.5
deg$g=ifelse(deg$P.Value>0.05,'stable',
 ifelse( deg$logFC > logFC_t,'UP',
 ifelse( deg$logFC < -logFC_t,'DOWN','stable') )
)
table(deg$g)
head(deg)
deg$symbol=rownames(deg)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
df <- bitr(unique(deg$symbol), fromType = "SYMBOL",
 toType = c( "ENTREZID"),
 OrgDb = org.Hs.eg.db)
head(df)
DEG=deg
head(DEG)

DEG=merge(DEG,df,by.y='SYMBOL',by.x='symbol')
head(DEG)
save(DEG,file = 'anno_DEG.Rdata')

gene_up= DEG[DEG$g == 'UP','ENTREZID'] 
gene_down=DEG[DEG$g == 'DOWN','ENTREZID']

或者你提供其它方法，拿到如下的数据，会R语言的都知道是简单的向量罢了。

image-20191210111843553

如果你确实差异分析都不会，就是拿到上下调基因集有困难，看：多分组差异分析续集 里面的教程哈。

最简单的超几何分布检验

这里就拿KEGG数据库举例吧，拿自己判定好的上调基因集进行超几何分布检验，如下：

gene_down
gene_up
enrichKK <- enrichKEGG(gene = gene_up,
 organism = 'hsa',
 #universe = gene_all,
 pvalueCutoff = 0.1,
 qvalueCutoff =0.1)
head(enrichKK)[,1:6] 
browseKEGG(enrichKK, 'hsa04512')
dotplot(enrichKK)
enrichKK=DOSE::setReadable(enrichKK, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
enrichKK

得到的enrichKK是Y叔定义好的对象，这个对象值得大家用心摸索一下，毕竟Y叔后续的可视化函数，都是基于此。

最基础的条形图和点图

其实就是超几何分布检验的结果，如何可视化，毕竟大家都喜欢看图

image-20191210112432726

上面的表格进行可视化，可以是最基础的条形图和点图，使用Y叔的可视化函数即可，作用于其定义好的对象。

#条带图
barplot(enrichKK,showCategory=20)
#气泡图
dotplot(enrichKK)

这个没有啥悬念，虽然也是可以调整一些参数，可视化上面的表格里面的其它结果

image-20191210112716630

一般来说，我们会喜欢点图一点。

通路与基因之间的关系可视化

这里有两种形式，同一个函数，不同的参数即可：

cnetplot(enrichKK, categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
cnetplot(enrichKK, foldChange=geneList, circular = TRUE, colorEdge = TRUE)

这个就不好说大家会更喜欢哪个了

image-20191210112902619

展现的是不同的通路，共有基因情况，可以看到有一些基因是可以属于不同的通路的。

通路与通路之间的连接展示

其实上面的通路与基因之间的关系可视化已经告诉我们了，有一些通路是可以共享一些基因的，这样的话，这些共享基因的通路的相似性就会高一点，它们可以看做是通路模块。

所以就有可视化函数来具体的查看通路与通路之间的连接，不需要展示基因了，代码如下：

emapplot(enrichKK)

可以看到， p53通路和细胞周期通路，连接起来了，说明它们两个共有基因数量多。

而 toll-like受体，IL-17等，都是免疫相关通路，所以连接起来了。

image-20191210113223688

当然了，Y叔自己的例子更漂亮，而且多元化，毕竟参数都是他定义的。

image-20191210113423925

热图展现通路与基因之间的关系

前面的cnetplot就是展现通路与基因之间的关系，但是它的特色的炫酷，其实看起来很烧脑，如果仅仅是为了展现通路与基因之间的关系，热图更靠谱，代码如下：

heatplot(enrichKK)

image-20191210113843661

如果你仔细看这个图，就会发现，它不仅仅是通路与基因之间的关系，还包含了通路与通路之间的连接关系，如果两个通路共有的基因比较多，当然在emapplot里面，就会连接在一起。

我个人认为，其实 heatplot是最强大的，但是呢， 没有cnetplot和emapplot炫酷，而barplot和dotplot就太朴素了。最后，如果你是做GO数据库呢，其实还有一个goplot可以试试看，当然是以Y叔的书为主啦。

图都可以可以修饰的

对ggplot语法有所了解的朋友都知道，其产出的图都是可以叠加各种属性的。同理，Y叔的这些可视化函数，也可以继续映射基因的其它属性，比如变化倍数等等。

如果你还不知道Y叔是谁

南方医科大学教授博导，看他的招博士通知吧：常年招硕博士、博士后一起搞生信或被生信搞

​

Comments are closed.