---
output:
  knitrBootstrap::bootstrap_document:
    theme.chooser: TRUE
    highlight.chooser: TRUE
  html_document:
    toc: true
    theme: united
---

Compare vs count-based methods
==============================


This report compares the original implementation of `derfinder` available at [alyssafrazee/derfinder](https://github.com/alyssafrazee/derfinder), `DESeq2`, `edgeR-robust` and our implementation of `derfinder` using the exonic segments described in [evaluate.html](../evaluate/evaluate.html).

# 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.


This first code chunk loads the necessary data.


```{r 'setup', bootstrap.show.code = FALSE, bootstrap.show.message = FALSE}
library('edgeR')
library('DESeq2')
library('GenomicRanges')

## Load data
load("../coverageToExon/covToEx-ucsc.Rdata")
load("../derAnalysis/run2-v1.0.10/groupInfo.Rdata")
if(file.exists("../derAnalysis/run2-v1.0.10/colsubset.Rdat")) {
    load("../derAnalysis/run2-v1.0.10/colsubset.Rdata")
} else {
    colsubset <- seq_len(length(groupInfo))
}

## 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')
```


## 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')
```


# Comparison

```{r 'setup-comp', bootstrap.show.code = FALSE, bootstrap.show.message = FALSE}
library('TxDb.Hsapiens.UCSC.hg19.knownGene')
#library('knitr')
library('derfinder')
library('derfinderHelper')
library('derfinderPlot')
library('qvalue')
#library('bumphunter')
load('../simulation_info.Rdata')
load('../derAnalysis/run2-v1.0.10/fullRegions.Rdata')
load('../derAnalysis/run2-v1.0.10/models.Rdata')
load('../derAnalysis/run2-v1.0.10/chr22/optionsStats.Rdata')
load('../CoverageInfo/fullCov.Rdata')
names(fullRegions) <- seq_len(length(fullRegions))
```

## Exonic segments

Just as in [evaluate.html](../evaluate/evaluate.html) we can compare the results against the exonic segments. From that document:

Next we can evaluate the simulation by classifying the exonic segments as whether they should be DE or not. Then, we can find if the DERs overlap such segments.

The following code is a subset of [evaluate.html](../evaluate/evaluate.html) and generates the exonic segments.


```{r 'segments', bootstrap.show.code = FALSE}
## Find exons
txdb <- keepSeqlevels(TxDb.Hsapiens.UCSC.hg19.knownGene, 'chr22')
txinfo <- select(txdb, keys = chosen$ucsckg_id, columns = columns(txdb), keytype = 'TXNAME')

## Buiild GRangesList with exons grouped by transcript
tx <- split(GRanges(seqnames = txinfo$EXONCHROM, IRanges(start = txinfo$EXONSTART, end = txinfo$EXONEND), strand = txinfo$EXONSTRAND), txinfo$TXNAME)
tx <- tx[match(chosen$ucsckg_id, names(tx))]

## Gene level: DE if at least one transcript is DE
gene <- data.frame(gene_id = unique(chosen$gene_id))
gene$DE <- sapply(gene$gene_id, function(x) { any(chosen$DE[chosen$gene_id == x])  })
gene$case <- sapply(gene$gene_id, function(x) { unique(chosen$case[chosen$gene_id == x])  })

## Identify exonic segments
segments <- GRangesList(lapply(gene$gene_id, function(x) {
    i <- chosen$ucsckg_id[ chosen$gene_id == x]
    
    ## Find segments
    segs <- disjoin(unlist(tx[i]))
    ov <- findOverlaps(segs, tx[i])
    
    ## Find DE status per segment
    segs$DE <- as.vector(tapply(subjectHits(ov), queryHits(ov), function(y) {
        any(chosen$DE[ chosen$gene_id == x])
    }))
    
    ## Finish
    return(segs)
}))
names(segments) <- gene$gene_id
segs <- unlist(segments)
```

## derfinder-original

The original derfinder implementation does not support multi-group comparisons. So we performed all pair comparisons for the three groups. We can extract those that are differentially expressed (states 3 and 4) with a FDR adjusted p-value $<$ 0.05.

```{r 'original', bootstrap.show.code = FALSE}
load('../deranalysis-original/chr22/regions-merged-chr22-AB.Rdata')
original <- list('AB' = regions.merged)
load('../deranalysis-original/chr22/regions-merged-chr22-AC.Rdata')
original <- c(original, list('AC' = regions.merged))
load('../deranalysis-original/chr22/regions-merged-chr22-BC.Rdata')
original <- c(original, list('BC' = regions.merged))
original <- lapply(original, function(x) { 
    GRanges(seqnames = x$chr, ranges = IRanges(x$start, x$end), strand = '*', state = x$state, mean.t = x$mean.t, mean.fold.change = x$mean.fold.change)
})

## Load p-value info
original <- lapply(names(original), function(x) {
    load(paste0('../deranalysis-original/chr22/pvals-chr22-', x, '.Rdata'))
    original[[x]]$pvalue <- pvals
    return(original[[x]])
})


## Identify DE regions
ori.de <- lapply(original, function(x) {
    state.idx <- x$state > 2
    x[state.idx & !is.na(x$pvalue)]
})

## Adjust p-values for FDR and select significant ones
ori.de.sig <- do.call(c, ori.de)
ori.de.sig$qvalue <- qvalue(ori.de.sig$pvalue)$qvalue
ori.de.sig <- ori.de.sig[ori.de.sig$qvalue < 0.05]
```

Due to the simulation setup, only one of the three groups is differentially expressed at a time. One option is to compare the exonic segments versus the resulting DERs and ask if the exonic segments overlap at least one DER as shown below.



```{r 'only1', bootstrap.show.code = FALSE}
## Overlaps at least 1 DER
addmargins(table('DE status' = segs$DE, 'Overlaps DER (sig) -- original, min 1' = countOverlaps(segs, ori.de.sig) > 0))
```

The results from the new implementation are shown below (extracted from [evaluate.html](../evaluate/evaluate.html)). The original implementation has 15 false positive (versus 0 in the new implemenation) and less false negatives (11 versus 30).

```{r 'newFDR', bootstrap.show.code = FALSE}
## Compare against new derfinder:
## Check result with FDR sig 
fullRegions$significantFDR <- factor(fullRegions$qvalues < 0.05, levels = c('TRUE', 'FALSE'))
fdr <-  fullRegions[fullRegions$significantFDR == 'TRUE']
addmargins(table('DE status' = segs$DE, 'Overlaps DER (sig FDR) -- new' = countOverlaps(segs, fdr) > 0))
```

However, because of the simulation setup, with the original implementation we would actually expect at least two DERs to overlap each exonic segment that was set to be DE. When doing this comparison, we get 3 false positives and 28 false negatives versus 0 and 30 in the new implementation.



```{r 'only2', bootstrap.show.code = FALSE}
## Identical results? to derfinder-original if we ask segments to overlap 
## two DERs given that we are doing 3 group comparisons with 1 group
## higher/lower than the other two
addmargins(table('DE status' = segs$DE, 'Overlaps DER (sig) -- original, min 2' = countOverlaps(segs, ori.de.sig) > 1))

identical(
    addmargins(table('DE status' = segs$DE, 'Overlaps DER (sig)' = countOverlaps(segs, fdr) > 0)),
    addmargins(table('DE status' = segs$DE, 'Overlaps DER (sig)' = countOverlaps(segs, ori.de.sig) > 1))
)
```

While the global numbers are similar, the agreement at the exonic segment level is not perfect between the original and new implementations as shown below. 

```{r 'newVsOnly2', bootstrap.show.code = FALSE}
addmargins(table('new' = countOverlaps(segs, fdr) > 0, 'original -- min 2' = countOverlaps(segs, ori.de.sig) > 1))
```

The following code was adapted from [evaluate.html](../evaluate/evaluate.html) and explores the false negatives when requiring at least 2 DERs to overlap each exonic segment. As with the new implementation, the original implementation struggles in genes where one of the transcripts was set to be DE while the other wasn't.



```{r 'seg.fn'}
## Explore false negative segments using sig DERs
seg.fn <- which(segs$DE & !countOverlaps(segs, ori.de.sig) > 1)

## Some segments are short
summary(width(segs[seg.fn]))
## new:
##
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     5.0    73.5   119.0   242.7   241.0  1844.0
chosen[chosen$gene_id == '6523', ]

## 6 / 28 of the segmeents are from gene with id 6523
## new is 9 / 30 for the same gene
tail(sort(table(names(segs[seg.fn]))))

## Cases of the genes with at least one FN segment
table(tapply(subset(chosen, gene_id %in% names(seg.fn))$case, subset(chosen, gene_id %in% names(seg.fn))$gene_id, unique))
## new:
## 
##   bothDE    oneDE singleDE 
##        1       10        5

## Type of gene where the segments come from. Mostly oneDE genes
table(sapply(names(segs[seg.fn]), function(x) { unique(chosen$case[chosen$gene_id == x]) }))

## new:
## 
##   bothDE    oneDE singleDE 
##        1       24        5
```


## Counts-based

Just as a verification step, the exonic segments used are disjoint as well as the exons used for the counts-based methods. Also, each exonic segment overlaps only one exon.

```{r 'verify'}
max(countOverlaps(segs)) - 1
max(countOverlaps(exons)) - 1
table(countOverlaps(segs, exons))
```

### DESeq2


Next, we can compare the exonic segments versus the significant DE exons as determined by `DESeq2`. Significance can be determined by regular p-values, FDR FDR adjusted p-values (BH method), and FWER adjusted p-values (Holm method). In all cases we use a cutoff of $<$ 0.05.

```{r 'comp-fun', bootstrap.show.code = FALSE}
count_comp <- function(info, ptype = 'padj', cut = 0.05) {
    if(ptype == 'padj') {
        idx <- info$padj < cut
    } else if (ptype == 'pvalue') {
        idx <- info$pvalue < cut
    } else {
        p <- p.adjust(info$pvalue, ptype)
        idx <- p < cut
    }
    idx[is.na(idx)] <- FALSE

    ## Overlaps at least 1 DE exon
    addmargins(table('DE status' = segs$DE, 'Overlaps DE exon' = countOverlaps(segs, info[idx]) > 0))
}
```


```{r 'deseq-comp'}
## Regular p-values
count_comp(deseq, ptype = 'pvalue')

## FDR adjusted p-values by method BH
count_comp(deseq)

## FWER adjusted p-values by method Holm
count_comp(deseq, ptype = 'holm')
identical(count_comp(deseq, ptype = 'holm'), count_comp(deseq, ptype = 'bonferroni'))
```

As expected the number of false positives is the highest with regular p-values and the lowest with FWER adjusted p-values. In all cases, there are fewer false negatives compared to `derfinder`.

### edgeR-robust

Similar results are shown below when using `edgeR-robust`.

```{r 'edger-comp'}
## Regular p-values
count_comp(edger, ptype = 'pvalue')

## FDR adjusted p-values by method BH
count_comp(edger)

## FWER adjusted p-values by method Holm
count_comp(edger, ptype = 'holm')
identical(count_comp(edger, ptype = 'holm'), count_comp(edger, ptype = 'bonferroni'))
```



## Overlap

```{r 'ov-comp-setup', bootstrap.show.code = FALSE}
## Some formatting and subsets
fullRegions$sigFDR <- as.logical(fullRegions$significantFDR)
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 (FDR)' = ders$sigFDR, 'Overlaps sig DE exon' = countOverlaps(ders, counts[counts$sig]) > 0))
        } else {
            res <- addmargins(table(ders$sigFDR, countOverlaps(ders, counts[counts$sig], minoverlap = minov) > 0, dnn = c('Significant DER (FDR)', 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 (FDR)' = countOverlaps(counts, ders[ders$sigFDR]) > 0))
        } else {
            res <- addmargins(table(counts$sig[width(counts) >= minov], countOverlaps(counts[width(counts) >= minov], ders[ders$sigFDR], minoverlap = minov) > 0, dnn = c('Significant DE exon', paste0('Overlaps sig DER (FDR, 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$sigFDR)
    } else if (case == 'TRUE-FALSE') {
        i <- which(!countOverlaps(ders, counts[counts$sig], minoverlap = minov) > 0 & ders$sigFDR)
    } 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$sigFDR], minoverlap = minov) > 0 & !counts$sig)
    } else if (case == 'TRUE-FALSE') {
        i <- which(!countOverlaps(counts, ders[ders$sigFDR], 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$sigFDR], minoverlap = minov)),
            width_exon = summary(width(counts[i])),
            exons_per_der_table = table(table(subjectHits(findOverlaps(counts[i], ders[ders$sigFDR], minoverlap = minov)))),
            exons_per_der = sort(table(subjectHits(findOverlaps(counts[i], ders[ders$sigFDR], 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$sigFDR], ignore.strand = TRUE))$distance),
            distance_nearest = distanceToNearest(counts[i], ders, ignore.strand = TRUE),
            distance_nearest_sig = distanceToNearest(counts[i], ders[ders$sigFDR], 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])
```

Most of the non-significant DERs that overlap a significant exon are shorter than 20bp with the longest being 37bp long. When restricting the analysis to 20bp or longer DERs, 6 the DERs overlap a single exon. It is an exon from a gene with two isoforms where only one of them was set to be DE. As can be seen in the figure below, the F-statistics curve oscillates and a lower cutoff could have resulted in significant DERs in this exon.

```{r 'vis-long-exon', fig.width = 10, fig.height = 7, bootstrap.show.code = FALSE}
## Identify long exon
long_exon <- deseq[which(deseq$sig)[as.integer(names(explore_ov(fullRegs20, deseq, minov = 20L)$ders_per_exon)[explore_ov(fullRegs20, deseq, minov = 20L)$ders_per_exon == 6])]]

## What was it set to be?
mcols(segs[subjectHits(findOverlaps(long_exon, segs))])$DE
gene[gene$gene_id %in% unique(names(segs[subjectHits(findOverlaps(long_exon, segs))])),]

## Calculate F-stats
range <- start(long_exon):end(long_exon)
dat <- fullCov$chr22[range, ]

## Log2 transform
for(i in 1:30) 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(long_exon, GenomicState.Hsapiens.UCSC.hg19.knownGene$fullGenome, verbose = FALSE)
## Remove symbol name because it gets chomped on the plot
mcols(annoReg$annotationList[[1]])$symbol <- NA

## Plot long exon
long_exon_plot <- function() {
    plotRegionCoverage(long_exon, getRegionCoverage(fullCov, long_exon, verbose = FALSE), groupInfo, data.frame(name = "sig DE exon via DESeq2, no sig DERs", distance = NA, region = 'SLC5A1'), annoReg, verbose = FALSE, ask = FALSE, txdb = txdb)

    ## 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), yaxt = 'n')
    axis(2, at = c(0, 10, 20, 30), c(0, 10, 20, NA), las = 2, tick = TRUE)
    abline(h = optionsStats$cutoffFstatUsed, col = 'red')
    abline(h = 0, col = 'grey')
}

## Save in pdf
pdf(file = 'long_exon.pdf', width = 10, height = 7)
long_exon_plot()
dev.off()
#system('open long_exon.pdf')

## Render in png
long_exon_plot()
```



#### 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])
```

We can further explore the cases where significant DE exons did not overlap a significant DER and compare against the exonic segments.


```{r 'ov-comp-deseq-counts-segs', bootstrap.show.code = FALSE}
## Further explore 
i <- explore_ov_counts(fullRegions, deseq, 'TRUE-FALSE')[['i']]

## Exclude cases where an exon overlaps exonic segments from both strands
i.simple <- i[countOverlaps(deseq[i], segs, ignore.strand = TRUE) == 1]

## Number of cases excluded because of strand overlaps
length(i) - length(i.simple)

## Check with exonic segments
table(mcols(segs)$DE[subjectHits(findOverlaps(deseq[i.simple], segs ))])

## Gene cases that were set to be DE
table(gene$case[gene$gene_id %in% unique(names(segs)[intersect(subjectHits(findOverlaps(deseq[i.simple], segs )), which(mcols(segs)$DE))])])

## Exon mistmatched cases (DE), what was the gene status?
table(sapply(names(segs)[intersect(subjectHits(findOverlaps(deseq[i.simple], segs )), which(mcols(segs)$DE))], function(x) { gene$case[gene$gene_id == x] }))

## Gene cases that were not set to be DE
table(gene$case[gene$gene_id %in% unique(names(segs)[intersect(subjectHits(findOverlaps(deseq[i.simple], segs )), which(!mcols(segs)$DE))])])

## Exon mistmatched cases (not DE), what was the gene status?
table(sapply(names(segs)[intersect(subjectHits(findOverlaps(deseq[i.simple], segs )), which(!mcols(segs)$DE))], function(x) { gene$case[gene$gene_id == x] }))
```

Out of the 22 cases in disagreement (min overlap 0 bp), 4 of them are due to overlapping exons from both strands. 15 of these exons are correctly detected by `DESeq2` whereas 3 are not. The ones incorrectly identified by `derfinder` to not be DE are mostly (12 / 15) from genes with two isoforms with only one of them DE.


```{r 'vis-long-deseq2-error', fig.width = 10, fig.height = 7, bootstrap.show.code = FALSE}
## Identify longest exon with the error
err_deseq2 <- segs[intersect(subjectHits(findOverlaps(deseq[i.simple], segs )), which(!mcols(segs)$DE))]
err_deseq2 <- err_deseq2[which.max(width(err_deseq2))]

## What was it set to be?
mcols(err_deseq2)$DE

## Calculate F-stats
range <- start(err_deseq2):end(err_deseq2)
dat <- fullCov$chr22[range, ]

## Log2 transform
for(i in 1:30) 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(err_deseq2, GenomicState.Hsapiens.UCSC.hg19.knownGene$fullGenome, verbose = FALSE)
## Remove symbol name because it gets chomped on the plot
mcols(annoReg$annotationList[[1]])$symbol <- NA

## Plot long exon
fp_plot <- function() {
    plotRegionCoverage(err_deseq2, getRegionCoverage(fullCov, err_deseq2, verbose = FALSE), groupInfo, data.frame(name = "False positive via DESeq2/edgeR-robust", distance = NA, region = 'SF3A1'), annoReg, verbose = FALSE, ask = FALSE, txdb = txdb)

    ## 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), yaxt = 'n')
    axis(2, at = c(0, 5, 10), c(0, 5, 10), las = 2, tick = TRUE)
    abline(h = optionsStats$cutoffFstatUsed, col = 'red')
    abline(h = 0, col = 'grey')
}

## Save in pdf
pdf(file = 'fp_deseq2.pdf', width = 10, height = 7)
fp_plot()
dev.off()
#system('open fp_deseq2.pdf')

## Render in png
fp_plot()
```





### 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])
```

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])
```

As before, we can further explore the mismatched cases where `edgeR-robust` finds DE signal but `derfinder` does not.

```{r 'ov-comp-edger-counts-segs', bootstrap.show.code = FALSE}
## Further explore 
j <- explore_ov_counts(fullRegions, edger, 'TRUE-FALSE')[['i']]

## Exclude cases where an exon overlaps exonic segments from both strands
j.simple <- j[countOverlaps(edger[j], segs, ignore.strand = TRUE) == 1]

## Number of cases excluded because of strand overlaps
length(j) - length(j.simple)

## Check with exonic segments
table(mcols(segs)$DE[subjectHits(findOverlaps(edger[j.simple], segs ))])

## Gene cases that were set to be DE
table(gene$case[gene$gene_id %in% unique(names(segs)[intersect(subjectHits(findOverlaps(edger[j.simple], segs )), which(mcols(segs)$DE))])])

## Exon mistmatched cases (DE), what was the gene status?
table(sapply(names(segs)[intersect(subjectHits(findOverlaps(edger[j.simple], segs )), which(mcols(segs)$DE))], function(x) { gene$case[gene$gene_id == x] }))

## Gene cases that were not set to be DE
table(gene$case[gene$gene_id %in% unique(names(segs)[intersect(subjectHits(findOverlaps(edger[j.simple], segs )), which(!mcols(segs)$DE))])])

## Exon mistmatched cases (not DE), what was the gene status?
table(sapply(names(segs)[intersect(subjectHits(findOverlaps(edger[j.simple], segs )), which(!mcols(segs)$DE))], function(x) { gene$case[gene$gene_id == x] }))
```

### overall

The results are fairly similar to those from using `DESeq2`. As shown below, they agree on all but 6 of the exons analyzed. In particular, the longest false positive exon by `DESeq2` is also a FP with `edgeR-robust`.

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

## Find FP errors
err_edger <- segs[intersect(subjectHits(findOverlaps(edger[j.simple], segs )), which(!mcols(segs)$DE))]

## Does it overlap FP from DESeq2? Yes
countOverlaps(err_deseq2, err_edger) == 1
```

The following is a coverage plot of the longest exon where all methods agree that its differentially expressed.

```{r 'vis-agree-de', fig.width = 10, fig.height = 7, bootstrap.show.code = FALSE}
## Find longest exon where edgeR and DEseq2 agree that it's DE
l <- which.max(width(edger)[which(edger$sig & deseq$sig)])
agree <- edger[which(edger$sig & deseq$sig)[l]]

## Width of exon
width(agree)

## DERs this exon overlaps
width(fullRegions[subjectHits(findOverlaps(agree, fullRegions))])
mcols(fullRegions[subjectHits(findOverlaps(agree, fullRegions))])$sigFDR

## Case
gene[gene$gene_id == unique(names(segs[subjectHits(findOverlaps(agree,segs))])), ]

## Calculate F-stats
range <- start(agree):end(agree)
dat <- fullCov$chr22[range, ]

## Log2 transform
for(i in 1:30) 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(agree, GenomicState.Hsapiens.UCSC.hg19.knownGene$fullGenome, verbose = FALSE)
## Remove symbol name because it gets chomped on the plot
mcols(annoReg$annotationList[[1]])$symbol <- NA

## Plot long exon
agree_plot <- function() {
    plotRegionCoverage(agree, getRegionCoverage(fullCov, agree, verbose = FALSE), groupInfo, data.frame(name = "Full agreement: differentially expressed", distance = NA, region = 'CECR6'), annoReg, verbose = FALSE, ask = FALSE, txdb = txdb)

    ## 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), yaxt = 'n')
    axis(2, at = c(0, 50, 100, 150, 200, 250), c(0, NA, 100, NA, 200, NA), las = 2, tick = TRUE)
    abline(h = optionsStats$cutoffFstatUsed, col = 'red')
    abline(h = 0, col = 'grey')
}

## Save in pdf
pdf(file = 'agree.pdf', width = 10, height = 7)
agree_plot()
dev.off()
#system('open agree.pdf')

## Render in png
agree_plot()
```

The following is a coverage plot of the longest exon where all methods agree that its not differentially expressed. For `derfinder`, we are requiring DERs to be at least 20 bp long.

```{r 'vis-agree-not-de', fig.width = 10, fig.height = 7, bootstrap.show.code = FALSE}
## Find longest exon where edgeR and DEseq2 agree that it's not DE
## Exclude the one from chr12
l <- which(width(edger)[which(!edger$sig & !deseq$sig)] == 5789)
agree_notde <- edger[which(!edger$sig & !deseq$sig)[l]]

## Width of exon
width(agree_notde)

## DERs this exon overlaps
width(fullRegions[subjectHits(findOverlaps(agree_notde, fullRegions))])
mcols(fullRegions[subjectHits(findOverlaps(agree, fullRegions))])$sigFDR

## Case
gene[gene$gene_id == unique(names(segs[subjectHits(findOverlaps(agree_notde, segs))])), ]

## Calculate F-stats
range <- start(agree_notde):end(agree_notde)
dat <- fullCov$chr22[range, ]

## Log2 transform
for(i in 1:30) 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(agree_notde, GenomicState.Hsapiens.UCSC.hg19.knownGene$fullGenome, verbose = FALSE)
## Remove symbol name because it gets chomped on the plot
mcols(annoReg$annotationList[[1]])$symbol <- NA

## Plot long exon
agree_notde_plot <- function() {
    plotRegionCoverage(agree_notde, getRegionCoverage(fullCov, agree_notde, verbose = FALSE), groupInfo, data.frame(name = "Full agreement: not differentially expressed", distance = NA, region = 'TCF20'), annoReg, verbose = FALSE, ask = FALSE, txdb = txdb)

    ## 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), yaxt = 'n')
    axis(2, at = c(0, 5, 10), c(0, 5, 10), las = 2, tick = TRUE)
    abline(h = optionsStats$cutoffFstatUsed, col = 'red')
    abline(h = 0, col = 'grey')
}

## Save in pdf
pdf(file = 'agree_notde.pdf', width = 10, height = 7)
agree_notde_plot()
dev.off()
#system('open agree_notde.pdf')

## Render in png
agree_notde_plot()
```




# Conclusions

Globally, the results between the original implementation and our implementation of `derfinder` are very similar, with both resulting in the same empirical power and absence of false positives. Both struggle in the scenario where one of the two isoforms of a gene was set to be differentially expressed.


```{r 'empirical', bootstrap.show.code = FALSE}
emp_power <- function(info, ptype) {
    m <- count_comp(info, ptype)
    round(m[2, 2] / m[2, 3] * 100, 2)
}
emp_fpr <- function(info, ptype) {
    m <- count_comp(info, ptype)
    round(m[1, 2] / m[1, 3] * 100, 2)
}
emp_fdr <- function(info, ptype) {
    m <- count_comp(info, ptype)
    round(m[1, 2] / m[3, 2] * 100, 2)
}
```

When adjusting p-values to control the FDR, both `DESeq2` and `edgeR-robust` result in higher empirical power than `derfinder`: `r emp_power(deseq, 'padj')` and `r emp_power(edger, 'padj')` respectively versus `r round(100 * 139 / 169, 2)`. However, both `DESeq2` and `edgeR-robust` have non-zero false positive rates: `r emp_fpr(deseq, 'padj')` and `r emp_fpr(edger, 'padj')` respectively. Their false discovery rate is: `r emp_fdr(deseq, 'padj')` and `r emp_fdr(edger, 'padj')` respectively. Note that the number of false positives is reduced to 1 in both cases when controlling the FWER.

We knew from [evaluate.html](../evaluate/evaluate.html) that `derfinder` struggled in the scenario where a two transcript gene had only one of them set to be differentially expressed. `DESeq2` and `edgeR-robust` are able to correctly identify such exons as differentially expressed in most cases, but they also introduce some occasional false positives whose coverage plots do not reveal apparent differences between the groups.


# Reproducibility

```{r 'reproduc', bootstrap.show.code = FALSE}
## Reproducibility info
Sys.time()
proc.time()
options(width = 120)
devtools::session_info()
```
