用limma包的voom方法来做RNA-seq 差异分析

来源：生信菜鸟团评论12,121

大家都知道，这十几年来最流行的差异分析软件就是R的limma包了，但是它以前只支持microarray的表达数据。

考虑到大家都熟悉了它，它又发了一个voom的方法，让它从此支持RNA-seq的count数据啦！大家都知道芯片数据跟RNA-seq数据的本质就是value的分布不一样，所以各种针对RNA-seq的差异分析包也就是提出来一个新的normalization方法而已。而我们limma本身就提出了一个voom的方法来对RNA-seq数据进行normalization

使用这个包也需要三个数据：

  • 表达矩阵
  • 分组矩阵
  • 差异比较矩阵

用法与limma一模一样的，只是多了一个normalization而已。

读取自己的数据

一般我们会自己读取表达矩阵和分组信息，下面是一个例子：

  1. options(warn=-1)
  2. suppressMessages(library(DESeq2))
  3. suppressMessages(library(limma))
  4. suppressMessages(library(pasilla))
  5. data(pasillaGenes)
  6. exprSet=counts(pasillaGenes)
  7. head(exprSet)  ##表达矩阵如下
  1. ##             treated1fb treated2fb treated3fb untreated1fb untreated2fb
  2. ## FBgn0000003          0          0          1            0            0
  3. ## FBgn0000008         78         46         43           47           89
  4. ## FBgn0000014          2          0          0            0            0
  5. ## FBgn0000015          1          0          1            0            1
  6. ## FBgn0000017       3187       1672       1859         2445         4615
  7. ## FBgn0000018        369        150        176          288          383
  8. ##             untreated3fb untreated4fb
  9. ## FBgn0000003            0            0
  10. ## FBgn0000008           53           27
  11. ## FBgn0000014            1            0
  12. ## FBgn0000015            1            2
  13. ## FBgn0000017         2063         1711
  14. ## FBgn0000018          135          174
  1. group_list=pasillaGenes$condition
  1. ## [1] treated   treated   treated   untreated untreated untreated untreated
  2. ## Levels: treated untreated
  3. ##这是分组信息，7个样本，3个处理的，4个未处理的对照！

另外一个例子是从airway里面得到表达矩阵和分组信息

  1. library(airway)
  2. data(airway)
  3. exprSet=assays(airway)$counts  ## 表达矩阵
  4. group_list=colData(airway)$dex ## 分组信息

第一步：构建分组矩阵

  1. suppressMessages(library(limma))
  2. design <- model.matrix(~factor(group_list))
  3. colnames(design)=levels(factor(group_list))
  4. rownames(design)=colnames(exprSet)
  5. design
  1. ##              treated untreated
  2. ## treated1fb         1         0
  3. ## treated2fb         1         0
  4. ## treated3fb         1         0
  5. ## untreated1fb       1         1
  6. ## untreated2fb       1         1
  7. ## untreated3fb       1         1
  8. ## untreated4fb       1         1
  9. ## attr(,"assign")
  10. ## [1] 0 1
  11. ## attr(,"contrasts")
  12. ## attr(,"contrasts")$`factor(group_list)`
  13. ## [1] "contr.treatment"

第二步：根据分组信息和表达矩阵进行normalization

这里构建了一个对象 An object of class “EList”

  1. v <- voom(exprSet,design,normalize="quantile")
  2. ## 下面的代码如果你不感兴趣不需要运行，免得误导你
  3. ## 就是看看normalization前面的数据分布差异
  4. #png("RAWvsNORM.png")
  5. exprSet_new=v$E
  6. par(cex = 0.7)
  7. n.sample=ncol(exprSet)
  8. if(n.sample>40) par(cex = 0.5)
  9. cols <- rainbow(n.sample*1.2)
  10. par(mfrow=c(2,2))
  11. boxplot(exprSet, col = cols,main="expression value",las=2)
  12. boxplot(exprSet_new, col = cols,main="expression value",las=2)
  13. hist(exprSet)
  14. hist(exprSet_new)
  15. #dev.off()

用limma包的voom方法来做RNA-seq 差异分析

可以很明显看到normalization前后数据分布差异非常大！

第三步：做差异分析，提取差异分析结果

这里不需要差异比较矩阵了，因为我的分组矩阵表明样本就分成两组，直接提取coef=2的数据即可

  1. fit <- lmFit(v,design)
  2. fit2 <- eBayes(fit)
  3. tempOutput = topTable(fit2, coef=2, n=Inf)
  4. DEG_voom = na.omit(tempOutput)
  5. head(DEG_voom)
  1. ##                  logFC  AveExpr         t      P.Value    adj.P.Val
  2. ## FBgn0029167  2.1608689 7.880589  41.19142 5.195806e-08 0.0007518331
  3. ## FBgn0003137 -0.9560776 8.421903 -29.97211 2.920473e-07 0.0021129620
  4. ## FBgn0003748 -1.0385933 6.395990 -23.78787 1.020682e-06 0.0049230892
  5. ## FBgn0035085  2.5190255 5.221183  21.01339 1.993995e-06 0.0058809216
  6. ## FBgn0050185 -0.4886279 5.487547 -19.95422 2.635106e-06 0.0058809216
  7. ## FBgn0015568 -1.1040979 3.922438 -19.96954 2.624231e-06 0.0058809216
  8. ##                    B
  9. ## FBgn0029167 9.065020
  10. ## FBgn0003137 7.800063
  11. ## FBgn0003748 6.652695
  12. ## FBgn0035085 5.870585
  13. ## FBgn0050185 5.734162
  14. ## FBgn0015568 5.55756

差异分析结果就不解释啦！

原文来自：www.bio-info-trainee.com

发表评论

匿名网友

拖动滑块以完成验证
加载失败