查看原文
其他

使用barcodeplot可视化你的基因排序

生信技能树 生信技能树 2022-06-06

gsea分析这方面教程我在《生信技能树》公众号写了不少了,不管是芯片还是测序的表达矩阵,都是一样的,把全部基因排序即可:

  • 比如你有2万个基因,你根据自己的条件分组后算差异情况,根据差异把基因排序,然后看缺氧相关200个基因组成的集合在全部的排好序的2万个基因是散乱分布,还是集中于头部和尾部。
  • 当然了,基因集肯定不仅仅是缺氧这个生物学功能啦,在msigdb数据库有几万基因集合,其实生物学背景更重要。
  • 另外,基因的排序也不仅仅是条件分组后算差异来排序,也可以仅仅是表达量高低排序。

msigdb数据库网页里面有着丰富的基因集,MSigDB(Molecular Signatures Database)数据库中定义了已知的基因集合:http://software.broadinstitute.org/gsea/msigdb 包括H和C1-C7八个系列(Collection),每个系列分别是:

  • H: hallmark gene sets (癌症)特征基因集合,共50组,最常用;
  • C1: positional gene sets 位置基因集合,根据染色体位置,共326个,用的很少;
  • C2: curated gene sets:(专家)校验基因集合,基于通路、文献等:
  • C3: motif gene sets:模式基因集合,主要包括microRNA和转录因子靶基因两部分
  • C4: computational gene sets:计算基因集合,通过挖掘癌症相关芯片数据定义的基因集合;
  • C5: GO gene sets:Gene Ontology 基因本体论,包括BP(生物学过程biological process,细胞原件cellular component和分子功能molecular function三部分)
  • C6: oncogenic signatures:癌症特征基因集合,大部分来源于NCBI GEO  发表芯片数据
  • C7: immunologic signatures: 免疫相关基因集合。

如下所示的一个经典的gsea图:

文献:Decoding Human Megakaryocyte Development

它就很好的说明了这个 innate immune response基因集,是在排序好的2万多个基因里面是左偏的,也就是说这个通路是激活的。

但是这个gsea分析对绝大部分初学者来说理解起来很困难,而且代码实现难度也不小。最近看教程,发现其实limma包就有一个barcodeplot函数,其示例代码很容易让人理解:

library(limma)
stat <- rnorm(100)
fivenum(stat)
sel <- 1:10
sel2 <- 11:20
stat[sel] <- stat[sel]+1
stat[sel2] <- stat[sel2]-1
stat
plot(stat)

首先呢,上面的代码创造了一个向量,是100个数值,它们最开始时是乱序,但是我们人为的把前面的10个数值加上了1,这样前面的10个数值就会名列前茅。然后呢,我们人为的把第11到20个数值减去1,这样的它们数值会偏小,但是并不会垫底。

下面的barcodeplot可视化你的基因排序 ,很容易看出来:

# One directional
barcodeplot(stat, index = sel)
# Two directional
barcodeplot(stat, index = sel, index2 = sel2)

可以看到,第一个基因集合, 就是前面的10个数值其实是遥遥领先,而第二个基因集合,就是第11到20个数值会比较小,但不会是绝对的垫底。

 

亲爱的读者,发挥你聪明的小脑瓜,思考一下,假如你在前面把 人为的把第11到20个数值减去10,这样的话,第二个基因集合,是不是就可以看到很明显的垫底情况了?

上面的代码大量涉及到R基础知识:

需要把R的知识点路线图搞定,如下:

  • 了解常量和变量概念
  • 加减乘除等运算(计算器)
  • 多种数据类型(数值,字符,逻辑,因子)
  • 多种数据结构(向量,矩阵,数组,数据框,列表)
  • 文件读取和写出
  • 简单统计可视化
  • 无限量函数学习

文末友情推荐

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存