###################################################
### chunk number 1: load_library
###################################################
#line 39 "vignettes/goseq/inst/doc/goseq.Rnw"
library(goseq)


###################################################
### chunk number 2: set_width
###################################################
#line 42 "vignettes/goseq/inst/doc/goseq.Rnw"
options(width=84)


###################################################
### chunk number 3: read.data eval=FALSE
###################################################
## #line 61 "vignettes/goseq/inst/doc/goseq.Rnw"
## gene.vector=as.integer(assayed.genes%in%de.genes)
## names(gene.vector)=assayed.genes
## head(gene.vector)


###################################################
### chunk number 4: supported_genomes eval=FALSE
###################################################
## #line 82 "vignettes/goseq/inst/doc/goseq.Rnw"
## supportedGenomes()


###################################################
### chunk number 5: supported_geneids eval=FALSE
###################################################
## #line 86 "vignettes/goseq/inst/doc/goseq.Rnw"
## supportedGeneIDs()


###################################################
### chunk number 6: getLengthDataFromUCSC eval=FALSE
###################################################
## #line 116 "vignettes/goseq/inst/doc/goseq.Rnw"
## txsByGene=transcriptsBy(txdb,"gene")
## lengthData=median(width(txsByGene))


###################################################
### chunk number 7: edger_1
###################################################
#line 140 "vignettes/goseq/inst/doc/goseq.Rnw"
library(edgeR)
table.summary=read.table(system.file("extdata","Li_sum.txt",package='goseq'),sep='\t',header=TRUE,stringsAsFactors=FALSE)
counts=table.summary[,-1]
rownames(counts)=table.summary[,1]
grp=factor(rep(c("Control","Treated"),times=c(4,3)))
summarized=DGEList(counts,lib.size=colSums(counts),group=grp)


###################################################
### chunk number 8: edger_2
###################################################
#line 151 "vignettes/goseq/inst/doc/goseq.Rnw"
disp=estimateCommonDisp(summarized)
disp$common.dispersion
tested=exactTest(disp)
topTags(tested)


###################################################
### chunk number 9: edger_3
###################################################
#line 160 "vignettes/goseq/inst/doc/goseq.Rnw"
genes=as.integer(p.adjust(tested$table$p.value[tested$table$logFC!=0],method="BH")<.05)
names(genes)=row.names(tested$table[tested$table$logFC!=0,])
table(genes)


###################################################
### chunk number 10: head_genomes
###################################################
#line 170 "vignettes/goseq/inst/doc/goseq.Rnw"
head(supportedGenomes())


###################################################
### chunk number 11: head_geneids
###################################################
#line 176 "vignettes/goseq/inst/doc/goseq.Rnw"
head(supportedGeneIDs(),n=12)


###################################################
### chunk number 12: pwf
###################################################
#line 188 "vignettes/goseq/inst/doc/goseq.Rnw"
pwf=nullp(genes,"hg18","ensGene")
head(pwf)


###################################################
### chunk number 13: GO.wall
###################################################
#line 204 "vignettes/goseq/inst/doc/goseq.Rnw"
GO.wall=goseq(pwf,"hg18","ensGene")
head(GO.wall)


###################################################
### chunk number 14: GO.samp
###################################################
#line 215 "vignettes/goseq/inst/doc/goseq.Rnw"
GO.samp=goseq(pwf,"hg18","ensGene",method="Sampling",repcnt=1000)


###################################################
### chunk number 15: head_samp
###################################################
#line 218 "vignettes/goseq/inst/doc/goseq.Rnw"
head(GO.samp)


###################################################
### chunk number 16: plot_wal_v_samp
###################################################
#line 224 "vignettes/goseq/inst/doc/goseq.Rnw"
plot(log10(GO.wall[,2]), log10(GO.samp[match(GO.samp[,1],GO.wall[,1]),2]), xlab="log10(Wallenius p-values)", ylab="log10(Sampling p-values)", xlim=c(-3,0))
abline(0,1,col=3,lty=2)


###################################################
### chunk number 17: GO.nobias
###################################################
#line 235 "vignettes/goseq/inst/doc/goseq.Rnw"
GO.nobias=goseq(pwf,"hg18","ensGene",method="Hypergeometric")
head(GO.nobias)


###################################################
### chunk number 18: plot_wal_v_hyper
###################################################
#line 242 "vignettes/goseq/inst/doc/goseq.Rnw"
plot(log10(GO.wall[,2]), log10(GO.nobias[match(GO.nobias[,1],GO.wall[,1]),2]), xlab="log10(Wallenius p-values)", ylab="log10(Hypergeometric p-values)", xlim=c(-3,0), ylim=c(-3,0))
abline(0,1,col=3,lty=2)


###################################################
### chunk number 19: GO.limited
###################################################
#line 251 "vignettes/goseq/inst/doc/goseq.Rnw"
GO.MF=goseq(pwf,"hg18","ensGene",test.cats=c("GO:MF"))
head(GO.MF)


###################################################
### chunk number 20: enriched_GO
###################################################
#line 262 "vignettes/goseq/inst/doc/goseq.Rnw"
enriched.GO=GO.wall$category[p.adjust(GO.wall$over_represented_pvalue,method="BH")<.05]
head(enriched.GO)


###################################################
### chunk number 21: GO_explained
###################################################
#line 269 "vignettes/goseq/inst/doc/goseq.Rnw"
library(GO.db)
for(go in enriched.GO[1:10]){
	print(GOTERM[[go]])
	cat("--------------------------------------\n")
}


###################################################
### chunk number 22: getlength
###################################################
#line 282 "vignettes/goseq/inst/doc/goseq.Rnw"
len=getlength(names(genes),"hg18","ensGene")
length(len)
length(genes)
head(len)


###################################################
### chunk number 23: getgo
###################################################
#line 291 "vignettes/goseq/inst/doc/goseq.Rnw"
go=getgo(names(genes),"hg18","ensGene")
length(go)
length(genes)
head(go)


###################################################
### chunk number 24: conv_table
###################################################
#line 300 "vignettes/goseq/inst/doc/goseq.Rnw"
goseq:::.ID_MAP
goseq:::.ORG_PACKAGES


###################################################
### chunk number 25: norm_analysis eval=FALSE
###################################################
## #line 307 "vignettes/goseq/inst/doc/goseq.Rnw"
## pwf=nullp(genes,"hg18","ensGene")
## go=goseq(pwf,"hg18","ensGene")


###################################################
### chunk number 26: verbose_analysis eval=FALSE
###################################################
## #line 312 "vignettes/goseq/inst/doc/goseq.Rnw"
## gene_lengths=getlength(names(genes),"hg18","ensGene")
## pwf=nullp(genes,bias.data=gene_lengths)
## go_map=getgo(names(genes),"hg18","ensGene")
## go=goseq(pwf,"hg18","ensGene",gene2cat=go_map)


###################################################
### chunk number 27: KEGG_mappings eval=FALSE
###################################################
## #line 323 "vignettes/goseq/inst/doc/goseq.Rnw"
## #Get the mapping from ENSEMBL 2 Entrez
## en2eg=as.list(org.Hs.egENSEMBL2EG)
## #Get the mapping from Entrez 2 KEGG
## eg2kegg=as.list(org.Hs.egPATH)
## #Define a function which gets all unique KEGG IDs associated with a set of Entrez IDs
## grepKEGG=function(id,mapkeys){unique(unlist(mapkeys[id],use.names=FALSE))}
## #Apply this function to every entry in the mapping from ENSEMBL 2 Entrez to combine the two maps
## kegg=lapply(en2eg,grepKEGG,eg2kegg)
## head(kegg)


###################################################
### chunk number 28: KEGG eval=FALSE
###################################################
## #line 339 "vignettes/goseq/inst/doc/goseq.Rnw"
## pwf=nullp(genes,"hg18","ensGene")
## KEGG=goseq(pwf,gene2cat=kegg)
## head(KEGG)


###################################################
### chunk number 29: KEGG_goseq
###################################################
#line 349 "vignettes/goseq/inst/doc/goseq.Rnw"
pwf=nullp(genes,'hg18','ensGene')
KEGG=goseq(pwf,'hg18','ensGene',test.cats="KEGG")
head(KEGG)


###################################################
### chunk number 30: KEGG_from_db
###################################################
#line 361 "vignettes/goseq/inst/doc/goseq.Rnw"
kegg=as.list(org.Hs.egPATH)
head(kegg)


###################################################
### chunk number 31: countbias
###################################################
#line 376 "vignettes/goseq/inst/doc/goseq.Rnw"
countbias=rowSums(counts)[rowSums(counts)!=0]
length(countbias)
length(genes)


###################################################
### chunk number 32: GO.counts
###################################################
#line 384 "vignettes/goseq/inst/doc/goseq.Rnw"
pwf.counts=nullp(genes,bias.data=countbias)
GO.counts=goseq(pwf.counts,"hg18","ensGene")
head(GO.counts)


###################################################
### chunk number 33: setup
###################################################
#line 395 "vignettes/goseq/inst/doc/goseq.Rnw"
sessionInfo()