如何使用GEOquery和limma完成芯片数据的差异表达分析

差异表达分析是一种常用的方法，用于比较不同条件下基因表达水平的变化。GEOquery和limma是R语言中广泛使用的两个包，可用于处理和分析芯片数据的差异表达。下面是使用GEOquery和limma进行差异表达分析的步骤：

1. 下载和导入芯片数据

使用GEOquery包中的getGEO函数下载芯片数据并导入到R中。例如，如果您要下载GSE12345数据集，可以使用以下代码：

library(GEOquery)
gset <- getGEO("GSE12345")

2. 数据质量控制

在进行差异表达分析之前，需要对数据进行质量控制。使用GEOquery包中的summary函数和plotPCA函数可以对芯片数据进行基本的质量控制。例如，可以使用以下代码绘制PCA图：

library(limma)
library(GEOquery)
gset <- getGEO("GSE12345")
edata <- exprs(gset[[1]])
edata <- t(edata)
edata <- na.omit(edata)
fit <- prcomp(edata)
plotPCA(fit)

3. 数据预处理

对芯片数据进行归一化和标准化处理，以消除不同芯片之间的差异，并确保数据符合正态分布。常用的预处理方法包括RMA、GCRMA、MAS5等。使用limma包中的normalizeBetweenArrays函数可以对芯片数据进行预处理。例如，可以使用以下代码对芯片数据进行RMA预处理：

library(limma)
library(GEOquery)
gset <- getGEO("GSE12345")
edata <- exprs(gset[[1]])
edata <- t(edata)
edata <- na.omit(edata)
edata <- backgroundCorrect.RMA(edata)
edata <- normalize.quantiles.RMA(edata)
edata <- log2(edata)

4. 差异表达分析

使用limma包中的lmFit函数和eBayes函数可以进行差异表达分析。lmFit函数用于拟合线性模型，eBayes函数用于对差异表达结果进行统计显著性检验。例如，可以使用以下代码进行差异表达分析：

library(limma)
library(GEOquery)
gset <- getGEO("GSE12345")
edata <- exprs(gset[[1]])
edata <- t(edata)
edata <- na.omit(edata)
edata <- backgroundCorrect.RMA(edata)
edata <- normalize.quantiles.RMA(edata)
edata <- log2(edata)
factors <- c(0,0,0,1,1,1)
design <- model.matrix(~factors)
fit <- lmFit(edata, design)
fit <- eBayes(fit)
results <- topTable(fit, adjust="BH", sort.by="P", n=1000)

上述代码中，factors表示芯片数据中不同样本的分组信息，design表示设计矩阵，fit表示拟合的线性模型，results表示差异表达结果。

5. 结果分析

根据差异表达结果，可以进行进一步的功能分析、通路分析等。常用的工具包括ClusterProfiler、GOstats、KEGGprofile等。例如，可以使用以下代码进行GO分析：

library(ClusterProfiler)
gset <- getGEO("GSE12345")
edata <- exprs(gset[[1]])
edata <- t(edata)
edata <- na.omit(edata)
edata <- backgroundCorrect.RMA(edata)
edata <- normalize.quantiles.RMA(edata)
edata <- log2(edata)
factors <- c(0,0,0,1,1,1)
design <- model.matrix(~factors)
fit <- lmFit(edata, design)
fit <- eBayes(fit)
results <- topTable(fit, adjust="BH", sort.by="P", n=1000)

genes <- rownames(results)
geneList <- names(genes)
geneList <- names(genes)[abs(results$logFC) > 1 & results$adj.P.Val < 0.05]
ego <- enrichGO(geneList, OrgDb="org.Hs.eg.db", ont="BP")
barplot(ego)

上述代码中，geneList表示差异表达基因列表，ego表示进行GO分析所得到的富集结果。