
---
title: "Lab 6: Answer Key"
output:
  BiocStyle::html_document:
    toc: true
vignette: >
  % \VignetteIndexEntry{Lab 6: BiocFileCache and Annotation Resources AnswerKey}
  % \VignetteEngine{knitr::rmarkdown}
---

```{r style, echo = FALSE, results = 'asis'}
BiocStyle::markdown()
```

```{r setup, echo=FALSE, warning=FALSE}
options(max.print=1000)
suppressPackageStartupMessages({
    library(org.Hs.eg.db)
    library(hugene20sttranscriptcluster.db)
    library(EnsDb.Hsapiens.v86)
    library(org.Hs.eg.db)
    library(TxDb.Hsapiens.UCSC.hg19.knownGene)
    library(Homo.sapiens)
    library(BSgenome)
    library(BSgenome.Hsapiens.UCSC.hg19)
    library(AnnotationHub)
    library(Organism.dplyr)
    library(biomaRt)
    })
```

# ChipDb/OrgDb questions

1. What is the central key for each of these AnnoDb objects? Can you print out the first 10 keys for each? 

ChipDb central key are the probe ids by default. The OrgDb `org.Hs.eg.db` is an
entrez orgDb (`eg` in name) and therefore the default is Entrez Gene Id. To
verify we print out the first 10 keys

```{r}
head(keys(hugene20sttranscriptcluster.db), 10)

head(keys(org.Hs.eg.db), 10)
```

2. Can all of the available columns be used as keytypes for each of these?   

```{r}
keytypes(org.Hs.eg.db)
columns(org.Hs.eg.db)

setdiff(columns(hugene20sttranscriptcluster.db), keytypes(hugene20sttranscriptcluster.db))
# so in this case yes
```

3. What gene symbol corresponds to Entrez Gene ID 1000? 

```{r}
select(org.Hs.eg.db, "1000","SYMBOL")
select(hugene20sttranscriptcluster.db, "1000","SYMBOL", "ENTREZID")
```

4. What is the Ensembl Gene ID for PPARG?

```{r}
select(org.Hs.eg.db, "PPARG", "ENSEMBL", "SYMBOL")
```

5. What is the UniProt ID for GAPDH?

```{r}
select(org.Hs.eg.db, "GAPDH", "UNIPROT","SYMBOL")
```

6. How many of the probeset ids from the ChipDb map to a single gene? How many don't map to a gene at all? 

```{r}
z <- mapIds(hugene20sttranscriptcluster.db, keys(hugene20sttranscriptcluster.db),
            "ENTREZID","PROBEID", multiVals = "list")
#or 

z <- mapIds(hugene20sttranscriptcluster.db, keys(hugene20sttranscriptcluster.db),
            "SYMBOL","PROBEID", multiVals = "list")

table(sapply(z, function(x) length(x[!is.na(x)])))
```

7. Get all the gene alias for a given Entrez Gene ID displayed as a CharacterList. What are the alias for Entrez Gene ID 1000?

```{r}
z <- mapIds(org.Hs.eg.db, keys(org.Hs.eg.db), "ALIAS", "ENTREZID", multiVals="CharacterList")

z[["1000"]]

# or

select(org.Hs.eg.db, "1000", "ALIAS") 
```

# TxDb/EnsDb

1.  How many Bioconductor TxDb packages are there? 

```{r}

# look on the website  or
length(BiocManager::available("TxDb"))
length(BiocManager::available("EnsDb"))
```

2.  What information can be inferred from the name? What about for  TxDb.Ggallus.UCSC.galGal4.refGene, TxDb.Rnorvegicus.BioMart.igis, EnsDb.Rnorvegicus.v79?

The name infers the species, source, build, and table used to generate the species. 

| Resource                          | Species          | Source | Build         | Table     |
|-----------------------------------|------------------|--------|---------------|-----------|
| TxDb.Hsapiens.UCSC.hg19.knownGene | Homo sapiens     | UCSC   | hg19 (GRCh37) | knownGene | 
| TxDb.Ggallus.UCSC.galGal4.refGene | Gallus gallus    | UCSC   | galGal4       | refGene   |
| TxDb.Rnorvegicus.BioMart.igis     | Rattus norvegicus| Biomart| igis          |           |
| EnsDb.Rnorvegicus.v79             | Rattus norvegicus| Ensembl| v79           |           |   


3. What do you think is the central key for TxDb objects? What other keytypes can be used? Can all of the available columns be used as keytypes?  

```{r}
head(keys(TxDb.Hsapiens.UCSC.hg19.knownGene))
keytypes(TxDb.Hsapiens.UCSC.hg19.knownGene)
# central key is TXID 

setdiff(columns(TxDb.Hsapiens.UCSC.hg19.knownGene), keytypes(TxDb.Hsapiens.UCSC.hg19.knownGene))
# These are columns available that are NOT available as keytypes... 
# so No not all columns can be used as keys
```

4. List all the genomic positions of the transcripts by gene. 

```{r}
transcriptsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, by="gene")
```

5. List all the promoter regions 1000bp upstream and 200bp download

```{r}
promoters(TxDb.Hsapiens.UCSC.hg19.knownGene, upstream=1000, downstream=200,  use.names=TRUE)
```

6. How many transcripts does PPARG have, according to UCSC? And does Ensembl agree? Is it right to compare these results? 

```{r}
select(org.Hs.eg.db, "PPARG", "ENTREZID", "SYMBOL")
# so entrez id is 5468

transcriptsBy(TxDb.Hsapiens.UCSC.hg19.knownGene)[["5468"]]

select(org.Hs.eg.db, "PPARG", "ENSEMBL", "SYMBOL")
# so ensembl id is ENSG00000132170

transcriptsBy(EnsDb.Hsapiens.v86)[["ENSG00000132170"]]

# No they don't agree 10 from entrez and 14 from ensembl
# No the builds are different 

library(TxDb.Hsapiens.UCSC.hg38.knownGene)
transcriptsBy(TxDb.Hsapiens.UCSC.hg38.knownGene)[["5468"]]

# The still don't agree 17 to 14
# There is a difference. Which is better? for you to decide
```

7. How many genes are between 2858473 and 3271812 on chr2 in the hg19 genome? 
    - Hint: you make a `GRanges` object
   
```{r}
gns <- genes(TxDb.Hsapiens.UCSC.hg19.knownGene) 
gns[gns %over%  GRanges("chr2:2858473-3271812")]

## OR

subsetByOverlaps(gns, GRanges("chr2:2858473-3271812"))

```

# Organism.Db

1. Get all the GO terms for BRCA1

```{r}
select(Homo.sapiens, "BRCA1", "GO", "SYMBOL")
```

2. What gene does the UCSC transcript ID uc002fai.3 map to?

```{r}
select(Homo.sapiens, "uc002fai.3", "SYMBOL", "TXNAME")
```


3. How many other transcripts does that gene have?

```{r}
transcriptsBy(Homo.sapiens)[[mapIds(Homo.sapiens, "uc002fai.3", "ENTREZID","TXNAME")]]
```


4. Get all the transcripts from the hg19 genome build, along with their Ensembl gene ID, UCSC transcript ID and gene symbol

```{r}
transcriptsBy(Homo.sapiens, columns = c("ENSEMBL", "SYMBOL"))
```

# Organism.dplyr


1. How many supported organisms are implemented in Organism.dplyr?

```{r}
library(Organism.dplyr)
supportedOrganisms()
dim(supportedOrganisms())[1]
```

2. Display the ensembl Id and genename description for symbol "NAT2".

```{r}
src <- src_organism(dbpath = hg38light())
select(src, "NAT2", c("ensembl", "genename"),"symbol")
```

3. Show all the alias for "NAT2" in the database.

```{r}
select(src, "NAT2", c("alias"),"symbol")
```

4. Get all the promoter regions.

```{r}
promoters(src)
```

5. Extract the "id" table. 

```{r}
tbl(src, "id")
```

6. Display Gene ontology (GO) information for gene symbol "Nat2".

```{r}
src_tbls(src)
colnames(tbl(src, "id_go"))
inner_join(tbl(src, "id"), tbl(src, "id_go")) %>%
		    filter(symbol == "NAT2") %>%
		    dplyr::select(entrez, ensembl, symbol, go, evidence, ontology)

```

# BSgenome

1. How many BSgenomes are there in Bioconductor?

```{r}
length(BiocManager::available("BSgenome"))
```

2. Get the sequence from UCSC hg19 builds for chromosome 1. And print the frequency of each letter.

```{r}
library(BSgenome.Hsapiens.UCSC.hg19)
seq = getSeq(BSgenome.Hsapiens.UCSC.hg19, "chr1")
alphabetFrequency(seq)
# only bases
alphabetFrequency(seq, baseOnly=TRUE)
```

3. Get the sequence corresponding to chr6 35310335-35395968. Get the complement and reverse complement of it.

```{r}
seq = getSeq(BSgenome.Hsapiens.UCSC.hg19, GRanges("chr6:35310335-35395968"))
complement(seq)
reverseComplement(seq)
```

4. Get the sequences for all transcripts of the TP53 gene

```{r}
tp53 <- transcriptsBy(Homo.sapiens)[[mapIds(Homo.sapiens, "TP53", "ENTREZID","SYMBOL")]]
getSeq(Hsapiens, tp53)
```

# AnnotationHub

1. How many resources are available in AnnotationHub?

```{r}
hub <- AnnotationHub()
length(hub)
```

2. How many resources are on AnnotationHub for Atlantic salmon (Salmo salar)?

```{r}
length(query(hub, "salmo salar"))
```

3. Get the most recent Ensembl build for domesticated dog (Canis familiaris) and make a TxDb

```{r}
dog <- query(hub, c("canis familiaris","ensembl","gtf"))
doggy <- dog[["AH69344"]]
TxDoggy <- makeTxDbFromGRanges(doggy)
```

4. Get information on the following ids "AH73986", "AH73881","AH64538", "AH69344".

```{r}
getInfoOnIds(hub, c("AH73986", "AH73881","AH64538", "AH69344"))
```

# biomaRt


1. List available marts

```{r}
listMarts()
```

2. Use mart corresponding to Ensembl Genes 97 and list available datasets.

```{r}
mart <- useMart("ENSEMBL_MART_ENSEMBL")
head(listDatasets(mart))
```

3. Create mart with Homo sapiens ensembl genes dataset

```{r}
mart <- useMart("ENSEMBL_MART_ENSEMBL","hsapiens_gene_ensembl")
```

4. Get the Ensembl gene IDs and HUGO symbol for Entrez Gene IDs 672, 5468 and 7157

```{r}
getBM(c("ensembl_gene_id","hgnc_symbol","entrezgene_id"), "entrezgene_id", c("672","5468","7157"), mart)
```

5. What do you get if you query for the 'gene_exon' for GAPDH?

```{r}
res <- getBM(c("hgnc_symbol","gene_exon"),"hgnc_symbol","GAPDH",mart)

```

6. What are the avaiable search terms?

```{r}
nms <- listAttributes(mart)$name
```

# BiocFileCache

1. Create a temporary cache using `tempdir()`. We use a temporary directory so it
will be cleaned up at the end.  If using to manage file would want a more permanent location.

```{r}
# the ask is not needed here
# because this is a vignette and not an interactive session when knitting, 
# this gives permission to create the cache on the system. 
bfc <- BiocFileCache(tempdir(), ask=FALSE)
```

2. What is the path to your cache? 

```{r}
bfccache(bfc)
bfc
```

3. What columns do we store in the database by default? 

```{r}
names(bfcinfo(bfc))
```


4. Get a file path to save an object so that it is tracked in the cache. 
Assume the object will be saved as a RData object. 

```{r}
savepath <- bfcnew(bfc, "NewResource", ext=".RData")
savepath

# we can see the info has been added
bfcinfo(bfc)

# use this path to save an object
m = matrix(1:12, nrow=3)
save(m, file=savepath)
```

5. Add a remote resource to the cache. The `httpbin` site has mock urls that can be used for testing
"http://httpbin.org/get". Give this resource the name "TestWeb"

```{r}
url <- "http://httpbin.org/get"
resource1 <- bfcadd(bfc, "TestWeb", fpath=url)
rid1 <- names(resource1)
```

6. Add another remote resource to the cache but do not automatically download. Call this resource "TestNoDownload".
You can use the same url as above or any of your choosing. 

```{r}
resource2 <- bfcadd(bfc, "TestNoDownload", fpath=url, download=FALSE)
rid2 <- names(resource2)
```

7. Check that the resources have been added to the cache.
Do a query for resources matching "Test". How many resources match? 

```{r}
bfcinfo(bfc)
bfccount(bfc)

res <- bfcquery(bfc, "Test")
res
bfccount(res)
```

8. Do resources need to be updated? 

```{r}
bfcneedsupdate(bfc)
```

9. Make a data.frame of metadata and add it to the cache. Do a query that would use the metadata.

```{r}
# the metdata data.frmae can be for a subset of data or the entire cache.
# the constrain is a column rid must match the rids in the cache.
rid = c("BFC1", "BFC3")
vls1 = c(1000, 2000)
vls2 = c("An experiment", "Online resource")

meta <- as.data.frame(list(rid=rid, vls1=vls1, vls2=vls2))
bfcmeta(bfc, name="resourceData") <- meta

names(bfcinfo(bfc))
bfcquerycols(bfc)

library(dplyr)
bfcinfo(bfc) %>% dplyr::select("rid", "vls2")

res <- bfcquery(bfc, field="vls2", "experiment")
res
# get id of resource matching query
# use that to load the resource
bfcrid(res)

load(bfcrpath(bfc, rids=bfcrid(res)))
```
