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

Gene-level evaluation
=========================

This document is based on [all-exons.html](../all-exons/all-exons.html). It uses the gene sets defined by [counts-gene.html](../counts-gene/counts-gene.html) and runs `DESeq2` and `edgeR` analyses. The results are compared against the exonic segments first. Then, overlap comparisons between the genes and the DERs are performed as in [all-exons.html](../all-exons/all-exons.html).

```{r 'setup', bootstrap.show.code = FALSE, bootstrap.show.message = FALSE}
library('edgeR')
library('DESeq2')
library('GenomicRanges')
library('GenomeInfoDb')
library('TxDb.Hsapiens.UCSC.hg19.knownGene')

## Load data
load("../derAnalysis/run2-v1.0.10/groupInfo.Rdata")
load('../simulation_info.Rdata')
load('../counts-gene/summOv_comp.Rdata')
load('../counts-gene/summOv_inc.Rdata')
load('../counts-gene/summOv_rand.Rdata')
load('../counts-gene/gene_sets.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))
}

## DESeq2 analysis
run_deseq <- function(counts, genes, 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 <- genes[nonzero]
    mcols(deseq) <- cbind(mcols(deseq), results(dse))

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

    ## Save results
    save(deseq, file = paste0(file, '-DESeq2.Rdata'))
    
    ## End
    return(deseq)
}

## edgeR analysis
run_edger <- function(counts, genes, 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 <- genes[nonzero]
    mcols(edger) <- cbind(mcols(edger), DataFrame(lrw$table))
    mcols(edger)$pvalue <-  lrw$table$PValue
    mcols(edger)$padj <- p.adjust(lrw$table$PValue, 'BH')

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

    ## Save results
    save(edger, file = paste0(file, '-edgeR.Rdata'))
    
    ## End
    return(edger)
}
```


# Gene analysis

## Sets

Define the gene sets.

```{r defineSets}
genes <- list('complete' = gene_comp, 'incomplete' = gene_inc,
    'incomplete_random' = gene_rand)
```




## Counting

Load counting data produced by `summarizeOverlaps()`. Also checks that the names of the genes are in the correct order and identifies the genes that have 0 counts in all samples.

```{r count}
## Count from summarizeOverlaps()
counts <- list('complete' = summOv_comp, 'incomplete' = summOv_inc,
    'incomplete_random' = summOv_rand)
counts <- lapply(counts, function(x) { assay(x)[, paste0('sample', 1:30)]})

## Check names
identical(rownames(counts[[1]]), names(gene_comp))
identical(rownames(counts[[2]]), names(gene_inc))

## Which genes have 0 counts?
lapply(counts, function(x) { 
    rownames(x)[rowSums(x) == 0]
})

## How did those genes look like?
x <- gene_comp[rownames(counts[[1]])[rowSums(counts[[1]]) == 0]]
x

## Do they overlap each other when ignoring the strand?
strand(x) <- RleList('100422998' = Rle(factor("*", levels = c("+", "-", "*")), 1), '100423034' = Rle(factor("*", levels = c("+", "-", "*")), 1))

## They do
table(countOverlaps(x) - 1)
```

The two genes that have 0 counts are from different strands but overlap each other when we ignore the strand. Given the recommended settings for `summarizeOverlaps()` and the ambiguity introduced by non-strand specific RNA-seq data, we have to exclude these genes.

## DE

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

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

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

### Agreement

The following code compares the `DESeq2` and `edgeR` results.

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

There are only 2 and 1 disagreements between `DESeq2` and `edgeR` for the complete and incomplete gene sets respectively.


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

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

## 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 checks if an exonic segment overlaps a DE gene for `DESeq2` and `edgeR`, first by controlling the FDR and then by controlling the FWER.

```{r 'segments_comp'}
count_comp <- function(info, ptype = 'padj', cut = 0.05) {
    if(ptype == 'padj') {
        idx <- mcols(info)$padj < cut
    } else if (ptype == 'pvalue') {
        idx <- mcols(info)$pvalue < cut
    } else {
        p <- p.adjust(mcols(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))
}

## Default: adjusting p-values by FDR
lapply(deseq, count_comp)
lapply(edger, count_comp)

## Adjusting p-values by Holm method
lapply(deseq, count_comp, ptype = 'holm')
lapply(edger, count_comp, ptype = 'holm')
```

#### Empirical Power

Empirical power for `DESeq2` and `edgeR`. First for controlling the FDR and next for controlling the FWER.

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

## Adjusting p-values by Holm method
sapply(deseq, emp_power, ptype = 'holm')
sapply(edger, emp_power, ptype = 'holm')
```


#### Empirical FPR

Empirical false positive rate (FPR) for `DESeq2` and `edgeR`. First for controlling the FDR and next for controlling the FWER.

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

## Adjusting p-values by Holm method
sapply(deseq, emp_fpr, ptype = 'holm')
sapply(edger, emp_fpr, ptype = 'holm')
```

#### Empirical FDR

The empirical False Discovery Rate (FDR) is shown below for `DESeq2` and `edgeR`. First for controlling the FDR and next for controlling the FWER.

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

## Adjusting p-values by Holm method
sapply(deseq, emp_fdr, ptype = 'holm')
sapply(edger, emp_fdr, ptype = 'holm')
```


### Overlap

As with previous reports ([all-exons.html](../all-exons/all-exons.html) and [counts.based.html](../counts-based/counts.based.html)) we can compare DERs versus the genes by overlapping them. 

```{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 gene' = countOverlaps(ders, counts[mcols(counts)$sig]) > 0))
        } else {
            res <- addmargins(table(ders$sigFDR, countOverlaps(ders, counts[mcols(counts)$sig], minoverlap = minov) > 0, dnn = c('Significant DER (FDR)', paste0('Overlaps sig DE gene (min ', minov, 'bp)'))))
        }
    } else if (query == 'counts') {
        if(minov == 0) {
            res <- addmargins(table('Significant DE gene' = mcols(counts)$sig, 'Overlaps sig DER (FWER)' = countOverlaps(counts, ders[ders$sigFDR]) > 0))
        } else {
            res <- addmargins(table(mcols(counts)$sig[sapply(width(counts), sum) >= minov], countOverlaps(counts[sapply(width(counts), sum) >= minov], ders[ders$sigFDR], minoverlap = minov) > 0, dnn = c('Significant DE gene', paste0('Overlaps sig DER (FWER, min ', minov, 'bp)'))))
        }
    }
    return(res)
}

## Explore mistmatched cases for DERs vs genes direction
explore_ov <- function(ders, counts, case = "FALSE-TRUE", minov = 0L) {
    if(case == 'FALSE-TRUE') {
        i <- which(countOverlaps(ders, counts[mcols(counts)$sig], minoverlap = minov) > 0 & !ders$sigFDR)
    } else if (case == 'TRUE-FALSE') {
        i <- which(!countOverlaps(ders, counts[mcols(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[mcols(counts)$sig], minoverlap = minov)),
            width_der = summary(width(ders[i])),
            ders_per_gene_table = table(table(subjectHits(findOverlaps(ders[i], counts[mcols(counts)$sig], minoverlap = minov)))),
            ders_per_gene = sort(table(subjectHits(findOverlaps(ders[i], counts[mcols(counts)$sig], minoverlap = minov)))),
            i = i
        )
    } else {
        res <- list(
            width_der = summary(width(ders[i])),
            distance_nearest_sum = summary(mcols(distanceToNearest(ders[i], unlist(counts), ignore.strand = TRUE))$distance),
            distance_nearest_sig_sum = summary(mcols(distanceToNearest(ders[i], unlist(counts[mcols(counts)$sig]), ignore.strand = TRUE))$distance),
            distance_nearest = distanceToNearest(ders[i], unlist(counts), ignore.strand = TRUE),
            distance_nearest_sig = distanceToNearest(ders[i], unlist(counts[mcols(counts)$sig]), ignore.strand = TRUE),
            i = i
        )
    }
    
    return(res)
}

## Explore mistmatched cases for genes vs DERs direction
explore_ov_counts <- function(ders, counts, case = "FALSE-TRUE", minov = 0L) {
    counts <- counts[sapply(width(counts), sum) >= minov]
    if(case == 'FALSE-TRUE') {
        i <- which(countOverlaps(counts, ders[ders$sigFDR], minoverlap = minov) > 0 & !mcols(counts)$sig)
    } else if (case == 'TRUE-FALSE') {
        i <- which(!countOverlaps(counts, ders[ders$sigFDR], minoverlap = minov) > 0 & mcols(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_gene = summary(sapply(width(counts[i]), sum)),
            genes_per_der_table = table(table(subjectHits(findOverlaps(counts[i], ders[ders$sigFDR], minoverlap = minov)))),
            genes_per_der = sort(table(subjectHits(findOverlaps(counts[i], ders[ders$sigFDR], minoverlap = minov)))),
            i = i
        )
    } else {
        res <- list(
            width_gene = summary(sapply(width(counts[i]), sum)),
            distance_nearest_sum = summary(mcols(distanceToNearest(unlist(counts[i]), ders, ignore.strand = TRUE))$distance),
             distance_nearest_sig_sum = summary(mcols(distanceToNearest(unlist(counts[i]), ders[ders$sigFDR], ignore.strand = TRUE))$distance),
            distance_nearest = distanceToNearest(unlist(counts[i]), ders, ignore.strand = TRUE),
            distance_nearest_sig = distanceToNearest(unlist(counts[i]), ders[ders$sigFDR], ignore.strand = TRUE),
            i = i
        )
    }
    
    return(res)
}

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

#### Query: DERs

We first use the DERs as the query as shown below for `DESeq2` and `edgeR` using all the 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)
})
```

The results are identical between `DESeq2` and `edgeR` with surprisingly slightly better agreements in with the incomplete gene set instead of the complete one when requiring a minimum overlap of 20bp.

The common disagreement scenario is when a non significant DER overlaps a significant DE gene.


The following code explores the mismatches (min 20 bp overlap) using `DESeq2` when the DERs are the query.

```{r exploreQueryDER}
lapply(deseq, function(x) {
    explore_ov(fullRegions, x, minov = 20L)[1:3]
})

lapply(deseq, function(x) {
    explore_ov(fullRegions, x, 'TRUE-FALSE', minov = 20L)[1:3]
})
```

#### Query: genes

Next we can use the genes as the query.

```{r 'queryGenes'}
## 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)
})
```

As opposed to when DERs were used as the query, this time there is a small disagreement between `DESeq2` and `edgeR` when requiring a minimum 20bp overlap. Overlap, the most common disagreement case against the DERs is when a significant DE gene does not overlap a significant DER.


The following code explores the mismatches (min 20 bp overlap) using `DESeq2` when the genes are the query.

```{r exploreQueryGene}
## No cases FALSE-TRUE
#lapply(deseq, function(x) {
#    explore_ov_counts(fullRegions, x, minov = 20L)[1:3]
#})

lapply(deseq, function(x) {
    explore_ov_counts(fullRegions, x, 'TRUE-FALSE', minov = 20L)[1:3]
})
```



# Conclusions

`DESeq2` and `edgeR` result in the same empirical power regardless of the gene set used, although both have high empirical FPR (minimum 18.92 for `DESeq2` in complete gene set).

In the overlap comparison with a minimum overlap of 20 bp, when DERs are used as the query, the common disagreement scenario is when a non significant DER overlaps a significant DE gene. When genes are used as the query, the most common disagreement case against the DERs is when a significant DE gene does not overlap a significant DER.


# Reproducibility

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