在《CS7:Genomic coordination的富集性分析(1)》说到了seq2pathway这个包,其实是两部曲,seq2gene->gene2pathway,无非是把测序片段用临近的基因注释,包括和TSS overlap的基因,宿主基因,上下游的基因等,然后拿这些基因跑ORA,做富集,仅此而已,这个包支持的物种极有限,《CS4:关于ChIPseq注释的几个问题》这一文中讲到ChIPseeker支持所有有基因组注释的物种,而《clusterProfiler for enrichment analysis》也支持所有物种(即使你自己跑的电子注释,也能支持),那么使用ChIPseeker来做基因注释,然后衔接clusterProfiler就可以支持所有物种的测序片段进行功能富集分析了。

CS3: peak注释》本身就支持几种注释,另外我写了一个seq2gene的函数,套用seq2pathway的思路,把一个基因位置上所有关联的基因全部返回来,我们可以使用它去把基因位置信息转换成基因列表,然后用于富集分析,还是熟悉的味道,还是熟悉的配方🦄

Continue reading

Reassortment is an important strategy for influenza A viruses to introduce a HA subtype that is new to human populations, which creates the possibilities of pandemic.

A diagram showed above (Figure 2 of doi:10.1038/srep25549) is widely used to illustrate the reassortment events. While such diagrams are mostly manually draw and edit without software tool to automatically generate. Here, I implemented the hybrid_plot function for producing publication quality figure of reassortment events.

library(tibble)
library(ggplot2)


n <- 8

virus_info <- tibble(
    id = 1:7,
    x = c(rep(1990, 4), rep(2000, 2), 2009),
    y = c(1,2,3,5, 1.5, 3, 4),
    segment_color = list(
        rep('purple', n),
        rep('red', n),
        rep('darkgreen', n),
        rep('lightgreen', n),
        c('darkgreen', 'darkgreen', 'red', 'darkgreen', 'red', 'purple', 'red', 'purple'),
        c('darkgreen', 'darkgreen', 'red', 'darkgreen', 'darkgreen', 'purple', 'red', 'purple'),
        c('darkgreen', 'lightgreen', 'lightgreen', 'darkgreen', 'darkgreen', 'purple', 'red', 'purple'))
)

flow_info <- tibble(from = c(1,2,3,3,4,5,6),
                    to = c(5,5,5,6,7,6,7))
hybrid_plot(virus_info, flow_info)

Continue reading

取子集对于进化树可视化来说是非常常见的,我们要区分内部节点和外部节点,我们也可能想针对某些特点的节点进行注释。

ggplot2现在所有图层都不支持直接取子集,所以呢ggtree就自己定义了一些修改的图层,包括geom_text2, geom_label2, geom_point2geom_segment2,这些图层和ggplot2的版本唯一差别就是支持取子集。这样对于我们做注释来说,就更方便了。

比如说我想给内部节点打点,可以用:

ggtree(tree) + geom_text2(aes(label=label, subset=!isTip), hjust=-.2) +
        geom_point2(aes(subset=!isTip), color="red", size=3)

Continue reading

之前有写过读文献的文章,也介绍了如何获取文献,那是我还在读硕士的时候写的,覆盖了找文献、文献管理、下载文献、订阅期刊等方方面面,说到下文献,最大的神器当然是sci-hub,但是sci-hub的域名sci-hub.cc, sci-hub.bz, sci-hub.io纷纷都沦陷,最近国内有一神器geenmedical,整合了影响因子、引用数等指标,然而这一神器的下文献功能只是提供外链,有一种说法是geenmedical=pubmed+scihub,当然geenmedical的资源还可能来自于百度学术、researchGate等,不都是来自于sci-hub,然而我们必须承认的是最大的文献提供来源必须还得是sci-hub,所以如果sci-hub挂了,基于sci-hub的神器也就神奇不起来了。如果有人跟你说忘记sci-hub吧,用geenmedical,那绝对是标题党!这里并没有说geenmedical不好,它整合了多方面的资源并且本地化得非常好,很适合国内使用。

比如我在geenmedical里搜clusterProfiler,给出的下载链接第一个就是sci-hub的:http://doi.org.sci-hub.cc/10.1089/omi.2011.0118,然而这个链接是不可用的!当然它给出了比如百度学术的链接是可以下载到文献的。然而还是那句话,论「海盗精神」我只服sci-hub,虽然我之前也介绍了《sci-hub如果挂了,你还有神器下文献吗?》,但任何下文献的神器,在sci-hub面前都是渣渣。

Continue reading

我这个公众号不像大多数所谓的生信从入门到精通的各种其实只是搬运点入门教程的群众喜闻乐见的公众号。正如我在《为什么要开这个公众号》里说的,这是小众的,有个人色彩的各种原创文。我不可能像其它公众号一样招两小弟当客服,很多人在公众号后台向我扔了许多问题,由于个人精力有限,只能优先解答「知识星球」的问题,上次写的《同一数据多变量分组的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

Author's picture

Guangchuang Yu

Bioinformatics Professor @ SMU

Bioinformatics Professor

Guangzhou