使用limma包分析GEO芯片数据

limma 包是Bioconductor中分析芯片数据的一个标准包。它具有分析芯片的一整套流程（数据读入，背景矫正，芯片内归一化，芯片间归一化，差异基因删选），并且能够分析单色和双色芯片，在bioconductor中是一个比较权威的芯片分析包。似乎最近也可以用来分析基因组数据，大家可以去尝试一下。

本文结合自己分析GSE21009芯片的经验，大致概览了limma的分析流程。这篇文章为方便只分析了第三四个芯片（M2-B）。这也是我第一次使用limma包进行分析，如有不足，敬请大家批评指正！

对于GEO上的数据，R中有一个GEOquery包能够方便的获取NCBI上的芯片数据，同时也能够实现代码的自动化。不足的是，使用过程中出现了一些数据类型的混乱，下面会有所提及。

代码如下（包含分析所用的包）：

library(GEOquery)
library(Biobase)
library(limma)        #主要用这个包分析
library(plyr)

m4.p <- getGEO('GSE21009', GSEMatrix=FALSE, destdir="./origin2/")  #如此获得是Matrix类的表达数据，默认GSEMatrix=TRUE  

如果有原始数据，getGEO()命令会自动下载。并保存在destdir中。

如果想查看下载的这个大类中各个数据的信息，可使用以下的命令：

                                        # 基本信息的获取
## 查看包含的GSM类
names(GSMList(m4.p))
## 查看包含的GPL
names(GPLList(m4.p))
## 查看GSM对应的平台
lapply(GSMList(m4.p),function(x) {Meta(x)$platform})
## 获取第一个GSM的类型及信息
class(GSMList(m4.p)[[1]])
Meta(GSMList(m4.p)[[1]])        # 获取芯片信息
## 查看GSM的列信息
Columns(GSMList(m4.p)[[1]])
## 获取ID信息
Table(GPLList(m4.p)[[1]])$ID        # 获取GPL，忽略平台信息等，直接生成一个data.frame结构

以上是针对非Matrix类的，对于获取的Matrix类可使用head()，names()，varLabels()，featureNames()，sampleNames()，phenoData及pData获取其基本信息，请使用?+函数名查询。

本文从GEO获取的数据是经过一些处理的，似乎不能直接用limma的函数直接读取，因此第一步是先生成能够被limma分析的数据结构————RGList。

这是关于RGList类的介绍：

## RGList-class

### Slots/List Components:

###

### ‘RGList’ objects can be created by ‘new("RGList",RG)’ where ‘RG’

### is a list. Objects of this class contains no slots (other than

### ‘.Data’), but objects should contain the following list

### components:

### ‘R’ numeric matrix containing the red (cy5) foreground intensities. Rows correspond to spots and columns to arrays.

### ‘G’ numeric matrix containing the green (cy3) foreground intensities. Rows correspond to spots and columns to arrays.

### Optional components include

### ‘Rb’ numeric matrix containing the red (cy5) background intensities ### ‘Gb’ numeric matrix containing the green (cy3) background intensities ### ‘weights’ numeric matrix of same dimension as ‘R’ containing relative spot quality weights. Elements should be non-negative. ### ‘other’ list containing other matrices, all of the same dimensions as ‘R’ and ‘G’. ### ‘genes’ data.frame containing probe information. Should have one row for each spot. May have any number of columns. ### ‘targets’ data.frame containing information on the target RNA samples. Rows correspond to arrays. May have any number of columns. Usually includes columns ‘Cy3’ and ‘Cy5’ specifying which RNA was hybridized to each array.

### ‘printer’ list containing information on the process used to print the spots on the arrays. See PrintLayout.

### Valid ‘RGList’ objects may contain other optional components, but

### all probe or array information should be contained in the above

### components.

以下代码是将我们从GEO获取的原始数据转换成limma可以直接使用的RGList类：

### 读取数据
gs.3 <- Table(GSMList(m4.p)[[3]])
gs.4 <- Table(GSMList(m4.p)[[4]])
### 计算RGList R G Rb Gb数值
#### 原数据的修饰，方便引用
library(plyr)
colnames(gs.3) <- revalue(colnames(gs.3),c("F635 Median"="R","F532 Median"="G","F635 Median - B635"="R.b","F532 Median - B532"="G.b"))        #更改名字，以方便引用
colnames(gs.4) <- revalue(colnames(gs.4),c("F635 Median"="R","F532 Median"="G","F635 Median - B635"="R.b","F532 Median - B532"="G.b"))
#### 红绿值
or <- data.frame(myb.a=gs.3$R,myb.b=gs.4$R)        # Red
og <- data.frame(myb.a=gs.3$G,myb.b=gs.4$G)        # green
colnames(or) <- c("myb.a","myb.b")        # colnames经过上步木有预期改变，似乎要调整
colnames(og) <- c("myb.a","myb.b")        # colnames经过上步木有预期改变，似乎要调整
#### 背景值
orb <- data.frame(myb.a=gs.3$R-gs.3$R.b,myb.b=gs.4$R-gs.4$R.b)        # red background
ogb <- data.frame(myb.a=gs.3$G-gs.3$G.b,myb.b=gs.4$G-gs.4$G.b)        # green background
colnames(orb) <- c("myb.a","myb.b")
colnames(ogb) <- c("myb.a","myb.b")
#### 对GPL的修饰，主要是数据类型的改变（使用GEOquery，不知道为啥读出来的类型完全是错的）
gp.b <- Table(GPLList(m4.p)[[2]])
gp.b$ID <- as.integer(gp.b$ID)
gp.b$Block <- as.integer(gp.b$Block)
gp.b$Column <- as.integer(gp.b$Column)
gp.b$Row <- as.integer(gp.b$Row)
gp.b$Name <- as.character(gp.b$Name)
gp.b$GB_LIST <- as.character(gp.b$GB_LIST)
gp.b$SEQUENCE <- as.character(gp.b$SEQUENCE)
gp.b$LocusID <- as.character(gp.b$LocusID)
gp.b$Annotation <- as.character(gp.b$Annotation)
gp.b$SPOT_ID <- as.character(gp.b$SPOT_ID)
#### RGList中的genes值，似乎可以全部贴上去？
genes <- gp.b[,c(2,4,3,1,5)]
#### 合并成RGList类
mRG <- list(R=as.matrix(or),G=as.matrix(og),Rb=as.matrix(orb),Gb=as.matrix(ogb),genes=genes)
RG <- new("RGList",mRG)

而对于原始数据，可以直接使用readTargets()，read.maimages(), readGAL()命令生成RGList类，十分的方便。不过这样的函数只能处理图片软件读取出来的数值型的数据，不能直接读图（似乎R中就没有可以直接读图的包。。。）。

读入单通道数据，可以设定green.only = TRUE即可，然后对应读入columns = list(G = "Col1", Gb = "Col2")（这是抄袭的，自己没有实践过）。读入的数据，如果是单通道，则成为EListRaw class；如果是双通道，则是RGList class。

接下来就进入重点，使用limma包分析我们的RGList类。

limma包可以查看芯片的信号情况：

### 获得输出
RG$printer <- getLayout(RG$genes)
### 绘制芯片图
pdf(file="./result/B-M2/p.image.pdf",width=6,height=6)
imageplot(log2(RG$R[,1]), RG$printer, low="white", high="red")
dev.off()

使用pdf()仅仅只是将图保存下来。

后续分析，最重要的是获取gene的Locus号，这个比较简单，我使用的是merge()命令³。

diffgene.B.M2 <- merge(diffgene.B.M2,gp.b)