接前一篇： 用R和BioConductor進行基因芯片數據分析(五)：芯片間歸一化
經過一系列的預處理，包括缺失值填充，中位數計算以及歸一化，我們的數據終於可以用啦。
下面我們就來分析一下new population和old population的個體是否有差異表達基因。
判斷一個基因是否差異表達有許多方法，最早使用的就是看log ratio的絕對值是否大於2,這種方法早已廢棄。
下一個想到的也許是t-test，誠然t-test可以統計地判斷一個基因是否差異表達，但是對於有數千數萬基因的芯片來說，會有很高的錯誤發現率（False Discovery Rate, FDR)，如果 p value < 0.05,則10000個基因里有500個基因實際沒有差異表達而被誤認為是差異表達。因此t-test方法需要改進。
於是 Westfall & Young (1993) 提出了Step-down maxT and minP multiple testing procedures，大意就是比較幾個group間有沒有差異基因表達，就通過隨機置換這些group的標記，相當於隨機互換group的成員，模擬一個空分布(null distribution），以此計算調整后的p value，這個方法可以極大的減小Family-wise Error Rate (FWER).
以下分析就使用Step-down maxT and minP multiple testing procedures，需要求助於Bioconductor的multtest package的mt.maxT()函數。具體用法可通過help(mt.maxT)查看。
source("http://bioconductor.org/biocLite.R")
biocLite("multtest")
library(multtest)
classlabel<-c(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1)
maxttt<-mt.maxT(norm_log_btw_array,classlabel,B=100000)
默認隨機置換次數B=10000，對於microarray來說B應該比10000大很多，所以這里取B=100000
以下是畫圖：
rawp<-maxttt$rawp[order(maxttt$index)]
plot(sort(rawp),type='p',col=1,ylim=c(-0.05,1.00),ylab='p value')
lines(sort(maxttt$adjp),type='p',col='red')
#min adj-p: sort(maxttt$adjp)[1] 0.0163
#rawp: >sort(rawppp)[170] [1] 0.0493 > sort(rawppp)[171] [1] 0.0502 170個raw p小於0.05
abline(h=0.05,col='blue')
text(1000,c(0.6,0.7),labels=c('raw p-value','adjusted p-value'),col=c('black','red'))
text(1000,0.08,labels='p=0.05',col='blue')
可見調整后只有一個基因的p value小於0.05，而未調整的有170個基因的p value小於0.05，可以說雖然此方法降低了錯誤發現率，但是也導致了很高的False negative.
此外可以考慮使用multtest package的mt.rawp2adjp()函數，這個函數可以通過”Bonferroni”, “Holm”, “Hochberg”, “SidakSS”, “SidakSD”, “BH”, “BY”等方法調整p value，不過對我們的數據來說都過於嚴格了。
procs<-c("Bonferroni","Holm","Hochberg","SidakSS","SidakSD","BH","BY")
adjps<-mt.rawp2adjp(rawp,procs)
plot(sort(adjps$adjp[,1]),ylab='p value')
for (i in 2:8){
points(sort(adjps$adjp[,i]),col=i)
}
abline(h=0.05,col='blue')
text(1000,0.08,labels='p=0.05',col='blue')
因此可以考慮不這么嚴格的SAM (Significance Analysis of Microarrays)分析。有興趣的請看參考資料。。
參考資料：
課堂講義：Differential expression. Identifying differentially expressed genes — notions on multiple testing and p-value adjustments.
Dudoit, S., Yang, Y.H., Speed, T.P., and Callow, M.J. (2002), Statistical methods for identifying differentially expressed genes in replicated cDNA microarray experiments, Statistica Sinica 12(1):111-139.
Dudoit S., Shaffer J.P., Boldrick J.C. (2003). Multiple hypothesis testing in microarray experiments, Statistical Science, 18(1): 71-103.
Efron B., Tibshirani, R., Storey J.D., and Tusher V. (2001), Empirical Bayes analysis of a microarray experiment, Journal of the American Statistical Association 96:1151-1160.
SAM (Significance Analysis of Microarrays)相關：Tusher, V.G., Tibshirani, R., and Chu, G. (2001) Significance analysis of microarrays applied to the ionizing radiation response, PNAS 98:5116-5121.
FDR相關：Storey J.D. (2002), A direct approach to false discovery rates, JRSS-B 64(3):479-498.
from : http://azaleasays.com/tag/r/