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

All exon counting methods
=========================

This document is based on [counts.based.html](../counts-based/counts.based.html) and expands the ideas there to using 5 different exon sets. These are exons from:

* the genomic state object produced by `derfinder`: _genomicState_
* `exons()`: _txdb_
* `exonsBy(by = 'gene')`: _byGene_
* `disjointExons()`: _disjoint_
* `featureCounts()` as included in the `Rsubread` package: _featureCounts_.

Not all exonic sets are disjoint and we wanted to know how they would perform using the simulated data.

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

## Load data
load("../coverageToExon/covToEx-ucsc.Rdata")
load("../CoverageInfo/fullCov.Rdata")
load("../derAnalysis/run2-v1.0.10/groupInfo.Rdata")
load('../simulation_info.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')
}

## TxDb setup
library('GenomeInfoDb')
library('TxDb.Hsapiens.UCSC.hg19.knownGene')
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))]

## Calculate coverage
calc_cov <- function(exons) {
    print(system.time(reg_cov <- getRegionCoverage(fullCov = fullCov, exons, verbose = FALSE)))
    cov <- t(sapply(reg_cov, function(x) { colSums(x) })) / 100
    round(cov, 0)
}

## DESeq2 analysis
run_deseq <- function(counts, exons, file, groupInfo) {
    nonzero <- sapply(rowSums(counts), function(x) {x > 0})
    
    ## 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 = paste0(file, '-DESeq2.Rdata'))
    
    ## End
    return(deseq)
}

## edgeR analysis
run_edger <- function(counts, exons, file, groupInfo) {
    nonzero <- sapply(rowSums(counts), function(x) {x > 0})
    
    ## 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 = paste0(file, '-edgeR.Rdata'))
    
    ## End
    return(edger)
}
```


# Exon analysis

## Sets

The following code defines the five exonic sets. It then checks how many exons overlap other exons. If it's zero for all of them, then they are disjoint.

```{r defineSets}
## Genomic state
exons <- GenomicState.Hsapiens.UCSC.hg19.knownGene$fullGenome
exons <- exons[seqnames(exons) == 'chr22']
exons <- list('genomicState' = exons[exons$theRegion == 'exon'])

## TxDb
exons <- c(exons, list(txdb = exons(txdb)))

## TxDb by gene
exons <- c(exons, list(byGene = unlist(exonsBy(txdb, 'gene'))))

## TxDb disjoint
exons <- c(exons, list(disjoint = disjointExons(txdb)))

## feature counts
hg19 <- read.table(system.file("annot", "hg19_RefSeq_exon.txt", package = "Rsubread"), header = TRUE)
hg19 <- subset(hg19, Chr == 'chr22')
exons <- c(exons, list(featureCounts = GRanges(hg19$Chr, IRanges(hg19$Start, hg19$End), strand = hg19$Strand, GeneID = hg19$GeneID)))
rm(hg19)

## Note that in some exon sets, some exons overlap
lapply(exons, function(x) { table(countOverlaps(x) - 1)})
```

As expected, only the _genomicState_ and _disjoint_ sets are fully disjoint since they were designed to be that way.



## Counting

The following code counts the number of reads per exon using tools from `derfinder`. Alternatively, other software could be used for counting.

```{r count}
## Count using derfinder
counts <- lapply(exons, calc_cov)
```

## DE

The following code performs the differential expression analysis using `DESeq2` and `edgeR`.

```{r runDE}
## DESeq2
system.time( deseq <- mapply(run_deseq, counts, exons, names(exons), MoreArgs = list(groupInfo = groupInfo)) )

## edgeR
system.time( edger <- mapply(run_edger, counts, exons, names(exons), MoreArgs = list(groupInfo = groupInfo)) )
```


## Agreement

The code below compares the results between `DESeq2` and `edgeR`.

```{r edgeR-DESeq2-agreement}
agree <- function(deseq, edger) {
    addmargins(table('Significant DE gene -- DESeq2' = mcols(deseq)$sig, 'Significant DE gene -- edgeR' = mcols(edger)$sig))
}
mapply(agree, deseq, edger, SIMPLIFY = FALSE)
```

Overall, `DESeq2` and `edgeR` mostly agree regardless of the exon set. Their best agreement is in the _disjoint_ set with only 5 disagreements (`r round(5/283 * 100, 1)` percent).

## Compare

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

The following code compares the exonic segments (280 of them) against the resulting exon DE calls from `DESeq2` and `edgeR` (in that order). 

```{r 'segments_comp', 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))
}

lapply(deseq, count_comp)
lapply(edger, count_comp)
```

#### Empirical Power

The empirical power is shown below for `DESeq2` and `edgeR`.

```{r empPower}
emp_power <- function(info, ptype = 'padj') {
    m <- count_comp(info, ptype)
    round(m[2, 2] / m[2, 3] * 100, 2)
}


sapply(deseq, emp_power)
sapply(edger, emp_power)
```


#### Empirical FPR

The empirical False Positive Rate (FPR) is shown below for `DESeq2` and `edgeR`.

```{r empFPR}
emp_fpr <- function(info, ptype = 'padj') {
    m <- count_comp(info, ptype)
    round(m[1, 2] / m[1, 3] * 100, 2)
}

sapply(deseq, emp_fpr)
sapply(edger, emp_fpr)
```

#### Empirical FDR

The empirical False Discovery Rate (FDR) is shown below for `DESeq2` and `edgeR`.

```{r empFDR}
emp_fdr <- function(info, ptype = 'padj') {
    m <- count_comp(info, ptype)
    round(m[1, 2] / m[3, 2] * 100, 2)
}

sapply(deseq, emp_fdr)
sapply(edger, emp_fdr)
```



### Overlap

As in [counts.based.html](../counts-based/counts.based.html) we can compare the DERs and exons directly. 

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

## Some formatting and subsets
fullRegions$significantFDR <- factor(fullRegions$qvalues < 0.05, levels = c('TRUE', 'FALSE'))
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 (FWER)' = 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 (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$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)]
}
```

#### Query: DERs

First we use the DERs as the query and check if they overlap a significantly differentially expressed exon. The results are shown for `DESeq2` and `edgeR` using all DERs and then requiring a minimum overlap of 20 bp.

```{r 'queryDERs'}
## DESeq2
lapply(deseq, function(x) {
    ov_table(fullRegions, x)
})
lapply(deseq, function(x) {
    ov_table(fullRegs20, x, minov = 20L)
})

## edgeR
lapply(edger, function(x) {
    ov_table(fullRegions, x)
})
lapply(edger, function(x) {
    ov_table(fullRegs20, x, minov = 20L)
})
```

#### Query: exons

Next, we can use the exons as the query and check if they overlap a significant DER. The results are shown for `DESeq2` and `edgeR` using all DERs and then restricting the overlap to minimum 20 bp.

```{r 'queryExons'}
## DESeq2
lapply(deseq, function(x) {
    ov_table(fullRegions, x, 'counts')
})
lapply(deseq, function(x) {
    ov_table(fullRegs20, x, 'counts', minov = 20L)
})

## edgeR
lapply(edger, function(x) {
    ov_table(fullRegions, x, 'counts')
})
lapply(edger, function(x) {
    ov_table(fullRegs20, x, 'counts', minov = 20L)
})
```

# Conclusions

Using `DESeq2` with the _genomicState_ or _txdb_ exonic sets resulted in the highest power and reasonable FPR. While the _disjoint_ set analyzed by `DESeq2` resulted in the best FPR it had the lowest empirical power.

When using the DERs as the query, the _disjoint_ set has the lowest disagreement for both `DESeq2` and `edgeR`. However, the _genomicState_, _txdb_, and _byGene_ sets are the only ones that do not have cases of a significant DER overlapping a non significant DE exon.

Regardless of the exon set, when using the exons as query, there are no cases of a non significantly DE exon overlap a significant DER when requiring a 20 bp overlap.

# Reproducibility

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