應(yīng)學(xué)生及個別博友的要求，盡管專業(yè)博文點(diǎn)擊率和反應(yīng)均很差，但在去San Diego參加PAG會議之前罗心，還是抽時(shí)間給出【R高級教程】的第二專題里伯。專題一給出了聚類分析的示例，本專題主要談在表達(dá)譜芯片分析中如何利用Bioconductor鑒定差異表達(dá)基因渤闷。
鑒定差異表達(dá)基因是表達(dá)譜芯片分析pipeline中必須的分析步驟疾瓮。差異表達(dá)基因分析是根據(jù)表型協(xié)變量（分類變量）鑒定組間差異表達(dá)，它屬于監(jiān)督性分類的一種肤晓。在鑒定差異表達(dá)基因以前爷贫，一般需要對表達(dá)值實(shí)施非特異性過濾（在機(jī)器學(xué)習(xí)框架下屬于非監(jiān)督性分類）认然，因?yàn)檫m當(dāng)?shù)姆翘禺愋赃^濾可以提高差異表達(dá)基因的檢出率补憾、甚至是功效。R分析差異表達(dá)基因的library有很多卷员，但目前運(yùn)用最廣泛的Bioconductor包是limma盈匾。
本專題示例依然來自GEO數(shù)據(jù)庫中檢索號為GSE11787 的Affymetrix芯片的數(shù)據(jù)，數(shù)據(jù)介紹參閱專題一毕骡。
>library(limma)
>design <- model.matrix(~ -1+factor(c(1,1,1, 2,2,2)))
這個是根據(jù)芯片試驗(yàn)設(shè)計(jì)削饵，對表型協(xié)變量的水平進(jìn)行design，比如本例中共有6張芯片未巫，前3張為control對照組窿撬，后3張芯片為實(shí)驗(yàn)處理組，用1表示對照組叙凡，用2表示處理組劈伴。其他試驗(yàn)設(shè)計(jì)同理，比如2*2的因子設(shè)計(jì)試驗(yàn)握爷，如果每個水平技術(shù)重復(fù)3次跛璧，那么可以表示為：design <- model.matrix(~ -1+factor(c(1,1,1, 2,2,2, 3,3,3, 4,4,4)))。接上面的程序語句繼續(xù)：
>colnames(design) <- c("control", "LPS")
>fit <- lmFit(eset2, design)
>contrast.matrix <- makeContrasts(control-LPS, levels=design)
>fit <- eBayes(fit)
>fit2 <- contrasts.fit(fit, contrast.matrix)
>fit2 <- eBayes(fit2)
>results<-decideTests(fit2, method="global", adjust.method="BH", p.value=0.01, lfc=1.5)
>summary(results)
>vennCounts(results)
>vennDiagram(results)
比較遺憾的是新啼，目前l(fā)imma自帶的venn作圖函數(shù)不能做超過3維的高維venn圖追城，只能畫出3個圓圈的venn圖，即只能同時(shí)對三個coef進(jìn)行venn作圖燥撞。上面的venn圖只有一個coef座柱，太簡單了。下面是一個由本實(shí)驗(yàn)室芯片數(shù)據(jù)得出的三個coef的venn圖例：
>heatDiagram(results,fit2$coef)
紅色為control中（與LPS相比）的高表達(dá)基因物舒，綠色為control中（與LPS相比）的低表達(dá)基因色洞，x軸的數(shù)字表示差異表達(dá)基因在eset2中所處的位置。
>x<-topTable(fit2, coef=1, number=10000, adjust.method="BH", sort.by="B", resort.by="M")
>write.table(x, file="limma.xls", row.names=F, sep="\t")
將結(jié)果寫入limma.xls文件中茶鉴，內(nèi)容包括AveExpr值（比較組間絕對值的平均差異值）锋玲、logFC值（差異倍數(shù)）、t值涵叮、P值惭蹂、q值（即adj.P.Val值）和B值伞插。一般logFC值、P值盾碗、q值和AveExpr值用來作為差異表達(dá)的判斷標(biāo)準(zhǔn)媚污，比如差異倍數(shù)在2倍以上、絕對差異表達(dá)在10以上廷雅、P值小于0.01等耗美。在Excel文件中，根據(jù)各項(xiàng)判斷標(biāo)準(zhǔn)排序航缀，可以很容易地得到差異表達(dá)基因列表商架，這個列表可以用來進(jìn)行后續(xù)的分析，如GO注釋芥玉、基因網(wǎng)絡(luò)繪制等蛇摸。
專題一中提到實(shí)際研究中，一般只用差異表達(dá)基因進(jìn)行聚類分析灿巧，在R中赶袄，根據(jù)差異表達(dá)結(jié)果過濾表達(dá)值很簡單（具體的值可以依據(jù)芯片數(shù)據(jù)的實(shí)際情況設(shè)定，比如P值可以設(shè)寬松點(diǎn)0.05抠藕、logFC的絕對值也可設(shè)為1或2饿肺、絕對表達(dá)差異也可設(shè)低一點(diǎn)，如6或8這樣的值）：
>y <- x[x$P.Value < 0.01 & (x$logFC > 1.5 | x$logFC < -1.5& x$AveExpr > 10),]
>length(y$ID)
>eset3<-eset2[y$ID,]
經(jīng)過上面P值盾似、表達(dá)倍數(shù)差異和絕對差異的過濾敬辣，eset3中就只包含差異表達(dá)基因了，這樣eset3可用來進(jìn)行聚類分析了颜说。
【除Bioconductor的logo购岗，所有圖片均由程序運(yùn)行所得】
本文引用地址：http://blog.sciencenet.cn/blog-295006-403640.html此文來自科學(xué)網(wǎng)朱猛進(jìn)博客，轉(zhuǎn)載請注明出處门粪。
