第4章 读取微阵列数据

4.1 本章范围

本章涵盖除Affymetrix以外的大多数微阵列类型。从Affymetrix GeneChips读取数据，请使用affy，gcrma或aroma.affymetrix软件包来读取和标准化数据。

4.2 推荐文件

我们假设已经使用一个或多个微阵列进行了实验全部利用相同的探针库打印。扫描每个阵列以产生TIFF图像，然后使用诸如ArrayVision，ImaGene，GenePix，QuantArray或SPOT的图像分析程序处理TIFF图像，获得每个点的红色和绿色的前景和背景强度。然后将点强度从图像分析程序导出到一系列文本文件中。每个阵列对应一个文件，如果使用Imagene软件，每个阵列对应两个文件。

你需要有图像分析的输出文件。在大多数情况下，这些文件将包含ID以及探针的名称以及可能的其他注释信息。有些图像分析程序，例如SPOT，不会将探针ID写入输出文件中。在这种情况下，你还需要一个描述探针的基因列表文件。大多数情况下，还需要一个目标文件（targets file）描述与每个阵列的每个通道杂交的是哪个RNA样品。另一个可选文件是识别特殊探针（如控制点）的点类型文件（spot types file）。

4.3 目标框架

准备limma输入数据的第一步通常是创建一个目标文件，其中列出了与每个阵列的每个通道杂交的RNA靶标。它通常是以制表符分隔的文本格式，并且应该包含实验中的每个微阵列行。该文件可以有任何名称，但是默认为Targets.txt。如果它具有默认名称，使用下列命令可以将它读入R会话

> targets <- readTargets()

一旦读入R，它就成为了目标框架（targets frame）。

目标框架通常包含一个FileName栏，给出图像分析输出文件的名称，Cy3列给出片子中用Cy3染料标记的RNA类型，Cy5列给出片子中用Cy5染料标记的RNA型。其他列是可选的。该目标文件可以使用任何文本编辑器，但电子表格程序如Microsoft Excel，比较方便。对于Swirl案例研究的目标文件包括可选的SlideNumber和Date列：

给关联的微阵列创建短的可读标签通常是很方便的，特别是在文件名过长或非直观的情况下用于后续的图文输出。包含这些标签的列可以包含在目标文件中，例如Apoa1案例研究中使用的Name列：

此列可用于为目标框架创建行名称

> targets <- readTargets("targets.txt", row.names="Name")

当读取这些数据对象时，行名称可以转化成数组名称。

对于ImaGene文件，文件名列拆分成一个FileNameCy3列和一个FileNameCy5列，因为ImaGene软件将红色和绿色的强度存储在单独的文件中。这里是一个简短的例子：

4.4 读取双色强度数据

假设files是一个包含图像分析软件输出文件名称的字符向量。前景和背景强度可以使用下面的命令读入一个RGList对象

> RG <- read.maimages(files, source="<imageanalysisprogram>", path="<directory>")

其中<imageanalysisprogram>是图像分析程序的名称，<directory>是包含文件目录的完整路径。如果文件在当前的R工作目录中，参数path可省略；参考setwd帮助条目了解如何设置当前工作目录。文件名通常是从目标文件中读取的。例如，目标文件Targets.txt连同SPOT输出文件都存在于当前工作目录中，那么可能使用下面的命令行

> targets <- readTargets()
> RG <- read.maimages(targets$FileName, source="spot")

或者更简单地说，可能将目标框架本身作为files的一个参数

> RG <- read.maimages(targets, source="spot")

在这种情况下，软件将自动在目标框架中查找FileName列。

如果是GenePix输出文件，那么可能会使用下面的命令读取

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

来给出适当的目标文件。查看read.maimages帮助条目获得支持的其它图像分析程序。默认情况下，文件默认为制表符分隔，尽管其他分隔符可以用在sep=参数中来指定。

从ImaGene软件读取数据与其他图像分析程序有所不同，因为红色和绿色强度分别存储在单独的文件中。这意味着目标框架应该包括两个文件名栏，比如说，FileNameCy3和FileNameCy5，分别表示包含绿色和红色强度的文件名称。第4.3节给出了一个例子。ImaGene数据的典型代码可能是

> targets <- readTargets()
> files <- targets[,c("FileNameCy3","FileNameCy5")]
> RG <- read.maimages(files, source="imagene")

对于ImaGene数据，files参数read.maimages()是文件名的2列矩阵，而不是一个向量。

下表给出了前景和背景强度的默认估计值：

来源	前景	背景
agilent	Median Signal	Median Signal
agilent.mean	Mean Signal	Median Signal
agilent.median	Median Signal	Median Signal
bluefuse	AMPCH	None
genepix	F Mean	B Median
genepix.median	F Median	B Median
genepix.custom	Mean	B
imagene	Signal Mean	Signal Median, or Signal Mean if auto segmentation has been used
quantarray	Intensity	Background
scanarrayexpress	Mean	Median
smd.old	I_MEAN	B_MEDIAN
smd	Intensity (Mean)	Background (Median)
spot	mean	morph
spot.close.open	mean	morph.close.open

默认估计值可以通过指定read.maimages()中的columns参数重写。例如，假设GenePix已经使用了自定义的背景方法，然后你希望使用中位数进行前景估计。虽然limma软件包中的预设选择没有提供这种前景和背景的组合，但你可以指定它

> RG <- read.maimages(files,source="genepix",
+ columns=list(R="F635 Median",G="F532 Median",Rb="B635",Gb="B532"))

如果你的图像分析程序不在上述列表中，该怎么办？如果图像输出文件是标准格式，那么你可以自己提供注释和强度列名。
例如，

> RG <- read.maimages(files,
+ columns=list(R="F635 Mean",G="F532 Mean",Rb="B635 Median",Gb="B532 Median"),
+ annotation=c("Block","Row","Column","ID","Name"))

完全等同于source="genepix"。“标准格式”是指每个感兴趣的列都有一个独特的列名称标识，并且该文件中最后一行数据后没有其它行。文件开始时的头文件信息是可以接受的，但是文件末尾额外的行将导致读取失败。
查看你的数据来确定它是否正确读取是一个好主意。键入

> show(RG)

查看前几行数据的打印。同时尝试

> summary(RG$R)

查看每个阵列红色强度的5个数字摘要，等等。

可以通过几个步骤读取数据。如果RG1和RG2是对应于不同阵列的两个数据集，那么

> RG <- cbind(RG1, RG2)

可以将它们合并成一个大的数据集。数据集也可以子集化。例如RG[,1]是第一阵列中的数据，而RG[1:100,]是其中前100个基因的数据。

4.5 读取单通道Agilent强度数据

读取单通道数据类似于双色数据，除了应该添加参数green.only=TRUE到read.maimages()来说明没有红色通道。由Agilent特征提取软件生成的单通道Agilent强度，可以通过下面的命令读取

> x <- read.maimages(files, source="agilent", green.only=TRUE)

或者

> x <- read.maimages(targets, source="agilent", green.only=TRUE)

如果是双色数据，将会使用path参数：

> x <- read.maimages(files, source="agilent", path="<directory>", green.only=TRUE)

如果数据文件不在当前工作目录。green.only参数告诉read.maimages()函数输出一个EList对象而非RGList对象。原始强度将被存储在数据对象的E组分中，并且可以通过下面的命令检查

> summary(x$E)

Agilent特征提取软件拥有估计每个点前景和背景信号的能力，使用前景像素和背景像素的平均值或中值。默认read.maimages要读取前景和背景的中值信号。
或者

x <- read.maimages(targets, source="agilent.mean", green.only=TRUE)

将会读取前景信号平均值同时采用背景信号中位数。source的可能值是：

来源	前景	背景
agilent	Median Signal	Median Signal
agilent.mean	Mean Signal	Median Signal
agilent.median	Median Signal	Median Signal

至于双色数据，前景和背景估计的默认选项可以通过在read.maimages()函数指定columns参数来重写。

Agilent特征提取输出文件包含探针注释列以及强度列。默认情况下，如果存在，read.maimages()函数将读取下面的注释列：Row，Col，Start，Sequence，SwissProt，GenBank，Primate，GenPept，ProbeUID，ControlType，ProbeName，GeneName，SystematicName，Description。

完整的单通道Agilent数据工作案例参见第17.4节。

4.6 读取Illumina BeadChip数据

Illumina全基因组BeadChip需要特殊处理。Illumina图像由Bead-Scan软件扫描，Illumina的BeadStudio或GenomeStudio软件可用于导出探针摘要文件。探针摘要文件是包含强度数据的制表符分隔文件。通常，同时处理的所有阵列都将写入一个文件，不同列对应相应数组。我们建议从GenomeStudio导出没有背景校正或标准化的强度，因为这些预处理步骤可以通过limma函数更好地完成。还可以要求GenomeStudio导出控制探针的配置文件，我们也建议这样做。

Illumina文件与其他平台不同，每个图像输出文件包含来自多个阵列的数据，并且控制探针的强度被写入到和常规探针分开的不同文件中。这些文件的其他功能可以可选地用于预处理和过滤。Illumina探针摘要文件可以通过read.ilmn函数读取。典型的用法是

> x <- read.ilmn("probe profile.txt", ctrlfiles="control probe profile.txt")

probe profile.txt是主要探针摘要配置文件，control probe profile.txt是包含控件探针的配置文件名称。

如果要读取多个探针摘要文件，并且样本被汇总到目标框架中，那么可以使用read.ilmn.targets函数。

读取控制探针配置文件是可选的，但是推荐这样做。如果控制探针文件可用，那么Illumina数据可以利用neqc或nec函数进行背景校正和标准化。否则Illumina数据将像其他单通道平台数据一样被背景校正和标准化。

有关Illumina微阵列数据的完整工作案例研究，请参见第17.3节。

4.7 图像点质量权重

图像分析程序通常输出大量的信息，除了前景和背景强度，还提供了每个点的质量信息。有时候我们希望使用这些信息为每个点生成质量指标可以用于随后的分析步骤。一种方法是删除所有不符合特定质量标准的点。更复杂的方法是产生定量的质量指标，然后依赖其感知可靠性来逐级升高或降低每个点的权重。limma提供了一种同时支持这两种计算点重量方法的计算方式。

limma方法是计算每个点的定量质量权重。重量在limma中的处理方式类似在R中的大多数回归函数，例如lm()。0重量表示在所有分析中该点因为不可靠应该被忽略。重量1表示正常质量。大于或小于1的点质量将导致该点在随后的分析中被给予相对或多或少的权重。点重量小于0没有意义。

当强度从图像分析输出文件读取的同时，可以读取质量信息并计算点的质量权重。质量权重的计算由read.maimages()函数的wt.fun参数定义。这个参数是一个函数，其定义了如何使用从图像分析文件中发现的信息计算权重。获得好的点质量权重并不是直截了当的，非常依赖于图像分析软件的使用。limma提供了一些被研究人员发现有用的例子。

一些图像分析程序产生质量指标作为输出的一部分。例如，GenePix生成一个名为Flags的列，对于“正常”点为0，不同类别的问题点的负值越来越大。如果你正在读取GenePix图像分析文件，命令行

> RG <- read.maimages(files,source="genepix",wt.fun=wtflags(weight=0,cutoff=-50))

将读取强度数据，并计算点权重矩阵，给出零重量和Flags值小于-50的任何点。权重存储在RGList数据对象的weight组件中。权重可以通过在RG列表上操作的诸如normalizeWithinArrays之类的函数自动使用。

从图像像素方面来说，理想尺寸是一个完美的圆形点。在这个情况下，对于比这个理想尺寸要大得多或小得多的点，降低它们的权重可能有用。如果SPOT图像分析输出正在读取，以下调用

> RG <- read.maimages(files,source="spot",wt.fun=wtarea(100))

将给予正好具有100个像素的区域完全的权重，降低更小和更大的点的权重，具有0面积或超过理想尺寸两倍的点给予0重量。

计算点质量权重的适当方法取决于使用的图像分析程序。请查阅QualityWeights帮助条目以了解可用的质量权重函数。wt.fun参数非常灵活，允许你构建自己的权重。wt.fun参数可以是将数据集作为参数并计算所需权重的任何函数。例如，如果你希望给予所有不到-50的GenePix标志零权重，你可以使用

> myfun <- function(x) as.numeric(x$Flags > -50.5)
> RG <- read.maimages(files, source="genepix", wt.fun=myfun)

wt.fun可用于基于图像分析文件中任意数量的列计算权重。例如，一些研究人员喜欢过滤掉某些点，如果给定点的GenePix前景平均值和中位数差异超过某个阈值，比如说50，可以通过下面的命令实现

> myfun <- function(x, threshold=50) {
+ okred <- abs(x[,"F635 Median"]-x[,"F635 Mean"]) < threshold
+ okgreen <- abs(x[,"F532 Median"]-x[,"F532 Mean"]) < threshold
+ as.numeric(okgreen & okred)
+}
> RG <- read.maimages(files, source="genepix", wt.fun=myfun)

那么所有的“坏”点都会得到零权重，这在limma中相当于将它们标记出来。这里的myfun的定义可以用任何其他代码来替代，利用GenePix输出文件中的列来计算权重。

4.8 读取探针注释

由read.maimages()读取的RGList几乎总是有一个称为genes的组件，包含探针相关联的ID和其他注释信息。唯一的例外是SPOT数据，source=“spot”，或者当读取通用数据，source=“generic”，而不设置注释参数时，annotation=NULL。尝试

> names(RG$genes)

看是否已经设定了genes组分。

如果未设置genes组分，则需要从单独的文件中读取探针ID。如果已经使用Axon扫描仪扫描了阵列，探针ID将存储在一个制表符分隔的GenePix数组列表（GAL）文件中。如果GAL文件的扩展名为“gal”，并且在当前的工作目录中，那么它可能通过下面的命令被读入一个数据框中

> RG$genes <- readGAL()

非GenePix基因列表可以使用R基础包中的read.delim函数读入R中。

4.9 打印布局

打印布局是阵列上的点和区块的排列。了解打印布局和用多打印头的机械机器人印刷的旧式学术点阵列密切相关。区块有时称为打印头组或针组或多行列。每个区块对应于打印头上的打印针，用于打印阵列；阵列上区块的布局对应于打印头上打印针的布局。每个区块中的点数是打印头下降到阵列上的次数。在可能的情况下，例如Agilent，GenePix或ImaGene数据，read.maimages函数将打印布局信息设置在组件printer中。尝试

> names(RG$printer)

查看打印布局信息是否已设置。

如果你已经使用readGAL来设置了genes组件，那么你也可以使用getLayout来设置printer信息

> RG$printer <- getLayout(RG$genes)

请注意，这仅适用于GenePix GAL文件，而不适用于一般基因列表。

4.10 点类型文件

点类型文件（STF）是另一个可选的制表符分隔文本文件，允许识别来自基因列表条目中的不同类型的探针。它在识别不同类型的控制探针方面特别有用。STF用于设置阵列上每个探针的控制状态，使得图像可以以适当的方式突出不同类型的点。它通常是用于区分控制探针与常规探针对应的基因，并区分阴性对照和阳性对照，校正对照和比例等。STF应该有一个SpotType列给出不同点类型的名称。一个或多个其他列应该与基因列表中的列具有相同的名称，并应包含模式或正则表达式足以识别点的类型。任何其他列应该包含绘图属性，如颜色或符号，与点的类型相关联。每个点类型有一行用于区分。

STF使用简化的正则表达式来匹配模式。例如，AA*表示任何以AA开头的字符串，*AA表示以AA结尾的任何代码，AA表示这两个字母，*AA*表示任何包含AA的字符串，AA.意思是AA后紧跟一个其他字符，AA\.表示AA后跟一个句点，没有其他字符。对于那些熟悉正则表达式的人，任何其他正则表达式都是允许的，但^用于字符串开始和$用于结尾的字符串代码应该排除。注意，模式从第一个顺序匹配到最后，应该首先包括更一般的模式。第一行应该指定默认的spot-type并且应该具有所有模式匹配列的pattern*。

以下是适用于ApoAI数据的简短STF：

在此示例中，列ID和Name位于基因列表中，并包含匹配的模式。星号是可以代表任何东西的通配符。小心使用适当的大写或小写，不要插入任何额外的空格。剩下的列给出与不同类型的点关联的颜色。该代码假定探针注释数据框包括列ID和Name。如果GenePix被用于图像分析，情况通常是这样，但其他图像分析软件可能使用其他列名。

下面是一个适当的用于Lucidea通用记分卡控制点阵列的STF。

如果STF具有默认名称SpotTypes.txt，那么可以通过下面的命令读取

> spottypes <- readSpotTypes()

它通常被用作参数传递给controlStatus()的函数，来设置阵列每个点的状态，例如

> RG$genes$Status <- controlStatus(spottypes, RG)