基因芯片技术的特点使用寡聚核苷酸探针检测基因。前一节使用ReadAffy函数读取CEL文件获得的数据是探针水平的（probe level），即杂交信号，而芯片数据预处理的目的是将杂交信号转成表达数据（即表达水平数据，expression level data）。存储探针水平数据的是AffyBatch类对象，而表达水平数据为ExpressionSet类对象。基因芯片探针水平数据处理的R软件包有affy, affyPLM, affycomp, gcrma等，这些软件包都很有用。如果没有安装可以通过运行下面R语句安装：

source("http://bioconductor.org/biocLite.R")

biocLite(c("affy","gcrma","affyPLM","affycomp"))

Affy芯片数据的预处理一般有三个步骤：背景处理（background adjustment），归一化处理（normalization，或称为“标准化处理”），汇总（summarization）。最后一步获取表达水平数据。需要说明的是，每个步骤都有很多不同的处理方法（算法），选择不同的处理方法对最终结果有非常大的影响。选择哪种方法是仁者见仁智者见智，不同档次的杂志或编辑可能有不同的偏好。

0、需要了解的一点Affy芯片基础知识

Affy基因芯片的探针长度为25个碱基，每个mRNA用11~20个探针去检测，检测同一个mRNA的一组探针称为probe sets。由于探针长度较短，为保证杂交的特异性，affy公司为每个基因设计了两类探针，一类探针的序列与基因完全匹配，称为perfect match（PM）probes，另一类为不匹配的探针，称为mismatch （MM）probes。PM和MM探针序列除第13个碱基外完全一样，在MM中把PM的第13个碱基换成了互补碱基。PM和MM探针成对出现。我们先使用前一节的方法载入数据并修改芯片名称：

library(affy)

the.filter<-matrix(c("CEL file (*.cel)","*.cel","All (*.*)","*.*"),ncol=2,byrow=T)

cel.files<-choose.files(caption="Select CEL files",multi=TRUE,filters=the.filter,index=1)

data.raw<-ReadAffy(filenames=cel.files)

sampleNames(data.raw)<-paste("CHIP",1:length(cel.files),sep="")

用pm和mm函数可查看每个探针的检测情况：

>pm.data.raw<-pm(data.raw)

>head(pm.data.raw,2)

CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8

501131127.0166.3112139.8111.385.5126.3102.8

251604118.5105.082101.594.081.3103.8103.0

>mm.data.raw<-mm(data.raw)

>head(mm.data.raw,2)

CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8

50184389.088.080.591.077.07579.072.0

252316134.377.377.0107.898.57599.571.3

上面显示的列名称就是探针的名称。而基因名称用probeset名称表示：

>head(geneNames(data.raw))

[1]"244901_at""244902_at""244903_at""244904_at""244905_at""244906_at"

probeset不一定和实际基因一一对应，这些后面对探针进行基因名称映射时会看到。

一、背景处理（background adjustment）

虽然说是背景处理，但是这一步既处理背景值，又处理噪声信号。注意背景和噪声是两个概念，比如说乡间夜晚的蛙叫声虽嘈杂但很稳定，算是背景，如果突然来一声狗叫，那就是噪声。这两者在统计上可以区分。

芯片的背景处理理论上很简单，因为Affy公司设计MM的目的就是检测非特异杂交信号，PM -MM 不就结了？但是研究发现居然有多达30%的MM探针获得的信号强度比相应PM探针的还强（嘿嘿），所以啊，研究者的饭碗就出来了，整些看起来还合理的方法吧。

R软件包affy用于芯片背景噪声消减的函数是bg.correct()，而MAS和RMA方法是最常用的两种方法。

MAS方法将芯片分为k（默认值为16）个网格区域，用每个区域使用信号强度最低的2%探针去计算背景值和噪声。RMA方法的原理比较复杂，可以参看文献：R. A. Irizarry, B. Hobbs, F. Collin, et al. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics, 4:249–64, 2003b. 11, 18, 27, 232, 241, 432, 443。

>data.rma<-bg.correct(data.raw,method="rma")

Loadingrequiredpackage:AnnotationDbi

>data.mas<-bg.correct(data.raw,method="mas")

>class(data.rma)

[1]"AffyBatch"

attr(,"package")

[1]"affy"

>class(data.mas)

[1]"AffyBatch"

attr(,"package")

[1]"affy"

>class(the.data)

[1]"AffyBatch"

attr(,"package")

[1]"affy"

可以看到：ReadAffy()读入的CEL芯片数据以AffyBatch类数据形式存储，而背景消减后得到的依然是AffyBatch类数据。

MAS方法应用后PM和MM的信号强度都被重新计算。RMA方法仅使用PM探针数据，背景调整后MM的信号值不变。

>head(pm(data.raw)-pm(data.mas),2)

CHIP1    CHIP2    CHIP3    CHIP4    CHIP5    CHIP6    CHIP7    CHIP8

50113179.3434462.9287863.2347172.8700362.4808264.4335962.9744760.85734

25160477.5727463.0641563.7300180.6878463.0687366.5250963.3347760.33844

>head(pm(data.raw)-pm(data.rma),2)

CHIP1     CHIP2    CHIP3     CHIP4    CHIP5    CHIP6     CHIP7    CHIP8

501131111.2499102.2057693.22705116.3628293.7457176.14978102.8209485.21383

251604103.879588.6881172.5728690.7481582.2269372.6360087.7493285.31779

>head(mm(data.raw)-mm(data.mas),2)

CHIP1    CHIP2    CHIP3    CHIP4    CHIP5    CHIP6    CHIP7    CHIP8

50184379.3370862.9294263.2363172.8975462.4842064.4415862.9721660.85470

25231677.5618163.0636263.7313980.6620563.0707266.5130063.3311360.34233

>head(mm(data.raw)-mm(data.rma),2)#差值为全部为0，说明rma方法没有对mm数值进行处理

CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8

50184300000000

25231600000000

二、归一化处理（normalization）

同一个RNA样品用相同类型的几块芯片进行杂交，获得的结果（信号强度等）都不可能完全相同，甚至差别很大。为了使不同芯片获得的结果具有可比性，必需进行归一化处理。这一步的方法也很多。归一化处理的affy函数为normalize( )，以Affybatch对象和处理方法为参数。

1、线性缩放方法：这是Affy公司在其软件（4.0和5.0版本）中使用的方法。这种方法先选择一个芯片作为参考，将其他芯片和参考芯片逐一做线性回归分析，用回归直线（没有截距）对其他芯片的信号值做缩放。Affy公司的软件做回归分析前去除了2%最强和最弱信号。使用很简单，mas方法做背景消减后做归一化分析的R语句是：

>data.mas.ln<-normalize(data.mas,method="constant")

>head(pm(data.mas.ln)/pm(data.mas),5)

CHIP1    CHIP2    CHIP3    CHIP4    CHIP5    CHIP6    CHIP7    CHIP8

50113110.87876461.0021781.1405490.99288791.0294440.75775370.8833788

25160410.87876461.0021781.1405490.99288791.0294440.75775370.8833788

26189110.87876461.0021781.1405490.99288791.0294440.75775370.8833788

23038710.87876461.0021781.1405490.99288791.0294440.75775370.8833788

21733410.87876461.0021781.1405490.99288791.0294440.75775370.8833788

>head(mm(data.mas.ln)/mm(data.mas),5)

CHIP1    CHIP2    CHIP3    CHIP4    CHIP5    CHIP6    CHIP7    CHIP8

50184310.87876461.0021781.1405490.99288791.0294440.75775370.8833788

25231610.87876461.0021781.1405490.99288791.0294440.75775370.8833788

26260310.87876461.0021781.1405490.99288791.0294440.75775370.8833788

23109910.87876461.0021781.1405490.99288791.0294440.75775370.8833788

21804610.87876461.0021781.1405490.99288791.0294440.75775370.8833788

可以看出，线性缩放方法以第一块芯片为参考，它的数值没有被处理，而其他芯片都被缩放了。对同一块芯片，不同探针的缩放倍数是一个常数。PM和MM的缩放方法完全一样。

2、非线性缩放方法：此方法获得的结果比线性方法要好，做非线性拟合时不是取整张芯片而仅取部分（一列）作为基线。

>data.mas.nl<-normalize(data.mas,method="invariantset")

>head(pm(data.mas.nl)/pm(data.mas),5)

CHIP1    CHIP2    CHIP3    CHIP4    CHIP5    CHIP6     CHIP7     CHIP8

50113111.0496001.4169841.3592821.3991642.0231091.00857861.3123854

25160411.3483262.2789291.8378301.5866652.4305501.11258041.3059096

26189111.5640441.3973761.6751831.4966641.5564771.31673551.5508943

23038711.2587541.6829281.3898361.3601181.5038931.01445911.2238827

21733411.0093941.1265551.2414241.2294021.2629490.84172530.9933603

可以看到，同一芯片不同探针的信号值的缩放倍数是不一样的。

3、分位数（quantile）方法：这种方法认为（或假设）每张芯片探针信号的经验分布函数应完全一样，使用任两张芯片的数据做QQ图应该得到一条斜率为1截距为0的直线。

>data.mas.qt<-normalize(data.mas,method="quantiles")

>head(pm(data.mas.qt)/pm(data.mas),5)

CHIP1     CHIP2    CHIP3    CHIP4    CHIP5    CHIP6     CHIP7     CHIP8

5011310.71763740.96021341.1402331.1807011.1115931.1868050.81452381.0155723

2516040.69840190.98141221.1985101.2091941.1117261.1722720.82872131.0143434

2618910.69387050.99556101.1367341.2050781.1106201.1856080.83887291.0579171

2303870.76529230.97767551.1711231.1830991.1148641.1850520.81528250.9921032

2173340.95075450.95468441.0633201.1620731.1304111.1729860.79450020.9265588

4、其他：如循环局部加权回归法（Cyclic loess）和 Contrasts方法。

data.rma.lo<-normalize(data.rma,method="loess")

data.rma.ct<-normalize(data.rma,method="contrasts")

三、汇总（Summarization）：

最后一步汇总是使用合适的统计方法通过probeset（包含多个探针）的杂交信号计算出计算表达量。affy包的函数为computeExprSet。需要注意的是computeExprSet函数除需要指定统计方法外还需要指定PM校正的方式：computeExprSet(x, pmcorrect.method, summary.method, ...)

两个参数可以设定的值可以通过下面函数获得：

>pmcorrect.methods()

[1]"mas""methods""pmonly""subtractmm"

>generateExprSet.methods()

[1]"avgdiff""liwong""mas""medianpolish""playerout"

常用的汇总方法是medianpolish, liwong和mas。liwong方法仅使用PM做背景校正（pmcorrect.method="pmonly"）。例如：

>eset.rma.liwong<-computeExprSet(data.rma.lo,pmcorrect.method="pmonly",summary.method="liwong")

22810ids to be processed

||

|####################|

共有29个警告(用warnings()来显示)

计算过程需要一定的时间，耐心等候完成。最后的结果 ExpressionSet 类型数据：

>class(eset.rma.liwong)

[1]"ExpressionSet"

attr(,"package")

[1]"Biobase"

>class(data.raw)

[1]"AffyBatch"

attr(,"package")

[1]"affy"

四、三合一处理和“自动化”函数：

了解芯片预处理的原理和步骤后你完全可以用一个R函数完成三步处理。affy包提供的三合一处理函数为 expresso( )，用法为：

expresso(

afbatch,

# background correction

bg.correct = TRUE,

bgcorrect.method = NULL,

bgcorrect.param = list(),

# normalize

normalize = TRUE,

normalize.method = NULL,

normalize.param = list(),

# pm correction

pmcorrect.method = NULL,

pmcorrect.param = list(),

# expression values

summary.method = NULL,

summary.param = list(),

summary.subset = NULL,

# misc.

verbose = TRUE,

widget = FALSE)

例如：

eset.mas<-expresso(data.raw,bgcorrect.method="mas",normalize.method="constant",

pmcorrect.method="mas",summary.method="mas")

使用的各类方法可用bgcorrect.methods，pmcorrect.methods 和 express.summary.stat.methods函数了解。

expresso速度较慢，一些R软件包提供速度较快的预编译三合一函数，如affyPLM软件包的threestep函数。affyPLM提供的处理方法很丰富，有兴趣可以自学。

下面介绍介3个“自动化”函数，这些函数已经预定义了预处理各步骤所采用的方法和参数，不再需要设定。

1、mas5函数，由affy包提供：

>eset.mas5<-mas5(data.raw)

background correction:mas

PM/MM correction:mas

expression values:mas

background correcting...done.

22810ids to be processed

||

|####################|

mas5据说是expresso函数和mas方法的封装。但使用下面语句获得的结果似乎与 mas5(data.raw)的结果不一样（自己去验证看看）：

expresso(data.raw, bgcorrect.method="mas", normalize=FALSE, pmcorrect.method="mas", summary.method="mas")

2、rma函数，也是由affy包提供，其背景处理方法为rma法，归一化处理使用分位数法，而汇总方法使用medianpolish：

eset.rma<-rma(data.raw)

它等价于：

eset.rma<-expresso(data.raw,bgcorrect.method="rma",normalize.method="quantiles",pmcorrect.method="pmonly",summary.method="medianpolish")

但rma函数是预编译好的C语言函数，速度非常快！

3、gcrma函数，由gcrma包提供：

Affy的软件（比如mas方法）使用MM数据做背景处理，但由于MM出现的问题（上面提到过），这些方法可能高估了背景值。而rma方法在做背景处理时没有使用MM数据，这可能又低估了背景值。MM序列公布后有人对其特异性进行了评估，并使用这些评估结果建立了新方法。gcrma就是这样的方法，也是封装好的三合一方法。

>eset.gcrma<-gcrma(data.raw)

Adjustingforoptical effect........Done.

Computingaffinities.Done.

Adjustingfornon-specific binding........Done.

Normalizing

CalculatingExpression

五、芯片处理方法的优劣评估

芯片预处理的方法这么多，哪个好？我选哪个？知道得越多越迷惑。幸好这些已经有人做了，牛人Rafael Irizarry 和 Leslie Cope专门写了一个R软件包affycomp用于方法评估。

评估方法的优劣必需有数据，而且是包含已知因素的数据。affycomp需要两个系列的数据，一个是RNA稀释系列芯片数据，称为Dilution data，另一个是使用了内参/外标RNA的芯片，称为Spike-in data。Spike-in RNA是在目标物种中不存在、但在芯片上含有相应检测探针的RNA，比如Affy的拟南芥芯片上有几个人或细菌的基因检测探针。由于稀释倍数已知，内参/外标的RNA量和杂交特异性也已知，所以结果可以预测，也就可以用做方法评估。对于严格的芯片实验来说，这些实验都是必须的。但是绝大多数人不做，因为成本很高，尤其是只做几张芯片的时候，一般直接使用别人认可的方法如RMA或MAS。下面是affycomp说明文档比较MAS5和RMA方法优劣的一个演示图，软件具体用法此处不做介绍：