我这个公众号不像大多数所谓的生信从入门到精通的各种其实只是搬运点入门教程的群众喜闻乐见的公众号。正如我在《为什么要开这个公众号》里说的,这是小众的,有个人色彩的各种原创文。我不可能像其它公众号一样招两小弟当客服,很多人在公众号后台向我扔了许多问题,由于个人精力有限,只能优先解答「知识星球」的问题,上次写的《同一数据多变量分组的boxplot?》,图虽然简单,却穿着好多件马甲,而我把它扒光了给你看🙈

上面这个图,你看着高大上吧,我都可以吐它一脸口水。每一个有灰色背景的图,在x轴上violin都够到边界了,其实所有的violin都够到了,这证明什么?每个violin之间其实不可比较!你能想像几个独立的数据,在统一的bin width情况下,画density curve,竟然最高点都一样高吗?显然可能性几乎为0。这个如果使用ggplot2的话,可以使用scale='width'强制拉成一样高,但我不推荐,正如我前面说的,不可比较了。默认参数scale='area',积分面积一样,和density curve一样解析,另外的参数scale="count",高度与计数同比例,和histogram一样解析,而scale='width'强制拉成一样高,如果没有在显眼处说明,误导性太强。

画这种图也可以手工拼,这样就简单了。在你需要的情况下,加个灰色背景嘛,最后拼图嘛。当然拼图不一定要在illustrator里拼,比如你用grid,先画好坐标轴,然后水平上定义几个一样大的viewport,每一个violin都画在相应的viewport里面,对于画图函数来说,viewport就是整个画布了(虽然只是画布里的一块区域),所以你要么画violin,要么在画之前先画个矩阵,一路画下来,代码可以直接生成这样的图,但这图每一个violin都是独立画的(当然也不是完全独立,每一个水平上的ylim是有统一的),就算代码一步生成,也跟illustrator里拼没两样。


这里我要教你用ggplot2自动生成,其实解决思路早已推送过,请看《facet_plot:加图层到特定分面,方法二》,也正如我在《什么!你的图上有一双看不见的手》里说的,你们以为我在教ggtree,其实同时在教ggplot2。

Continue reading

这个问题,其实答案就存在于《听说你还不会画热图》,我们先生成一个矩阵:

set.seed(2017-11-12)
d = data.frame(matrix(rnorm(100), ncol=10))
colnames(d) = paste0('t', 1:10)
rownames(d) = paste0('g', 1:10)

这个矩阵是rnorm生成的随机数,有正有负,我们再生成一个矩阵,只有正数,并且数值上比第一个矩阵要大:

d2 = abs(d) * 1.2

Continue reading

掐架的额外收获

你昨天才做的分析,可能是几年前的结果!》这篇文章给大家敲了警钟,各种各样的web-server,要小心看有没有维护更新,有些是五年十年都不更新的,十分可怕。文章虽然讲的是富集分析,但其它分析工具你同样需要小心。

当然并不是说独立的软件/软件包就一定靠谱,如果软件自己打包了数据,同样要注意数据是否有更新,而如果数据不打包在软件里,而是在线获取,你同样也该留一下心。这也是clusterProfiler做富集分析的优势所在,KEGG数据是在线的,永远是最新的,而GO的数据不在软件包里,而依赖于别的数据包,而这些数据包是社区维护的(相对而言,个人的维护比较难以为继),就确保了数据一直在有更新维护的。

Bioconductor每半年发行一次,注释包同样每半年更新一次,所以你用clusterProfiler做GO分析,你用的GO数据库不会说超过半年没更新,而不像有些公司给出的结果,落后于这个世界不是一年两年这么简单。

Continue reading

听说你想打水印?

正如我在《用R画meme ;)》里说的,meme包就可以拿来画,当然这个包是设计来画meme的,所以文本只会打在中间,如果你想在右下角之类的,不好意思,我没提供选项让你调文本位置,氮素,meme依然是你的好帮手,怎么说呢,meme不单能帮你读图片,还能用grid画出来,这意味着你可以用grid.text随便加个文本就OK了。

实例演示一发:

library(meme)
u <- system.file("angry8.jpg", package="meme")
x <- meme(u)
print(x)
grid.text("biobabble", x=.9, y=.05, 
	gp=gpar(col='white', cex=1.2))

Continue reading

正如我在「知识星球」里说的,scholar这个包现在已经死了,因为作者不维护了,除非出现接盘侠。我自己fork了一个,打了补丁,现在是可用的。这个包应该说我用得蛮多的,你看我简历里的引用数柱状图,我不会无聊到过一段时间就去查,然后重新画图的,太low。都是自动化的,我到简历目录下打make指令,就会自动去爬数据、重新画图并生成简历,而我啥都不用干。

Continue reading

github上的问题,问了两个问题,这是其中第二个:

Meanwhile, when I fed the ENTREZID to enrichKEGG, it show me two unreasonable results:

kegg_enrich = enrichKEGG(
	gene = new_ids$ENTREZID,
	organism = 'hsa',
	keyType = 'ncbi-geneid',
	pvalueCutoff = 0.05,
	#pAdjustMethod = 'BH',
	#qvalueCutoff = 0.2,
	use_internal_data = FALSE
)
head(kegg_enrich)
              ID                           Description GeneRatio  BgRatio       pvalue   p.adjust     qvalue
hsa04380 hsa04380            Osteoclast differentiation    14/281 128/7299 0.0003817040 0.04767338 0.04508266
hsa04070 hsa04070 Phosphatidylinositol signaling system    12/281  99/7299 0.0003875885 0.04767338 0.04508266

I noted that only 281 genes are remained(there are 700+ genes in my list). In case that there is something wrong wtihin my gene list, I also tried my list with DAVID. It gave me reasonable results. So this is my second question, why enrichKEGG cannot recognize my geneids?

Continue reading

在很久很久以前,当我还在用perl的时候,有一个TGen-EUtils的perl模块,可以去爬PubMed,我喜欢用它来画所谓的pubmed趋势,比如说我在2012年的时候,就画了各个组学技术的pubmed趋势,通过每一年的文章发表量来看某个技术到底火不火,走势怎么样。

当年microarray还特别猛,但显然可以看到二代测序必须要火,也可以看到proteomics在吹了几年牛B之后,有回落的趋势。

Continue reading

Author's picture

Guangchuang Yu

a senior-in-age-but-not-senior-in-knowledge bioinformatician

Postdoc researcher

Hong Kong