基因芯片技术的特点使用寡聚核苷酸探针检测基因。前一节使用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方法优劣的一个演示图，软件具体用法此处不做介绍：
网友评论