3.1 R语言简介

R是一种统计计算程序。它是一种命令驱动语言，也就是说，你必须在其中键入命令，而不是使用鼠标指向并单击。本指南假设你已经从http://www.r-project.org成功下载并安装了R。不妨在R提示符下输入以下命令行开始使用：

> help.start()

或者如果你在使用Windows操作系统下的R，请按照下拉菜单项Help≻html help获得。从html帮助页面中的Packages≻limma链接将引导你找到limma软件包函数的帮助内容。

在使用任何limma命令之前，你必须在R提示符键入下述命令来加载软件包：

> library(limma)

你可以在R提示符下通过键入?和函数名获得任何载入软件包的帮助，例如：

> ?read.maimages

或等效地

> help("read.maimages")

来获得有关read.maimages函数的详细帮助。单独的函数帮助页面特别重要，列出了函数接受的所有参数以及参数的取值。

理解R语言的一个关键是理解你在R中创建的任何东西都是“对象”。对象可能包括数据集，变量，函数或者任何东西。例如：

> x <- 2

将创建一个变量x，并为它赋值2。在R对话的任何阶段，你都可以键入

> objects()

来获取你创建的所有对象的列表。你可以通过在提示符下键入对象名称来查看任何对象的内容，例如以下命令任一都将打印对象x的内容：

> show(x)
> x

我们希望你无需花费大量时间了解R语言本身就可以使用limma，但在这个方向上的一点认知将是非常有帮助的，特别是当你想要的做一些没有在limma或其他Bioconductor包中明确规定的内容。对于有关R语言的更多详细信息，请参阅在线帮助中提供的R导论。关于使用R进行统计分析的更多背景见[7]。

3.2 limma样本会话

这是一个针对分析可能看起来是什么样子的快速概述。第一个例子假设有4个双色微阵列的重复，其中第二和第四个是染色交换重复实验。我们假设图像已经经过GenePix处理以产生用于每个阵列的.gpr文件和一个目标文件targets.txt，目标文件中包含了.gpr文件的名称列。

> library(limma)
> targets <- readTargets("targets.txt")

设置一个过滤器，使任何-99或更小的标记点都获得零权重。

> f <- function(x) as.numeric(x$Flags > -99)

读取数据。

> RG <- read.maimages(targets, source="genepix", wt.fun=f)

以下命令用于实现一种自适应背景校正。这一步骤是可选的，推荐用于GenePix数据。

> RG <- backgroundCorrect(RG, method="normexp", offset=50)

打印头勒斯标准化：

> MA <- normalizeWithinArrays(RG)

通过拟合每个基因的线性模型估计倍数变化和标准误差。该设计矩阵标明哪些阵列是染色交换重复实验。

> fit <- lmFit(MA, design=c(-1,1,-1,1))

应用经验贝叶斯平滑标准误差。

> fit <- eBayes(fit)

显示前10名基因的统计资料。

> topTable(fit)

第二个例子假设Affymetrix微阵列与野生型（wt）或突变型（mu）RNA杂交。总共应有3个或3个以上的微阵列，以确保一些重复实验。现在目标文件被假设为具有另一Genotype列，指示杂交在每个微阵列上的RNA源。

> library(gcrma)
> library(limma)
> targets <- readTargets("targets.txt")

读取并预处理Affymetrix CEL文件数据。

> ab <- ReadAffy(filenames=targets$FileName)
> eset <- gcrma(ab)

为两个RNA源形成适当的设计矩阵并拟合线性模型。设计矩阵有两列。第一列表示野生型中的对数表达式，第二列表示突变型和野生型样品之间的对数比。有关设计矩阵的更多细节，请参见第9.2节。

> design <- cbind(WT=1, MUvsWT=targets$Genotype=="mu")
> fit <- lmFit(eset, design)
> fit <- eBayes(fit)
> topTable(fit, coef="MUvsWT")

该代码适合线性模型，平滑标准错误并显示突变型与野生型相比的前10个差异基因。

运行eBayes时，选项trend=TRUE和robust=TRUE常常也是有帮助的，可以增加某些类型数据的权重。例如：

> fit <- eBayes(fit, trend=TRUE, robust=TRUE)
> topTable(fit, coef="MUvsWT")

3.3 数据对象

在limma中创建和使用的主要类型的数据对象有6种：

EListRaw：原始表达列表。一种用于存储标准化之前单通道原始强度的类。强度未取对数。此类对象包含一行代表每个探针，一列代表每个微阵列。例如函数read.ilmn()创建了这样的一个这类对象。

EList：表达列表。包含背景校正和标准化的强度对数。通常使用normalizeBetweenArrays()或neqc()从 EListRaw对象创建。

RGList：红-绿色列表。一个用于存储原始双色强度的类，通常由read.maimages()当它们从图像分析输出文件被读入时创建。

MAList：双色强度转换为M值和A值，即转换为在对数尺度上的点内和全点对比度。一般使用MA.RG()或normalizeWithinArrays()从RGList创建。这类对象每个点都包含一行。也可能有多个点，因此每个探针有可能多于一行。

MArrayLM：微阵列线性模型。存储标准化强度或对数比拟合的基因型线性模型的结果。通常由lmFit()创建。这类对象通常是每个唯一的探针包含一行。

TestResults：存储每个探针测试一组对比度为零的结果。通常通过decideTests()创建。这类对象通常每个唯一的探针包含一行。

所有这些对象可以被认为是R中的任何一个列表，例如MA$M提取了M值的矩阵如果MA是一个MAList对象，或者fit$coef提取估计系数如果fit是一个MArrayLM对象。names(MA)显示了包含在该对象的所有内容。对于那些熟悉R中矩阵的人，所有这些对象也被设计为遵循与矩阵的许多类比。在RGList和MAList的情况下，行对应点，列对应阵列。在MarrayLM的情况下，行对应探针，列对应参数或对比度。函数summary，dim，length，ncol，nrow，dimnames，rownames，colnames具有用于这些类的方法。例如：

> dim(RG)
[1] 11088 4

说明RGList对象RG包含数据为11088点和4个阵列。

> colnames(RG)

给出该对象中的文件名或阵列名称，如果fit是一个MArrayLM对象，然后

> colnames(fit)

将给出线性模型拟合中系数的名称。

任何这类对象都可以子集化，所以RG[,j]用于描述阵列j，RG [i,]用于描述由索引i所指示的探针数据 。多个数据对象可以使用cbind，rbind或merge合并。于是

> RG1 <- read.maimages(files[1:2], source="genepix")
> RG2 <- read.maimages(files[3:5], source="genepix")
> RG <- cbind(RG1, RG2)

相当于

> RG <- read.maimages(files[1:5], source="genepix")

另外，如果控制状态已经在MAList对象被设置，那么

> i <- MA$genes$Status=="Gene"
> MA[i,]

可以用于在拟合线性模型之前从数据对象中消除控制点。