---
output:
  knitrBootstrap::bootstrap_document:
    theme.chooser: TRUE
    highlight.chooser: TRUE
---

Compare vs PNAS
===============


# Counts-based analysis

This section has the code for running `edgeR-robust` and `DESeq2` on the simulation data set using the known exons as features. It is based on [counts-based.html](../../simulation/counts-based/counts-based.html).


This first code chunk loads the necessary data.


```{r 'setup', bootstrap.show.code = FALSE, bootstrap.show.message = FALSE}
## Track time spent on making the report
startTime <- Sys.time()

library('edgeR')
library('DESeq2')
library('GenomicRanges')

## Load data
load("../coverageToExon/covToEx-ucsc.Rdata")
load("../derAnalysis/run3-v1.0.10/groupInfo.Rdata")
load("../derAnalysis/run3-v1.0.10/colsubset.Rdata")

## GenomicState object
if(file.exists('/home/epi/ajaffe/Lieber/Projects/RNAseq/derannotator/rdas/GenomicState.Hsapiens.UCSC.hg19.knownGene.rda')) {
    load('/home/epi/ajaffe/Lieber/Projects/RNAseq/derannotator/rdas/GenomicState.Hsapiens.UCSC.hg19.knownGene.rda')
} else if(file.exists('../../GenomicState.Hsapiens.UCSC.hg19.knownGene.rda')) {
    load('../../GenomicState.Hsapiens.UCSC.hg19.knownGene.rda')
} else {
    stop('Missing UCSC hg19 genomic state object')
}

## Annotation used
exons <- GenomicState.Hsapiens.UCSC.hg19.knownGene$fullGenome
exons <- exons[exons$theRegion == 'exon']

## Round matrix and remove exons with 0s
counts <- round(covToEx[, colsubset])
nonzero <- sapply(rowSums(counts), function(x) {x > 0})
```


## DESeq2

The following code performs the DESeq2 analysis. Code is based on [edgeR_Robust supplementary code](http://imlspenticton.uzh.ch/robinson_lab/edgeR_robust/). The main change is that it has been modified for the multi-group scenario.

```{r 'deseq2', bootstrap.show.code = FALSE}
## Round matrix and specify design
dse <- DESeqDataSetFromMatrix(counts[nonzero, ], data.frame(group = groupInfo), ~ group)

## Perform DE analysis
system.time( dse <- DESeq(dse, test = 'LRT', reduced = ~ 1) )

## Extract results
deseq <- exons[nonzero]
mcols(deseq) <- cbind(mcols(deseq), results(dse))

## Which are significant?
deseq$sig <- deseq$padj < 0.05
deseq$sig[is.na(deseq$sig)] <- FALSE

## Save results
save(deseq, file = 'deseq.Rdata')

## Adjust by Holm
deseq_holm <- deseq
deseq_holm$sig <- p.adjust(deseq_holm$pvalue, 'holm') < 0.05
```


## edgeR-robust

The following code performs the DESeq2 analysis. Code is based on [edgeR_Robust supplementary code](http://imlspenticton.uzh.ch/robinson_lab/edgeR_robust/). The main change is that it has been modified for the multi-group scenario.


```{r 'edgeR', bootstrap.show.code = FALSE}
## Determine design matrix
design <- model.matrix(~ groupInfo)

## Perform DE analysis
d <- DGEList(counts = counts[nonzero, ], group = groupInfo)
d <- calcNormFactors(d)
system.time(dw <- estimateGLMRobustDisp(d, design = design, prior.df = 10, maxit = 6))
fw <- glmFit(dw, design = design, coef = 2:3)
lrw <- glmLRT(fw, coef = 2:3)

## Extract results
edger <- exons[nonzero]
mcols(edger) <- cbind(mcols(edger), DataFrame(lrw$table))
edger$pvalue <-  lrw$table$PValue
edger$padj <- p.adjust(lrw$table$PValue, 'BH')

## Which are significant?
edger$sig <- edger$padj < 0.05
edger$sig[is.na(edger$sig)] <- FALSE

## Save results
save(edger, file = 'edger.Rdata')

## Adjust by Holm
edger_holm <- edger
edger_holm$sig <- p.adjust(edger_holm$pvalue, 'holm') < 0.05
```








## Overlap

```{r 'ov-comp-setup', bootstrap.show.code = FALSE}
## Load data
load('../derAnalysis/run3-v1.0.10/fullRegions.Rdata')

## Some formatting and subsets
names(fullRegions) <- seq_len(length(fullRegions))
fullRegions$sigFWER <- as.logical(fullRegions$significantFWER)
fullRegs20 <- fullRegions[width(fullRegions) >= 20]

## Overlap table for all 4 cases
ov_table <- function(ders, counts, query = 'der', minov = 0) {
    if(query == 'der') {
        if(minov == 0) {
            res <- addmargins(table('Significant DER (FWER)' = ders$sigFWER, 'Overlaps sig DE exon' = countOverlaps(ders, counts[counts$sig]) > 0))
        } else {
            res <- addmargins(table(ders$sigFWER, countOverlaps(ders, counts[counts$sig], minoverlap = minov) > 0, dnn = c('Significant DER (FWER)', paste0('Overlaps sig DE exon (min ', minov, 'bp)'))))
        }
    } else if (query == 'counts') {
        if(minov == 0) {
            res <- addmargins(table('Significant DE exon' = counts$sig, 'Overlaps sig DER (FWER)' = countOverlaps(counts, ders[ders$sigFWER]) > 0))
        } else {
            res <- addmargins(table(counts$sig[width(counts) >= minov], countOverlaps(counts[width(counts) >= minov], ders[ders$sigFWER], minoverlap = minov) > 0, dnn = c('Significant DE exon', paste0('Overlaps sig DER (FWER, min ', minov, 'bp)'))))
        }
    }
    return(res)
}

## Explore mistmatched cases for DERs vs Exons direction
explore_ov <- function(ders, counts, case = "FALSE-TRUE", minov = 0L) {
    if(case == 'FALSE-TRUE') {
        i <- which(countOverlaps(ders, counts[counts$sig], minoverlap = minov) > 0 & !ders$sigFWER)
    } else if (case == 'TRUE-FALSE') {
        i <- which(!countOverlaps(ders, counts[counts$sig], minoverlap = minov) > 0 & ders$sigFWER)
    } else{
        stop('invalid case')
    }
    if(length(i) == 0) return("No such cases")
    
    if(case == 'FALSE-TRUE') {
        res <- list(
            n_overlaps = table(countOverlaps(ders[i], counts[counts$sig], minoverlap = minov)),
            width_der = summary(width(ders[i])),
            ders_per_exon_table = table(table(subjectHits(findOverlaps(ders[i], counts[counts$sig], minoverlap = minov)))),
            ders_per_exon = sort(table(subjectHits(findOverlaps(ders[i], counts[counts$sig], minoverlap = minov)))),
            i = i
        )
    } else {
        res <- list(
            width_der = summary(width(ders[i])),
            distance_nearest_sum = summary(mcols(distanceToNearest(ders[i], counts, ignore.strand = TRUE))$distance),
            distance_nearest_sig_sum = summary(mcols(distanceToNearest(ders[i], counts[counts$sig], ignore.strand = TRUE))$distance),
            distance_nearest = distanceToNearest(ders[i], counts, ignore.strand = TRUE),
            distance_nearest_sig = distanceToNearest(ders[i], counts[counts$sig], ignore.strand = TRUE),
            i = i
        )
    }
    
    return(res)
}

## Explore mistmatched cases for Exons vs DERs direction
explore_ov_counts <- function(ders, counts, case = "FALSE-TRUE", minov = 0L) {
    counts <- counts[width(counts) >= minov]
    if(case == 'FALSE-TRUE') {
        i <- which(countOverlaps(counts, ders[ders$sigFWER], minoverlap = minov) > 0 & !counts$sig)
    } else if (case == 'TRUE-FALSE') {
        i <- which(!countOverlaps(counts, ders[ders$sigFWER], minoverlap = minov) > 0 & counts$sig)
    } else{
        stop('invalid case')
    }
    if(length(i) == 0) return("No such cases")
    
    if(case == 'FALSE-TRUE') {
        res <- list(
            n_overlaps = table(countOverlaps(counts[i], ders[ders$sigFWER], minoverlap = minov)),
            width_exon = summary(width(counts[i])),
            exons_per_der_table = table(table(subjectHits(findOverlaps(counts[i], ders[ders$sigFWER], minoverlap = minov)))),
            exons_per_der = sort(table(subjectHits(findOverlaps(counts[i], ders[ders$sigFWER], minoverlap = minov)))),
            i = i
        )
    } else {
        res <- list(
            width_exon = summary(width(counts[i])),
            distance_nearest_sum = summary(mcols(distanceToNearest(counts[i], ders, ignore.strand = TRUE))$distance),
             distance_nearest_sig_sum = summary(mcols(distanceToNearest(counts[i], ders[ders$sigFWER], ignore.strand = TRUE))$distance),
            distance_nearest = distanceToNearest(counts[i], ders, ignore.strand = TRUE),
            distance_nearest_sig = distanceToNearest(counts[i], ders[ders$sigFWER], ignore.strand = TRUE),
            i = i
        )
    }
    
    return(res)
}

noNA <- function(x) {
    x[!is.na(x)]
}
```


### DESeq2


#### Query: DERs

We can first compare the results by using the DERs as the query and the exons as the subject. The following output shows the comparison using all DERs and exploring the mismatched cases. Then its repeated using the DERs $\geq$ 20 bp and a minimum overlap of 20bp.

For the mismatched cases of non-significant DERs overlapping a significant exon, we check:

* how many exons each DER overlaps,
* the width of the DERs
* the frequency table of how many DERs overlap the same exon

For the other mismatched case, we check:

* the width of the DERs
* distance to nearest exon (regardless of exon size)
* distance to nearest significant DE exon (ibidem)

```{r 'ov-comp-deseq', bootstrap.show.code = FALSE}
## Overlap between DERs and significant DE exons
ov_table(fullRegions, deseq)

## Explore mismatched cases
#noNA(explore_ov(fullRegions, deseq)[1:3])
#noNA(explore_ov(fullRegions, deseq, 'TRUE-FALSE')[1:3])

## Min 20 bp overlap, using only DERs 20 bp long
ov_table(fullRegs20, deseq, minov = 20L)

## Explore mismatched cases, min 20bp overlap
noNA(explore_ov(fullRegs20, deseq, minov = 20L)[1:3])
noNA(explore_ov(fullRegs20, deseq, 'TRUE-FALSE', minov = 20L)[1:3])

## Holm vs BH
addmargins(table('DESeq2 Holm' = deseq_holm$sig, 'DESeq2 BH' = deseq$sig))

## Use Holm and min 20 bp ov
ov_table(fullRegs20, deseq_holm, minov = 20L)
```

Most of the DERs are shorter than 20bp (`r round(sum(width(fullRegions) < 20) / length(fullRegions) * 100, 2)` percent), so we'll focus on the longer ones. The majority of the mismatches are from non significant DERs that overlap a significant exon. 

As expected, when controlling the FWER instead of the FDR, most of the DE exons are no longer significant. Using FWER-controlled DE exons, most of the DERs 20bp or longer agree with the exons as not being significantly DE.


#### Query: exons

We can now repeat the comparison using the exons as the query and the DERs as the subject.

For the mismatched cases of non-significant exons overlapping a significant DER, we check:

* how many DERs each exon overlaps,
* the width of the exons
* the frequency table of how many exons overlap the same DER

For the other mismatched case, we check:

* the width of the exons
* distance to nearest DER (regardless of DER size)
* distance to nearest significant DER (ibidem)

```{r 'ov-comp-deseq-counts', bootstrap.show.code = FALSE}
## Overlap between exons and significant DERs
#ov_table(fullRegions, deseq, 'counts')

## Explore mismatched cases
#noNA(explore_ov_counts(fullRegions, deseq)[1:3])
#noNA(explore_ov_counts(fullRegions, deseq, 'TRUE-FALSE')[1:3])

## Overlap between exons and significant DERs, min 20 bp
ov_table(fullRegions, deseq, 'counts', 20L)

## Explore mismatched cases
noNA(explore_ov_counts(fullRegions, deseq, minov = 20L)[1:3])
noNA(explore_ov_counts(fullRegions, deseq, 'TRUE-FALSE', minov = 20L)[1:3])

## Now with Holm
ov_table(fullRegions, deseq_holm, 'counts', 20L)
```

From these results, we can see that `derfinder` is more conservative.




### edgeR-robust

#### Query: DERs

Similar comparison using DERs as query and exons as subject with `edgeR-robust` results.

```{r 'ov-comp-edger', bootstrap.show.code = FALSE}
## Overlap between DERs and significant DE exons
#ov_table(fullRegions, edger)

## Explore mismatched cases
#noNA(explore_ov(fullRegions, edger)[1:3])
#noNA(explore_ov(fullRegions, edger, 'TRUE-FALSE')[1:3])

## Min 20 bp overlap, using only DERs 20 bp long
ov_table(fullRegs20, edger, minov = 20L)

## Explore mismatched cases, min 20bp overlap
noNA(explore_ov(fullRegs20, edger, minov = 20L)[1:3])
noNA(explore_ov(fullRegs20, edger, 'TRUE-FALSE', minov = 20L)[1:3])

## Holm vs BH
addmargins(table('edgeR Holm' = edger_holm$sig, 'edger BH' = edger$sig))

## With Holm, 20bp
ov_table(fullRegs20, edger_holm, minov = 20L)
```

The results are fairly similar to those from using `DESeq2`.

#### Query: exons

Similar comparison using exons as query and DERs as subject with `edgeR-robust` results.

```{r 'ov-comp-edger-counts', bootstrap.show.code = FALSE}
## Overlap between exons and significant DERs
#ov_table(fullRegions, edger, 'counts')

## Explore mismatched cases
#noNA(explore_ov_counts(fullRegions, edger)[1:3])
#noNA(explore_ov_counts(fullRegions, edger, 'TRUE-FALSE')[1:3])

## Overlap between exons and significant DERs, min 20 bp
ov_table(fullRegions, edger, 'counts', 20L)

## Explore mismatched cases
noNA(explore_ov_counts(fullRegions, edger, minov = 20L)[1:3])
noNA(explore_ov_counts(fullRegions, edger, 'TRUE-FALSE', minov = 20L)[1:3])

## With Holm, 20 bp
ov_table(fullRegions, edger_holm, 'counts', 20L)
```


### overall

While the DERs vs exons results are fairly similar between `edgeR-robust` and `DESeq2`, as shown below the number of mismatched cases is high compared to the number of cases both counts-based methods agree. This is also true when controlling the FWER to determine significance.

```{r 'deseq-vs-edger'}
## edgeR vs DESeq2
addmargins(table('edgeR-robust (FDR)' = edger$sig, 'DESeq2 (FDR)' = deseq$sig))

## Control FWER
addmargins(table('edgeR-robust (FWER)' = edger_holm$sig, 'DESeq2 (FWER)' = deseq_holm$sig))

## Only sig if both edgeR and DEseq2 say it is
both <- deseq
both$sig <- both$sig & edger$sig

## Same, for holm
both_holm <- deseq_holm
both_holm$sig <- both_holm$sig & edger_holm$sig
```

We can consider an exon to be DE only if both `edgeR-robust` and `DESeq2` find that its significantly DE. The next sections use this information.

#### Query: DERs

```{r 'ov-comp-both', bootstrap.show.code = FALSE}
## Overlap between DERs and significant DE exons
#ov_table(fullRegions, both)

## Explore mismatched cases
#noNA(explore_ov(fullRegions, both)[1:3])
#noNA(explore_ov(fullRegions, both, 'TRUE-FALSE')[1:3])

## Min 20 bp overlap, using only DERs 20 bp long
ov_table(fullRegs20, both, minov = 20L)

## Explore mismatched cases, min 20bp overlap
noNA(explore_ov(fullRegs20, both, minov = 20L)[1:3])
noNA(explore_ov(fullRegs20, both, 'TRUE-FALSE', minov = 20L)[1:3])

## Holm vs BH
addmargins(table('Both Holm' = both_holm$sig, 'Both BH' = both$sig))

## Use Holm and min 20 bp ov
ov_table(fullRegs20, both_holm, minov = 20L)
```

The trends observed previously are maintained in this comparison with a reduction of cases where the exon is DE. This is expected due to the non-perfect agreement between `DESeq2` and `edgeR-robust`.


```{r 'regionPlot-setup', bootstrap.show.code = FALSE, bootstrap.show.message = FALSE}
library('TxDb.Hsapiens.UCSC.hg19.knownGene')
library('derfinder')
library('derfinderHelper')
library('derfinderPlot')
load('../derAnalysis/run3-v1.0.10/models.Rdata')
load('../derAnalysis/run3-v1.0.10/chr22/optionsStats.Rdata')
load("../CoverageInfo/fullCov.Rdata")

def.par <- par()
def.par <- def.par[-which(names(def.par) %in% c('cin', 'cra', 'csi', 'cxy', 'din', 'page'))]

regPlot <- function(region, title) {
    ## Calculate F-stats
    range <- start(region):end(region)
    dat <- fullCov[[as.character(seqnames(region))]][range, colsubset]

    ## Log2 transform
    for(i in seq_len(length(groupInfo))) dat[[i]] <- log2(dat[[i]] + 32) 

    ## Calculate f-stats
    fstats <- as.numeric(fstats.apply(data = dat, mod = models$mod, mod0 = models$mod0))

    ## Find annotation
    annoReg <- annotateRegions(region, GenomicState.Hsapiens.UCSC.hg19.knownGene$fullGenome, verbose = FALSE)
    symbol <- mcols(annoReg$annotationList[[1]])$symbol
    symbol <- as.character(noNA(symbol)[[1]])
    if(length(symbol) > 1) symbol <- symbol[1]
    symbol <- ifelse(is.null(symbol), NA, symbol)
    ## Remove symbol name because it gets chomped on the plot
    mcols(annoReg$annotationList[[1]])$symbol <- NA
    
    par(def.par)

    ## Plot long exon
    plotRegionCoverage(region, getRegionCoverage(fullCov, region, verbose = FALSE), groupInfo, data.frame(name = title, distance = NA, region = symbol), annoReg, verbose = FALSE, ask = FALSE, txdb = TxDb.Hsapiens.UCSC.hg19.knownGene)

    ## Add F-stat track
    par(fig = c(0, 1, 0.065, 0.125), new = TRUE, xaxt = 'n', oma = c(0, 0, 0, 0), mar = c(0, 4.5, 0, 1.1))
    plot(y = fstats, x = range, ylab = 'F-stat', type = 'l', xlab = '', bty = 'n', ylim = c(0, max(fstats[is.finite(fstats)], optionsStats$cutoffFstatUsed) * 1.1), las = 2, yaxt = 'n')
    y.max <- round(max(c(optionsStats$cutoffFstatUsed, fstats[is.finite(fstats)]), na.rm = TRUE), 0)
    axis(2, at = c(0, round(y.max / 2, 0), y.max), c(0, round(y.max / 2, 0), y.max), las = 2, tick = TRUE)
    abline(h = optionsStats$cutoffFstatUsed, col = 'red')
    abline(h = 0, col = 'grey')
}

sortWidth <- function(regions) {
    regions[order(width(regions), decreasing = TRUE)]
}
```

We can now make plots to explore some DERs for each of the cases.

```{r 'query-der-plots', fig.width = 10, fig.height = 7, bootstrap.show.code = FALSE, dev = 'CairoPNG'}
query_der_plots <- function() {
    sapply(sortWidth(fullRegs20[countOverlaps(fullRegs20, both[both$sig], minoverlap = 20L) > 0 & fullRegs20$sigFWER])[1:10], function(reg) {
        regPlot(reg, 'DER query: DE agreement')
    })

    sapply(sortWidth(fullRegs20[countOverlaps(fullRegs20, both[both$sig], minoverlap = 20L) == 0 & !fullRegs20$sigFWER])[1:10],  function(reg) {
        regPlot(reg, 'DER query: not DE agreement')
    })

    sapply(sortWidth(fullRegs20[countOverlaps(fullRegs20, both[both$sig], minoverlap = 20L) == 0 & fullRegs20$sigFWER])[1:10], function(reg) {
        regPlot(reg, 'DER query: only exon not DE')
    })

    sapply(sortWidth(fullRegs20[countOverlaps(fullRegs20, both[both$sig], minoverlap = 20L) > 0 & !fullRegs20$sigFWER])[1:10], function(reg) {
        regPlot(reg, 'DER query: only exon DE')
    })
}
pdf(file = 'query_der_plots.pdf', width = 10, height = 7)
query_der_plots()
dev.off()

query_der_plots()
```


#### Query: exons

As was shown with either `DESeq2` or `edgeR-robust` results, `derfinder` is more conservative than the counts-based methods.

```{r 'ov-comp-both-counts', bootstrap.show.code = FALSE}
## Overlap between exons and significant DERs, min 20 bp
ov_table(fullRegions, both, 'counts', 20L)

## Explore mismatched cases
noNA(explore_ov_counts(fullRegions, both, minov = 20L)[1:3])
noNA(explore_ov_counts(fullRegions, both, 'TRUE-FALSE', minov = 20L)[1:3])

## With Holm, 20 bp
ov_table(fullRegions, both_holm, 'counts', 20L)
```

We can now visually explore some exons for each of the four cases.


```{r 'query-exon-plots', fig.width = 10, fig.height = 7, bootstrap.show.code = FALSE, dev = 'CairoPNG'}
query_exon_plots <- function() {
    sapply(sortWidth(both[width(both) >= 20 & both$sig & countOverlaps(both, fullRegions[fullRegions$sigFWER], minoverlap = 20L) > 0])[1:10], function(reg) {
        regPlot(reg, 'Exon query: DE agreement')
    })

    sapply(sortWidth(both[width(both) >= 20 & !both$sig & countOverlaps(both, fullRegions[fullRegions$sigFWER], minoverlap = 20L) == 0])[1:10], function(reg) {
        regPlot(reg, 'Exon query: not DE agreement')
    })

    sapply(sortWidth(both[width(both) >= 20 & !both$sig & countOverlaps(both, fullRegions[fullRegions$sigFWER], minoverlap = 20L) > 0])[1:10], function(reg) {
        regPlot(reg, 'Exon query: only exon not DE')
    })

    sapply(sortWidth(both[width(both) >= 20 & both$sig & countOverlaps(both, fullRegions[fullRegions$sigFWER], minoverlap = 20L) == 0])[1:10], function(reg) {
        regPlot(reg, 'Exon query: only exon DE')
    })
}
pdf(file = 'query_exon_plots.pdf', width = 10, height = 7)
query_exon_plots()
dev.off()

query_exon_plots()
```





# Finding regions of interest


The code in this section is partially based on `/home/epi/ajaffe/Lieber/Projects/RNAseq/HippoPublic/clean_previous_hits.R`.


First we find the regions of the genome corresponding to the genes of interest. That is, the genes from the original paper that were differentially expressed in the conditions we are analyzing.

```{r pkgs, bootstrap.show.message=FALSE, bootstrap.show.code = FALSE}
## Required pkgs
library("GenomicRanges")
library("ggbio")
library("reshape2")
library("plyr")
library("scales")
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
library("org.Hs.eg.db")
```


```{r identify, bootstrap.show.code = FALSE}
## Gene symbol names of interest
#### Original names
## symbols <- c("HIST1H4E", "RN7SK", "CDR1", "SNORD89", "SNORA73A", "SCARNA17", "PAPD1", "CACNB2", "LRCH4", "SNORD42A", "SNORA47", "LENG8", "FAM123A", "HIVEP3", "HNRPH1", "ZGPAT", "ERF", "SNORD116-29", "C9orf139", "C9orf3", "KCNA2", "EXOC6B", "CENTB5", "TAOK2", "TNRC6C", "ADAMTS4", "MSH4", "C16orf72", "CCR5")

## http://www.genenames.org/data/hgnc_data.php?hgnc_id=25532
## http://www.genenames.org/data/hgnc_data.php?hgnc_id=26360
## http://www.genenames.org/data/hgnc_data.php?hgnc_id=5041
## http://www.genenames.org/data/hgnc_data.php?hgnc_id=16754

## Updated symbols
symbols <- c("HIST1H4E", "RN7SK", "CDR1", "SNORD89", "SNORA73A", "SCARNA17", "MTPAP", "CACNB2", "LRCH4", "SNORD42A", "SNORA47", "LENG8", "AMER2", "HIVEP3", "HNRNPH1", "ZGPAT", "ERF", "SNORD116-29", "C9orf139", "C9orf3", "KCNA2", "EXOC6B", "ACAP3", "TAOK2", "TNRC6C", "ADAMTS4", "MSH4", "C16orf72", "CCR5")


## Map gene symbol names to entrezid's

keys <- keys(org.Hs.eg.db, keytype = "ENTREZID")
columns <- c("SYMBOL")
map <- select(org.Hs.eg.db, keys, columns, keytype = "ENTREZID")
idx <- sapply(symbols, function(x) { 
	res <- which(map$SYMBOL == x)
	ifelse(length(res) > 0, res, NA)
})
ids <- map$ENTREZID[idx]
names(ids) <- names(idx)

## Remove those not-found
ids <- ids[!is.na(ids)]

## Find the exons

txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
exonsUCSC <- exons(txdb, vals = list(gene_id =  ids), columns = c("gene_id", "exon_id", "tx_name", "tx_id"))

## some gene ids have multiples, straighten out
gids <- as.list(exonsUCSC$gene_id)
for(i in seq(along=gids)) gids[[i]] <- gids[[i]][which(gids[[i]] %in% ids)]
gids <- unlist(gids)


## split into list by EID
exonListUCSC <- split(exonsUCSC, gids)
exonListUCSC <- exonListUCSC[ ids[ids %in% gids] ]
### drop duplicated exons
# exonListUCSC = lapply(exonListUCSC, function(x) x[!duplicated(x)])
# identical(names(exonListUCSC), ids[ids %in% gids]) # TRUE

## Not found
ids[which(!ids %in% gids)]

## Find them manually
# http://www.ncbi.nlm.nih.gov/gene/125050
# http://www.ncbi.nlm.nih.gov/gene/6080
missing <- GRanges(seqnames=c("chr6", "chr1"), ranges=IRanges(start=c(52860418, 28833877), end=c(52860749, 28834083)))
toAdd <- split(missing, 1:2)
names(toAdd) <- ids[which(!ids %in% gids)]

## Reduce to min/max per gene
windows <- c(GRangesList(lapply(exonListUCSC, range)), toAdd)

## Save for later use
save(windows, ids, idx, file="windows.Rdata")
```

# Original genes

In this section, we make a plot for each gene showing the coverage data and whether `derfinder` identified candidate DERs as described in the main text.

```{r plots, message=FALSE, fig.width=20, fig.height=10, dev="CairoPNG", bootstrap.show.code = FALSE}
## Find chrs used
chrs <- as.character(unique(unlist(seqnames(windows), use.names=FALSE)))

## Build ideograms
data(hg19IdeogramCyto, package = "biovizBase")
p.ideos <- lapply(chrs, function(xx) { 
	plotIdeogram(hg19IdeogramCyto, xx)
})
names(p.ideos) <- chrs


## Filter data
fullCovSmall <- lapply(chrs, function(chr) {
	fullCov[[chr]][, colsubset]
})
names(fullCovSmall) <- chrs
rm(fullCov)

## Main plotting function
plotClusterCustom <- function(cluster, regions, titleName, coverageInfo, groupInfo, titleUse="fwer", txdb=NULL, p.ideogram=NULL, maxExtend=300L, colsubset=NULL, forceLarge=FALSE) {

	stopifnot(is.factor(groupInfo))
	if(is.null(colsubset)) colsubset <- seq_len(length(groupInfo))
	
	## Window length
	l <-  width(cluster) + 2 * min(maxExtend, width(cluster))
	
	if(l > 1e5 & !forceLarge) {
		message(paste("No plot will be made because the data is too large. The window size exceeds 100 kb."))
		return(invisible(l))
	}
	
	wh <- resize(cluster, l, fix="center")
	title <- paste("Window view for ENTREZ Symbol", titleName)
	
	## Plot the ideogram if not supplied
	if(is.null(p.ideogram)) {
		chr <- as.character(seqnames(wh))
		## Now load the ideogram info
		hg19IdeogramCyto <- NULL
		load(system.file("data", "hg19IdeogramCyto.rda", package="biovizBase", mustWork=TRUE))
		p.ideogram <- plotIdeogram(hg19IdeogramCyto, chr)
	}
	
	## Regions found (from the view)
	neighbors <- regions[queryHits(findOverlaps(regions, wh))]
	if(length(neighbors) == 0) {
		neighbors <- wh
		neighbors$significant <- NA
		neighbors$significantQval <- NA
        neighbors$significantFWER <- NA
	} 
	if(titleUse == "pval") {
		p.region <- autoplot(neighbors, aes(fill=significant)) + 
		scale_fill_manual(values=c("chartreuse4", "wheat2"), limits=c("TRUE", "FALSE")) 
	} else if (titleUse == "qval" ){
		p.region <- autoplot(neighbors, aes(fill=significantQval)) +
		scale_fill_manual(values=c("chartreuse4", "wheat2"), limits=c("TRUE", "FALSE")) 
	} else if (titleUse == "fwer" ){
		p.region <- autoplot(neighbors, aes(fill=significantFWER)) +
		scale_fill_manual(values=c("chartreuse4", "wheat2"), limits=c("TRUE", "FALSE")) 
	} else {
		p.region <- autoplot(neighbors)
	}

	## Graphical parameters
	nGroups <- length(levels(groupInfo))
	
	## Construct the coverage plot
	pos <- start(wh):end(wh)
	rawData <- as.data.frame(coverageInfo[pos, colsubset])
	rawData$position <- pos
	covData <- melt(rawData, id.vars="position")
	covData$group <- rep(groupInfo, each=nrow(rawData))
	p.coverage <- ggplot(covData, aes(x=position, y=value, group=variable, colour=group)) + geom_line(alpha=1/nGroups) + scale_y_continuous(trans=log2_trans())
	
	## Construct mean by group coverage plot
	meanCoverage <- ddply(covData, c("position", "group"), summarise, meanCov=mean(value))
	p.meanCov <- ggplot(meanCoverage, aes(x=position, y=meanCov, colour=group)) + geom_line(alpha=1/max(1, 1/2 * nGroups)) + scale_y_continuous(trans=log2_trans())
	
	## Annotation info and final plot
	if(is.null(txdb)) {
		p.transcripts <- FALSE
	} else {
		## The tryCatch is needed because not all regions overlap a transcript
		p.transcripts <- tryCatch(autoplot(txdb, which = wh, names.expr = "tx_name(gene_id)"), error = function(e) { FALSE })
	}	
	if(!is.logical(p.transcripts)) {
		result <- tracks(p.ideogram, "Coverage" = p.coverage, "Mean coverage" = p.meanCov, "Regions" = p.region, "tx_name\n(gene_id)" = p.transcripts, heights = c(2, 4, 4, 1.5, 3), xlim=wh, title=title) + ylab("") + theme_tracks_sunset()		
	} else {
		result <- tracks(p.ideogram, "Coverage" = p.coverage, "Mean coverage" = p.meanCov, "Regions" = p.region, heights = c(2, 5, 5, 2), xlim=wh, title=title) + ylab("") + theme_tracks_sunset()
	}
	return(result)	
}

## Plotting function
regionClusterPlot <- function(i, tUse="fwer") {
	## Chr specific selections
	chr <- as.character(seqnames(windows[[i]]))
	p.ideo <- p.ideos[[chr]]
	covInfo <- fullCovSmall[[chr]]
	
	## Make the plot
	p <- plotClusterCustom(windows[[i]], regions=fullRegions, titleName=names(ids)[ids == names(windows)[i]], coverageInfo=covInfo, groupInfo=groupInfo, titleUse=tUse, txdb=txdb, p.ideogram=p.ideo, forceLarge=TRUE)
	print(p)
	rm(p.ideo, covInfo)
	
	return(invisible(TRUE))	
}

## Make plots
for(i in seq_len(length(windows))) {
	regionClusterPlot(i)
}
```




# Reproducibility

Date the report was generated.

```{r reproducibility1, echo=FALSE, bootstrap.show.code=FALSE}
## Date the report was generated
Sys.time()
```

Wallclock time spent generating the report.

```{r "reproducibility2", echo=FALSE, bootstrap.show.code=FALSE}
## Processing time in seconds
totalTime <- diff(c(startTime, Sys.time()))
round(totalTime, digits=3)
```

`R` session information.

```{r "reproducibility3", echo=FALSE, bootstrap.show.code=FALSE, bootstrap.show.message=FALSE}
## Session info
options(width=120)
devtools::session_info()
```
