本文补充些芯片数据分析注意事项。
GEOquery 与 limma
前面文章的分析案例都是从原始数据开始分析，如果想采用上传作者处理好的表达矩阵开始分析，一般用 GEOquery 下载数据并用 limma 进行差异分析。
因为 GEOquery 返回 ExpressionSet 对象，在《ExpressionSet 数据结构》文章已详细介绍 ExpressionSet 相关操作；limma 差异分析代码也在前 2 篇文章有。这里不再赘述代码操作，只讲一讲注意事项。
注意 getGEO 返回对象的类型，使用 GEO 参数时返回为列表，此时需要取列表的元素再操作。
> geo1 <- getGEO(GEO = "GSE45636", destdir = "GSE45636", getGPL = FALSE)
> class(geo1)
[1] "list"
# 列表，而不是 ExpressionSet
> exprs(geo1)
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘exprs’ for signature ‘"list"’
# 取列表第一个，是需要的 ExpressionSet
> exprs(geo1[[1]])[1:3, 1:3]
GSM1111159 GSM1111160 GSM1111161
1007_s_at 7.450834 6.962312 7.216552
1053_at 6.788506 6.122643 6.597648
117_at 3.557593 3.657426 3.361082
如果自己下载好文件使用 filename 参数，将返回 ExpressionSet 对象。
> geo2 <- getGEO(filename = "GSE45636/GSE45636_series_matrix.txt.gz", destdir = "GSE45636", getGPL = FALSE)
> class(geo2)
[1] "ExpressionSet"
attr(,"package")
[1] "Biobase"
如果使用作者处理好的表达数据，要注意是不是 log2 后的，因为 limma 默认输入数据为 log2 后的。这个没有特别的方法，打开数据集的某个样品详情页，查看数据处理过程。譬如里面提到了用 RMA 算法，那就基本可以确定是 log2 后的了。也可以看看数据大小，譬如有上千的数值，基本可以确定没有 log2 处理。
芯片注释
如果数据集的芯片平台在网站 Bioconductor - 3.13 AnnotationData Packages 有注释 R 包，推荐优先使用。注释 R 包使用方法见《AnnotationDbi 教程》文章。
GEOquery 也能下载平台注释信息，但往往由于网络问题不成功。
> gpl1 <- getGEO(GEO = "GPL570")
> class(gpl1)
[1] "GPL"
attr(,"package")
[1] "GEOquery"
> Table(gpl1)[1:3, 1:5]
ID GB_ACC SPOT_ID Species Scientific Name Annotation Date
1 1007_s_at U48705 NA Homo sapiens Oct 6, 2014
2 1053_at M87338 NA Homo sapiens Oct 6, 2014
3 117_at X51757 NA Homo sapiens Oct 6, 2014
另外，手动从 GEO 下载注释文件并导入 R 也行，避免了 GEOquery 的网络问题。
另外，大部分时候都会遇到多个探针对应一个基因的情况，而我们的分析往往需要基因层面。这时候可以选择将多个探针合并，取平均值或中位数或其他。也可以根据分析目的，不进行合并，比如说并不需要基因的表达只需要基因差异分析，可以直接对探针差异分析，然后看对应了哪些基因。
基因 ID 转换
基因 ID 转换推荐使用基因组注释 R 包，或者 biomaRt 包在线转换，具体操作见《批量转换基因ID》文章。另外，芯片数据参杂着各种平台的基因 ID 才是最噩梦的地方，只有耐心手动整理才行。
