limma包的使用技巧

limmar package是一个功能比较全的包，既含有cDNA芯片的RAW data输入、前处理（归一化）功能，同时也有差异化基因分析的“线性”算法（limma: Linear Models for Microarray Data），特别是对于“多因素实验（multifactor designed experiment）”。limmar包的可扩展性非常强，单通道（one channel）或者双通道（tow channel）数据都可以分析差异基因，甚至也包括了定量PCR和RNA-seq（第一次见分析microarray的包也能处理RNA-seq）。

limmar package是一个集大成的包，对载入数据、数据前处理（背景矫正、组内归一化和组间归一化都有很多种选择方法）、差异基因分析都有很多的选择。而且，所设计的线性回归和经验贝叶斯方法找差异基因非常值得学习。

1. 读入样本信息

使用函数读readTargets(file="Targets.txt", path=NULL, sep=" ", row.names=NULL, quote=""",...)。这个函数其实是一个包装了的read.table()，读入的是样本的信息，创建的对象类似于marray包的marrayInfos和Biobase包的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=" ", quote=NULL, ...)

参数说明：files需要通过函数dir(pattern = "Mypattern")配合正则表达式筛选（规范命名很重要），同时该函数可以读入符合格式的压缩过的文件，比如*.txt.gz的文件，这极大的减小的数据储存大小；source的取值分为两类，一类是“高富帅”，比如“agilent”、“spot”等等（下表），它们是特定扫描仪器的特定输出格式；如果不幸是“屌丝”，即格式是自己规定的，可以选定source="generic"，这时需要规定columns；任何cDNA文件都要有R/G/Rb/Gb四列（Mean或者Median）；annotation可以规定哪些是注释列；wt.fun用于对点样点进行质量评估，取值为0表示这些点将在后续的分析中被剔除，取值位1表示需要保留，对点样点的评估依赖于图像扫描软件的程序设定，比如SPOT和GenePix软件，查看QualityWeights（现成函数或者自己写函数）。

读入单通道数据：读入单通道数据，可以设定green.only = TRUE即可，然后对应读入columns = list(G = "Col1", Gb = "Col2")。

读入的数据，如果是单通道，则成为EListRaw class；如果是双通道，则是RGList class。

数据操作：

cbind()：合并数据；

“[”：分割数据；

RGList class有的names是 "R"，"G"，"Rb"，"Gb"，"weights"，"printer"，"genes"，"targets"，"notes"：  R/G/Rb/Gb分别红和绿的前景和背景噪音；weight是扫描软件的质量评估；printer是点样规则（printer layout）；genes是基因注释；target是样本注释；notes是一般注释。可以通过myRGList$names进行相应的取值和赋值。

3. 前处理

cDNA芯片前处理的流程是：背景去除（background correction）--> 组内归一化（with-array normalization）--> 组间归一化（between-array normalization）。最后一步组间归一化可选，注意trade-off原则。

3.1 背景去除： 

backgroundCorrect(RG, method="auto", offset=0, printer=RG$printer, normexp.method="saddle", verbose=TRUE)

RG：可以是EListRaw，也可以是RGList class；

method：取“auto”，“none”，“subtract”，"half"，"minimum"，"movingmin"，"edwards"或者"normexp"；

normexp.method：在method取“normexp”时，可以取"saddle"，"mle"，"rma"或者"rma75"。

作者建议使用“mle”或者“saddle”，两个运行速度也差不多。如果数据中含有“形态背景估计（morphological background estimates）”，比如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=0)

method："none"，"median"，"loess"，"printtiploess"，"composite"，"control"和"robustspline"。初始值设定为“printtiploess”，但是对于Agilent芯片或者小样本芯片（每个“点样块”少于150个探针）。可用选用的loess或者robustspline。“loess归一化方法假设了有相当大的一部分探针没有发生差异化表达，而不是上调或者下调的基因数差不多或者差异化程度围绕着0波动。”（参考文献1）。需要注意，组内归一化只是对单张芯片的M值进行了规整（不影响A值），但是对与组间各个通道没有进行比较。

weight：是图像处理软件对探针权重的标记，如果使用weight，那么在归一化过程中weight为0的探针不会影响其他探针，这并不意味着这些探针会被剔除，它们同样会被归一化，也会出现在归一化的结果中。如果不想使用，那么weight设为NULL。

iterations：设定loess的循环数，循环数越多，结果越强健（robust）。

3.3 组间归一化

normalizeBetweenArrays(object, method=NULL, targets=NULL, cyclic.method="fast", ...)

method：“none"，"scale"，"quantile"，"Aquantile"，"Gquantile"，“Rquantile”，“Tquantile”，“cyclicloess”。“scale”对log-ratio数值进行归一化；“quantile”对密度（intensity）进行归一化；“Aquantile”对A数值进行归一化，不调正M数值；"Gquantile"和“Rquantile”则分别对Green和Red通道进行归一化，适用于Green或者Red为“参考标准（common reference）”的样品；“Tquantile”表示样品分开归一化或者Green和Red一先一后进行归一化。

以下是一个cDNA芯片的密度图，分别为原始、去背景（method = "normexp", normexp.method="mle"）、组内归一化（method = "robustspline"）和组间归一化（method = “quantile”），使用plotDensities()绘制。

4. 线性模型设计（linear model desigh）

“线性模型”主要分析差异化表达基因，使用这种方法需要准备两个矩阵：“实验设计矩阵（design matrix）”和“对比矩阵（contrast matrix）”。“实验设计矩阵”主要用于每张芯片上用的RNA样品种类，“对比矩阵”主要用于规定实验组和对照组。

“实验设计矩阵”：列表示RNA样品种类，对于oligo芯片或者使用common reference的cDNA芯片（简称第一类），列数与不同RNA样品数恰好相同；行数等于芯片数。对于不同染料对应不同样品的cDNA芯片（简称第二类），列数比RNA样品数恰好少一，行数等于芯片数（如要要衡量染料效应（dye effect），则列数与第一类相同）。对于第一类芯片，使用model.matrix()函数创建“实验设计矩阵”：比如oligo芯片model.matrix(~ 0+factor(c(1,1,1,2,2,3,3,3)))表示三类RNA样品；对于cDNA芯片，modelMatrix(targets,ref="EGFP")创建以“EGFP”为common reference的矩阵。对于第二类芯片，需要手动创建。

“对比矩阵”通过函数makeContrasts()函数创建。

使用步骤：创建“实验设计矩阵” --> lmFit() --> 创建“对比矩阵”（此步可选，如正交实验） -->contrasts.fit()（接上步，可选）--> eBayes() --> topTable()。文档详细列出了各种各样的例子。

 比较有用的例子：

4.1 交换染料的技术重复

===========================================================================

> beta7pq$targets

      FileNames SubjectID  Cy3  Cy5 Date of Blood Draw Date of Scan

1 6Hs.195.1.gpr         1 b7 - b7 +         2002.10.11   2003.07.25

2   6Hs.168.gpr         3 b7 + b7 -         2003.01.16   2003.08.07

3   6Hs.166.gpr         4 b7 + b7 -         2003.01.16   2003.08.07

4 6Hs.187.1.gpr         6 b7 - b7 +         2002.09.16   2003.07.18

5   6Hs.194.gpr         8 b7 - b7 +         2002.09.18   2003.07.25

6 6Hs.243.1.gpr        11 b7 + b7 -         2003.01.13   2003.08.06

         Labels

1 6Hs.195.1.gpr

2   6Hs.168.gpr

3   6Hs.166.gpr

4 6Hs.187.1.gpr

5   6Hs.194.gpr

6 6Hs.243.1.gpr

# 假定这六张芯片每两两做技术重复

> design <- c(1, -1, -1, 1, 1, -1)

> corfit <- duplicateCorrelation(beta7pq, design, ndups=1, block = 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.Value  adj.P.Val        B

6647  4.870266e-07 0.01129122 5.341680

11025 1.507433e-06 0.01747417 4.663784

4910  9.223463e-06 0.04706320 3.434271

11470 1.054914e-05 0.04706320 3.336518

7832  1.182200e-05 0.04706320 3.252921

22582 1.217992e-05 0.04706320 3.230931

20812 1.939031e-05 0.05662524 2.882772

1755  2.042581e-05 0.05662524 2.843204

21405 2.743542e-05 0.05662524 2.616544

11717 2.833754e-05 0.05662524 2.591458

# 也可以采用以下方法

> design <- modelMatrix(beta7pq$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 <- eBayes(fit2)

> topTable(fit2, adjust = "fdr")

        P.Value  adj.P.Val        B

6647  8.351055e-07 0.01936109 4.430939

11025 3.555781e-06 0.02976803 3.700994

11431 3.851970e-06 0.02976803 3.656632

6211  6.431916e-06 0.03727939 3.362305

11470 2.201258e-05 0.09075900 2.586146

4910  2.495165e-05 0.09075900 2.501700

1755  3.124884e-05 0.09075900 2.347645

7832  3.131780e-05 0.09075900 2.346120

11987 5.008082e-05 0.12764486 2.014907

21859 5.505731e-05 0.12764486 1.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 = 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")

> 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, array=1, xlab="A", ylab="M", main=colnames(MA)[array], xlim=NULL, ylim=NULL, status, values, pch, col, cex, legend=TRUE, zero.weights=FALSE, ...)画单个MA图。但这个函数有些问题，有时画不出。所以，完全可以自己来画，比如：

==========================================================

# plot the MAplot befor normalization

> M <- log2((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)

# we can also plot the MAplot above using marray package

> beta7ma <- as(beta7, "marrayRaw")

> plot(maA(beta7ma)[,1], maM(beta7ma)[, 1])

# plot the MAplot after normalization

> plot(beta7q$A[, 1], beta7q$M[, 1], cex = 0.5)

===========================================================

使用plotPrintTipLoess(object,layout,array=1,span=0.4,...)画print-tip的MA图，下图分别是归一化前和后的print-tip的MA图。

使用boxplot()画盒箱图，使用plotDensity()画密度曲线图。下图为原始数据、组内归一化（method = "normexp", normexp.method="mle"）和组建归一化（method = "scale"）的盒箱图。

最后，还可以使用 volcanoplot(fit, coef=1, highlight=0, names=fit$genes$ID, xlab="Log Fold Change", ylab="Log Odds", pch=16, cex=0.35, ...)画“火山图”（highlight标记前N个差异基因）；

使用vennDiagram(object, include="both", names, mar=rep(1,4), cex=1.5, lwd=1, circle.col, counts.col, show.include, ...)画“韦恩图”；

使用plotMDS()绘制多纬标度图（multidimensional scaling plot, MDS）。

参考文献：limma: Linear Models for Microarray Data User’s Guide