使用camera进行基因集分析

2021-05-27 15:51:18 浏览数 (3)

发表这个算法的文章是2012年:《Camera: a competitive gene set test accounting for inter-gene correlation》,链接是:https://academic.oup.com/nar/article/40/17/e133/2411151 如果是第一次接触生物信息学软件的小伙伴可能是会觉得有点老。

而且呢,里面的公式一大把,看起来会比较头疼,其实可以把它理解为 gsea的类似方法即可,gsea分析这方面教程我在《生信技能树》公众号写了不少了,不管是芯片还是测序的表达矩阵,都是一样的,把全部基因排序即可:

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

目前这个 camera 功能实现起来非常的方便,因为它仅仅是limma包的一个函数:camera: Competitive Gene Set Test Accounting for Inter-gene Correlation

首先呢,我们认为的创造一个表达量矩阵,1000行代表1000个基因,然后6个样品区分成为2个组,这样的话可以很方便做差异分析:

代码语言:javascript复制
library(limma)
y <- matrix(rnorm(1000*6),1000,6)
design <- cbind(Intercept=1,Group=c(0,0,0,1,1,1))
design
# First set of 20 genes are genuinely differentially expressed
index1 <- 1:20
y[index1,4:6] <- y[index1,4:6] 1
# Second set of 20 genes are not DE
index2 <- 21:40

可以看到,我们把前面的20个基因整体呢在右边的3个样品表达量增加了1,所以前面的20个基因组成的集合理论上应该是功能富集啦。但是index2 <- 21:40代表的基因集我们并没有处理它,它就是一个随机数,所以理论上不应该是被富集。

接下来就使用limma包的一个函数:camera: Competitive Gene Set Test Accounting for Inter-gene Correlation 对这两个基因集进行统计学检验吧:

代码语言:javascript复制
> camera(y, index1, design)
     NGenes Direction       PValue
set1     20        Up 3.110234e-10
> camera(y, index2, design)
     NGenes Direction    PValue
set1     20        Up 0.9888623
> 
> 
> camera(y, list(set1=index1,set2=index2), 
         design, inter.gene.cor=NA)
     NGenes Correlation Direction      PValue         FDR
set1     20 -0.01269341        Up 0.002283103 0.004566206
set2     20 -0.02253970        Up 0.988594938 0.988594938
> camera(y, list(set1=index1,set2=index2), 
         design, inter.gene.cor=0.01)
     NGenes Direction       PValue          FDR
set1     20        Up 3.110234e-10 6.220468e-10
set2     20        Up 9.888623e-01 9.888623e-01

可以看到,只有第一个被我们认为干扰的基因集,被认为是统计学显著的上升。

同样的,我们也可以很方便的可视化一下:

代码语言:javascript复制
 stat=rowMeans(y[,4:6]) - rowMeans(y[,1:3])
# One directional
barcodeplot(stat, index = index1)
# Two directional
barcodeplot(stat, index = index1,
            index2 = index2)

如下所示:

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

  • 生信分析人员如何系统入门R(2019更新版)

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

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

0 人点赞