### R code from vignette source 'AnnotationExercisesBioc2011.Rnw'

###################################################
### code chunk number 1: setup
###################################################
library(Annotations)


###################################################
### code chunk number 2: install.packages (eval = FALSE)
###################################################
## install.packages("Annotations_1.0.4.tar.gz", repos=NULL, 
##                  type="source")


###################################################
### code chunk number 3: loadData
###################################################
  library(Annotations)
  load(system.file("data","tt.Rda",package="Annotations"))
  tt


###################################################
### code chunk number 4: get
###################################################
  library(org.Hs.eg.db)
  get("100127974", org.Hs.egSYMBOL)
  get("100127974", org.Hs.egPMID)
  otherGenes <- get("14702039", revmap(org.Hs.egPMID))
  head(otherGenes)


###################################################
### code chunk number 5: merge
###################################################
tt
toTable(head(org.Hs.egSYMBOL))
toTable(head(org.Hs.egCHR))
tt1 <- merge(tt, toTable(org.Hs.egSYMBOL), by.x="ID", by.y="gene_id")
merge(tt1, toTable(org.Hs.egCHR), by.x="ID", by.y="gene_id")


###################################################
### code chunk number 6: merge
###################################################
## One way
goIDs <- names(get("10013",org.Hs.egGO))
head(Term(goIDs))
## another way
goTab <- toTable(org.Hs.egGO["10013"])
head(cbind(goTab, Term(goTab$go_id)))


###################################################
### code chunk number 7: loadTxDb
###################################################
library("TxDb.Hsapiens.UCSC.hg18.knownGene")
txdb <- Hsapiens_UCSC_hg18_knownGene_TxDb
## Set ALL of the chromsomed to be inactive
isActiveSeq(txdb)[seqlevels(txdb)] <- FALSE
## Now set only chr7 to be active
isActiveSeq(txdb)["chr7"] <- TRUE
## Extract transcripts
txs <- transcriptsBy(txdb, by="gene")


###################################################
### code chunk number 8: readTable
###################################################
  load(system.file("data",
                   "ttDEX.Rda",
                   package="Annotations"))


###################################################
### code chunk number 9: readTableAns
###################################################
## find the Gene Symbols
library(org.Hs.eg.db)
mget(as.character(unique(ttDEX$geneID)), org.Hs.egSYMBOL, ifnotfound=NA)
## find the ranges of the transcripts
tdx <- names(txs) %in% as.character(unique(ttDEX$geneID)) 
txs[tdx]
## find the ranges of the exons too
exs <- exonsBy(txdb, by="gene")
tdx <- names(exs) %in% as.character(unique(ttDEX$geneID))
exs[as.character(unique(ttDEX$geneID))]


###################################################
### code chunk number 10: readGappedAligns
###################################################
  ga <- readGappedAlignments(system.file("extdata",
                                         "chr7Cont1.bam",
                                         package="Annotations"))


###################################################
### code chunk number 11: readGappedAlignAns
###################################################
library(org.Hs.eg.db)
get("CFTR", revmap(org.Hs.egSYMBOL))
txs["1080"]
## then you can count
countOverlaps(txs["1080"], ga)


###################################################
### code chunk number 12: fdb
###################################################
## get oreganno data:
fdb <- makeFeatureDbFromUCSC(genome="mm9",track="oreganno",
                             tablename="oreganno")

## filter out c-somes for txdb
isActiveSeq(txdb)[seqlevels(txdb)] <- FALSE
isActiveSeq(txdb)[c(paste("chr",1:22,sep=""),"chrM","chrX","chrY")] <- TRUE

## get elements
txs <- transcriptsBy(txdb, by="gene")
f <- features(fdb)

## use findOverlaps
res <- findOverlaps(txs,f)
head(matchMatrix(res))

## get data for 1st match
txs[73]
f[8287]

## now get the gene with the most matches
mtchs <- matchMatrix(res)
sort(table(mtchs[,1]))

## now look that gene up.
txs[6307]
## get the oreganno elements:
mtchs[mtchs[,1]==6307,]
f[mtchs[mtchs[,1]==6307,2]]
## NOW get the gene symbols
get("6307", org.Hs.egSYMBOL)