```{r setup, include=FALSE}
library(edgeR)
library(org.Mm.eg.db)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(ggbio)
library(GenomicRanges)
library(GenomicAlignments)
```

```{r prepareData}
# load data
load("../Robjects/DE.Rdata")
# make results table
results <- as.data.frame(topTags(lrt_BvsL, n = Inf))
# identify differentially expressed genes
de <- decideTestsDGE(lrt_BvsL)
# create a vector of differentially expressed genes
detags <- rownames(dgeObj)[as.logical(de)]
# fetch annotations
ann <- select(org.Mm.eg.db,
              keys = rownames(results),
              columns = c("ENTREZID","SYMBOL","GENENAME"))
# annotate the results
results.annotated <- cbind(results, ann)
# load the txDB object
tx <- TxDb.Mmusculus.UCSC.mm10.knownGene
exo <- exonsBy(tx, "gene")
myreg <- reduce(exo[["110308"]])
# make gene model
geneMod <- autoplot(tx, which = myreg)
```


> ## Challenge 1 {.challenge}
>
> Re-visit the `plotSmear` plot from above and use the `text` function to add labels for the names of the top 200 most DE genes
>

```{r solutionChallenge1, echo=FALSE, fig.height=5, fig.width=10}
n_genes <- 20
plotSmear(lrt_BvsL, de.tags=detags)
text(results.annotated$logCPM[seq(n_genes)],
     results.annotated$logFC[seq(n_genes)],
     labels = results.annotated$SYMBOL[seq(n_genes)],
     col = "blue")
```

> ## Challenge 2 {.challenge}
>
> Use the TxDb.Mmusculus.UCSC.mm10.knownGene package to retrieve the exon coordinates for the genes `50916`, `110308`, `12293` 
>

```{r solutionChallenge2, echo=FALSE, warning=FALSE, message=FALSE}

keys <- c("50916", "110308", "12293")
select(tx, 
       keys=keys,
       keytype = "GENEID",
       columns=c("EXONCHROM","EXONSTART","EXONEND")
      )

```

> ## Challenge 3a {.challenge}
>
> Create tracks to compare the coverage of the gene Krt5 for the samples MCL1.DG, MCL1.DH, MCL1.LA and MCL1.LB
>

```{r solutionChallenge3a, echo=FALSE, fig.height=5, fig.width=10}
bam <- readGAlignments(file="../data/bams_small/MCL1.DG.small.bam",
                       param=ScanBamParam(which=myreg),use.names = TRUE)
reads.MCL1.DG <- autoplot(bam, stat = "coverage") + 
    xlim(GRanges("chr15", IRanges(101707000, 101713000)))

bam <- readGAlignments(file="../data/bams_small/MCL1.DH.small.bam",
                       param=ScanBamParam(which=myreg),use.names = TRUE)
reads.MCL1.DH <- autoplot(bam, stat = "coverage") + 
    xlim(GRanges("chr15", IRanges(101707000, 101713000)))


bam <- readGAlignments(file="../data/bams_small/MCL1.LA.small.bam",
                       param=ScanBamParam(which=myreg),use.names = TRUE)
reads.MCL1.LA <- autoplot(bam, stat = "coverage") + 
    xlim(GRanges("chr15", IRanges(101707000, 101713000)))

bam <- readGAlignments(file="../data/bams_small/MCL1.LB.small.bam",
                       param=ScanBamParam(which=myreg),use.names = TRUE)
reads.MCL1.LB <- autoplot(bam, stat = "coverage") + 
    xlim(GRanges("chr15", IRanges(101707000, 101713000)))


tracks(mm10=geneMod, MCL1.DG=reads.MCL1.DG, MCL1.Dh=reads.MCL1.DH, MCL1.LA=reads.MCL1.LA, MCL1.LB=reads.MCL1.LB) 

```

