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或其他存檔函式儲存。