搜索
查看: 5059|回复: 0

[software] limma常用的预处理差异表达一体包

[复制链接]

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-7-10 09:31:57 | 显示全部楼层 |阅读模式
## Dependencies

This document has the following dependencies:

```{r dependencies, warning=FALSE, message=FALSE}
library(limma)
library(leukemiasEset)
```

Use the following commands to install these packages in R.

```{r biocLite, eval=FALSE}
source("http://www.bioconductor.org/biocLite.R")
biocLite(c("limma", "leukemiasEset"))
```

## Overview

`r Biocpkg(“limma”)`是用于分析微阵列和RNA-seq数据的非常受欢迎的软件包。

LIMMA代表“微阵列数据的线性模型”。 也许不出意料的是,limma包含适合广泛类别的统计模型(称为“线性模型”)的功能。 这些模型的例子包括线性回归和方差分析。 虽然limma的大多数功能已经被开发用于微阵列数据,但是limma的模型拟合程序对于许多类型的数据是有用的,并且不限于微阵列。 例如,我在我的DNA甲基化分析研究中使用了limma。

### Obtaining data
我们使用`bioCxPlus(“leukemiasEset”)数据集的`leukemiasEset`数据集 这是一个“ExpressionSet”。
```{r load}
library(leukemiasEset)
data(leukemiasEset)
leukemiasEset
table(leukemiasEset$LeukemiaType)#计数
```

这是不同类型白血病的数据。 代码“NoL”不是白血病,也就是正常对照.

让我们询问哪些基因在“ALL”类型和正常对照之间有差异表达。 首先我们对数据进行子集并进行清理
```{r subset}
ourData <- leukemiasEset[, leukemiasEset$LeukemiaType %in% c("ALL", "NoL")]
ourData$LeukemiaType <- factor(ourData$LeukemiaType)
```

### A linear model
现在我们做一个标准的limma模型拟合
```{r limma}
design <- model.matrix(~ ourData$LeukemiaType)
fit <- lmFit(ourData, design)
fit <- eBayes(fit)
topTable(fit)
```

这里发生的是一个常见的limma(和朋友)工作流程。首先,通过所谓的“设计矩阵”来定义兴趣的比较(和实验的设计)。该矩阵基本上涵盖了我们对设计知识的一切;在这种情况下,有两个组(我们在下面的设计中有更多的说明)。接下来,该模型是适合的。其次是使用所谓的经验贝叶斯程序(这是真正有益于奇迹的limma中的一步),在基因之间借用力量。因为这个设计只有两个组,所以只有一个可能的比较:两组之间哪些基因不同。这个问题由'topTable()'函数检查,它列出了最差异表达基因。在一个比较复杂的设计中,`topTable()'函数需要被告知哪些比较的兴趣来进行总结。


输出的一个重要部分是“logFC”,它是foldchange。取log   要解释此数量的符号,您需要知道这是否为“ALL-NOL”(在这种情况下,正值在ALL中上调)或相反。在这种情况下,这是由参考水平确定的,该参考水平是因子的第一级。
```{r level}
ourData$LeukemiaType
```

正的数值意味着它在癌症中被下调。 您可以使用`relevel()`函数来更改因子的参考级别。 您也可以通过手动计算“logFC”来确认这一点,这对于知道是有用的。 我们来计算前面的差异表达基因的倍数变化:
```{r FCbyHand}
topTable(fit, n = 1)
genename <- rownames(topTable(fit, n=1))
typeMean <- tapply(exprs(ourData)[genename,], ourData$LeukemiaType, mean)
typeMean
typeMean["NoL"] - typeMean["ALL"]
```

人工检查有时是有用的,以确保您有正确的解释。最后,请注意,在计算“logFC”时,limma不会与方法的差异做任何不同的事情;所有的统计改进都集中于计算更好的t统计和p值。

具有统计经验的读者将会注意到,我们正在做的是比较两组;这与经典的t统计量是一样的。我们在这里计算的确是一个t统计量,但是方差估计(t统计的分母)通过借助基因的强度来调节*(eBayes())是这样计算的;这被称为t统计。

`topTable()`的输出包括

- “logFC”:案例和控件之间的日志折叠更改。
- t:用于评估差异表达的t统计量。
- “P.Value”:差分表达式的p值;此值不会针对多次测试进行调整。
- `adj.P.Val`:p值调整为多次测试。不同的调整方法是可用的,默认是Benjamini-Horchberg。

## More on the design

在前面的分析中,我们设置了以下的模型
```{r design2}
design <- model.matrix(~ ourData$LeukemiaType)
```

我们来看一下
```{r headDesign}
head(design)
```
这个矩阵有两列,因为这个概念设计中有两个参数:两个组中的每个表达式。 在这个**参数化中**列1表示“ALL”组的表达式,第2列表示从“NoL”组到“ALL”组的表达式级别的差异。 通过测试第二个参数(等于两组之间的表达差异)是否等于零,来测试两个组具有相同的表达水平。

不同的参数是
```{r design3}
design2 <- model.matrix(~ ourData$LeukemiaType - 1)
head(design2)
colnames(design2) <- c("ALL", "NoL")
```

在本设计中,设计矩阵的两列对应的两个参数代表两组中的表达水平。 科学问题被翻译成询问这两个参数是否相同。 让我们看看我们如何形成对比矩阵
```{r design4}
fit2 <- lmFit(ourData, design2)
contrast.matrix <- makeContrasts("ALL-NoL", levels = design2)
contrast.matrix
```

这里我们说我们对'ALL-NOL'感兴趣,这与上面我们正在做的事情相反(那里是“NoL-ALL”;因为“NoL”是自然的参考组,这更有意义了)
```{r cont.fit}
fit2C <- contrasts.fit(fit2, contrast.matrix)
fit2C <- eBayes(fit2C)
topTable(fit2C)
```

请注意,除了“logFC”列之外,这与上面的“topTable()”完全相同。
## Background: Data representation in limma
  It has a slot called `genes` which is basically equivalent to `featureData` for `ExpressionSet`s (ie. information about which genes are measured on the microarray) as well as a `targets` slot which is basically the `pData` information from `ExpressionSet`.
如上所述,`r Biocpkg(“limma”)`直接在`ExpressionSet`上工作。 它也可以直接用于矩阵。 但是,limma还有一个“RGList”类,它代表一个双色微阵列。 存储在这个类中的基本数据是“ExpressionSet”的,但它至少有两个表达式测量矩阵R(红色)和“G”(绿色),另外还有两个额外的背景估计矩阵 和“Gb”)。 它有一个称为“基因”的slot,其基本上等同于“ExpressionSet”的“featureData”(即关于在微阵列上测量哪些基因的信息)以及基本上是来自“ExpressionSet”中“pData”信息的targets”slot。

## Background: The targets file

limma引入了所谓的“目标”文件的概念。 这只是一个简单的文本文件(通常是TAB或逗号分隔),它保存表型数据。 这个想法是,许多用户在电子表格程序中创建此文本文件更容易,然后将其读入R并将信息存储在数据对象中。


回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-11-21 10:04 , Processed in 0.027175 second(s), 23 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.