这次以R语言处理表达芯片教程中的GSE号为例子，具体可以参看http://www.biotrainee.com/thread-566-1-1.html 首先我们先进GSE42872网站看一下这个芯片数据是什么平台的，如：Affymetrix Human Gene 1.0 ST Array [transcript (gene) version]，这说明这个芯片是Affymetrix的一个human表达谱芯片，用来检测基因表达情况。
接下来就之前文章回顾基因表达芯片分析（数据处理）中同样的方法先对数据进行处理
用getGEO下载表达矩阵，其实这里不下载也没事，因为后续数据分析还是会从raw data开始，再到标准化后的表达矩阵。用slotNames等函数主要还是为了看下gset[[1]]这个对象的一些信息
gset <- getGEO("GSE42872",getGPL = FALSE)
class(gset[[1]])
gset[[1]]
slotNames(gset[[1]])
annotation(gset[[1]])
phenoData(gset[[1]])
featureData(gset[[1]])
assayData(gset[[1]])然后获取表达矩阵和比较感兴趣的pdata
head(featureNames(gset[[1]]))
pdata <- pData(gset[[1]])
exprSet_ori <- exprs(gset[[1]])
dim(exprSet_ori)下面开始则是比较重要的步骤，从芯片的raw data数据开始到最后的差异表达分析
数据标准化
由于网速不太好，我还是选择人工下载再放置在data文件中，然后用read.celfiles读入
file_CELs <- list.celfiles("data", listGzipped = TRUE, full.name = TRUE)
rawAffy <- read.celfiles(filenames = file_CELs, phenoData = phenoData(gset[[1]]), sampleNames = rownames(pData(gset[[1]])))标准化这步跟之前的略有不同，区别在于我们要获取转录本的表达量，而不是probeset；所以在使用rma函数时略有不同（其实rma默认参数就是基于转录本的）
eset <- rma(rawAffy, target = "core")
show(eset)
dim(eset)
featureData(eset) <- getNetAffx(eset, "transcript")Select one probeset per gene
由于通常多个探针会对应一个基因，所以我们需要在后续分析中选其中最有代表性的一个探针的表达量来代表基因的表达水平
如Differential gene expression analysis I (microarray data)教程中所介绍是使用IQR值最大的方法，也就是选择在所有样本中表达量变化最大的那个探针。还有一种方法（平时听得较多）是选用表达量最大的那个探针。不管哪种标准，其处理方法其实都是类似的，比如以IQR最大的方法为例，先计算每个探针的IQR值
iqrs <- apply(exprs(eset), 1, IQR)然后使用genefilter包的findLargest函数，其实这步自己来写R代码来实现也可以的，嫌麻烦就用这个函数，并且这时还需要一个芯片对应的注释R包。比如这个芯片是Affymetrix Human Gene 1.0 ST Array，那么需要下载的注释包则是hugene10sttranscriptcluster.db；如果是小鼠的 Affymetrix Mouse Gene 1.0 ST Array，则对应的注释包为mogene10sttranscriptcluster.db
library(hugene10sttranscriptcluster.db)
library(genefilter)
prbs <- findLargest(featureNames(eset), testStat = iqrs, data = "hugene10sttranscriptcluster.db")
eset2 <- eset[prbs, ]
dim(eset2)然后获取表达矩阵
exprSet <- exprs(eset)如果想选择取探针中表达量最大的方法，只需要将iqrs向量变成每个探针在所有样本中表达量的平均值（或者总和）即可
Probe ID转化为symbol
这步是分析芯片数据中必不可少的，比如我们可以看下exprSet表达矩阵的探针ID
head(exprSet, 2)
##GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620
##8144866 5.429810 5.628468 5.376701 6.054502 6.106508 5.975622
##8066431 9.543157 9.406544 9.500876 10.150883 10.251234 10.006384这个8144866数字就代表着一个探针，这对后续下游分析是无法识别的，所以要将其转化为更有识别性的ID，比如gene symbol，这里也需要用到芯片对应的注释R包
我们可以看下hugene10sttranscriptcluster.db这个包有哪些调用数据的函数
columns(hugene10sttranscriptcluster.db)
##[1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
##[6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
##[11] "GO" "GOALL" "IPI" "MAP" "OMIM"
##[16] "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID"
##[21] "PROBEID" "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG"
##[26] "UNIGENE" "UNIPROT" 使用SYMBOL从注释包中提取出探针ID与symbol对应关系的数据框
probe2symbol <- toTable(hugene10sttranscriptclusterSYMBOL)然后将exprSet中的探针ID转化为symbol
exprSet <- as.data.frame(exprSet)
exprSet$probe_id <- rownames(exprSet)
exprSet_symbol <- merge(probe2symbol, exprSet, by = "probe_id")
dim(exprSet_symbol)再把exprSet中symbol列转化为列名，并只保留各个样本表达量的列（去除掉其他不必要的列）
rownames(exprSet_symbol) <- exprSet_symbol$symbol
exprSet_symbol <- exprSet_symbol[, -c(1,2)]差异表达分析
芯片数据的差异表达分析一般首选的就是limma包。当然limma现在不仅可以处理芯片数据，而且可以处理RNA-seq数据。如本文的例子，是两组样本单因素实验设计，代码就比较简单了
先设置分组信息（分组信息可以看pdata）
library(limma)
condition <- factor(c(rep("control", 3), rep("vemurafenib", 3)), levels = c("control", "vemurafenib"))
design <- model.matrix(~condition)然后常规的差异分析，这里默认是使用limma-trend方法（如果是RNA-seq数据还可以使用voom方法）
fit <- lmFit(log2(exprSet_symbol), design)
fit=eBayes(fit)
output <- topTable(fit, coef=2,n=Inf)最后统计下差异基因数目
sum(output$adj.P.Val<0.01)
##[1] 5274GO/KEGG富集分析以及PPI
后续的GO/KEGG分析的方法跟RNA-seq数据，具体可参照之前的RNA-seq分析的博文转录组差异表达分析小实战（二）
PPI分析可以使用在线的STRING网站，也可以使用stringDB包，具体可以参照之前的一篇博文STRINGDB包的简单使用
参考资料：
Differential gene expression analysis I (microarray data)
http://www.bio-info-trainee.com/1586.html
http://www.biotrainee.com/thread-566-1-1.html
本文出自于http://www.bioinfo-scrounger.com转载请注明出处