科研有捷徑，輸入代碼，一鍵獲取科研成果！就是這么省事，來具體看下有多方便！

搜索http://985.so/a9kb查看全部代碼（目前共計50+持續(xù)新增中），也可以點擊右側(cè)【目錄】，可以看到更多有趣的代碼

真香提示：文末可以知道如何獲取代碼~?

在處理高通量數(shù)據(jù)過程中，我們往往需要借助一些表達譜信息進行比對分析，通過不同類型樣本之間的比較，獲得與疾病表型存在顯著相關(guān)性的基因標識。

GEO數(shù)據(jù)庫http://www.ncbi.nlm.nih.gov/geo/為我們提供了大量和疾病相關(guān)的表達譜信息，其包含物種豐富（人，鼠，?植物等），表達譜類型多樣（芯片數(shù)據(jù)，甲基化數(shù)據(jù)， RNA-seq數(shù)據(jù)等），并且為使用者提供了各種可下載的數(shù)據(jù)資源。其中比較常用的一種下載資源是CEL格式文件，即原始數(shù)據(jù)，cel文件就是Affymetrix掃描儀產(chǎn)出的原始數(shù)據(jù)文件，記錄了每個probe 的signal intensity。

這里我們通過從GEO數(shù)據(jù)庫下載的一個樣例，逐步介紹RMA進行CEL表達譜數(shù)據(jù)標準化，以及l(fā)imma算法提取差異表達基因的過程。

首先我們從GEO下載瘢痕瘤的表達譜數(shù)據(jù)GSE7890，下載地址http://www.ncbi.nlm.nih.gov/geo/query/acc.cgiacc=GSE7890，下載CEL文件到本地，解壓，這里默認目錄為F:\。之后調(diào)用R語言的RMA包，標準化過程如下：

直方圖為標準化前后數(shù)據(jù)分布的比較，可以觀察出基因的以2為底對數(shù)轉(zhuǎn)換后的信號值集中在什么區(qū)間

箱式圖為RMA方法消除樣本間差異前后的數(shù)據(jù)比較，紅色圖代表處理前數(shù)據(jù)分布，藍色為處理后數(shù)據(jù)分布，可見RMA方法可以有效去除樣本間差異。

實現(xiàn)代碼如下：

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

biocLite("affy")

biocLite("Biobase")

biocLite("tkWidgets")

library(affy)

AffyData<-ReadAffy(widget=TRUE)?????

exprsSet.RMA<-rma(AffyData)???##用rma方法處理數(shù)據(jù)

exp.RMA<-exprs(exprsSet.RMA)??##?exprsSet.MAS5<-mas5(AffyData)??##用mas5方法處理數(shù)據(jù)

write.exprs(exprsSet.RMA,?file="mydata.txt")??##結(jié)果保存

hist(log2(intensity(AffyData[,4])),breaks=100,col="blue")????##畫芯片數(shù)據(jù)的直方圖，因為實驗數(shù)據(jù)變化范圍大所以取log2，結(jié)果集中在6~8，大于10?的很少，說明樣本僅少部分高表達，與實際吻合。

par(mfrow=c(1,2))

boxplot(AffyData,col="red")

boxplot(data.frame(exp.RMA),col="blue")?##?畫箱式圖，比較數(shù)據(jù)分布情況

plot(exp.RMA[,1:2],log="xy",pch=".",main="all")

plot(pm(AffyData)[,1:2],log="xy",pch=".",main="pm")

plot(mm(AffyData)[,1:2],log="xy",pch=".",main="mm")

plot(exp.RMA[,1:2],pch=".",main="exp.RMA")??##?散點圖評價兩個樣本之間的數(shù)據(jù)差異

library(genefilter)???##?函數(shù)rowSds在genenfliter中

rsd<-rowSds(exp.RMA)???##?計算行的標準差

sel<-order(rsd,decreasing=TRUE)[1:20]??##?選擇表達差異高的top20，order返回行號

heatmap(exp.RMA[sel,],col=rainbow(256))?##?畫熱圖

之后我們獲得了標準化后的表達譜數(shù)據(jù)，接下來我們可以利用limma package的算法對表達譜數(shù)據(jù)進行差異基因提取。

需要設(shè)置的參數(shù)包括表達譜數(shù)據(jù)的文件目錄，case組和control組的數(shù)據(jù)，以及輸出文件的文件名，這里默認為case_control_limma.txt。那么對于我們的樣例文件，我們將其放置于F：\盤下，case組（treated with**）為9個，control組（treated without**）為10個，實現(xiàn)的代碼如下：

library(limma)

setwd("F:\\")

exp<-read.table("profile_normalized.txt",header=T,sep="\t")

rownames(exp)<-exp[,1]

exp<-exp[2:20]

exp<-as.matrix(exp)

samps<-factor(c(rep("case",9),rep("control",10)))

design?<-?model.matrix(~0+samps)?;

colnames(design)?<-?c("case","control")?

fit?<-?lmFit(exp,?design)?

cont.matrix<-makeContrasts(case-control,levels=design)?

fit2?<-?contrasts.fit(fit,?cont.matrix)?

fit2?<-?eBayes(fit2)?

final<-topTable(fit2,?coef=1,?number=dim(exp)[1],?adjust.method="BH",?sort.by="B",?resort.by="M")

write.table(final,paste("case_control","_limma.txt"),quote=FALSE,sep="\t")

最后輸出的結(jié)果如下圖所示：

表中共包含8列，分別代表探針I(yè)D，校正后P值，校正前P值，t/B統(tǒng)計量，對數(shù)轉(zhuǎn)換的fold change值，gene symbol以及基因名稱。接下來使用者可根據(jù)需要，設(shè)定P值和logFC的閾值，從而提取出顯著差異表達的基因。

當然數(shù)據(jù)標準化和提取差異基因的方法很多，例如Z-score標準化，通過計算均值和標準差進行零均值校正，提出基因固有表達差異的影響也是一種常用的標準化方法，另外還有類似Python語言中preprocessing包中的scale，normalize等方法。

差異基因的提取方法還包括student T test，秩和檢驗，SAM package，甚至更高級一些的機器學習中的feature selection也可以實現(xiàn)這一目的。使用者可根據(jù)需要或數(shù)據(jù)類型進行選擇。