library('GenomicRanges')

## Loading required package: methods
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## 
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## 
## The following object is masked from 'package:stats':
## 
##     xtabs
## 
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, as.vector, cbind,
##     colnames, do.call, duplicated, eval, evalq, Filter, Find, get,
##     intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rep.int, rownames, sapply, setdiff, sort,
##     table, tapply, union, unique, unlist, unsplit
## 
## Loading required package: S4Vectors
## Loading required package: stats4
## Loading required package: IRanges
## Loading required package: GenomeInfoDb

library('data.table')
library('rCharts')
library('knitr')
library('ggplot2')
library('devtools')
opts_chunk$set(bootstrap.show.code = FALSE, dev = 'CairoPNG')

## Check that data was loaded
stopifnot(all(c('fullRegions', 'regionMat', 'analysisPath') %in% ls()))

## Fix DERs
ders <- fullRegions
names(ders) <- NULL

## Fix region matrices
regs <- unlist(GRangesList(lapply(regionMat, '[[', 'regions')))
names(regs) <- NULL

## Assign seqlengths
data(hg19Ideogram, package = 'biovizBase', envir = environment())
seqlengths(ders) <- seqlengths(hg19Ideogram)[names(seqlengths(ders))]
seqlengths(regs) <- seqlengths(hg19Ideogram)[names(seqlengths(regs))]

## Sort
ders <- sort(ders)
regs <- sort(regs)

## Check names match
identical(seqlengths(regs), seqlengths(ders))

## [1] TRUE

## Save the raw data
save(ders, file = file.path(analysisPath, 'dersOriginal.Rdata'))
save(regs, file = file.path(analysisPath, 'regsOriginal.Rdata'))

## Filter out those with less than 3bp
c('# ders under 3bp' = sum(width(ders) < 3), '# regs under 3bp' = sum(width(regs) < 3))

## # ders under 3bp # regs under 3bp 
##            16078            36576

ders <- ders[width(ders) >= 3]
regs <- regs[width(regs) >= 3]

## Save data used
save(ders, file = file.path(analysisPath, 'ders.Rdata'))
save(regs, file = file.path(analysisPath, 'regs.Rdata'))

## Construct logical Rle indexes for bases with some region
build_index <- function(gr) {
    res <- lapply(names(seqlengths(gr)), function(chr) {
        chr.len <- seqlengths(gr)[chr]
        ir <- sort(ranges(gr[seqnames(gr) == chr]))
        log <- c(rep(c(FALSE, TRUE), length(ir)), FALSE)
    
        starts <- ends <- rep(NA, length(ir) * 2)
        i <- rep(c(TRUE, FALSE), length(ir))
        starts[i] <- start(ir)
        ends[i] <- end(ir)
    
        starts[!i] <- ends[i] + 1
    
        if(max(ends, na.rm = TRUE) < chr.len) {
            ends[!i] <- c(starts[i] - 1, chr.len)[-1]
        } else {
            ends[!i] <- c(starts[i] - 1, NULL)[-1]
            starts <- starts[- length(starts)]
            log <- log[- length(log)]
        }
   
        if(starts[1] != 1) {
            ends <- c(starts[1] - 1, ends)
            starts <- c(1, starts)
        } else {
            log <- log[-1]
        }
   
        widths <- mapply(function(s, e) { e - s + 1}, starts, ends)
    
        Rle(log, widths)
    })
    names(res) <- names(seqlengths(gr))
    return(res)
}
index.ders <- build_index(ders)
index.regs <- build_index(regs)

library('epivizr')
mgr <- startEpiviz()
ders_dev <- mgr$addDevice(ders[!as.logical(ders$significantFWER)], "DERs no sig FWER")
ders_sig_dev <- mgr$addDevice(ders[as.logical(ders$significantFWER)], "DERs sig FWER")
regs_dev <- mgr$addDevice(regs, "Region Matrix")

## SOX11
mgr$navigate("chr2", 5810000, 5850000)

## MEX3A
mgr$navigate("chr1", 156040000, 156090000)

## VASH2
mgr$navigate("chr1", 213120000, 213170000)

## TG:
mgr$navigate("chr8", 134040000, 134120000)

## IGF2BP2
mgr$navigate("chr3", 185350000, 185410000)

## FBN3
mgr$navigate("chr19", 8130000, 8180000)

## End
mgr$stopServer()

## Number of regions
c('ders' = length(ders), 'regs' = length(regs))

##   ders   regs 
##  15900 180397

## Size of regions
c('ders' = summary(width(ders)), 'regs' = summary(width(regs)))

##    ders.Min. ders.1st Qu.  ders.Median    ders.Mean ders.3rd Qu. 
##         3.00         4.00         7.00        12.34        14.00 
##    ders.Max.    regs.Min. regs.1st Qu.  regs.Median    regs.Mean 
##       182.00         3.00        12.00        36.00        56.87 
## regs.3rd Qu.    regs.Max. 
##        70.00      5572.00

## Merge all the indexes
index.all <- mapply(function(der, reg) {
    both <- der & reg
    only.der <- der & !reg
    only.reg <- !der & reg
    none <- !der & !reg
    
    res <- list('both' = both, 'only.der' = only.der, 'only.reg' = only.reg,
        'none' = none, 'all.der' = der, 'all.reg' = reg)
    return(list(res))
}, index.ders, index.regs)

## Find number of base-pairs in each index
index.num <- data.frame(do.call(rbind, lapply(index.all, function(x) { sapply(x, sum)})))
index.num$chrLen <- seqlengths(ders)
index.num$chr <- rownames(index.num)
rownames(index.num) <- NULL

## Print info
summary(index.num)

##       both          only.der         only.reg            none          
##  Min.   :  471   Min.   :  10.0   Min.   :  17706   Min.   : 48039170  
##  1st Qu.: 4646   1st Qu.: 663.8   1st Qu.: 300546   1st Qu.: 79841586  
##  Median : 6343   Median : 940.0   Median : 374345   Median :133877898  
##  Mean   : 7251   Mean   : 921.5   Mean   : 420243   Mean   :128558143  
##  3rd Qu.: 9428   3rd Qu.:1097.8   3rd Qu.: 534129   3rd Qu.:161639197  
##  Max.   :27552   Max.   :2347.0   Max.   :1040506   Max.   :248194137  
##     all.der         all.reg            chrLen              chr           
##  Min.   :  481   Min.   :  18177   Min.   : 48129895   Length:24         
##  1st Qu.: 5413   1st Qu.: 304750   1st Qu.: 80415720   Class :character  
##  Median : 7281   Median : 381242   Median :134429206   Mode  :character  
##  Mean   : 8172   Mean   : 427494   Mean   :128986559                     
##  3rd Qu.:10670   3rd Qu.: 541522   3rd Qu.:162132764                     
##  Max.   :28918   Max.   :1054137   Max.   :249250621

## Overall info
overallInfo <- colSums(index.num[, -ncol(index.num)])
overallInfo

##       both   only.der   only.reg       none    all.der    all.reg 
##     174014      22117   10085842 3085395439     196131   10259856 
##     chrLen 
## 3095677412

## Overall info in percent
overallInfo / sum(as.numeric(index.num$chrLen)) * 100

##         both     only.der     only.reg         none      all.der 
## 5.621193e-03 7.144478e-04 3.258040e-01 9.966786e+01 6.335641e-03 
##      all.reg       chrLen 
## 3.314252e-01 1.000000e+02

d1 <- data.table(data.frame(row = seq_len(nrow(index.num)), index.num, check.names=FALSE))
t1 <- dTable(d1, sPaginationType=  'full_numbers', iDisplayLength=25,
    sScrollX='100%')
t1$print("bases", cdn=TRUE)

## Find number of segments in each index
index.seg <- data.frame(do.call(rbind, lapply(index.all, function(x) { 
    sapply(x, function(y) {
        sum(runValue(y))
    })
})))
index.seg$chr <- rownames(index.seg)
rownames(index.seg) <- NULL

## Print info
summary(index.seg)

##       both           only.der        only.reg          none      
##  Min.   :  26.0   Min.   :  2.0   Min.   :  391   Min.   :  372  
##  1st Qu.: 357.5   1st Qu.:113.2   1st Qu.: 5697   1st Qu.: 5486  
##  Median : 508.0   Median :136.5   Median : 7558   Median : 7276  
##  Mean   : 546.9   Mean   :143.9   Mean   : 7978   Mean   : 7627  
##  3rd Qu.: 694.2   3rd Qu.:171.2   3rd Qu.:10105   3rd Qu.: 9745  
##  Max.   :1726.0   Max.   :372.0   Max.   :19596   Max.   :18908  
##     all.der          all.reg          chr           
##  Min.   :  28.0   Min.   :  369   Length:24         
##  1st Qu.: 456.0   1st Qu.: 5410   Class :character  
##  Median : 623.5   Median : 7158   Mode  :character  
##  Mean   : 662.5   Mean   : 7517                     
##  3rd Qu.: 808.8   3rd Qu.: 9615                     
##  Max.   :1885.0   Max.   :18604

## Overall info
colSums(index.seg[, -ncol(index.seg)])

##     both only.der only.reg     none  all.der  all.reg 
##    13125     3454   191474   183045    15900   180397

d2 <- data.table(data.frame(row = seq_len(nrow(index.seg)), index.seg, check.names=FALSE))
t2 <- dTable(d2, sPaginationType=  'full_numbers', iDisplayLength=25,
    sScrollX='100%')
t2$print("segments")

## Get an idea of the width of the segments in each index
index.width <- data.frame(do.call(rbind, lapply(index.all, function(x) {
    tmp <- data.frame(do.call(rbind, lapply(x, function(y) {
        summary(runLength(y)[runValue(y)])
    })), check.names = FALSE)
    tmp$index <- names(x)
    rownames(tmp) <- NULL
    return(tmp)
})), check.names = FALSE)
index.width$chr <- rep(names(seqlengths(ders)), each = 6)
rownames(index.width) <- NULL

## Print info
summary(index.width)

##       Min.          1st Qu.          Median            Mean          
##  Min.   :1.000   Min.   : 3.00   Min.   :  4.00   Min.   :     5.00  
##  1st Qu.:1.000   1st Qu.: 4.00   1st Qu.:  6.00   1st Qu.:    11.63  
##  Median :1.000   Median : 8.00   Median : 23.00   Median :    31.70  
##  Mean   :1.694   Mean   : 8.26   Mean   : 33.74   Mean   :  4198.87  
##  3rd Qu.:3.000   3rd Qu.:12.00   3rd Qu.: 36.00   3rd Qu.:    56.19  
##  Max.   :3.000   Max.   :39.75   Max.   :259.00   Max.   :159600.00  
##     3rd Qu.             Max.             index          
##  Min.   :   6.00   Min.   :       7   Length:144        
##  1st Qu.:  13.00   1st Qu.:      98   Class :character  
##  Median :  44.25   Median :     307   Mode  :character  
##  Mean   : 276.39   Mean   : 1609538                     
##  3rd Qu.:  70.00   3rd Qu.:    2746                     
##  Max.   :3379.00   Max.   :30260000                     
##      chr           
##  Length:144        
##  Class :character  
##  Mode  :character  
##                    
##                    
## 

d3 <- data.table(data.frame(row = seq_len(nrow(index.width)), index.width, check.names=FALSE))
t3 <- dTable(d3, sPaginationType=  'full_numbers', iDisplayLength=25,
    sScrollX='100%')
t3$print("widths")

ov20 <- findOverlaps(ders, regs, minoverlap = 20L)

counts <- list()
for(type in c("any", "within", "equal")) {
    ct.ders <- countOverlaps(ders, regs, minoverlap = 20L, type = type)
    plot(log10(table(ct.ders)), main = paste("DERs in regs for type", type))
    ct.regs <- countOverlaps(regs, ders, minoverlap = 20L, type = type)
    plot(log10(table(ct.regs)), main = paste("Regs in DERs for type", type))
    counts <- c(counts, list(table(ct.ders), table(ct.regs)))
}

cts <- as.integer(unique(unlist(lapply(counts, names))))

nOverlap20 <- do.call(rbind, lapply(counts, function(x) { 
    df <- data.frame(nOverlap = cts, freq = x[match(cts, names(x))],
        row.names = seq_len(length(cts)))
    df$observed <- !is.na(df$freq)
    df$freq[is.na(df$freq)] <- 0
    df$cumFreq <- cumsum(df$freq)
    df$cumPerc <- df$cumFreq / max(df$cumFreq) * 100
    return(df)
}))
nOverlap20$type <- factor(rep(c('any', 'within', 'equal'), each = length(cts) * 2), levels = c('any', 'within', 'equal'))
nOverlap20$match <- rep(rep(c('DERs-in-regs', 'regs-in-DERs'), each = length(cts)), 3)
#nOverlap <- nOverlap[complete.cases(nOverlap), ]
rownames(nOverlap20) <- NULL

nOverlap20$alpha <- ifelse(nOverlap20$observed, 1, 1/3)

## Make a nice plot
ggplot(data = nOverlap20, aes(x = nOverlap, y = cumFreq, colour = match, alpha = alpha)) + geom_point() + facet_grid( . ~ type )# + geom_smooth(se=FALSE)

## Show cumulative percents
ggplot(data = nOverlap20, aes(x = nOverlap, y = cumPerc, colour = match, linetype = match)) + geom_line(lwd=1) + facet_grid( . ~ type )

## Percent with widths < 20L
small <- c('ders' = sum(width(ders) < 20) / length(ders), 'regs' = sum(width(regs) < 20) / length(regs)) * 100
data.frame('under-20' = small, '20-and-above' = 100 - small, check.names = FALSE)

##      under-20 20-and-above
## ders 82.43396     17.56604
## regs 34.13582     65.86418

## Base level agreement
c('regs' = overallInfo['both'] / (overallInfo['both'] + overallInfo['only.reg']) * 100, 'ders' = overallInfo['both'] / (overallInfo['both'] + overallInfo['only.der']) * 100)

## regs.both ders.both 
##  1.696067 88.723353

## Overlap (min 20) agreement
c('regs' = 100 - subset(nOverlap20, match == 'regs-in-DERs' & nOverlap == 0 & type == 'any')$cumPerc, 'ders' = 100 - subset(nOverlap20, match == 'DERs-in-regs' & nOverlap == 0 & type == 'any')$cumPerc)

##      regs      ders 
##  1.139709 16.578616

ov1 <- findOverlaps(ders, regs, minoverlap = 1L)

counts <- list()
for(type in c("any", "within", "equal")) {
    ct.ders <- countOverlaps(ders, regs, minoverlap = 1L, type = type)
    plot(log10(table(ct.ders)), main = paste("DERs in regs for type", type))
    ct.regs <- countOverlaps(regs, ders, minoverlap = 1L, type = type)
    plot(log10(table(ct.regs)), main = paste("Regs in DERs for type", type))
    counts <- c(counts, list(table(ct.ders), table(ct.regs)))
}

cts <- as.integer(unique(unlist(lapply(counts, names))))

nOverlap1 <- do.call(rbind, lapply(counts, function(x) { 
    df <- data.frame(nOverlap = cts, freq = x[match(cts, names(x))], 
        row.names = seq_len(length(cts)))
    df$observed <- !is.na(df$freq)
    df$freq[is.na(df$freq)] <- 0
    df$cumFreq <- cumsum(df$freq)
    df$cumPerc <- df$cumFreq / max(df$cumFreq) * 100
    return(df)
}))
nOverlap1$type <- factor(rep(c('any', 'within', 'equal'), each = length(cts) * 2), levels = c('any', 'within', 'equal'))
nOverlap1$match <- rep(rep(c('DERs-in-regs', 'regs-in-DERs'), each = length(cts)), 3)
#nOverlap <- nOverlap[complete.cases(nOverlap), ]
rownames(nOverlap1) <- NULL

nOverlap1$alpha <- ifelse(nOverlap1$observed, 1, 1/3)

## Overlap (min 1bp) agreement
c('regs' = 100 - subset(nOverlap1, match == 'regs-in-DERs' & nOverlap == 0 & type == 'any')$cumPerc, 'ders' = 100 - subset(nOverlap1, match == 'DERs-in-regs' & nOverlap == 0 & type == 'any')$cumPerc)

##      regs      ders 
##  4.069912 82.176101

## Make a nice plot
ggplot(data = nOverlap1, aes(x = nOverlap, y = cumFreq, colour = match, alpha = alpha)) + geom_point() + facet_grid( . ~ type )# + geom_smooth(se=FALSE)

## Show cumulative percents
ggplot(data = nOverlap1, aes(x = nOverlap, y = cumPerc, colour = match, linetype = match)) + geom_line(lwd=1) + facet_grid( . ~ type )

save(index.all, index.num, index.seg, index.width, nOverlap20, ov20, nOverlap1, ov1, overallInfo, file = file.path(analysisPath, "comparison-results.Rdata"))

# Load fullRegions.Rdata and regionMat.Rdata before this step
library('rmarkdown')
library('knitrBootstrap')
render('step7-regMatVsDERs.Rmd')

## [1] "2014-11-23 18:06:46 EST"

## Session info-----------------------------------------------------------------------------

##  setting  value                                      
##  version  R version 3.1.1 Patched (2014-10-16 r66782)
##  system   x86_64, linux-gnu                          
##  ui       X11                                        
##  language (EN)                                       
##  collate  en_US.UTF-8                                
##  tz       

## Packages---------------------------------------------------------------------------------

##  package        * version date       source                                   
##  BiocGenerics   * 0.12.1  2014-11-15 Bioconductor                             
##  Cairo            1.5.6   2014-06-26 CRAN (R 3.1.0)                           
##  chron            2.3.45  2014-02-11 CRAN (R 3.1.1)                           
##  colorspace       1.2.4   2013-09-30 CRAN (R 3.1.0)                           
##  data.table     * 1.9.4   2014-10-02 CRAN (R 3.1.1)                           
##  devtools       * 1.6.1   2014-10-07 CRAN (R 3.1.1)                           
##  digest           0.6.4   2013-12-03 CRAN (R 3.1.0)                           
##  evaluate         0.5.5   2014-04-29 CRAN (R 3.1.0)                           
##  formatR          1.0     2014-08-25 CRAN (R 3.1.1)                           
##  GenomeInfoDb   * 1.2.3   2014-11-15 Bioconductor                             
##  GenomicRanges  * 1.18.3  2014-11-20 Bioconductor                             
##  ggplot2        * 1.0.0   2014-05-21 CRAN (R 3.1.0)                           
##  gtable           0.1.2   2012-12-05 CRAN (R 3.1.0)                           
##  htmltools        0.2.6   2014-09-08 CRAN (R 3.1.1)                           
##  IRanges        * 2.0.0   2014-10-15 Bioconductor                             
##  knitr          * 1.8     2014-11-11 CRAN (R 3.1.1)                           
##  knitrBootstrap * 1.0.0   2014-11-19 Github (jimhester/knitrBootstrap@76c41f0)
##  labeling         0.3     2014-08-23 CRAN (R 3.1.1)                           
##  lattice          0.20.29 2014-04-04 CRAN (R 3.1.1)                           
##  markdown         0.7.4   2014-08-24 CRAN (R 3.1.1)                           
##  MASS             7.3.35  2014-09-30 CRAN (R 3.1.1)                           
##  mime             0.2     2014-09-26 CRAN (R 3.1.1)                           
##  munsell          0.4.2   2013-07-11 CRAN (R 3.1.0)                           
##  plyr             1.8.1   2014-02-26 CRAN (R 3.1.0)                           
##  proto            0.3.10  2012-12-22 CRAN (R 3.1.0)                           
##  rCharts        * 0.4.5   2014-11-17 Github (ramnathv/rCharts@929875d)        
##  Rcpp             0.11.3  2014-09-29 CRAN (R 3.1.1)                           
##  reshape2         1.4     2014-04-23 CRAN (R 3.1.0)                           
##  RJSONIO          1.3.0   2014-07-28 CRAN (R 3.1.1)                           
##  rmarkdown      * 0.3.3   2014-09-17 CRAN (R 3.1.1)                           
##  rstudioapi       0.1     2014-03-27 CRAN (R 3.1.1)                           
##  S4Vectors      * 0.4.0   2014-10-15 Bioconductor                             
##  scales           0.2.4   2014-04-22 CRAN (R 3.1.0)                           
##  stringr          0.6.2   2012-12-06 CRAN (R 3.1.0)                           
##  whisker          0.3.2   2013-04-28 CRAN (R 3.1.0)                           
##  XVector          0.6.0   2014-10-15 Bioconductor                             
##  yaml             2.1.13  2014-06-12 CRAN (R 3.1.1)