8.1 背景介绍

limma软件包使用称为线性模型的方法来分析设计的微阵列实验。这种方法允许分析非常一般的实验，就像分析简单的重复实验一样容易。这种方法在[42, 55]中进行了概述。该方法需要指定一个或两个矩阵。第一个是设计矩阵，表明有每个阵列使用了哪些RNA样本。第二个是对比矩阵，它规定了想要在RNA样品之间进行哪些比较。对于非常简单的实验，可能不需要指定对比矩阵。

方法的理念如下，你必须从将数据拟合到一个线性模型开始，完全建模数据的系统部分。该模型由设计矩阵指定。设计矩阵的每一行对应于实验的阵列，每一列对应于描述实验中RNA源的系数。Affymetrix或单通道数据，或者是有共同参照的双色数据，将需要许多系数，和RNA源的数量一致。直接设计的双色数据，需要的系数比RNA源数目少一个，除非你想要估计每个基因的染料效应，在这种情况下RNA源的数量和系数的数量相同。任何一套独立的系数都可以，只要它们足以描述实验。该步骤的主要目的是估计数据的变异性，因此系统部分需要进行建模，从而可以和随机变异进行区分。

在实践中，从可能想要回答的问题角度来说，要求系数与RNA源数目完全相同太过苛刻。你可能比较感兴趣的是更多或更少的RNA源之间的比较。因此，提供对比步骤，以便可以利用初始系数，并以想要回答任何问题的方式进行比较，不管这些可能有多少。

如果有来自Affymetrix实验的数据，比如来自单通道点微阵列或来自使用共同参照的点微阵列，那么线性建模与普通方差分析或多重回归分析相同，只是需要为每个基因拟合一个模型。具有这种类型的数据，可以像单变量数据普通建模一样创建设计矩阵。如果你使用直接设计的点状微阵列的数据，即没有共同参照的连接设计，那么线性建模方法非常强大，但是设计矩阵的创建可能需要更多的统计知识。

对于统计分析和评估差异表达，limma使用经验贝叶斯方法以减轻估计的对数倍变化的标准误差。这导致更稳定的推导和改进的效力，特别是对于具有少量阵列的实验[42, 26]。对于阵列内部有重复点的阵列，limma使用合并的相关方法充分利用重复点[39]。

8.2单通道设计

Affymetrix数据通常使用affy包进行标准化。我们在这里假设数据可用作名为eset的ExpressionSet对象。这样的对象将包含一个插槽，表示每个阵列上每个基因的对数表达值，可以使用exprs(eset)提取。Affymetrix和其他单通道微阵列数据非常可能用普通线性模型或anova模型分析。与微阵列数据的差异在于它几乎总是需要提取特定的感兴趣的对比，因此在R中为因子提供的标准参数通常不足。

limma中有很多方法可以分析一个复杂实验。一个直接的策略是设置最简单的设计矩阵，然后从拟合中提取感兴趣的对比。

假设有三个RNA源需要进行比较。假设前三个阵列与RNA1杂交，接下来两个与RNA2杂交，其次三个与RNA3杂交。假设RNA源之间的所有两两比较都是值得关注的。我们假设数据已经被标准化并存储在ExpressionSet对象中，例如通过

> data <- ReadAffy()
> eset <- rma(data)

下面的命令可以创建适当的设计矩阵，并使用线性模型进行拟合

> design <- model.matrix(~ 0+factor(c(1,1,1,2,2,3,3,3)))
> colnames(design) <- c("group1", "group2", "group3")
> fit <- lmFit(eset, design)

为了三组之间的两两比较，可以创建合适的对比矩阵

> contrast.matrix <- makeContrasts(group2-group1, group3-group2, group3-group1, levels=design)
> fit2 <- contrasts.fit(fit, contrast.matrix)
> fit2 <- eBayes(fit2)

第2组与第1组中差异表达的基因列表排名可以通过下面的命令获得

> topTable(fit2, coef=1, adjust="BH")

每个假设检验的结果可以使用

> results <- decideTests(fit2)

代表每个比较中显著基因数目的维恩图可以通过下面的命令获得

> vennDiagram(results)

8.3 共同参照设计

现在考虑双色微阵列实验，其中阵列已经使用了共同参照。这样的实验可以非常类似于Affymetrix实验进行分析，只是必须进行染料交换。最简单的方法是使用modelMatrix()函数和目标文件设计矩阵。例如，我们考虑Joëlle Michaud，Catherine Carmichael和Hamish Scott博士在沃尔特伊丽莎豪尔研究所比较转录因子在人类细胞系中的作用实验的一部分。目标文件如下：

> targets <- readTargets("runxtargets.txt")
> targets
      SlideNumber        Cy3        Cy5
1            2144       EGFP       AML1
2            2145       EGFP       AML1
3            2146       AML1       EGFP
4            2147       EGFP  AML1.CBFb
5            2148       EGFP  AML1.CBFb
6            2149  AML1.CBFb       EGFP
7            2158       EGFP       CBFb
8            2159       CBFb       EGFP
9            2160       EGFP  AML1.CBFb
10           2161  AML1.CBFb       EGFP
11           2162       EGFP  AML1.CBFb
12           2163  AML1.CBFb       EGFP
13           2166       EGFP       CBFb
14           2167       CBFb       EGFP

在实验中，绿色荧光蛋白（EGFP）被用作共同参照。利用一个腺病毒系统将各种转录因子转运到HeLa细胞的细胞核中。在这里我们考虑转录因子AML1，CBFbeta或两者都考虑。生成一个简单的设计矩阵和线性模型：

> design <- modelMatrix(targets,ref="EGFP")
> design
  AML1  AML1.CBFb  CBFb
1    1          0     0
2    1          0     0
3   -1          0     0
4    0          1     0
5    0          1     0
6    0         -1     0
7    0          0     1
8    0          0    -1
9    0          1     0
10   0         -1     0
11   0          1     0
12   0         -1     0
13   0          0     1
14   0          0    -1
> fit <- lmFit(MA, design)

将每个转录因子与EGFP进行比较，并且分别比较组合转录因子与AML1和CBFb是有意义的。形成的适当的对比矩阵如下：

> contrast.matrix <- makeContrasts(AML1,CBFb,AML1.CBFb,AML1.CBFb-AML1,AML1.CBFb-CBFb,
+ levels=design)
> contrast.matrix
          AML1  CBFb  AML1.CBFb  AML1.CBFb - AML1  AML1.CBFb - CBFb
AML1         1     0          0                -1                 0
AML1.CBFb    0     0          1                 1                 1
CBFb         0     1          0                 0                -1

现在可以扩展线性模型拟合，并计算经验贝叶斯统计量：

> fit2 <- contrasts.fit(fit, contrasts.matrix)
> fit2 <- eBayes(fit2)

8.4 双色直接设计

没有共同参照的双色设计需要最多的统计知识来选择适当的设计矩阵。直接设计是没有单一RNA源与每个阵列杂交的设计。例如，我们考虑由沃尔特伊丽莎豪尔研究所的Mireille Lahoud博士进行的比较三种不同树突细胞群体（DC）的基因表达实验。

该实验涉及三个染料交换实验中的6个cDNA微阵列，每对比较两种DC类型。设计如上图所示。目标文件如下：

> targets
          SlideNumber     FileName   Cy3  Cy5
ml12med            12 ml12med.spot   CD4  CD8
ml13med            13 ml13med.spot   CD8  CD4
ml14med            14 ml14med.spot    DN  CD8
ml15med            15 ml15med.spot   CD8   DN
ml16med            16 ml16med.spot   CD4   DN
ml17med            17 ml17med.spot    DN  CD4

这种实验的设计矩阵有很多有效的选择，没有一个绝对正确的选择。我们选择的设计矩阵如下：

> design <- modelMatrix(targets, ref="CD4")
Found unique target names:
CD4 CD8 DN
> design
        CD8 DN
ml12med   1  0
ml13med  -1  0
ml14med   1 -1
ml15med  -1  1
ml16med   0  1
ml17med   0 -1

在这个设计矩阵中，CD8和DN群体已经与CD4群体进行了比较。通过线性模型估计的系数将对应于CD8与CD4的对数比（第一列）和DN与CD4的对数比（第二列）。

在表达数据适当标准化后，使用下面命令拟合线性模型

> fit <- lmFit(MA, design)

线性模型现在可以被用来回答任何感兴趣的问题。对于这个实验，目的是在三个DC群体之间进行两两比较。使用对比矩阵完成

> contrast.matrix <- cbind("CD8-CD4"=c(1,0),"DN-CD4"=c(0,1),"CD8-DN"=c(1,-1))
> rownames(contrast.matrix) <- colnames(design)
> contrast.matrix
    CD8-CD4 DN-CD4 CD8-DN
CD8       1      0      1
DN        0      1     -1

对比矩阵可用于扩展线性模型拟合，然后计算经验贝叶斯统计：

> fit2 <- contrasts.fit(fit, contrast.matrix)
> fit2 <- eBayes(fit2)