limma包的使用技巧

1、limma包的使用技巧原文地址：limma包的使用技巧作者：牛牛龙limmarpackage是一个功能比较全的包，既含有cDNA芯片的RAWdata输入、前处理(归一化)功能，同时也有差异化基因分析的“线性”算法(limma:LinearModelsforMicroarrayData)，特别是对于“多因素实验(multifactordesignedexperiment)”。limmar包的可扩展性非常强，单通道(onechannel)或者双通道(towchannel)数据都可以分析差异基因，甚至也包括了定量PCR和RNA-seq(第一次见分析microarray的包也能处理RNA-seq)。l

2、immarpackage是一个集大成的包，对载入数据、数据前处理(背景矫正、组内归一化和组间归一化都有很多种选择方法)、差异基因分析都有很多的选择。而且，所设计的线性回归和经验贝叶斯方法找差异基因非常值得学习。1.读入样本信息使用函数读readTargets(file="Targets.txt",path=NULL,sep="t",s=NULL,quote=""",.)。这个函数其实是一个包装了的read.table()，读入的是样本的信息，创建的对象类似于marray包的marrayInfos和Biobas

3、e包的AnnotatedDataFrame。2.读入探针密度数据与marray包一致，Bioconductor不能读入原始的TIFF图像文件，只能输出扫描仪输入的、转换成数字信号的文本文件。使用函数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、.）参数说明：files需要通过函数dir（pattern="Mypattern"）配合正则表达式筛选（规范命名很重要），同时该函数可以读入符合格式的压缩过的文件，比如*.txt.gz的文件，这极大的减小的数据储存大小；source的取值分为两类，一类是“高富帅”，比如“agilent"、“spot”等等（下表），它们是特定扫描仪器的特定输由格式；如果不幸是“展丝”，即格式是自己规定的，可以选定source="generic"，这时需要规定columns；任何cDNA文件都要有R/G/Rb/Gb四列（Mean或者Median）；annotati

5、on可以规定哪些是注释列；wt.fun用于对点样点进行质量评估，取值为0表示这些点将在后续的分析中被剔除，取值位1表示需要保留，对点样点的评估依赖于图像扫描软件的程序设定，比如SPOT和GenePix软件，查看QualityWeights(现成函数或者自己写函数)。读入单通道数据：读入单通道数据，可以设定green.only=TRUE即可，然后对应读入columns=list(G="Col1",Gb="Col2")。读入的数据，如果是单通道，则成为EListRawclass；如果是双通道，则是RGListclass。数据操作：cbind()：合并数据；“

6、”：分割数据；RGListclass有的names是"R","G","Rb","Gb"weights”printer","genes"，"targets"，"notes"：R/G/Rb/Gb分别红和绿的前景和背景噪音；weight是扫描软件的质量评估；printer是点样规则(printerlayout)；genes是基因注释；target是样本注释；notes是一般注释。可以通过myRGList$names进行相应的取值和赋值。3. 前处理cD

7、NA芯片前处理的流程是：背景去除(backgroundcorrection)->组内归一化(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”时，可以取"saddle"，"mle"，"rma"或者"rma75"。作者建议使用“mle”或者“saddle"，两个运行速度也差不多。如果数据中含有“形态背景

9、估计(morphologicalbackgroundestimates)”，比如SPOT和GenePix图像处理软件，那么method="subtract"也就较好的表现。3.2组内归一化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"。初始值设定为“printtiploess”，但是对于Agilent芯片或者小样本芯片（每个“点样块”少于150个探针）。可用选用的loess或者robustspline。"loess归一化方法假设了有相当大的一部分探针没有发生差异化表达，而不是上调或者下调的基因数差不多或者差异

11、化程度围绕着0波动。”（参考文献1）。需要注意，组内归一化只是对单张芯片的M值进行了规整（不影响A值），但是对与组间各个通道没有进行比较。weight：是图像处理软件对探针权重的标记，如果使用weight，那么在归一化过程中weight为0的探针不会影响其他探针，这并不意味着这些探针会被剔除，它们同样会被归一化，也会出现在归一化的结果中。如果不想使用，那么weight设为NULL。iterations：设定loess的循环数，循环数越多，结果越强健（robust）。3.3 组间归一化normalizeBetweenArrays(object,method=NULL,targets=NULL,c

12、yclic.method="fast",.)method：“none"，"scale"，"quantile"，"Aquantile"，"Gquantile"，“Rquantile"，“Tquantile"，“cyclicloess"。“scale”对log-ratio数值进行归一化；“quantile”对密度(intensity)进行归一化；“Aquantile”对A数值进行归一化，不调正M数值；"Gquantile"和“Rquanti

13、le”则分别对Green和Red通道进行归一化，适用于Green或者Red为“参考标准(commonreference)”的样品；“Tquantile”表示样品分开归一化或者Green和Red一先一后进行归一化。以下是一个cDNA芯片的密度图，分别为原始、去背景(method="normexp",normexp.method="mle")、组内归一化(method="robustspline")和组间归一化(method=“quantile”)，使用plotDensities()绘制。4.线性模型设计(linearmodeldesig

14、h)“线性模型”主要分析差异化表达基因，使用这种方法需要准备两个矩阵：“实验设计矩阵(designmatrix)”和“对比矩阵(contrastmatrix)”。“实验设计矩阵”主要用于每张芯片上用的RNA样品种类，“对比矩阵”主要用于规定实验组和对照组。“实验设计矩阵”：列表示RNA样品种类，对于oligo芯片或者使用commonreference的cDNA芯片(简称第一类)，列数与不同RNA样品数恰好相同；行数等于芯片数。对于不同染料对应不同样品的cDNA芯片(简称第二类)，列数比RNA样品数恰好少一，行数等于芯片数(如要要衡量染料效应(dyeeffect)，则列数与第一类相同)。对于第一

15、类芯片，使用model.matrix()函数创建“实验设计矩阵”：比如oligo芯片model.matrix(0+factor(c(1,1,1,2,2,3,3,3)表示三类RNA样品；对于cDNA芯片，modelMatrix(targets,ref="EGFP")创建以“EGFP”为commonreference的矩阵。对于第二类芯片，需要手动创建。“对比矩阵”通过函数makeContrasts()函数创建。使用步骤：创建“实验设计矩阵”->lmFit()->创建“对比矩阵”(此步可选，如正交实验)->contrasts.fit()

16、(接上步，可选)->eBayes()->topTable()。文档详细列出了各种各样的例子。比较有用的例子：>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#假定这六张芯片每两两做技术重复>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芯片设计，而第二种是为非对称设计。4.2 类型2的芯片参考文献8.4和8.5节。其中用到了model.matrix()这个函数，什么意思？4.3 正交实验参考文献8.6和8.7节。4.4 将两个通道分开分析参考文献8.8节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,也可以是导数；low/high:规定颜色，比如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的格式六个一组保存成图片。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,)画单个MA图。但这个函数有些问题，有时画不由。所以,完全可以自己来画，比如：# 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