limma包的使用技巧

1、limma包的使用技巧原文地址：limma包的使用技巧作者：牛牛龍limmarpackage是一個(gè)功能比較全的包，既含有cDNA芯片的RAWdata輸入、前處理(歸一化)功能，同時(shí)也有差異化基因分析的“線性”算法(limma:LinearModelsforMicroarrayData)，特別是對于“多因素實(shí)驗(yàn)(multifactordesignedexperiment)”。limmar包的可擴(kuò)展性非常強(qiáng)，單通道(onechannel)或者雙通道(towchannel)數(shù)據(jù)都可以分析差異基因，甚至也包括了定量PCR和RNA-seq(第一次見分析microarray的包也能處理RNA-seq)。l

2、immarpackage是一個(gè)集大成的包，對載入數(shù)據(jù)、數(shù)據(jù)前處理(背景矯正、組內(nèi)歸一化和組間歸一化都有很多種選擇方法)、差異基因分析都有很多的選擇。而且，所設(shè)計(jì)的線性回歸和經(jīng)驗(yàn)貝葉斯方法找差異基因非常值得學(xué)習(xí)。1.讀入樣本信息使用函數(shù)讀readTargets(file="Targets.txt",path=NULL,sep="t",s=NULL,quote=""",.)。這個(gè)函數(shù)其實(shí)是一個(gè)包裝了的read.table()，讀入的是樣本的信息，創(chuàng)建的對象類似于marray包的marrayInfos和Biobas

3、e包的AnnotatedDataFrame。2.讀入探針密度數(shù)據(jù)與marray包一致，Bioconductor不能讀入原始的TIFF圖像文件，只能輸出掃描儀輸入的、轉(zhuǎn)換成數(shù)字信號的文本文件。使用函數(shù)read.maimages（files=NULL,source="generic",path=NULL,ext=NULL,names=NULL,columns=NULL,other.columns=NULL,annotation=NULL,green.only=FALSE,wt.fun=NULL,verbose=TRUE,sep="t",quote=NULL,

4、.）參數(shù)說明：files需要通過函數(shù)dir（pattern="Mypattern"）配合正則表達(dá)式篩選（規(guī)范命名很重要），同時(shí)該函數(shù)可以讀入符合格式的壓縮過的文件，比如*.txt.gz的文件，這極大的減小的數(shù)據(jù)儲(chǔ)存大?。籹ource的取值分為兩類，一類是“高富帥”，比如“agilent"、“spot”等等（下表），它們是特定掃描儀器的特定輸由格式；如果不幸是“展絲”，即格式是自己規(guī)定的，可以選定source="generic"，這時(shí)需要規(guī)定columns；任何cDNA文件都要有R/G/Rb/Gb四列（Mean或者M(jìn)edian）；annotati

5、on可以規(guī)定哪些是注釋列；wt.fun用于對點(diǎn)樣點(diǎn)進(jìn)行質(zhì)量評估，取值為0表示這些點(diǎn)將在后續(xù)的分析中被剔除，取值位1表示需要保留，對點(diǎn)樣點(diǎn)的評估依賴于圖像掃描軟件的程序設(shè)定，比如SPOT和GenePix軟件，查看QualityWeights(現(xiàn)成函數(shù)或者自己寫函數(shù))。讀入單通道數(shù)據(jù)：讀入單通道數(shù)據(jù)，可以設(shè)定green.only=TRUE即可，然后對應(yīng)讀入columns=list(G="Col1",Gb="Col2")。讀入的數(shù)據(jù)，如果是單通道，則成為EListRawclass；如果是雙通道，則是RGListclass。數(shù)據(jù)操作：cbind()：合并數(shù)據(jù)；“

6、”：分割數(shù)據(jù)；RGListclass有的names是"R","G","Rb","Gb"weights”printer","genes"，"targets"，"notes"：R/G/Rb/Gb分別紅和綠的前景和背景噪音；weight是掃描軟件的質(zhì)量評估；printer是點(diǎn)樣規(guī)則(printerlayout)；genes是基因注釋；target是樣本注釋；notes是一般注釋。可以通過myRGList$names進(jìn)行相應(yīng)的取值和賦值。3. 前處理cD

7、NA芯片前處理的流程是：背景去除(backgroundcorrection)->組內(nèi)歸一化(with-arraynormalization)->組間歸一化(between-arraynormalization)。最后一步組間歸一化可選，注意trade-off原則。3.1 背景去除：backgroundCorrect(RG,method="auto",offset=0,printer=RG$printer,normexp.method="saddle",verbose=TRUE)RG：可以是EListRaw，也可以是RGList

8、class；method:取“auto”，“none”，“subtract"，"half","minimum","movingmin","edwards"或者"normexp"normexp.method：在method取“normexp”時(shí)，可以取"saddle"，"mle"，"rma"或者"rma75"。作者建議使用“mle”或者“saddle"，兩個(gè)運(yùn)行速度也差不多。如果數(shù)據(jù)中含有“形態(tài)背景

9、估計(jì)(morphologicalbackgroundestimates)”，比如SPOT和GenePix圖像處理軟件，那么method="subtract"也就較好的表現(xiàn)。3.2組內(nèi)歸一化normalizeWithinArrays(object,layout,method="printtiploess",weights=object$weights,span=0.3,iterations=4,controlspots=NULL,df=5,robust="M",bc.method="subtract",offset=

10、0)method："none"，"median"，"loess"，"printtiploess"，"composite"，"control"和"robustspline"。初始值設(shè)定為“printtiploess”，但是對于Agilent芯片或者小樣本芯片（每個(gè)“點(diǎn)樣塊”少于150個(gè)探針）?？捎眠x用的loess或者robustspline。"loess歸一化方法假設(shè)了有相當(dāng)大的一部分探針沒有發(fā)生差異化表達(dá)，而不是上調(diào)或者下調(diào)的基因數(shù)差不多或者差異

11、化程度圍繞著0波動(dòng)?！保▍⒖嘉墨I(xiàn)1）。需要注意，組內(nèi)歸一化只是對單張芯片的M值進(jìn)行了規(guī)整（不影響A值），但是對與組間各個(gè)通道沒有進(jìn)行比較。weight：是圖像處理軟件對探針權(quán)重的標(biāo)記，如果使用weight，那么在歸一化過程中weight為0的探針不會(huì)影響其他探針，這并不意味著這些探針會(huì)被剔除，它們同樣會(huì)被歸一化，也會(huì)出現(xiàn)在歸一化的結(jié)果中。如果不想使用，那么weight設(shè)為NULL。iterations：設(shè)定loess的循環(huán)數(shù)，循環(huán)數(shù)越多，結(jié)果越強(qiáng)?。╮obust）。3.3 組間歸一化normalizeBetweenArrays(object,method=NULL,targets=NULL,c

12、yclic.method="fast",.)method：“none"，"scale"，"quantile"，"Aquantile"，"Gquantile"，“Rquantile"，“Tquantile"，“cyclicloess"?！皊cale”對log-ratio數(shù)值進(jìn)行歸一化；“quantile”對密度(intensity)進(jìn)行歸一化；“Aquantile”對A數(shù)值進(jìn)行歸一化，不調(diào)正M數(shù)值；"Gquantile"和“Rquanti

13、le”則分別對Green和Red通道進(jìn)行歸一化，適用于Green或者Red為“參考標(biāo)準(zhǔn)(commonreference)”的樣品；“Tquantile”表示樣品分開歸一化或者Green和Red一先一后進(jìn)行歸一化。以下是一個(gè)cDNA芯片的密度圖，分別為原始、去背景(method="normexp",normexp.method="mle")、組內(nèi)歸一化(method="robustspline")和組間歸一化(method=“quantile”)，使用plotDensities()繪制。4.線性模型設(shè)計(jì)(linearmodeldesig

14、h)“線性模型”主要分析差異化表達(dá)基因，使用這種方法需要準(zhǔn)備兩個(gè)矩陣：“實(shí)驗(yàn)設(shè)計(jì)矩陣(designmatrix)”和“對比矩陣(contrastmatrix)”?！皩?shí)驗(yàn)設(shè)計(jì)矩陣”主要用于每張芯片上用的RNA樣品種類，“對比矩陣”主要用于規(guī)定實(shí)驗(yàn)組和對照組?！皩?shí)驗(yàn)設(shè)計(jì)矩陣”：列表示RNA樣品種類，對于oligo芯片或者使用commonreference的cDNA芯片(簡稱第一類)，列數(shù)與不同RNA樣品數(shù)恰好相同；行數(shù)等于芯片數(shù)。對于不同染料對應(yīng)不同樣品的cDNA芯片(簡稱第二類)，列數(shù)比RNA樣品數(shù)恰好少一，行數(shù)等于芯片數(shù)(如要要衡量染料效應(yīng)(dyeeffect)，則列數(shù)與第一類相同)。對于第一

15、類芯片，使用model.matrix()函數(shù)創(chuàng)建“實(shí)驗(yàn)設(shè)計(jì)矩陣”：比如oligo芯片model.matrix(0+factor(c(1,1,1,2,2,3,3,3)表示三類RNA樣品；對于cDNA芯片，modelMatrix(targets,ref="EGFP")創(chuàng)建以“EGFP”為commonreference的矩陣。對于第二類芯片，需要手動(dòng)創(chuàng)建?！皩Ρ染仃嚒蓖ㄟ^函數(shù)makeContrasts()函數(shù)創(chuàng)建。使用步驟：創(chuàng)建“實(shí)驗(yàn)設(shè)計(jì)矩陣”->lmFit()->創(chuàng)建“對比矩陣”(此步可選，如正交實(shí)驗(yàn))->contrasts.fit()

16、(接上步，可選)->eBayes()->topTable()。文檔詳細(xì)列出了各種各樣的例子。比較有用的例子：>beta7Pq$targetsFileNamesSubjectIDCy3Cy5DateofBloodDrawDateofScan16Hs.195.1.gpr1 b7-b7+2002.10.112003.07.252 6Hs.168.gpr3 b7+b7-2003.01.162003.08.0736Hs.166.gpr4 b7+b7-2003.01.162003.08.0746Hs.187.1.gpr6b7-b7+2002.09.162003.0

17、7.1856Hs.194.gpr8b7-b7+2002.09.182003.07.2566Hs.243.1.gpr11b7+b7-2003.01.132003.08.06Labels1 6Hs.195.1.gpr2 6Hs.168.gpr3 6Hs.166.gpr4 6Hs.187.1.gpr5 6Hs.194.gpr6 6Hs.243.1.gpr#假定這六張芯片每兩兩做技術(shù)重復(fù)>design<-c(1,-1,-1,1,1,-1)>corfit<-duplicateCorrelation(beta7pq,design,ndups=1,bl

18、ock=c(1,1,2,2,3,3),weight=NULL)>fit<-lmFit(beta7pq,design,block=c(1,1,2,2,3,3),cor=corfit$consensus)>fit<-eBayes(fit)>topTable(fit,adjust="dfr")P.Valueadj.P.ValB66474.870266e-070.011291225.341680110251.507433e-060.017474174.66378449109.223463e-060.0470632

19、03.434271114701.054914e-050.047063203.33651878321.182200e-050.047063203.252921225821.217992e-050.047063203.230931208121.939031e-050.056625242.88277217552.042581e-050.056625242.843204214052.743542e-050.056625242.616544117172.833754e-050.056625242.591458#也可以采用以下方法>design<-modelMatrix(bet

20、a7pq$targets,ref="b7-")>design<-cbind(Dye=1,desigh)>fit<-lmFit(beta7pq,design)>colnames(design)2<-"b7">cont.matrix<-makeContrasts(MUvsWT=b7/2,levels=design)>fit2<-contrasts.fit(fit,cont.matrix)>fit2&

21、lt;-eBayes(fit2)>topTable(fit2,adjust="fdr")P.Valueadj.P.ValB66478.351055e-070.019361094.430939110253.555781e-060.029768033.700994114313.851970e-060.029768033.65663262116.431916e-060.037279393.362305114702.201258e-050.090759002.58614649102.495165e-050.090759002.50170017553.124884e-0

22、50.090759002.34764578323.131780e-050.090759002.346120119875.008082e-050.127644862.014907218595.505731e-050.127644861.946501#兩者有些不同，原因可能是第一種方法是專門為“對稱”的cDNA芯片設(shè)計(jì)，而第二種是為非對稱設(shè)計(jì)。4.2 類型2的芯片參考文獻(xiàn)8.4和8.5節(jié)。其中用到了model.matrix()這個(gè)函數(shù)，什么意思？4.3 正交實(shí)驗(yàn)參考文獻(xiàn)8.6和8.7節(jié)。4.4 將兩個(gè)通道分開分析參考文獻(xiàn)8.8節(jié)5.做圖5.1 背景和前景imageplot(z,layout,low

23、=NULL,high=NULL,ncolors=123,zerocenter=NULL,zlim=NULL,mar=c(2,1,1,1),legend=TRUE,.)z:可以定義R/Rb/G/Gb,也可以是導(dǎo)數(shù)；low/high:規(guī)定顏色，比如low="white",high="red"下圖是一張芯片的綠色前景圖、紅色背景圖和log-ratio圖(M值圖)：>imageplot(log2(beta7$G,1),beta7$printer,low="white",high="green")&gt

24、;imageplot(log2(beta7$Rb,1),beta7$printer,low="white",high="red")>imageplot(beta7q$M,1,beta7$printer)=也可以使用imageplot3by2(RG,z="Gb",prefix=paste("image",z,sep="-"),path=NULL,zlim=NULL,common.lim=TRUE,)，將圖以3*2的格式六個(gè)一組保存成圖片。5.2 M-A圖使用plotMA(MA,arr

25、ay=1,xlab="A",ylab="M",main=colnames(MA)array,xlim=NULL,ylim=NULL,status,values,pch,col,cex,legend=TRUE,zero.weights=FALSE,)畫單個(gè)MA圖。但這個(gè)函數(shù)有些問題，有時(shí)畫不由。所以,完全可以自己來畫，比如：# plottheMAplotbefornormalization>M<-10g2(beta7$R,1-beta7$Rb,1)/(beta7$G,1-beta7$Gb,1)>A<-(log2(beta7$R,1-beta7$Rb,1)+log2(beta7$G,1-beta7$Gb,1)/2>plot(A,M,cex=0.5)# wecana