基因芯片（Affymetrix）分析4：GO和KEGG分析
基因列表的分析一般都会涉及GO和KEGG分析，Bioconductor提供了很多这方面的R工具包。
选择工作目录，读入上一次分析和保存的数据：
results.sig <- read.csv("results.lim.7d.csv", header=TRUE, as.is=TRUE) head(results.sig)
## X logFC AveExpr t P.Value adj.P.Val B ## 1 254818_at 6.215 10.363 41.38 7.304e-10 1.169e-05 11.538 ## 2 254805_at 6.844 7.280 30.81 6.095e-09 4.878e-05 10.431 ## 3 245998_at 2.778 10.011 25.44 2.411e-08 9.545e-05 9.528 ## 4 265119_at 4.380 8.282 24.07 3.588e-08 9.545e-05 9.240 ## 5 256114_at 4.461 7.668 23.92 3.745e-08 9.545e-05 9.208 ## 6 265722_at -2.913 9.276 -23.91 3.760e-08 9.545e-05 9.205
colnames(results.sig)[1] <- "probe_id" ; genes.sig <- results.sig[, 1]
1 获取AGI、GO和KEGG注释
library(ath1121501.db) ls("package:ath1121501.db")
## [1] "ath1121501" "ath1121501ACCNUM" ## [3] "ath1121501ARACYC" "ath1121501ARACYCENZYME" ## [5] "ath1121501CHR" "ath1121501CHRLENGTHS" ## [7] "ath1121501CHRLOC" "ath1121501CHRLOCEND" ## [9] "ath1121501.db" "ath1121501_dbconn" ## [11] "ath1121501_dbfile" "ath1121501_dbInfo" ## [13] "ath1121501_dbschema" "ath1121501ENZYME" ## [15] "ath1121501ENZYME2PROBE" "ath1121501GENENAME" ## [17] "ath1121501GO" "ath1121501GO2ALLPROBES" ## [19] "ath1121501GO2PROBE" "ath1121501MAPCOUNTS" ## [21] "ath1121501ORGANISM" "ath1121501ORGPKG" ## [23] "ath1121501PATH" "ath1121501PATH2PROBE" ## [25] "ath1121501PMID" "ath1121501PMID2PROBE" ## [27] "ath1121501SYMBOL"
ath1121501GO为拟南芥基因的GO数据库，ath1121501PATH为KEGG pathway数据库。但不是每一个基因（probeset）都有GO或KEGG注释，哪些基因有注释可以用mappedkeys函数获得：
length(mappedkeys(ath1121501PATH))
## [1] 3018
length(mappedkeys(ath1121501GO))
## [1] 20299
head(mappedkeys(ath1121501PATH))
## [1] "261579_at" "261569_at" "261583_at" "261574_at" "261043_at" "261044_at"
有PATH注释的probesets只有3018个，而有GO注释的有2万多个。
通过ath1121501XXXX获得的数据是AnnotationDbi软件包定义的ProbeAnnDbBimap类型数据，它们可以用as.list转成列表形式。列表内每一个基因的注释内容也是列表形式：
all.path <- ath1121501PATH[mappedkeys(ath1121501PATH)] class(all.path)
## [1] "ProbeAnnDbBimap" ## attr(,"package") ## [1] "AnnotationDbi"
as.list(all.path)[1]
## $`261579_at` ## [1] "00190"
转换成列表类型的ProbeAnnDbBimap数据仍然是列表，但PATH和ACCNUM数据是二级列表（列表下只有一级列表），而GO数据是三级列表（列表下还有两级的列表）。所以得先编写get.GO函数，它把as.list产生的GO三级列表转成二级结构，和AGI和KEGG的列表类似，方便后面的统一处理：
get.GO <- function(the.keys, goList){ results <- NULL for (i in 1:length(the.keys)){ n <- length(goList[[i]]) info <- NULL for(j in 1:n){info <- c(info, goList[[i]][[j]]$GOID)} info <- list(info) names(info) <- the.keys[i] results <- c(results, info) } results }
使用这个函数和下列代码就可以获得AGI、GO和KEGG注释：
library(plyr) results.anno <- results.sig[,1:2] for(i in 1:3){ anno <- switch(i, ath1121501ACCNUM, ath1121501GO, ath1121501PATH) anno.label <- switch(i, "AGI", "GO", "PATH") mapped.probes <- mappedkeys(anno) mapped.present <- intersect(genes.sig, mapped.probes) mapped.anno <- as.list(anno[mapped.present]) if(anno.label=="GO") mapped.anno <- get.GO(mapped.present, mapped.anno) mapped.anno <- llply(mapped.anno, unique) mapped.anno <- ldply(mapped.anno, paste, collapse="; ") colnames(mapped.anno) <- c("probe_id", anno.label) results.anno <- merge(results.anno, mapped.anno, by.x = "probe_id", all = TRUE) }
上面代码有两点要注意：
- switch()函数使用。switch()是非常神奇的条件转向开关函数，它的参数（列表）可以是各种类型，变量、表达式、函数等都可以使用。
- 列表到数据框类型数据的转换，我们使用了plyr软件包的llply和ldply函数。plyr是很著名的软件包，用于数据糅合。这不属于本节的讨论范围，先不介绍，请自行学习使用。
由于探针id是唯一的，上面的代码用它作为关键字糅合数据。得到的结果是数据框：
str(results.anno)
## 'data.frame': 740 obs. of 5 variables: ## $ probe_id: chr "245015_at" "245042_at" "245088_at" "245196_at" ... ## $ logFC : num 1.37 -1.13 -1.62 -1.71 1.44 ... ## $ AGI : chr "ATCG00490" "AT2G26540" "AT2G39850" "AT1G67750" ... ## $ GO : chr "GO:0006091; GO:0006354; GO:0009737; GO:0015977; GO:0015979; GO:0018119; GO:0046686; GO:0016020; GO:0005618; GO:0009507; GO:0009"| __truncated__ "GO:0006779; GO:0006780; GO:0033014; GO:0009507; GO:0004852" "GO:0008152; GO:0006508; GO:0043086; GO:0005576; GO:0009505; GO:0004252; GO:0042802" "GO:0008150; GO:0005576; GO:0030570" ... ## $ PATH : chr NA NA NA "00040" ...
这样每一个探针都得到了对应的AGI、GO和KEGG途径注释（如果有）。其他类型数据如Pubmed ID可以使用类似方法获得，但编程之前得先了解它们的数据结构，最直接的方法就是使用head，summary和str等函数查看。
得到的结果用write.csv或其他存盘函数保存。
2 GO和KEGG富集分析
Bioconductor中有不少软件包可以进行GO和KEGG统计分析和作图，如GOstats和KEGGgraph，但我不建议使用它们。
对于这类分析，我推荐使用另外一个不是R软件包的网络分析软件：Cytoscape。它是免费的开源软件，由多所大学和几个公司联合开发和维护，已经逐渐成为网络分析的标准工具，软件网址为：http://www.cytoscape.org/
R软件包RCytoscape可以把R的分析结果推送到Cytoscape，充分利用R的统计功能和Cytoscape的可视化能力。但是RCytoscape不一定能跟上Cytoscape的软件更新步伐，所以最好还是把R分析的结果保存成文件，再用Cytoscape直接分析。
事实上，如果使用Cytoscape进行GO和KEGG富集或网络分析，我们只需要获得AGI列表，而不需要获得GO和KEGG注释。
3 Session Info
sessionInfo()
## R version 3.1.0 (2014-04-10) ## Platform: x86_64-pc-linux-gnu (64-bit) ## ## locale: ## [1] LC_CTYPE=zh_CN.UTF-8 LC_NUMERIC=C ## [3] LC_TIME=zh_CN.UTF-8 LC_COLLATE=zh_CN.UTF-8 ## [5] LC_MONETARY=zh_CN.UTF-8 LC_MESSAGES=zh_CN.UTF-8 ## [7] LC_PAPER=zh_CN.UTF-8 LC_NAME=C ## [9] LC_ADDRESS=C LC_TELEPHONE=C ## [11] LC_MEASUREMENT=zh_CN.UTF-8 LC_IDENTIFICATION=C ## ## attached base packages: ## [1] parallel tcltk stats graphics grDevices utils datasets ## [8] methods base ## ## other attached packages: ## [1] plyr_1.8.1 ath1121501.db_2.14.0 org.At.tair.db_2.14.0 ## [4] RSQLite_0.11.4 DBI_0.2-7 AnnotationDbi_1.27.3 ## [7] GenomeInfoDb_1.1.2 Biobase_2.25.0 BiocGenerics_0.11.0 ## [10] zblog_0.1.0 knitr_1.5 ## ## loaded via a namespace (and not attached): ## [1] evaluate_0.5.3 formatR_0.10 highr_0.3 IRanges_1.99.2 ## [5] Rcpp_0.11.1 S4Vectors_0.0.2 stats4_3.1.0 stringr_0.6.2 ## [9] tools_3.1.0
更多精彩内容
