1、下载安装R相关包
> source("http://bioconductor.org/biocLite.R") # 由于我之前已经安装了，这里提示升级
# 网上找了下解决方法
> remove.packages("BiocInstaller")
> install.packages("BiocInstaller") # 更新安装包
# 安装、更新bioconductor核心安装包
> biocLite()
# 安装 GEO 包
> biocLite("GEOquery")
# 出现如下问题
installation path not writeable, unable to update packages: cluster, foreign, MASS, Matrix, nlme, survival
# 解决方法：windows以管理员身份打开，linux建议的解决方法为使用管理者权限启动R，即 sudo R
GEOquery 是 NCBI 存储标准化的转录组数据的基因表达综合数据库 GEO 的接口程序。
2、下载芯片数据
教程中我们使用 Dr Andrew Browning 发表的数据集 GSE20986。HUVECs（人脐静脉内皮细胞）是从人胎儿脐带血中提取出来的，通常用来研究内皮细胞的病理生理机制。这个实验设计中，从捐献者的眼组织中提取的虹膜、视网膜、脉络膜微血管内皮细胞用来和 HUVEC 细胞进行比对，以便考查 HUVEC 细胞能否替代其他细胞作为研究眼科疾病的样本。 GEO 页面同时也包含了关于本实验研究的其他信息。对于每组样本有三次测量，样品分成四组 iris, retina, HUVEC, choroidal 。
实验平台是 GPL570 -- GEO 数据库对人类转录组芯片 Affymetrix Human Genome U133 Plus 2.0 Array 的缩写。通过 GSM 链接 GSM524662，我们可以得到各个芯片的更多实验条件信息。同一个实验设计中的不同芯片的实验条件应该是相同的。不同细胞系的提取条件，细胞生长条件，RNA 提取程序，RNA 处理程序，RNA 与探针杂交条件，用哪种仪器扫描芯片，这些数据都以符合 MIAME 标准的形式存储着。
对于每一个芯片，数据表中存储着探针组和对应探针组标准化之后的基因表达量值。表头中提供了表达量值标准化所用的算法。本教程中使用 GC-RMA 算法进行数据标准化，关于 GC-RMA 算法的详细细节可以参考这篇文献。
2.1 下载原始数据
首先，我们从 GEO 数据库下载原始数据，导入 GEOquery 包，用它下载原始数据，下载的原始数据约 53MB。
# 先设置我们的工作路径
> dir.create('D:/TCGA/microarray_analysis',recursive = T)
> setwd('D:/TCGA/microarray_analysis')
> getwd()
[1] "D:/TCGA/microarray_analysis"
# 下载数据
> library(GEOquery)
> getGEOSuppFiles("GSE20986")
# 查看下载好的文件
> list.files('GSE20986')
[1] "GSE20986_RAW.tar" # 可以看到下载了一个打包文件
> setwd('GSE20986')
> getwd()
[1] "D:/TCGA/microarray_analysis/GSE20986"
> untar('GSE20986_RAW.tar',exdir = 'data') # 解包文件到data目录下
> setwd('data')
> getwd()
[1] "D:/TCGA/microarray_analysis/GSE20986/data" # 当前的工作目录
# 查看解包后的文件，可以看到共12个压缩文件，每个样本一个
> list.files()
[1] "GSM524662.CEL.gz" "GSM524663.CEL.gz" "GSM524664.CEL.gz" "GSM524665.CEL.gz" "GSM524666.CEL.gz" "GSM524667.CEL.gz"
[7] "GSM524668.CEL.gz" "GSM524669.CEL.gz" "GSM524670.CEL.gz" "GSM524671.CEL.gz" "GSM524672.CEL.gz" "GSM524673.CEL.gz"
# 将这些文件进行解压缩
> cels <- list.files() ; sapply(cels, gunzip)
GSM524662.CEL.gz GSM524663.CEL.gz GSM524664.CEL.gz GSM524665.CEL.gz GSM524666.CEL.gz GSM524667.CEL.gz GSM524668.CEL.gz
13555726 13555055 13555639 13560122 13555663 13557614 13556090
GSM524669.CEL.gz GSM524670.CEL.gz GSM524671.CEL.gz GSM524672.CEL.gz GSM524673.CEL.gz
13560054 13555971 13554926 13555042 13555290
# 查看解压后的文件
> list.files()
[1] "GSM524662.CEL" "GSM524663.CEL" "GSM524664.CEL" "GSM524665.CEL" "GSM524666.CEL" "GSM524667.CEL" "GSM524668.CEL"
[8] "GSM524669.CEL" "GSM524670.CEL" "GSM524671.CEL" "GSM524672.CEL" "GSM524673.CEL"
Affymetrix的探针（proble）一般是长为25碱基的寡聚核苷酸；探针总是以perfect match 和mismatch成对出现，其信号值称为PM和MM，成对的perfect match 和mismatch有一个共同的affyID。
CEL文件：信号值和定位信息。
CDF文件：探针对在芯片上的定位信息。
当然，如果没有特殊要求，也可直接使用GEOquery包中的getGEO方法直接下载处理好的表达矩阵：
> gse20986 <- getGEO("GSE20986", destdir=".") # 会下载处理好的表达矩阵
3、芯片实验信息整理
在对数据进行分析之前，我们需要先整理好实验设计信息。这其实就是一个文本文件，包含芯片名字、此芯片上杂交的样本名字。为了方便在 R 中 使用 simpleaffy 包读取实验信息文本文件，需要先编辑好格式：
# windows下Rstudio里面自带的terminal具有部分linux的终端功能，可以借用一下
sxuan@DESKTOP-GNA2EAF /d/TCGA/microarray_analysis/GSE20986/data
$ ls -1 *.CEL > finame
sxuan@DESKTOP-GNA2EAF /d/TCGA/microarray_analysis/GSE20986/data
$ awk 'BEGIN{OFS="\t";print "Name","Filename","Target"}{print $1,$1}' finame > phenodata.txt
sxuan@DESKTOP-GNA2EAF /d/TCGA/microarray_analysis/GSE20986/data
$ cat phenodata.txt
Name Filename Target
GSM524662.CEL GSM524662.CEL
GSM524663.CEL GSM524663.CEL
GSM524664.CEL GSM524664.CEL
GSM524665.CEL GSM524665.CEL
GSM524666.CEL GSM524666.CEL
GSM524667.CEL GSM524667.CEL
GSM524668.CEL GSM524668.CEL
GSM524669.CEL GSM524669.CEL
GSM524670.CEL GSM524670.CEL
GSM524671.CEL GSM524671.CEL
GSM524672.CEL GSM524672.CEL
GSM524673.CEL GSM524673.CEL
最终的实验信息文本需要包含三列数据（用 tab 分隔），分别是 Name, FileName, Target。本教程中 Name 和 FileName 这两栏是相同的，当然 Name 这一栏可以用更加容易理解的名字代替。
Target 这一栏数据是芯片上的样本标签，例如 iris, retina, HUVEC, choroidal，这些标签定义了那些数据是生物学重复以便后面的分析。具体的信息可以在GEO数据库样本信息中找到，如下图。
image.png
最后的实验信息文本文件如下：
sxuan@DESKTOP-GNA2EAF /d/TCGA/microarray_analysis/GSE20986/data
$ cat phenodata.txt
Name Filename Target
GSM524662.CEL GSM524662.CEL iris
GSM524663.CEL GSM524663.CEL retina
GSM524664.CEL GSM524664.CEL retina
GSM524665.CEL GSM524665.CEL iris
GSM524666.CEL GSM524666.CEL retina
GSM524667.CEL GSM524667.CEL iris
GSM524668.CEL GSM524668.CEL choroid
GSM524669.CEL GSM524669.CEL choroid
GSM524670.CEL GSM524670.CEL choroid
GSM524671.CEL GSM524671.CEL huvec
GSM524672.CEL GSM524672.CEL huvec
GSM524673.CEL GSM524673.CEL huvec
4、载入数据并对其进行标准化
需要先安装 simpleaffy 包，simpleaffy 包提供了处理 CEL 数据的程序，可以对 CEL 数据进行标准化同时导入实验信息（即前一步中整理好的实验信息文本文件内容），导入数据到 R 变量 celfiles 中：
# 中间可能会有报错，按照报错提示解决即可
> biocLite("simpleaffy")
> library(simpleaffy)
> celfiles <- read.affy(covdesc='phenodata.txt', path='.')
你可以通过输入 ‘celfiles’ 来确定数据导入成功并添加芯片注释（第一次输入 ‘celfiles’ 的时候会进行注释）：
> celfiles
installing the source package ‘hgu133plus2cdf’
trying URL 'https://bioconductor.org/packages/3.6/data/annotation/src/contrib/hgu133plus2cdf_2.18.0.tar.gz'
Content type 'application/x-gzip' length 4353300 bytes (4.2 MB)
downloaded 4.2 MB
* installing *source* package 'hgu133plus2cdf' ...
** R
** data
** preparing package for lazy loading
Warning: replacing previous import 'AnnotationDbi::tail' by 'utils::tail' when loading 'hgu133plus2cdf'
Warning: replacing previous import 'AnnotationDbi::head' by 'utils::head' when loading 'hgu133plus2cdf'
** help
*** installing help indices
converting help for package 'hgu133plus2cdf'
finding HTML links ... 好了
geometry html
hgu133plus2cdf html
hgu133plus2dim html
** building package indices
** testing if installed package can be loaded
*** arch - i386
Warning: replacing previous import 'AnnotationDbi::head' by 'utils::head' when loading 'hgu133plus2cdf'
Warning: replacing previous import 'AnnotationDbi::tail' by 'utils::tail' when loading 'hgu133plus2cdf'
*** arch - x64
Warning: replacing previous import 'AnnotationDbi::tail' by 'utils::tail' when loading 'hgu133plus2cdf'
Warning: replacing previous import 'AnnotationDbi::head' by 'utils::head' when loading 'hgu133plus2cdf'
* DONE (hgu133plus2cdf)
In R CMD INSTALL
The downloaded source packages are in
‘C:\Users\sxuan\AppData\Local\Temp\RtmpItQIXU\downloaded_packages’
AffyBatch object
size of arrays=1164x1164 features (21 kb)
cdf=HG-U133_Plus_2 (54675 affyids)
number of samples=12
number of genes=54675
annotation=hgu133plus2
notes=
Warning messages:
1: replacing previous import ‘AnnotationDbi::tail’ by ‘utils::tail’ when loading ‘hgu133plus2cdf’
2: replacing previous import ‘AnnotationDbi::head’ by ‘utils::head’ when loading ‘hgu133plus2cdf’
现在我们需要对数据进行标准化，使用 GC-RMA 算法对 GEO 数据库中的数据进行标准化，第一次运行的时候需要下载一些其他的必要文件：
> celfiles.gcrma <- gcrma(celfiles)
Adjusting for optical effect............Done.
Computing affinities[1] "Checking to see if your internet connection works..."
installing the source package ‘hgu133plus2probe’
trying URL 'https://bioconductor.org/packages/3.6/data/annotation/src/contrib/hgu133plus2probe_2.18.0.tar.gz'
... ... # 这里省略了一些下载信息
.Done.
Adjusting for non-specific binding............Done.
Normalizing
Calculating Expression
如果你想看标准化之后的数据，输入 celfiles.gcrma， 你会发现提示已经不是 AffyBatch object 了，而是 ExpressionSet object，是已经标准化了的数据：
> celfiles.gcrma
ExpressionSet (storageMode: lockedEnvironment)
assayData: 54675 features, 12 samples
element names: exprs
protocolData
sampleNames: GSM524662.CEL GSM524663.CEL ... GSM524673.CEL (12 total)
varLabels: ScanDate
varMetadata: labelDescription
phenoData
sampleNames: GSM524662.CEL GSM524663.CEL ... GSM524673.CEL (12 total)
varLabels: sample Filename Target
varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu133plus2
5、芯片数据可视化
再进行下一步的数据分析之前，我们有必要对数据质量进行检查，确保没有其他的问题。首先，可以通过对标准化之前和之后的数据画箱线图来检查 GC-RMA 标准化的效果：
> # 载入色彩包
> library(RColorBrewer)
> # 设置调色板
> cols <- brewer.pal(8, "Set1")
> # 对标准化之前的探针信号强度做箱线图
> boxplot(celfiles, col=cols)
> # 对标准化之后的探针信号强度做箱线图，需要先安装 affyPLM 包，以便解析 celfiles.gcrma 数据
> biocLite("affyPLM")
> library(affyPLM)
> boxplot(celfiles.gcrma, col=cols)
> # 标准化前后的箱线图会有些变化
> # 但是密度曲线图看起来更直观一些
> # 对标准化之前的数据做密度曲线图
> hist(celfiles, col=cols)
> # 对标准化之后的数据做密度曲线图
> hist(celfiles.gcrma, col=cols)
数据标准化之前的箱线图
数据标准化之后的箱线图
数据标准化之前的密度图
数据标准化之后的密度图
通过这些图我们可以看出这12张芯片数据之间差异不大，标准化处理将所有芯片信号强度标准化到具有类似分布特征的区间内。为了更详细地了解芯片探针信号强度，我们可以使用 affyPLM 对单个芯片 CEL 数据进行可视化。
> # 从 CEL 文件读取探针信号强度:
> celfiles.qc <- fitPLM(celfiles)
> # 对芯片 GSM24662.CEL 信号进行可视化:
> image(celfiles.qc, which=1, add.legend=TRUE)
> # 对芯片 GSM524665.CEL 进行可视化
> # 这张芯片数据有些人为误差
> image(celfiles.qc, which=4, add.legend=TRUE)
> # affyPLM 包还可以画箱线图
> # RLE (Relative Log Expression 相对表达量取对数) 图中
> # 所有的值都应该接近于零。 GSM524665.CEL 芯片数据由于有人为误差而例外。
> RLE(celfiles.qc, main="RLE")
> # 也可以用 NUSE (Normalised Unscaled Standard Errors)作图比较.
> # 对于绝大部分基因，标准差的中位数应该是1。
> # 芯片 GSM524665.CEL 在这个图中，同样是一个例外
> NUSE(celfiles.qc, main="NUSE")
单个芯片AffyPLM信号图
存在人为误差的单个芯片AffyPLM信号图
affyPLM包画箱线图
cel数据的NUSE图
我们还可以通过层次聚类来查看样本之间的关系：
> eset <- exprs(celfiles.gcrma)
> distance <- dist(t(eset),method="maximum")
> clusters <- hclust(distance)
> plot(clusters)
CEL数据的层次聚类图
图形显示，与其他眼组织相比 HUVEC 样品是单独的一组，表现出组织类型聚集的一些特征，另外 GSM524665.CEL 数据在此图中并不显示为异常值。
6、数据质量控制
现在我们可以对数据进行分析了，分析的第一步就是要过滤掉数据中的无用数据，例如作为内参的探针数据，基因表达无明显变化的数据（在差异表达统计时也会被过滤掉），信号值与背景信号差不多的探针数据。
下面的 nsFilter 参数是为了不删除没有 Entrez Gene ID 的位点，保留有重复 Entrez Gene ID 的位点：
> celfiles.filtered <- nsFilter(celfiles.gcrma, require.entrez=FALSE, remove.dupEntrez=FALSE)
> # 哪些位点被过滤掉了？为什么？
> celfiles.filtered$filter.log
$numLowVar
[1] 27307
$feature.exclude
[1] 62
我们可以看出有 27307 个探针位点因为无明显表达差异（LowVar）被过滤掉，有 62 个探针位点因为是内参而被过滤掉。
7、查找有表达差异的探针位点
现在有了过滤之后的数据，我们就可以用 limma 包进行差异表达分析了。首先，我们要提取样本的信息：
> samples <- celfiles.gcrma$Target
> # 检查数据的分组信息
> samples
[1] "iris" "retina" "retina" "iris" "retina" "iris" "choroid"
[8] "choroid" "choroid" "huvec" "huvec" "huvec"
> # 将分组数据转换为因子类型变量
> samples <- as.factor(samples)
> # 检查转换的因子变量
> samples
[1] iris retina retina iris retina iris choroid choroid choroid
[10] huvec huvec huvec
Levels: choroid huvec iris retina
> # 设置实验分组
> design <- model.matrix(~0 + samples)
> colnames(design) <- c("choroid", "huvec", "iris", "retina")
> # 检查实验分组
> design
choroid huvec iris retina
1 0 0 1 0
2 0 0 0 1
3 0 0 0 1
4 0 0 1 0
5 0 0 0 1
6 0 0 1 0
7 1 0 0 0
8 1 0 0 0
9 1 0 0 0
10 0 1 0 0
11 0 1 0 0
12 0 1 0 0
attr(,"assign")
[1] 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$samples
[1] "contr.treatment"
现在我们将芯片数据进行了标准化和过滤，也有样品分组和实验分组信息，可以将数据导入 limma 包进行差异表达分析了：
> library(limma)
> # 将线性模型拟合到过滤之后的表达数据集上
> fit <- lmFit(exprs(celfiles.filtered$eset), design)
> # 建立对比矩阵以比较组织和细胞系
> contrast.matrix <- makeContrasts(huvec_choroid = huvec - choroid, huvec_retina = huvec - retina, huvec_iris <- huvec - iris, levels=design)
> # 检查对比矩阵
> contrast.matrix
Contrasts
Levels huvec_choroid huvec_retina huvec_iris
choroid -1 0 0
huvec 1 1 1
iris 0 0 -1
retina 0 -1 0
> # 现在将对比矩阵与线性模型拟合，比较不同细胞系的表达数据
> huvec_fits <- contrasts.fit(fit, contrast.matrix)
> # 使用经验贝叶斯算法计算差异表达基因的显著性
> huvec_ebFit <- eBayes(huvec_fits)
> # 返回对应比对矩阵 top 10 的结果
> # coef=1 是 huvec_choroid 比对矩阵, coef=2 是 huvec_retina 比对矩阵
> topTable(huvec_ebFit, number=10, coef=1)
ID logFC AveExpr t P.Value adj.P.Val
6147 204779_s_at 7.367947 4.171874 72.77016 3.290669e-15 8.985500e-11
7292 207016_s_at 6.934796 4.027229 57.37259 3.712060e-14 5.068076e-10
8741 209631_s_at 5.193313 4.003541 51.22423 1.177337e-13 1.071612e-09
26309 242809_at 6.433514 4.167462 48.52518 2.042404e-13 1.394247e-09
6828 205893_at 4.480463 3.544350 40.59376 1.253534e-12 6.845801e-09
20232 227377_at 3.670688 3.209217 36.03427 4.200854e-12 1.911809e-08
6222 204882_at -5.351976 6.512018 -34.70239 6.154318e-12 2.400711e-08
27051 38149_at -5.051906 6.482418 -31.44098 1.672263e-11 5.282592e-08
6663 205576_at 6.586372 4.139236 31.31584 1.741131e-11 5.282592e-08
6589 205453_at 3.623706 3.210306 30.72793 2.109110e-11 5.759136e-08
B
6147 20.25563
7292 19.44681
8741 18.96282
26309 18.70767
6828 17.75331
20232 17.01960
6222 16.77196
27051 16.08854
6663 16.05990
6589 15.92277
分析结果的各列数据含义：
第一列是探针组在表达矩阵中的行号； # 有的没有此列
第二列“ID” 是探针组的 AffymatrixID；
第三列“logFC”是两组表达值之间以2为底对数化的的变化倍数（Fold change, FC），由于基因表达矩阵本身已经取了对数，这里实际上只是两组基因表达值均值之差；
第四列“AveExpr”是该探针组所在所有样品中的平均表达值；
第五列“t”是贝叶斯调整后的两组表达值间 T 检验中的 t 值；
第六列“P.Value”是贝叶斯检验得到的 P 值；
第七列“adj.P.Val”是调整后的 P 值；
第八列“B”是经验贝叶斯得到的标准差的对数化值。
如果要设置一个倍数变化阈值，并查看不同阈值返回了多少基因，可以使用 topTable 的 lfc 参数，参数设置为 5,4,3,2 时返回的基因个数：
> nrow(topTable(huvec_ebFit,coef=1, number=10000, lfc=5))
[1] 88
> nrow(topTable(huvec_ebFit,coef=1, number=10000, lfc=4))
[1] 194
> nrow(topTable(huvec_ebFit,coef=1, number=10000, lfc=3))
[1] 386
> nrow(topTable(huvec_ebFit,coef=1, number=10000, lfc=2))
[1] 1016
> nrow(topTable(huvec_ebFit,coef=1, number=10000, lfc=1))
[1] 4025
> probeset.list <- topTable(huvec_ebFit, coef=1, number=10000, lfc=4)
> head(probeset.list)
logFC AveExpr t P.Value adj.P.Val B
204779_s_at 7.367790 4.171707 72.77347 3.284937e-15 8.969850e-11 20.25762
207016_s_at 6.936667 4.027733 57.39252 3.694641e-14 5.044293e-10 19.44987
209631_s_at 5.192949 4.003992 51.24892 1.170273e-13 1.065182e-09 18.96660
242809_at 6.433238 4.168870 48.51842 2.043082e-13 1.394710e-09 18.70852
205893_at 4.480331 3.543714 40.56477 1.261400e-12 6.888757e-09 17.75050
204882_at -5.353547 6.513690 -34.70646 6.141288e-12 2.395629e-08 16.77396
9、注释差异分析结果的基因ID
为了将探针集注释上基因 ID （转换为geneID）我们需要先安装一些数据库的包和注释的包，之后可以提取 topTable 中的探针 ID 并注释上基因 ID：
> biocLite("hgu133plus2.db")
> library(hgu133plus2.db)
> library(annotate)
> gene.symbols <- getSYMBOL(rownames(probeset.list), "hgu133plus2")
> results <- cbind(gene.symbols, probeset.list)
> head(results)
gene.symbols logFC AveExpr t P.Value adj.P.Val B
204779_s_at HOXB7 7.367790 4.171707 72.77347 3.284937e-15 8.969850e-11 20.25762
207016_s_at <NA> 6.936667 4.027733 57.39252 3.694641e-14 5.044293e-10 19.44987
209631_s_at <NA> 5.192949 4.003992 51.24892 1.170273e-13 1.065182e-09 18.96660
242809_at IL1RL1 6.433238 4.168870 48.51842 2.043082e-13 1.394710e-09 18.70852
205893_at NLGN1 4.480331 3.543714 40.56477 1.261400e-12 6.888757e-09 17.75050
204882_at ARHGAP25 -5.353547 6.513690 -34.70646 6.141288e-12 2.395629e-08 16.77396
# 可以看到这里有一些芯片ID已经找不到对应的基因ID了，我们将这些失败的提出来先
> trans.failed <- results[is.na(results$gene.symbols),]
> head(trans.failed)
gene.symbols logFC AveExpr t P.Value adj.P.Val B
207016_s_at <NA> 6.936667 4.027733 57.39252 3.694641e-14 5.044293e-10 19.449872
209631_s_at <NA> 5.192949 4.003992 51.24892 1.170273e-13 1.065182e-09 18.966602
243154_at <NA> 4.213968 3.476199 28.24818 4.932212e-11 1.035992e-07 15.294507
209994_s_at <NA> -5.949092 7.766907 -19.86274 1.695876e-09 1.403261e-06 12.349999
1556314_a_at <NA> -6.200352 7.680901 -18.81114 2.916279e-09 2.041844e-06 11.859483
240466_at <NA> 4.671497 4.946474 14.78163 3.156679e-08 1.134161e-05 9.610741
# 看下有多少行失败的
> nrow(trans.failed)
[1] 16
下面我们使用另一种方法来进行注释探针ID，虽然也是使用annotate包，这里虽然探针都注释到的比上面多，但有的一个探针注释到了多个ID:
# 下载芯片平台对应的数据包
> biocLite("hgu133plus2.db")
> library("annotate")
> probes <- as.character(rownames(probeset.list))
> out <- select(hgu133plus2.db, probes,"SYMBOL")
'select()' returned 1:many mapping between keys and columns # 表示有的探针注释到了多个基因ID
> head(out)
PROBEID SYMBOL
1 204779_s_at HOXB7
2 207016_s_at ALDH1A2
3 207016_s_at LOC101928635
4 209631_s_at GPR37
5 209631_s_at SEL1L2
6 242809_at IL1RL1
> probeset.list$PROBEID <- rownames(probeset.list)
> results.final <- merge(out, probeset.list, by="PROBEID", sort=F)
> head(results.final)
PROBEID SYMBOL logFC AveExpr t P.Value adj.P.Val B
1 204779_s_at HOXB7 7.367790 4.171707 72.77347 3.284937e-15 8.969850e-11 20.25762
2 207016_s_at ALDH1A2 6.936667 4.027733 57.39252 3.694641e-14 5.044293e-10 19.44987
3 207016_s_at LOC101928635 6.936667 4.027733 57.39252 3.694641e-14 5.044293e-10 19.44987
4 209631_s_at GPR37 5.192949 4.003992 51.24892 1.170273e-13 1.065182e-09 18.96660
5 209631_s_at SEL1L2 5.192949 4.003992 51.24892 1.170273e-13 1.065182e-09 18.96660
6 242809_at IL1RL1 6.433238 4.168870 48.51842 2.043082e-13 1.394710e-09 18.70852
> nrow(probeset.list);nrow(results.final)
[1] 194
[1] 222
查看一下两种方法注释到结果的区别：
> dupID <- results.final[duplicated(results.final$PROBEID),]; nrow(dupID)
[1] 28
> b<- results.final[is.na(results.final$SYMBOL),]; b
PROBEID SYMBOL logFC AveExpr t P.Value adj.P.Val B
13 243154_at <NA> 4.213968 3.476199 28.24818 4.932212e-11 1.035992e-07 15.294507
32 1556314_a_at <NA> -6.200352 7.680901 -18.81114 2.916279e-09 2.041844e-06 11.859483
49 240466_at <NA> 4.671497 4.946474 14.78163 3.156679e-08 1.134161e-05 9.610741
# 从上面可以看到仍然是有三个探针没有注释到
# 将第二种方法一对多注释的和未注释到的合并起来
> a<-as.factor(dupID$PROBEID); f <- as.factor(dupID$PROBEID)
> c <- c(f, as.character(b$PROBEID))
> d<-data.frame(c);d
c
1 204670_x_at
2 207016_s_at
3 208306_x_at
4 209312_x_at
5 209631_s_at
6 209994_s_at
7 211796_s_at
8 215193_x_at
9 218805_at
10 221974_at
11 226558_at
12 229748_x_at
13 64064_at
14 243154_at
15 1556314_a_at
16 240466_at
比较后发现这两种方法出现问题的16个探针是相同的，因此直接使用第一种方法getSYMBOL()更加简单。
最后，将注释好的差异表达基因写入文件进行后续分析。
# 将为注释到的探针移除
> nrow(results);results<-na.omit(results);nrow(results)
[1] 194
[1] 178
> write.table(results, file="resluts.txt", sep="\t",row.names = F)
参考：
http://bioinformatics.knowledgeblog.org/2011/06/20/analysing-microarray-data-in-bioconductor/
Bioconductor 分析芯片数据教程
