接前一篇： 用R和BioConductor進行基因芯片數據分析(三)：計算median
歸一化是從normalization翻譯過來的。歸一化的目的是使各次/組測量或各種實驗條件下的測量可以相互比較，消除測量間的非實驗差異。非實驗差異可能來源於樣品制備，點樣，雜交過程，雜交信號處理等。
歸一化的方法有很多，對於寡聚核苷酸芯片（單通道，以Affymetrix為代表）和cDNA芯片（雙通道，紅綠染料）也有所不同。以下討論針對雙通道芯片進行，當然也可能適用於單通道，請讀者自辨。
歸一化通常分為bulk normalization和control-based normalization. 前者假定僅有一小部分基因表達值在不同條件下有差異，因此使用所有的基因作為標准進行歸一化;而后者使用表達值被認為是不變的那些對照基因作為標准進行歸一化。這2種方法的假設前提都未必時時成立，因此需要具體情況具體分析。
“bulk” normalization又細分為許多方法，最簡單的就是global normalization,這種方法假設紅色染料的信號強度與綠色染料的信號強度是正相關的，即R=kG (R：紅色信號強度，G：綠色信號強度，k：常數），因此表達信號強度的以2為底的對數比log ratio在歸一化之后相當於平移了一個常量c=logk：log(R/G) → log(R/G) – c = log(R/G) – logk = log[R/(kG)].
c的計算方法通常是用所有紅色信號強度和除以所有綠色信號強度和然后取以2為底的對數，即c=log[total(R)/total(G)]
Intensity-dependent normalization通常比global normalization更好，因為后者的假設通常並不完全正確，通常情況下，log ratio是和信號強度值有關的，即log(R/G) → log(R/G) – c(A)，這里A=1/2*log(R*G)，是log product intensity或簡稱log intensity。
對於我們的數據mediandata，我們可以直觀的看到log ratio和信號強度的關系：
上圖叫做MA-plot(MA圖)，縱坐標是M=log(R/G)=log(new/old),橫坐標是A=1/2*log(R*G)=1/2*log(new*old)。
其中的藍色曲線是lowess回歸函數(什么是lowess)。（注：由於原始數據有5行有0值，導致有些M或A=Inf或-Inf，無法進行Lowess回歸，因此人工刪除了這5行，處理后的中位數數據mediandata在這里下載。當然你也可以用原始數據求出M和A值，自己刪除Inf值對應的mediandata中的行。
在R中畫出上圖的代碼：
mediandata<-read.table('mediandata.txt',header=TRUE)
mediandata<-mediandata[,-1] #移除第一列ID列
MA<-matrix(data = NA, nrow =dim(mediandata)[1], ncol = dim(mediandata)[2], byrow = TRUE, dimnames = NULL)
new<-0
old<-0
for (i in 1:dim(mediandata)[1]){
for (j in 1:(dim(mediandata)[2]/2)){
new<-mediandata[i,2*j-1]
old<-mediandata[i,2*j]
MA[i,2*j]<-log(new/old)/log(2) #M=log(new/old)/log2
MA[i,2*j-1]<-1/2*log(new*old)/log(2) #A=1/2*log(new*old)/log2
}
}
plot(MA[,1],MA[,2],xlab='A',ylab='M')
lines(lowess(MA[,1],MA[,2],f=0.2,iter=2),lwd=2,col='blue')
#僅畫出第一張芯片上2個雜交樣本的MA圖，使用MA[,3],MA[,4]可以畫出第2張芯片的圖，依此類推。
可以看出原始數據的log ratio受到log intensity的影響，因此需要intensity-based normalization.
R的lowess函數返回一個$y對象，儲存每個A值(lowess返回的$x對象）對應的~M值，而歸一化后的M’=M-~M=M-$y($x)
這樣歸一化后得到10個芯片的log ratio即10列數據，但是為了后面分析的方便，我想得到10個new的強度值和10個old的強度值共20列數據，那怎么辦呢？
答案很簡單，就是假設new的強度值在歸一化后不變，而僅僅改變old的強度值，得到old’=old*2^($y($x)) 注：$y($x)是2的指數，推導很簡單
下面是R代碼：
normed<-matrix(data = NA, nrow =dim(MA)[1], ncol = dim(MA)[2], byrow = TRUE, dimnames = NULL) #new—奇數列; old—偶數列
for (j in 1:(dim(MA)[2]/2)){
out_lowess<-lowess(MA[,2*j-1],MA[,2*j],f=0.2,iter=2)
#A=MA[,1],M=MA[,2]
loc_lowess<-cbind(out_lowess$x,out_lowess$y)
for (i in 1:dim(MA)[1]){
normed[i,2*j-1]<-mediandata[i,2*j-1] #歸一化后的new’強度值
normed[i,2*j]<-mediandata[i,2*j]*2^(loc_lowess[,2][loc_lowess[,1]==MA[i,2*j-1]][1]) #歸一化后的old’強度值
}
}
搞定，看看效果吧：
MAnormed<-matrix(data = NA, nrow =dim(MA)[1], ncol = 2, byrow = TRUE, dimnames = NULL)
MAnormed[,2]<-log(normed[,1]/normed[,2])/log(2) #M=log(new/old)/log2
MAnormed[,1]<-1/2*log(normed[,1]*normed[,2])/log(2) #A=1/2*log(new*old)/log2
plot(MAnormed[,1],MAnormed[,2],xlab='A',ylab='M')
lines(lowess(MAnormed[,1],MAnormed[,2],f=0.2,iter=2),lwd=2,col='blue')
現在的lowess回歸曲線是一條直線了，說明歸一化后的log ratio與強度值無關了。
可以從另一個角度看歸一化的效果：
plot(density(normed[,1]),type='line',col='red',xlab='intensity')
points(density(normed[,2]),type='line',col='green')
points(density(mediandata[,1]),type='line',col='blue')
points(density(mediandata[,2]),type='line',col='black')
text(2.2,c(0.09,0.11,0.13,0.15),labels=c('BEFORE normalization black','BEFORE normalization blue','AFTER normalization green','AFTER normalization red'),col=c('black','blue','green','red'))
關於芯片內歸一化，可以參考以下資料：
cDNA雙通道芯片：
Yang Y.H., Dudoit S., Luu P., Speed T.P. (2001) Normalization for cDNA microarray data, SPIE BiOS 2001, San Jose CA;
Yang Y.H., Dudoit S., Luu P., Lin D.M., Peng V., Ngai J., Speed T.P. (2002), Normalization for cDNA microarray data: a robust composite method addressing single and multiple slide systematic variation, Nucleic Acids Research 30(4);
寡聚核苷酸單通道芯片：
Bolstad B.M., Irizarry R. A., Astrand M., Speed T.P. (2003) A comparison of normalization methods for high density oligonucleotide array data based on bias and variance, Bioinformatics 19(2): 185-193