# export
write.table(dxr1, 'exon_results.txt', sep = '\t', row.names = F)

# explore results
sig <- dxr1[!is.na(dxr1$padj),]
sig <- sig[sig$padj <= 0.05,]
sig_df <- arrange(as.data.frame(sig), log2fold_shCtr_sh2)
write.table(sig_df, 'exon_results_sig.txt', sep = '\t', row.names = F)
DEXSeqHTML(dxr1)

# plots
plotDEXSeq(dxr1, "PLAUR", legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2)
plotDEXSeq(dxr1, "PLAUR", displayTranscripts=TRUE, legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2)
plotDEXSeq(dxr1, "PLAUR", expression=FALSE, norCounts=TRUE, legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2)

# re-orient exons and compute p values (first exon versus last exon)
x <- data.frame(dxr1)
s <- split(x, x$groupID)

output <- NULL
for (elmt in s) {
  line <- NULL
  if (elmt$genomicData.strand[1] == '-') {
    line <- c(elmt$groupID[1],
              elmt$featureID[length(elmt$featureID)],
              elmt$padj[length(elmt$padj)],
              elmt$featureID[1],
              elmt$padj[1])
  }
  else {
    line <- c(elmt$groupID[1],
              elmt$featureID[1],
              elmt$padj[1],
              elmt$featureID[length(elmt$featureID)],
              elmt$padj[length(elmt$padj)])
  }
  output <- rbind(output, line)
}

o <- data.frame(output, row.names = NULL, stringsAsFactors = F)
colnames(o) <- c('Group', 'FirstExon', 'pval_first', 'LastExon', 'pval_last')
o$pval_first <- as.numeric(o$pval_first)
o$pval_last <- as.numeric(o$pval_last)

p_df <- NULL
for (p in seq(0,1,0.01)) {
  p_df <- rbind(p_df, c(as.numeric(table(o[!is.na(o[,3]),][,3]<p)[2]) / nrow(o[!is.na(o[,3]),]),
                        as.numeric(table(o[!is.na(o[,5]),][,5]<p)[2]) / nrow(o[!is.na(o[,5]),])))
}

plot(p_df, type = 'l')
lines(seq(0,1,0.01), seq(0,1,0.01))