---
title: "Evaluate results"
author: "L Collado-Torres"
date: "`r doc_date()`"
output: 
  BiocStyle::html_document
---


Here we evaluate the performance of the different pipelines we used to analyze the simulated data. Evaluating the performance is tricky and depends on how you define the reference set. Here we show the results by using all exons and exons that overlap only 1 transcript. The second set removes ambiguity in determining the truth and is thus the most useful. Other possible reference sets can be considered and you can find the results in previous versions of this document. We removed them from the final version to avoid confusion.

# Setup

First, we load the required libraries.

```{r 'loadLibs', message = FALSE, warning = FALSE}
## Load required libraries
library('GenomicRanges')
library('TxDb.Hsapiens.UCSC.hg19.knownGene')
library('knitr')
library('devtools')
library('ggplot2')
```

## Load data

Next we load the required data from the simulation setup as well as the statistical results from the different pipelines. We then format the results in a way that we can easily access them later on.

```{r 'loadData'}
## Load simulation information
load(file.path('..', 'generateReads', 'simulation_info.Rdata'))

## Load transcripts used
txdb <- keepSeqlevels(TxDb.Hsapiens.UCSC.hg19.knownGene, 'chr17')
txinfo <- select(txdb, keys = chosen$ucsckg_id, columns = columns(txdb), keytype = 'TXNAME')

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


if(!identical(nrow(chosen), length(tx))) {
    message(paste(Sys.time(), 'Subsetting list of transcripts to those that were used when simulating the reads'))
    tx <- tx[names(tx) %in% chosen$gene_id]
    stopifnot(unique(sum(width(tx)) - chosen$width) == 0)
}

## Load stat results
inputs <- c('featureCounts', 'railMatrix', 'regionMatrix', 'ballgown')
reps <- 1:3
softwares <- c('DESeq2', 'edgeR', 'limma')
statuses <- c('complete', 'incomplete')

## Function for loading results from deseq2-edger folder
loadDE <- function(software, file) {
    load(file.path('..', 'deseq2-edger', paste0('stats-', file, '-', software, '.Rdata')))
    if(software == 'DESeq2') {
        DEres <- deseq
    } else if (software == 'edgeR') {
        DEres <- edger
    } else if (software == 'limma') {
        DEres <- limma
    }
    return(DEres)
}

## Function for loading ballgown results
loadBG <- function(level, file) {
    load(file.path('..', 'ballgown', paste0('bgres-', level, '-', file, '.Rdata')))
    return(bgres)
}

## Actually load the stats results
stats <- lapply(inputs, function(input) {
    inputres <- lapply(reps, function(rep) {
        if(input %in% c('featureCounts', 'ballgown')) {
            statusres <- lapply(statuses, function(status) {
                if(input == 'ballgown') {
                    res <- lapply(c('trans', 'exon'), loadBG, file = 
                        paste0('R', rep, ifelse(status == 'complete', '-comp', '-inc')))
                    names(res) <- c('trans', 'exon')
                } else if (input == 'featureCounts') {
                    res <- lapply(softwares, loadDE, file = paste0('featureCounts-R', rep, 
                        ifelse(status == 'complete', '-comp', '-inc')))
                    names(res) <- softwares
                }
                return(res)      
            })
            names(statusres) <- statuses
            return(statusres)
        } else {
            softres <- lapply(softwares, loadDE, file = paste0(input, '-R', rep))
            names(softres) <- softwares
            return(softres)
        }
    })
    names(inputres) <- reps
    return(inputres)
})
names(stats) <- inputs

## Group the stats results into a list without nested levels
stats_GR <- lapply(inputs, function(input) {
    if(input %in% c('featureCounts', 'ballgown')) {
        res <- lapply(stats[[input]], unlist)
    } else {
        res <- unlist(stats[[input]])
    }
    res <- unlist(res)
    return(res)
})
names(stats_GR) <- inputs
stats_GR <- unlist(stats_GR)
```

## Evaluation functions

The following code defines the functions we will use to evaluate the statistical results.

```{r 'evalFuns'}
## count_comp compares the information at hand versus the reference set
count_comp <- function(info, rep = 1, type = 'padj', reference, cut = 0.05) {
    if(type == 'padj') {
        idx <- mcols(info)$padj < cut
    } else if (type == 'qval') {
        idx <- mcols(info)$qval < cut
    } else if (type == 'pbonf') {
        idx <- mcols(info)$pbonf < cut
    }
    
    idx[is.na(idx)] <- FALSE

    ## Overlaps at least 1 DE 'reference' unit
    addmargins(table('DE truth' = mcols(reference)[[paste0('DEr', rep)]],
        'Overlaps DE' = countOverlaps(reference, info[idx]) > 0))
}


## Functions for evaluating empirical power, FPR and FDR, plus summarizing 
## the info
emp_power <- function(m, digits = 1) {
    round(m[2, 2] / m[2, 3] * 100, digits)
}
emp_fpr <- function(m, digits = 1) {
    round(m[1, 2] / m[1, 3] * 100, digits)
}
emp_fdr <- function(m, digits = 1) {
    round(m[1, 2] / m[3, 2] * 100, digits)
}

## Detailed table with the results for each replicate of the simulation
emp <- function(tables) {
    empirical <- data.frame(
        power = sapply(tables, emp_power),
        FPR = sapply(tables, emp_fpr),
        FDR = sapply(tables, emp_fdr)
    )
    empirical$replicate <- as.integer(sapply(strsplit(names(tables), '\\.'), '[[', 2))
    empirical$AnnotationComplete <- NA
    empirical$AnnotationComplete[grepl('\\.complete', names(tables))] <- TRUE
    empirical$AnnotationComplete[grepl('\\.incomplete', names(tables))] <- FALSE
    empirical$Aligner <- 'HISAT'
    empirical$Aligner[grepl('rail', names(tables))] <- 'Rail-RNA'
    empirical$SummaryMethod <- 'StringTie'
    empirical$SummaryMethod[grepl('Matrix', names(tables))] <- 'derfinder'
    empirical$SummaryMethod[grepl('featureCounts', names(tables))] <- 'featureCounts'
    empirical$StatMethod <- 'DESeq2'
    empirical$StatMethod[grepl('edgeR', names(tables))] <- 'edgeR'
    empirical$StatMethod[grepl('limma', names(tables))] <- 'limma'
    empirical$StatMethod[grepl('ballgown.*exon', names(tables))] <- 'ballgown-exon'
    empirical$StatMethod[grepl('ballgown.*trans', names(tables))] <- 'ballgown-trans'
    rownames(empirical) <- NULL
    return(empirical)
}

## This function takes the result from emp() and summarizes it by showing
## the minimum and maximum value per replicate
emp_sum_paste <- function(x) { paste0('(', x[1], '-', x[2], ')') }
emp_sum <- function(empinfo) {
    empinfo$situation <- paste(empinfo$AnnotationComplete, empinfo$Aligner, 
        empinfo$SummaryMethod, empinfo$StatMethod, sep=';')
    empsum <- data.frame(
        'Power' = sapply(tapply(empinfo$power, empinfo$situation, range),
            emp_sum_paste),
        'FPR' = sapply(tapply(empinfo$FPR, empinfo$situation, range),
            emp_sum_paste),
        'FDR' = sapply(tapply(empinfo$FDR, empinfo$situation, range),
            emp_sum_paste)
    )
    empsum$AnnotationComplete <- as.logical(sapply(strsplit(rownames(empsum), ';'), '[[', 1))
    empsum$Aligner <- sapply(strsplit(rownames(empsum), ';'), '[[', 2)
    empsum$SummaryMethod <- sapply(strsplit(rownames(empsum), ';'), '[[', 3)
    empsum$StatMethod <- sapply(strsplit(rownames(empsum), ';'), '[[', 4)
    empsum <- empsum[c(1, 12, 2, 13, 3, 14, 4, 15, 5, 16, 6, 9, 7, 10, 8, 11), ]
    rownames(empsum) <- NULL
    empsum 
}

## Summarize by taking the mean
emp_sum_mean <- function(empinfo) {
    empinfo$situation <- paste(empinfo$AnnotationComplete, empinfo$Aligner, 
        empinfo$SummaryMethod, empinfo$StatMethod, sep=';')
    empsum <- data.frame(
        'Power' = tapply(empinfo$power, empinfo$situation, mean, na.rm = TRUE),
        'FPR' = tapply(empinfo$FPR, empinfo$situation, mean, na.rm = TRUE),
        'FDR' = tapply(empinfo$FDR, empinfo$situation, mean, na.rm = TRUE)
    )
    empsum$AnnotationComplete <- as.logical(sapply(strsplit(rownames(empsum), ';'), '[[', 1))
    empsum$Aligner <- sapply(strsplit(rownames(empsum), ';'), '[[', 2)
    empsum$SummaryMethod <- sapply(strsplit(rownames(empsum), ';'), '[[', 3)
    empsum$StatMethod <- sapply(strsplit(rownames(empsum), ';'), '[[', 4)
    empsum <- empsum[c(1, 12, 2, 13, 3, 14, 4, 15, 5, 16, 6, 9, 7, 10, 8, 11), ]
    rownames(empsum) <- NULL
    empsum 
}

## index_comp is similar to count_comp but it returns the actual indices for
## the false positives, true positives, etc
index_comp <- function(info, rep = 1, type = 'padj', reference, cut = 0.05) {
    if(type == 'padj') {
        idx <- mcols(info)$padj < cut
    } else if (type == 'qval') {
        idx <- mcols(info)$qval < cut
    } else if (type == 'pbonf') {
        idx <- mcols(info)$pbonf < cut
    }
    
    idx[is.na(idx)] <- FALSE

    ## Overlaps at least 1 DE 'reference' unit
    ov <- countOverlaps(reference, info[idx]) > 0
    TP <- mcols(reference)[[paste0('DEr', rep)]] & ov
    TN <- !mcols(reference)[[paste0('DEr', rep)]] & !ov
    FP <- !mcols(reference)[[paste0('DEr', rep)]] & ov
    FN <- mcols(reference)[[paste0('DEr', rep)]] & !ov
    return(list(TruePositive = TP, TrueNegative = TN, FalsePositive = FP, FalseNegative = FN))
}

## case_result uses the information from index_comp to make a nice summary table
case_result <- function(idx, r, reference) { 
    res <- lapply(names(idx), function(i) {
        cases <- mcols(reference)[[paste0('case', r)]][idx[[i]]]
        if(length(cases) == 0) return(NULL)
        data.frame(case = cases, result = i, stringsAsFactors = FALSE)
    })
    res <- do.call(rbind, res)
    addmargins(table("Case" = res$case, "Result" = res$result))
}

## Chose the most frequent case, used for choosing the case in the exons 
## reference sets.exon
mostFreq <- function(cases) {
    if(length(cases) == 1) return(cases)
    names(sort(table(cases), decreasing = TRUE))[1]
}

## Process exons set to create a reference set
processExons <- function(exons_set) {
    ov <- findOverlaps(exons_set, tx)
    ## Find DE status per exon
    sHits <- subjectHits(ov)
    qHits <- queryHits(ov)
    exons_set$DEr1 <- as.vector(tapply(sHits, qHits, function(y) {
        ## Note that we could change this from any() to all() and that
        ## would impact change whether the exon is labeled as DE or not.
        any(chosen$rep1[y] != 'normal')
    }))
    exons_set$DEr2 <- as.vector(tapply(sHits, qHits, function(y) {
        any(chosen$rep2[y] != 'normal')
    }))
    exons_set$DEr3 <- as.vector(tapply(sHits, qHits, function(y) {
        any(chosen$rep3[y] != 'normal')
    }))

    exons_set$cases1 <- as.vector(tapply(sHits, qHits, function(y) {
        chosen$case1[y]
    }))
    exons_set$cases2 <- as.vector(tapply(sHits, qHits, function(y) {
        chosen$case2[y]
    }))
    exons_set$cases3 <- as.vector(tapply(sHits, qHits, function(y) {
        chosen$case3[y]
    }))
    exons_set$case1 <- unname(sapply(exons_set$cases1, mostFreq))
    exons_set$case2 <- unname(sapply(exons_set$cases2, mostFreq))
    exons_set$case3 <- unname(sapply(exons_set$cases3, mostFreq))
    return(exons_set)
}
```

## Reference sets

We can now proceed to defining our reference sets. They are:

1. Exon level with DE status given by whether the exon overlaps any transcript with low or high expression. Note that if two transcripts share the same exon, that exon will appear twice in this reference set.
1. Exons that overlap only 1 transcript with DE status given by the transcript they overlap. These exons are not ambiguous in the truth assignment and is thus the most reliable reference set.

```{r 'defineRefs'}
## Transcript level info
#trans <- tx
#mcols(trans)$DEr1 <- chosen$rep1 != 'normal'
#mcols(trans)$DEr2 <- chosen$rep2 != 'normal'
#mcols(trans)$DEr3 <- chosen$rep3 != 'normal'
#mcols(trans)$case1 <- chosen$case1
#mcols(trans)$case2 <- chosen$case2
#mcols(trans)$case3 <- chosen$case3

## Transcript level information with DE status given by the gene level
#trans_case <- trans
#mcols(trans_case)$DEr1 <- chosen$case1 %in% c('allDE', 'singleDE', 'someDE')
#mcols(trans_case)$DEr2 <- chosen$case2 %in% c('allDE', 'singleDE', 'someDE')
#mcols(trans_case)$DEr3 <- chosen$case3 %in% c('allDE', 'singleDE', 'someDE')


## Create exon level reference
exons <- processExons(unlist(tx))

## Create disjoint exon level reference
#exons_disjoin <- processExons(disjoin(unlist(tx)))

## Exons overlapping only 1 transcript
exons_one <- exons[countOverlaps(exons, tx) == 1]

## Disjoint exons overlapping only 1 transcript
#exons_disjoin_one <- exons_disjoin[countOverlaps(exons_disjoin, trans) == 1]


## Explore the number and percent of DE cases
percDE <- function(reference) {
    res <- lapply(1:3, function(r) {
        round(table(mcols(reference)[[paste0('DEr', r)]]) / length(reference) * 100, 2)
    })
    names(res) <- paste0('rep', 1:3)
    unlist(res)
}

## First at the transcript level, which matches exactly the 2/6 of transcripts 
## set to be DE (1/6 low, 1/6 high)
#percDE(trans)

## Next at the transcript level with DE status given by the gene level
## The percent DE is highly increased given that genes with more than 1
## transcript are likely to have at least 1 DE.
#percDE(trans_case)

## At the exon level, the percent with the truth set as DE is similar than a
## the transcript (DE by gene) level and is highly increased.
percDE(exons)

## Exons overlapping only 1 transcript. The percent is close to the one used in
## the simulation setup. That is: 2/6 of transcripts set to be DE (1/6 low, 1/6 
## high)
percDE(exons_one)


## For each reference set, see how many other units of that set each piece overlaps
table(countOverlaps(exons) - 1)
table(countOverlaps(exons_one) - 1)
```

Note how the overlap between the different units of the same reference set is considerable for the complete exon reference set.


Number of units per reference level:

* Exon level: `r length(exons)`
* Exon level subsetted to those that overlap only 1 transcript: `r length(exons_one)`


# Evaluation (FDR controlled)


## Evaluate results

The following code actually runs the evaluation functions for the different reference sets.

```{r 'runEval'}
## Get type and replicate info
types <- ifelse(grepl('ballgown', names(stats_GR)), 'qval', 'padj') ## For FDR 
#types <- 'pbonf' ## For FWER
replicates <- as.integer(sapply(strsplit(names(stats_GR), '\\.'), '[[', 2))


## Evaluate at the exons level
## DE if any of the transcripts it overlaps is DE
tables_exons <- mapply(count_comp, stats_GR, replicates, types,
    SIMPLIFY = FALSE, MoreArgs = list(reference = exons, cut = 0.05))
empirical_exons <- emp(tables_exons)
index_exons <- mapply(index_comp, stats_GR, replicates, types,
    SIMPLIFY = FALSE, MoreArgs = list(reference = exons, cut = 0.05))
case_result_exons <- mapply(case_result, index_exons, replicates,
    SIMPLIFY = FALSE, MoreArgs = list(reference = exons))
    
## Evaluate at the exons level (overlapping only 1 transcript)
## DE if any of the transcripts it overlaps is DE
tables_exons_one <- mapply(count_comp, stats_GR, replicates, types,
    SIMPLIFY = FALSE, MoreArgs = list(reference = exons_one, cut = 0.05))
empirical_exons_one <- emp(tables_exons_one)
index_exons_one <- mapply(index_comp, stats_GR, replicates, types,
    SIMPLIFY = FALSE, MoreArgs = list(reference = exons_one, cut = 0.05))
case_result_exons_one <- mapply(case_result, index_exons_one, replicates,
    SIMPLIFY = FALSE, MoreArgs = list(reference = exons_one))
```

## Summaries

Below are summaries of showing the minimum and maximum empirical power, false positive rate (FPR) and false discovery rate (FDR) per replicate for each of the analysis pipelines by the different reference sets.

```{r 'displayEmpiricalSummaries', results = 'asis'}
## Exon level
kable(emp_sum(empirical_exons), format = 'html')
## Exon level (overlapping only 1 transcript)
kable(emp_sum(empirical_exons_one), format = 'html')
```

From these summaries we can see that in most cases the minimum and maximum values are nearly the same, with more variability for the `HISAT` -> `StringTie` -> `ballgown` pipeline results. Note that for this pipeline we conducted tests at the transcript level (StatMethod = `ballgown-trans`) and at the exon level (`ballgown-exon`). `edgeR` achieves slightly greater empirical power than `DESeq2` at the cost of a higher empirical FPR and FDR.

The sets where we consider exons overlapping only 1 transcript shows a clear difference in power for the methods that rely on annotation between the analyses that used the complete annotation and those that used the incomplete one. The `derfinder` analyses do have higher FPR and FDR than other pipelines, but achieve nearly the same power as the other methods when the annotation is complete.


## Details

The following tables show the results for all the replicates.

```{r 'displayEmpirical', results = 'asis'}
## Exon level
kable(empirical_exons, format = 'html', row.names = TRUE)
## Exon level (overlapping only 1 transcript)
kable(empirical_exons_one, format = 'html', row.names = TRUE)
```


## Specific case

In this section we focus on the specific case of using `HISAT` -> `derfinder` -> `edgeR` from the first replicate because it has the highest empirical FPR from all the `derfinder` analyses when using the exon reference set (row 20 in the first table).


When using the exon level reference set, only about a tenth [(`r case_result_exons[['regionMatrix.1.edgeR']][2, 1]` + `r case_result_exons[['regionMatrix.1.edgeR']][2, 2]`) / `r case_result_exons[['regionMatrix.1.edgeR']][2, 5]` = `r round(100 * (case_result_exons[['regionMatrix.1.edgeR']][2, 1] + case_result_exons[['regionMatrix.1.edgeR']][2, 2]) / case_result_exons[['regionMatrix.1.edgeR']][2, 5], 2)`%] of the `noneDE` cases are incorrectly labeled. The numbers decrease in the subset of exons that overlap only 1 transcript: [(`r case_result_exons_one[['regionMatrix.1.edgeR']][2, 1]` + `r case_result_exons_one[['regionMatrix.1.edgeR']][2, 2]`) / `r case_result_exons_one[['regionMatrix.1.edgeR']][2, 5]` = `r round(100 * (case_result_exons_one[['regionMatrix.1.edgeR']][2, 1] + case_result_exons_one[['regionMatrix.1.edgeR']][2, 2]) / case_result_exons_one[['regionMatrix.1.edgeR']][2, 5], 2)`%]. 

```{r 'regionMatrix.1.edgeR'}
## Compare all references for one simulation scenario
case_result_exons[['regionMatrix.1.edgeR']]
case_result_exons_one[['regionMatrix.1.edgeR']]
```

Next we show the same tables but in percent of the total reference set units.

```{r 'regionMatrix.1.edgeR.percent'}
## Now in percent
round(case_result_exons[['regionMatrix.1.edgeR']] /
    case_result_exons[['regionMatrix.1.edgeR']][6, 5] * 100, 2)
round(case_result_exons_one[['regionMatrix.1.edgeR']] /
    case_result_exons_one[['regionMatrix.1.edgeR']][6, 5] * 100, 2)
```

In both the exons reference set, most of the false negatives are from the `someDE` scenario. Those are ambiguous cases and it's the main reason why the power is decreased when using this reference set. Note how the percent of total reference units from the `someDE` case that result in false negatives decreases from `r round(case_result_exons[['regionMatrix.1.edgeR']] / case_result_exons[['regionMatrix.1.edgeR']][6, 5] * 100, 2)[5, 1]` to `r round(case_result_exons_one[['regionMatrix.1.edgeR']] / case_result_exons_one[['regionMatrix.1.edgeR']][6, 5] * 100, 2)[5, 1]` between the exons set that does not take into account multiplicity and the subset of exons overlapping only 1 transcript.

```{r 'saveResults', echo = FALSE}
## Calculate summaries
empirical_exons_sum <- emp_sum(empirical_exons)
empirical_exons_one_sum <- emp_sum(empirical_exons_one)

## Save results
save(tables_exons, empirical_exons, empirical_exons_sum, case_result_exons, file = 'results-exons.Rdata')
save(tables_exons_one, empirical_exons_one, empirical_exons_one_sum, case_result_exons_one, file = 'results-exons_one.Rdata')
```


## At different cutoffs

The following plots shows the observed FDR and observed FPR for controlling the FDR at 0.01, 0.05, 0.1, 0.15 and 0.2.

```{r curve-fdr}
## Evaluate at the exons level (overlapping only 1 transcript)
## DE if any of the transcripts it overlaps is DE at different cutoffs
cuts <- c(0.01, 0.05, 0.1, 0.15, 0.2)
to_char <- function(x) { x[is.na(x)] <- 'NA'; return(x) }
emp_exons_one_cuts <- do.call(rbind, lapply(cuts, function(cut) {
    res <- emp_sum_mean(emp(mapply(count_comp, stats_GR, replicates, types,
        SIMPLIFY = FALSE, MoreArgs = list(reference = exons_one, cut = cut))))
    res$cut <- cut
    return(res)
}))
assign_cluster <- function(sum_met, stat_met, ann) {
    if(sum_met == 'featureCounts') {
        res <- ifelse(ann, 'fC-complete', 'fC-incomplete')
    } else if(sum_met == 'derfinder') {
        res <- 'derfinder'
    } else if(sum_met == 'StringTie') {
        res <- ifelse(ann, 'ST-complete', 'ST-incomplete')
    }
    return(res)
}
emp_exons_one_cuts$cluster <- with(emp_exons_one_cuts, mapply(assign_cluster, SummaryMethod, StatMethod, AnnotationComplete))

## Use only HISAT-aligned for simplicity
g1 <- ggplot(data = subset(emp_exons_one_cuts, Aligner == 'HISAT'),
    aes(x = FDR, y = Power, shape = StatMethod, color = cluster)) +
    geom_point(size = 3) + geom_line() + ylab('Empirical power') +
    xlab('Observed FDR (in percent)') +
    theme_linedraw(base_size = 16) +
    scale_color_brewer(palette = 'Set1', name = 'Group') +
    scale_shape_discrete(name = 'Statistical\nmethod')
g1


## Use only HISAT-aligned for simplicity
g2 <- ggplot(data = subset(emp_exons_one_cuts, Aligner == 'HISAT'),
    aes(x = FPR, y = Power, shape = StatMethod, color = cluster)) +
    geom_point(size = 3) + geom_line() + ylab('Empirical power') +
    xlab('Observed FPR (in percent)') +
    theme_linedraw(base_size = 16) +
    scale_color_brewer(palette = 'Set1', name = 'Group') +
    scale_shape_discrete(name = 'Statistical\nmethod')
g2

pdf('curve_controlFDR.pdf', width = 9)
g1
g2
dev.off()

#case_res_exons_one_cuts <- lapply(cuts, function(cut) {
#    i <- mapply(index_comp, stats_GR, replicates, types,
#        SIMPLIFY = FALSE, MoreArgs = list(reference = exons_one, cut = cut))
#    res <- mapply(case_result, i, replicates,
#        SIMPLIFY = FALSE, MoreArgs = list(reference = exons_one))
#    return(res)
#})
```

The following table shows the actual numbers used in the previous two plots. The empirical power, FPR and FDR are the mean values across the 3 simulation replicates.

```{r 'display-cuts-FDR', results = 'asis'}
## Exon level (overlapping only 1 transcript), showing mean by cutoff
kable(emp_exons_one_cuts, format = 'html', row.names = TRUE)
```




# Evaluation (FWER controlled)


## Evaluate results

The following code actually runs the evaluation functions for the different reference sets.

```{r 'runEval-fwer', eval = FALSE}
## Get type and replicate info
types <- 'pbonf' ## For FWER
replicates <- as.integer(sapply(strsplit(names(stats_GR), '\\.'), '[[', 2))


## Evaluate at the exons level
## DE if any of the transcripts it overlaps is DE
tables_exons <- mapply(count_comp, stats_GR, replicates, types,
    SIMPLIFY = FALSE, MoreArgs = list(reference = exons, cut = 0.05))
empirical_exons <- emp(tables_exons)
index_exons <- mapply(index_comp, stats_GR, replicates, types,
    SIMPLIFY = FALSE, MoreArgs = list(reference = exons, cut = 0.05))
case_result_exons <- mapply(case_result, index_exons, replicates,
    SIMPLIFY = FALSE, MoreArgs = list(reference = exons))
    
## Evaluate at the exons level (overlapping only 1 transcript)
## DE if any of the transcripts it overlaps is DE
tables_exons_one <- mapply(count_comp, stats_GR, replicates, types,
    SIMPLIFY = FALSE, MoreArgs = list(reference = exons_one, cut = 0.05))
empirical_exons_one <- emp(tables_exons_one)
index_exons_one <- mapply(index_comp, stats_GR, replicates, types,
    SIMPLIFY = FALSE, MoreArgs = list(reference = exons_one, cut = 0.05))
case_result_exons_one <- mapply(case_result, index_exons_one, replicates,
    SIMPLIFY = FALSE, MoreArgs = list(reference = exons_one))
```

## Summaries

Below are summaries of showing the minimum and maximum empirical power, false positive rate (FPR) and false discovery rate (FDR) per replicate for each of the analysis pipelines by the different reference sets.

```{r 'displayEmpiricalSummaries-fwer', results = 'asis', eval = FALSE}
## Exon level
kable(emp_sum(empirical_exons), format = 'html')
## Exon level (overlapping only 1 transcript)
kable(emp_sum(empirical_exons_one), format = 'html')
```


## Details

The following tables show the results for all the replicates.

```{r 'displayEmpirical-fwer', results = 'asis', eval = FALSE}
## Exon level
kable(empirical_exons, format = 'html', row.names = TRUE)
## Exon level (overlapping only 1 transcript)
kable(empirical_exons_one, format = 'html', row.names = TRUE)
```

```{r 'saveResults-fwer', echo = FALSE}
## Calculate summaries
empirical_exons_sum <- emp_sum(empirical_exons)
empirical_exons_one_sum <- emp_sum(empirical_exons_one)

## Save results
save(tables_exons, empirical_exons, empirical_exons_sum, case_result_exons, file = 'results-exons-FWER.Rdata')
save(tables_exons_one, empirical_exons_one, empirical_exons_one_sum, case_result_exons_one, file = 'results-exons_one-FWER.Rdata')
```


## At different cutoffs

The following plots shows the observed FDR and observed FPR for controlling the FWER at 0.01, 0.05, 0.1, 0.15 and 0.2.

```{r curve-fwer, eval = FALSE}
## Evaluate at the exons level (overlapping only 1 transcript)
## DE if any of the transcripts it overlaps is DE at different cutoffs
emp_exons_one_cuts <- do.call(rbind, lapply(cuts, function(cut) {
    res <- emp_sum_mean(emp(mapply(count_comp, stats_GR, replicates, types,
        SIMPLIFY = FALSE, MoreArgs = list(reference = exons_one, cut = cut))))
    res$cut <- cut
    return(res)
}))

## Use only HISAT-aligned for simplicity
g3 <- ggplot(data = subset(emp_exons_one_cuts, Aligner == 'HISAT'),
    aes(x = FDR, y = Power, shape = StatMethod, color = cluster)) +
    geom_point(size = 3) + geom_line() + ylab('Empirical power') +
    xlab('Observed FDR (in percent)') +
    theme_linedraw(base_size = 16) +
    scale_color_brewer(palette = 'Set1', name = 'Group') +
    scale_shape_discrete(name = 'Statistical\nmethod')
g3

## Use only HISAT-aligned for simplicity
g4 <- ggplot(data = subset(emp_exons_one_cuts, Aligner == 'HISAT'),
    aes(x = FPR, y = Power, shape = StatMethod, color = cluster)) +
    geom_point(size = 3) + geom_line() + ylab('Empirical power') +
    xlab('Observed FPR (in percent)') +
    theme_linedraw(base_size = 16) +
    scale_color_brewer(palette = 'Set1', name = 'Group') +
    scale_shape_discrete(name = 'Statistical\nmethod')
g4

pdf('curve_controlFWER.pdf', width = 9)
g3
g4
dev.off()
```

The following table shows the actual numbers used in the previous two plots. The empirical power, FPR and FDR are the mean values across the 3 simulation replicates.

```{r 'display-cuts-FWER', results = 'asis', eval = FALSE}
## Exon level (overlapping only 1 transcript), showing mean by cutoff
kable(emp_exons_one_cuts, format = 'html', row.names = TRUE)
```



# Reproducibility

```{r 'reproducibility'}
## Reproducibility info
Sys.time()
proc.time()
options(width = 120)
session_info()
```
