library(DMRcaller)
library(RUnit)

# create synthetic data
step <- 25
numberOfCytosines <- 200
startPos <- seq(from = 1, by=step, length.out=numberOfCytosines)

syntheticMethylationData1 <- GRanges(
  seqnames              = factor(rep("Chr1", length(startPos))),
  ranges                = IRanges(startPos, startPos),
  strand                = rep("+", length(startPos)),
  context               = rep("CG", length(startPos)),
  readsM                = rep(0, length(startPos)),
  readsN                = rep(0, length(startPos)),
  trinucleotide_context = rep("CGH", length(startPos)))

syntheticMethylationData1$readsN <- sample(50:150, 
                                           length(syntheticMethylationData1$readsN), 
                                           replace = TRUE)
syntheticMethylationData1$readsN <- c(rep(50,50),rep(100, 100), rep(50,50))


proportion1 <- c(rep(0.1, 10), rep(0.5, 30), rep(0.2, 3), rep(0.45, 7), 
                 rep(0.2, 50), rep(0.8, 20), rep(0.3, 10), rep(0.045, 20), 
                 rep(0.3, 20), rep(0.7, 20), rep(0.3, 10))
syntheticMethylationData1$readsM <- syntheticMethylationData1$readsN*proportion1

syntheticMethylationData2 <- syntheticMethylationData1
proportion2 <- c(rep(0.045, 130), rep(0.6, 20), rep(0.045, 50))
syntheticMethylationData2$readsM <- round(proportion2*syntheticMethylationData2$readsN)


#DMRs
vals <- rep(0, length(proportion1))
vals[(syntheticMethylationData1$readsM/syntheticMethylationData1$readsN - 
        syntheticMethylationData2$readsM/syntheticMethylationData2$readsN)>0.4] <- 1
vals[(syntheticMethylationData2$readsM/syntheticMethylationData2$readsN - 
        syntheticMethylationData1$readsM/syntheticMethylationData1$readsN)>0.4] <- -1



#join the differentially methylated cytosines into regions
rle <- rle(vals)
rle$cumulative <- cumsum(rle$lengths)
endOfRuns <- rle$cumulative 

DMRs <- GRanges(
  seqnames    = "Chr1",
  ranges      = IRanges(((endOfRuns - rle$lengths + 1)*step), (endOfRuns*step)),
  strand      = "+",
  direction   = rle$values
)
DMRs <- DMRs[DMRs$direction!=0]


genes <- GRanges(
  seqnames    = "Chr1",
  ranges      = IRanges(c(300, 1800),c(1200,2100)),
  strand      = "+"
)



# check if compute coverage works
test_computeMethylationDataCoverage <- function() {
  breaks <- c(1,10,40,90,100)
  answer <- c(1,1,1,0.5,0)
  checkEqualsNumeric(computeMethylationDataCoverage(syntheticMethylationData1, 
                                                    context="CG", 
                                                    breaks=breaks), 
                     answer, tolerance=1.0e-8)
  checkEqualsNumeric(computeMethylationDataCoverage(syntheticMethylationData2, 
                                                    context="CG", 
                                                    breaks=breaks), 
                     answer, tolerance=1.0e-8)
}


# check if compute coverage works
test_computeDMRs <- function() {
  DMRs_CG_noise_filter <-  computeDMRs(syntheticMethylationData1, 
                                       syntheticMethylationData2, 
                                       NULL, 
                                       method = "noise_filter", 
                                       context = "CG", 
                                       pValueThreshold = 0.01, 
                                       minReadsPerCytosine = 3, 
                                       minProportionDifference = 0.4, 
                                       minGap = 50, 
                                       minSize = 50, 
                                       test = "score", 
                                       windowSize=50)

  checkTrue(all(overlapsAny(DMRs_CG_noise_filter, DMRs)))
  checkTrue(all(overlapsAny(DMRs, DMRs_CG_noise_filter)))
  
  DMRs_CG_bins <-  computeDMRs(syntheticMethylationData1, 
                               syntheticMethylationData2, 
                               NULL, 
                               method = "bins",
                               context = "CG", 
                               binSize = 100, 
                               test = "score", 
                               pValueThreshold = 0.01,
                               minCytosinesCount = 3, 
                               minProportionDifference = 0.4, 
                               minGap = 100, 
                               minReadsPerCytosine = 3, 
                               cores = 1)
  checkTrue(all(overlapsAny(DMRs_CG_bins, DMRs)))
  checkTrue(all(overlapsAny(DMRs, DMRs_CG_bins)))
}
 
test_computeDMRs <- function() {
}


test_filterDMRs <- function() {
  DMRs_CG_genes <-  filterDMRs(syntheticMethylationData1, 
                               syntheticMethylationData2,
                               genes, 
                               context = "CG", 
                               test = "score", 
                               pValueThreshold=0.01, 
                               minCytosinesCount = 1,
                               minProportionDifference = 0.4, 
                               minReadsPerCytosine = 3)  
  checkTrue(DMRs_CG_genes == genes[overlapsAny(genes, DMRs)])
}

test_filterDMRs <- function() {
  DMRs_CG_noise_filter <-  computeDMRs(syntheticMethylationData1, 
                                       syntheticMethylationData2, 
                                       NULL, 
                                       method = "noise_filter",
                                       context="CG", 
                                       pValueThreshold = 0.01, 
                                       minReadsPerCytosine = 3, 
                                       minProportionDifference = 0.4, 
                                       minGap = 50, 
                                       minSize = 50,
                                       test = "score", 
                                       windowSize=50)
  
  DMRs_CG_noise_filter_filtered <- filterDMRs(syntheticMethylationData1,
                                              syntheticMethylationData2, 
                                              DMRs_CG_noise_filter, 
                                              context = "CG", 
                                              minCytosinesCount = 10, 
                                              minProportionDifference = 0.55,
                                              minReadsPerCytosine = 20,
                                              pValueThreshold = 0.01)
  
  checkTrue(all(overlapsAny(DMRs_CG_noise_filter_filtered, DMRs[3:5])))
  checkTrue(all(overlapsAny(DMRs[3:5], DMRs_CG_noise_filter_filtered)))
}

test_mergeDMRsIteratively <- function() {
  DMRs_CG_bins <-  computeDMRs(syntheticMethylationData1, 
                               syntheticMethylationData2, 
                               NULL, 
                               context = "CG", 
                               binSize = 100, 
                               method = "bins",
                               test = "score", 
                               pValueThreshold = 0.01, 
                               minCytosinesCount = 3, 
                               minProportionDifference = 0.4, 
                               minGap = 0, 
                               minReadsPerCytosine = 3, 
                               cores = 1)
  
  DMRs_CG_bins_merged <- mergeDMRsIteratively(DMRs_CG_bins, 
                                              minGap = 200,
                                              respectSigns = TRUE, 
                                              syntheticMethylationData1, 
                                              syntheticMethylationData2,
                                              context = "CG", 
                                              minProportionDifference = 0.4, 
                                              minReadsPerCytosine = 3, 
                                              pValueThreshold = 0.01,
                                              test = "score",
                                              alternative = "two.sided")
  checkTrue(length(DMRs_CG_bins_merged) < length(DMRs_CG_bins))
  checkTrue(all(overlapsAny(DMRs_CG_bins_merged, DMRs_CG_bins)))
  checkTrue(all(overlapsAny(DMRs_CG_bins, DMRs_CG_bins_merged)))
}

