搜索
查看: 1870|回复: 0

[annotation] GO.db 这个包里面记录着GO数据库的各种信息

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-3-1 17:07:56 | 显示全部楼层 |阅读模式
把下面代码运行一遍就明白了,我能做的就这些了
[AppleScript] 纯文本查看 复制代码
###################################
## Basic usage of GO information ##
###################################
library(GOstats); library(GO.db); library(ath1121501.db); library(annotate) # Loads the required libraries.
goann <- as.list(GOTERM) # Retrieves full set of GO annotations.
zz <- eapply(GOTERM, function(x) x@Ontology); table(unlist(zz)) # Calculates the number of annotations for each ontology category.
?GOTERM # To find out, how to access the different GO components.
GOTERM$"GO:0003700"; GOMFPARENTS$"GO:0003700"; GOMFCHILDREN$"GO:0003700" 
   # Shows how to print out the GO annotations for one entry and how to retrieve its direct parents and children.
GOMFANCESTOR$"GO:0003700"; GOMFOFFSPRING$"GO:0003700" # Prints out complete lineages of parents and children for a GO ID.
goterms <- unlist(eapply(GOTERM, function(x) x@Term)); goterms[grep("molecular_function", goterms)] 
   # Retrieves all GO terms and prints out only those matching a search string given in the grep function. The same can 
   # be done for the definition field with 'x@Definition'. A set of GO IDs can be provided as well: goterms[GOMFANCESTOR$"GO:0005507"]
go_df <- data.frame(GOID=unlist(eapply(GOTERM, function(x) x@GOID)), Term=unlist(eapply(GOTERM, function(x) x@Term)), Ont=unlist(eapply(GOTERM, function(x) x@Ontology))) 
   # Generates data frame of the commonly used GO components: GOID, GO Term and Ontology Type.
affyGO <- eapply(ath1121501GO, getOntology, "MF"); table(sapply(affyGO, length)) 
   # Retrieves MF GO terms for all probe IDs of a chosen Affy chip and calculates how many probes have multiple GO terms 
   # associated. Use "BP" and "CC" arguments to retrieve BP/CC GO terms.
affyGOdf <- data.frame(unlist(affyGO)); affyGOdf <- data.frame(AffyID=row.names(affyGOdf), GOID=affyGOdf[,1]); affyGOdf <- merge(affyGOdf, go_df, by.x="GOID", by.y="GOID", all.x=T) 
   # Converts above MF list object into a data frame. The AffyID occurence counts are appended to AffyIDs. The last step 
   # merges the two data frames: 'affyGOdf' and 'go_df'.
unique(lookUp("GO:0004713", "ath1121501", "GO2ALLPROBES")) # Retrieves all Affy IDs that are associated with a GO node.
z <- affyGO[c("254759_at", "260744_at")]; as.list(GOTERM)[z[[1]]] 
   # Retrieves GO IDs for set of Affy IDs and then the corresponding GO term for first Affy ID.
a <- data.frame(unlist(z)); a <- data.frame(ID=row.names(a), a); b <- data.frame(goterms[as.vector(unlist(z))]); b <- data.frame(ID=row.names(b), b); merge(b, a, by.x = "ID", by.y="unlist.z.") 
   # Merges Affy ID, GO ID and GO annotation information.
affyEv <- eapply(ath1121501GO, getEvidence); table(unlist(affyEv, use.names = FALSE)) 
   # Provides evidence code information for each gene and summarizes the result.
test1 <- eapply(ath1121501GO, dropECode, c("IEA", "NR")); table(unlist(sapply(test1, getEvidence), use.names = FALSE)) 
   # This example shows how one can remove certain evidence codes (e.g. IEA, IEP) from the analysis.

##############################################
## GO term enrichment analysis with GOstats ##
##############################################
## Example of how to test a sample set of probe set keys for over-representation of GO terms using a hypergeometric distribution
## test with the function hyperGTest(). For more information, read the GOstatsHyperG manual.
library(ath1121501.db); library(ath1121501cdf)
affySample <- c("266592_at", "266703_at", "266199_at", "246949_at", "267370_at", "267115_s_at", "266489_at", "259845_at", "266295_at", "262632_at")
geneSample <- as.vector(unlist(mget(affySample, ath1121501ACCNUM, ifnotfound=NA)))
affyUniverse <- ls(ath1121501cdf)
geneUniverse <- as.vector(unlist(mget(affyUniverse, ath1121501ACCNUM, ifnotfound=NA)))
params <- new("GOHyperGParams", geneIds = geneSample, universeGeneIds = geneUniverse, annotation="ath1121501", ontology = "MF", pvalueCutoff = 0.5, conditional = FALSE, testDirection = "over")
hgOver <- hyperGTest(params)
summary(hgOver)
htmlReport(hgOver, file = "MyhyperGresult.html") 

请仔细观察!
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-16 09:12 , Processed in 0.026806 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.